A GIS-Based Approach for Shallow Landslides Risk Assessment in the Giampilieri and Briga Catchments Areas (Sicily, Italy)

Abstract

Shallow landslides pose a widely growing hazard and risk, globally and particularly in Mediterranean areas. The implementation of adequate adaptation and mitigation measures necessarily requires the development of practical and affordable methodologies and technologies for assessing the shallow landslides hazard and its territorial impact. The assessment of shallow landslide hazard maps involves two different and sequential steps: the susceptibility and the runout analysis, respectively, aimed at the identification of the initiation and the propagation areas. This paper describes the application in the Giampilieri and Briga Villages area (Sicily, Italy) of a shallow landslide risk process at a basin scale with an innovative approach in the runout assessment segment. The runout analysis was conducted using specific GIS tools employing an empirical–geometric approach at a basin scale. The exposure and vulnerability values of the elements at risk were assigned using a qualitative and semi-quantitative approach, respectively. The results highlight the effectiveness of the procedure in producing consistent runout hazard and risk assessments in the valley areas where the more important and vulnerable exposed elements are located. This study contributes to addressing the public administration demand for valuable and user-friendly tools to manage and drive regional planning.

1. Introduction

Landslides have significant adverse economic and societal impacts, leading to loss of life and damage to structures and infrastructures. In the past 50 years, landslide-related disasters have increased tenfold, with the risk set to rise due to two key factors: climate change and urbanisation [1]. Shallow landslides, such as debris flows, mudflows, and debris avalanches [2], are particularly influenced by the climate change, given the close relationship between intense rainfall and landslide triggering [3,4,5]. Heavy rainfall events are intensifying due to recent decades’ warming [6] aligning with increased atmospheric moisture [7]. Consequently, climate change is expected to impact shallow landslides activity differently based on local climatic features [8,9,10]. Climate projections anticipate an escalation of shallow landslide hazards in tropical areas, where the intensity of cyclones and monsoons is on the rise. Similarly, the Mediterranean region may experience comparable effects due to intense cyclones, known as Medicanes, acquiring tropical characteristics accompanied by extreme winds and rainfall [11]. This phenomenon results in both larger and more frequent occurrences [12]. In mountainous and continental climate areas, drier summers may lead to a decrease in debris supply from the cliffs driven by thermoclastism, subsequently reducing the material necessary to trigger shallow landslides [13]. Conversely, in temperate and dry climate areas, warming scenarios could be associated with landslide-driven destructive impacts owing to the increase in extreme fire and post-wildfire extreme rainfall events [14].

Similar to the correlation observed with the climate, the connection between urbanisation and landslide risk is robust, and also, the incidences of shallow landslides triggered by human activity are on the rise [15]. The urbanisation processes impact both high-income countries, where industrial settlements are expanding, and low- and medium-income countries, where unplanned urban settlements continue to sprawl [16]. The abandonment of hilly and mountainous settlements resulting from urbanisation also instigates processes that directly influence slope stability [17,18,19,20]. Quantitatively assessing these changes is a critical aspect for fostering landslide risk consideration into urban planning tools [21].

In this context, the development of practical and cost-effective methodologies and technologies for monitoring and assessing the hazard of shallow landslides and its territorial impact is imperative to implementing necessary adaptation and mitigation measures [22]. Quantifying shallow landslides risk as interactions among environmental, climate, and human factors is a challenging task, but it is increasingly feasible thanks to improved satellite and remote sensing observations and open data initiatives. These initiatives provide accurate geospatial information about local rainfall, land cover, digital elevation models, and landslides [23,24,25]. Despite technological advances and innovative approaches [26,27], shallow landslide risk assessment tools often face high levels of uncertainty, connected primarily with the difficulty of modelling a highly complex physical process. Anticipated future climatic and socioeconomic changes, including increased urbanisation and resulting land use changes, contribute to further complicating uncertainties in predictive scenarios [28].

An open challenge in shallow landslides hazard analysis on the basin and regional scales is the runout assessment. Debris flows, mudflows, and debris avalanches typically originate from steep slopes, often in sparsely populated areas. The displaced material can travel downstream over considerable distances from the initiation areas, frequently encroaching upon human settlements, roads, and critical infrastructures [29]. Nevertheless, landslide hazard and risk maps currently implemented by regional and local authorities often feature solely an inventory of past landslides and/or results of a landslide susceptibility analysis, focusing primarily on potential source areas, where landslide failures may occur. Regional hazard maps seldom offer results of a shallow landslides runout analysis, and the identification of the runout areas is frequently overlooked. Integrating the susceptibility analysis with runout-focused approaches could facilitate a more detailed hazard and risk assessment and contribute to strengthening the efficacy of the mitigation measures [30].

Two main types of procedures are widely used for the shallow landslide runout investigation: empirical–statistical expressions based on geometric parameters correlations and analytical/numerical models that simulate physical process parameters [31,32,33]. The choice between the two approaches depends on the study area’s extension and objective, involving different details, scales, and efforts. Although recent studies have demonstrated its applicability in more extensive scenarios [34], the numerical modelling approach is particularly well suited for characterising a single event at the slope scale and designing slope stabilisation works. Empirical/geometric methods rely on simple correlations between certain landslide parameters and runout and can be applied at the basin and regional scales to analyse wider areas [35]. Despite the extreme simplification of the flow dynamic, empirical/geometric methods can reliably predict the extent of shallow landslide propagation [36]. Furthermore, according to more recent guidelines, they can effectively support decision makers in landslide risk assessment and land planning [37].

The significant parameters of the other two controlling factors of the risk, i.e., exposure and vulnerability, are also heavily contingent on the scale considered.

The elements at risk are typically categorised into four main asset types [38]: physical (buildings, cultural heritage, and transportation networks); social (population); economic (productive sites and agricultural areas); and environmental (land cover and protected areas). The exposure assessment process of delineating the quantity, quality, and spatial distribution of the elements at risk within a study area heavily varies at the local [39], basin/regional [40], national [41], and global scales [42]. Official large-scale maps serve as the primary source of spatial information for structures and infrastructures, and social characteristics are usually acquired by demographic databases. Where specific and detailed data on the elements at risk are lacking, land cover types may serve as proxies [43]. Future landslide risk scenarios, subject to variations due to social and climate changes, can be assessed by projecting the distribution of the elements at risk calibrated using a historical dataset [44].

Vulnerability refers to the conditions influenced by the physical, social, economic, and environmental factors that shape a community’s susceptibility to the impact of hazards [45,46,47]. It is widely agreed on that the current risk has increased to higher physical exposure, as well as to lower adaptive capabilities of existing settlements and communities. Nevertheless, in contrast to other natural risks, there is no common approach for assessing vulnerability in communities prone to landslides. This poses a significant obstacle in establishing the necessary tools for effective emergency and disaster management of these phenomena. Various methods are employed to analyse vulnerability, including empirical approaches that provide discrete values for a range of landslide intensity [48,49]. All approaches heavily rely on historical data specific for the region and the type of landslide [50]. Several factors contribute to the landslide impact on buildings, including material, age, height/size, and foundations [51], and the different types of damage, combined with landslide intensity, can be described using Boolean, qualitative, semiquantitative, and, with greater difficulties, quantitative variables [52]. Even when information on past landslides is available, details of the damage response of specific elements at risk to a particular type and magnitude of the process are often missing [53,54]. When detailed measurements are accessible, the development of an extensive database allows the extrapolation of significant fragility curves and the vulnerability assessment for buildings of different typologies [55,56,57,58].

In this study, a GIS-based process is applied for shallow landslide risk analysis in two small river basins on the Ionian east side of the Peloritani Mountains (Sicily, Italy). An innovative GIS tool for shallow landslides runout assessment and intensity evaluation in propagation areas was trialled. The innovative GIS tool aims to bridge a gap in the shallow landslide hazard assessment process at the basin scale, enhancing risk assessments in valley areas where urban settlements are concentrated. The reliability of the runout assessment process is tested matching estimated and observed data and evaluating the relative differences. Additionally, a comparison between the study results and current hazard and risk maps was performed to highlight the opportunities provided by the novel procedure.

2. Study Area

The study area encompasses the Giampilieri and Briga basins, spanning 20 km² between the Peloritani mountain ridge and the coast, situated 15 km south of Messina (Figure 1). Over recent years, both basins have experienced significant landslide events, notably on 25 March 2007, 1 October 2009, and 1 March 2011. The intense rainfall on 1 October 2009 led to widespread slope failures and triggered numerous shallow landslides, as documented by several studies [59,60]. Small villages and the transport system in the southern municipalities of Messina, as well as in the Itala and Scaletta Zanclea areas, were severely damaged by these landslides (Figure 2). Regrettably, the event, known as the “Messina flood”, also resulted in multiple fatalities and injuries.

The Peloritani Mountain Belt constitutes the southern extremity of the Calabride–Peloritan Arc in Nord-East Sicily. It is comprised of southeast-verging layers of Pre-Alpine metamorphic unit nappes overlapping sedimentary Maghrebid units [61]. Outcropping formations are mainly composed of micashists of various metamorphic grade and secondly by the deposition of sedimentary covers. The rapid crustal uplift, starting from the upper Pliocene and lower Pleistocene, has significantly shaped the overall high energy relief geomorphology of the region. The mountain sectors exhibit “V”-shaped valleys with steep slopes, eroded by torrent-like straight watercourses, creating narrow and deeply embanked features within high rock walls. Wide and over-flooded valleys appear in the coastal parts with extensive alluvial fans or dejection cones, ruins of old cliffs, and fluvial and marine terraces. The low geomorphological evolution, typical of a recent uplift, fosters intensive erosion processes and causes a widespread instability condition along the slopes. The alteration and degradation of crystalline lithotypes, coupled with the high erodibility of the sedimentary deposits, determine a continuous sediment supply from hillslopes, promoting soil formation. Brown soils, leptic cambisols, and eutric cambisols are present in the area, usually sandy with low clay content (approximately 10%), and clearly show tixotrophic characteristics [62]. The soil thickness on the slopes ranges from 50 to 80 cm, influenced by the preservation status of the terraces.

The Ionian side of the Peloritani Mountains has undergone significant urban development since the post-Second World War period (Figure 3). The study area now features a continuous urban landscape along the coastline, where the population migrated in past decades, leaving the outback original settlements. The number of buildings along the coast has significantly increased over the years due to the presence of the railway, established since 1866, and the Messina–Catania motorway (A19), inaugurated in 1971. Inland and higher areas, used in the past for grazing and agriculture, are now covered by woods, shrubs, and meadows. In the hills and valleys, areas devoted to various agriculture activities prevail. Currently, the east coast is characterised by the presence of several critical linear infrastructures, including roads (with S.S. 114 being the major one), the A19 motorway, the Messina–Catania railway, water (e.g., the aqueduct serving Messina and managed by AMAM), gas pipelines, and power and communication lines. These infrastructures are located side by side on a narrow strip of land between the hills and the sea, often less than 100 m wide.

An analysis of the land use maps of the Corine Land Cover project [63,64] indicates evolving trends in land cover and land use that have occurred in recent times, from 1990 to 2018, characterised mainly by (

Table 1. Changes in land use classes within the study area based on Corine Land Cover data from 1990, 2006, and 2018. The last column highlights changes relating to land use between 2018 and 1990. The cells are coloured red and green to, respectively, indicate a decrease or increase in land use in percentage for each category. The code classes 313 (mixed forest) and 323 (sclerophyllous vegetation) were not mapped in the 1990 dataset; thus, the variations of those classes are related to the time interval 2006–2018.

Code	Description	1990		2006		2018		Δ 2018–1990
		Km2	%	Km2	%	Km2	%	%
111	Continuous urban fabric	0.43	2.08	0.41	1.97	0.68	3.29	57%
112	Discontinuous urban fabric	0.42	2.03	0.42	2.04	0	0.00	−100%
222	Fruit trees and berry plantations	2.88	13.88	2.88	13.91	2.92	14.15	1.2%
223	Olive groves	0.48	2.31	0.48	2.31	0	0.00	−100%
243	Land principally occupied by agriculture	3.09	14.87	6.75	32.56	2.71	13.16	−12%
311	Broad-leaved forest	3.24	15.62	3.24	15.65	1.95	9.47	−40%
313	Mixed forest	0	0.00	1.33	6.43	0	0.00	−100%
321	Natural grasslands	0.16	0.76	1.16	5.57	0	0.00	−100%
322	Moors and heathland	5.72	27.55	0	0.00	0	0.00	−100%
323	Sclerophyllous vegetation	0	0.00	0.65	3.14	6.63	32.14	919%
324	Transitional woodland-shrub	3.58	17.24	3.41	16.43	0	0.00	−100%
333	Sparsely vegetated areas	0.74	3.56	0	0.00	5.73	27.79	674%

; Figure 4):

• a defragmentation and homogenisation process of the territory through a reduction of categories from 11 in 1990 to 6 in 2018;
• an increase of continuous urban settlements (57%) and the disappearance of the discontinuous urban fabric;
• the overall stability of the areas with fruit trees and berry plantations; in contrast, a noticeable disappearance of olive groves, which were previously visible in the hills along the eastern coast;
• the reduction of the territories covered by broad-leaved forests (−40%) and occupied by agriculture (−12%);
• the disappearance of natural grasslands, moors, shrubs, and mixed forests and heaths and bushes in favour of sclerophyllous vegetation (919%) and areas with sparse vegetation (674%).

3. Materials and Methods

This study implemented and tested a landslide risk assessment process, with innovative elements in the propagation analysis, a focal aspect of the shallow landslides hazard evaluation. The suitable thematic data and a landslide inventory were extracted from a geodatabase implemented by ENEA in the years following the “Messina flood” for the northeastern sector of Sicily. The empirical/geometrical approach was applied to identify the areas prone to landslide propagation. These areas were then discretised for intensity distribution based on the assessment of velocity and kinetic energy. The risk map of the study area was produced through a GIS overlay among areas with different hazard values and exposure and vulnerability maps. The exposure and vulnerability assessments were based on qualitative approaches derived from the literature.

3.1. Dataset

The landslide inventory of the study area contains several hundreds of rapid shallow landslides triggered by the heavy rainfall that occurred on 1 October 2009, identified through aerial images observation and field surveys [65]. One hundred elements were extracted from the inventory to constitute a representative subset of all known rapid shallow landslides of the area.

High-scale aerial images, acquired ad hoc after the event of 1 October 2009 by the Civil Protection Department, were used for geomorphological interpretation.

A 2 m cell size Digital Terrain Model (DTM; ATA flight 2012–2013), freely provided by Sicilian Regional Geospatial Portal (https://www.sitr.regione.sicilia.it/geoportale/it/Metadata/Details/946, accessed on 15 February 2024), constitutes the basis for the morphometric analysis.

In the absence of a detailed soil map, the 1:50,000 geological map of the Messina Province [61] was used as a proxy to geotechnically discretise the shallow material covering the bedrock. The fourteen geological formations outcropping in the study area were distinguished into six lithotechnical classes: granular A and B, stratified A and B, massive, and crystalline. The unit weight values, among other lithotechnical parameters, have been attributed to the six classes according to the geotechnical report of the General Master Plan of the Messina Municipality (

Table 2. Unit weight values for each of the lithotechnical units identified in the study area.

Lithotechnical Unit	Unit Weight Value
Granular A	16.5
Granular B	19.1
Stratified A	20.5
Stratified B	20.0
Massive	21.0
Crystalline	22.0

).

Level 3 of the Corine Land Cover (CLC) dataset from 1990 to 2018 was acquired and formed the basis for the discretisation of the exposure. Twelve distinct classes were identified among the artificial surfaces (continuous and discontinuous fabric); the agricultural areas (fruit trees, olive groves, and land principally occupied by agriculture); and the forest and semi-natural areas (broad-leaved and mixed forest, grassland, moors, sclerophyllous vegetation, transitional wood and shrub, and sparsely vegetated areas). Further detailed information was taken from the 2012 Regional Technical Map (1:10,000) of the Sicilian Regional Administration, focusing on buildings of public interest and the electricity grid, as well as from Open Street Maps, particularly concerning the road network and railway [66].

A soil thickness map (Figure 5) was produced using the GIST model approach [67], a semi-objective complex process based on the definition of the four indices of the equation:

h_i = −K_c ∗ C ∗ η ∗ ψ⁻¹

(1)

where h_i represents the soil thickness in each cell of the DTM, K_c is a calibration parameter that adjusts the normalised values of the other indices to real thickness values, C is an index linked to the slope curvature, η is an index linked to the position along the slope profile, and ψ is linked to the critical slope threshold. The soil thickness map production is based on the elaboration of

• the DTM through the “r.slope.aspect” and “r.flow” tools, the GRASS algorithms available in the QGIS free open source suite;
• the slope values;
• the lithotechnical characteristics in correspondence with the inventoried debris flows.

A 114 field measurements dataset, located with “Garmin GPS MAP 65”, was used to identify the K_c parameter and calibrate the equation.

3.2. Hazard Assessment

The selected shallow landslide dataset was divided into two sub-datasets (Figure 6) for the propagation assessment phase: 50 phenomena were used to calibrate the site-specific parameters of a runout equation. The same specific parameters were applied to the remaining 50 initiation areas used as the validation dataset. These latter were selected with care to avoid bias in the validation results. Particularly, the landslides with the scarp, the initiation area, the path, or the toe not always evidently shown in the aerial images, as well as phenomena with particularly lobed shapes, a long path clearly channelled into the drainage network, and evident interactions with anthropic structures, were excluded.

The 50 training landslides dataset was used to extract the parameters of the simply Legros runout equation [68]:

L = aV^b

(2)

where L is the runout planar distance, and V is the volume of the initiation area. Through a regression analysis of the runout distance and the volume (R² = 0.7514), the site-specific parameters of the Legros equation were defined: a = 9.3935 and b = 0.4592 [69] (Figure 7).

Subsequently, a GIS tools package specifically implemented for the assessment of the shallow landslide runout was applied to the validation dataset of 50 landslide initiation areas. This GIS tools process, mainly based on the DTM elaboration and on a supplementary dataset, aimed to rapidly evaluate the runout of shallow landslides at the basin scale in a very short time. The propagation dynamics were greatly simplified by introducing necessary approximations. The tools were implemented using several algorithms available in the processing instruments library of QGIS.

The first step of the process (Max Value Points script) allows to find the points with the highest altitude within each initiation area. These points have been considered as the landslide initiation points and the origin of the flow paths. Successively, the three Water Drop Paths scripts generated, firstly as a raster and then as a new vector layer, show the paths that a drop of water starting from the Max Value Points would take over the land surface until reaching the lower part of the DTM area. Applying the site-specific parameters of the Legros equation, the fifth script (Legros runout) cuts each line of the drop water path layer at the specific estimated runout distance. Once the runout path of each landslide has been delineated, the sixth script (Kinetic features) assigns a velocity and a kinetic energy distribution along the paths. Based on field measurements, the tool assumes that, starting from 0 m/s in the initiation areas, the mobilised masses accelerate up to the point of peak velocity (8.4 m/s), reached at one-quarter of the path, and then decelerate. Assuming a constant volume along the path and a homogenous rheologic characteristic of the mobilised masses, the kinetic energy variation along each shallow landslide’s path was estimated. A 10 m buffer of the drawn path of each flow was produced, attributing different kinetic energy values, discretised into four classes (

Table 3. Classification into four intensity classes of the twenty sectors of the polygons of each landslide path according to the kinetic energy and potential level of damage.

Kinetic Energy Value (kJ)	Potential Damage Level	Intensity Class
0–500	Low	1
500–2000	Medium	2
2000–6000	High	3
>6000	Very high	4

), to twenty different sectors of the polygons of each landslide path. Using this as a proxy of the impact force, the analysed territory was discretised by distinguishing areas with different susceptibilities to be reached by the mobilised material with varying destructive capacities (Figure 8).

The automatised process executes the six scripts in less than an hour.

3.3. Exposure Assessment

The description of the exposed elements started from Level 3 of the 2018 CLC map, in which the study area is exclusively divided into 6 classes. Based on the observation of Google Maps images (updated to 2021), some categories of exposed elements have been integrated or modified in accordance with the CLC nomenclature guidelines [70]. In particular, the historic centres of the hilly villages were integrated in the “urban areas with continuous fabric” CLC class due to the urban structures and the network transport occupying more than 80% of the territories. Furthermore, some areas occupied between 30% and 80% by artificial structures were drawn as “urban areas with a discontinuous fabric”. Integrating CLC and other available data, sixteen different types of exposed elements were identified: continuous and discontinuous urban fabric; four typology of public interest buildings (schools, churches, cemeteries, and other community buildings); cultivated areas; wooded areas; pastures and shrubs; tracks; roads; motorway; railway; gas station; transmission; and communication tower (Figure 9;

Table 4. Extension, in km² and percentage of occurrence, of the exposed elements classified according to the sixteen categories identified in the study area.

Exposed Element Categories	Area (km2)	Area (%)
Continuous urban fabric	0.7455	3.616
Discontinuous urban fabric	0.4301	2.086
Schools	0.0030	0.015
Churches	0.0040	0.019
Cemeteries	0.0105	0.051
Other community buildings	0.0005	0.002
Cultivated areas	4.9694	24.100
Wooded areas	1.9309	9.364
Pastures and shrubs	11.9026	57.724
Tracks	0.2364	1.146
Roads	0.2703	1.311
Motorway	0.0658	0.319
Railway	0.0424	0.206
Gas station	0.0015	0.007
Transmission tower	0.0067	0.033
Communication tower	0.0004	0.002

). A relative value from 1 to 4 (lower to higher importance) was assigned to each class of exposed elements based on qualitative evaluations about the social and economic importance (

Table 5. Classification of each class of the exposed elements in the study area with exposure values from 1 to 4 (lower to higher importance) considering qualitative evaluations of the social and economic relevance.

Data Source	Former Category	Exposed Elements	Exposure Value
CLC 2018	Continuous urban fabric	Continuous urban fabric	4
CLC 2018	Discontinuous urban fabric	Discontinuous urban fabric	3
RTM 2012	School complexes	Schools	4
RTM 2012	Religious complex	Churches	4
RTM 2012	Cemeterial complex	Cemeteries	3
RTM 2012	Social complexes	Other community buildings	4
CLC 2018	Fruit trees; land principally occupied by agriculture.	Cultivated areas	2
CLC 2018	Broad-leaved forest	Wooded areas	1
CLC 2018	Sparsely vegetated areas; sclerophyllous vegetation.	Pastures and shrubs	1
OSM	Track	Tracks	1
OSM	Unclassified; tertiary; residential; primary; pedestrian.	Roads	2
OSM	Motorway	Motorway	3
OSM	Railway	Railway	3
RTM 2012	Fuel Distribution supply	Gas station	3
RTM 2012	Electricity pylon	Transmission tower	3
RTM 2012	Antenna for telecommunication and broadcasting structures	Communication tower	3

).

3.4. Vulnerability Assessment

Unfortunately, a detailed description of the buildings located in the study area, in terms of construction materials, building height, and/or foundations depth, was not available in vector format. Similarly, a database with the surveyed damage data from shallow landslides events that occurred in the past in the area was lacking. Consequently, a simplified qualitative approach was used by attributing a degree of loss of the exposed element for different phenomena intensities according to some literature data and considerations.

Slightly modifying some empirical qualitative approaches, four different types of damage [48,50] and four corresponding classes of loss degree [47] were distinguished (

Table 6. Vulnerability classification according to different typologies of damage and the relative loss degree.

Damage Typology	Loss Degree	Vulnerability Class
Superficial	0–0.25	1
Functional	0.25–0.5	2
Structural	0.5–0.75	3
Total	0.75–1	4

). An estimated resistance category was assigned to the exposed elements consisting of buildings [50] (

Table 7. Estimated resistance classification of the exposed buildings according to their constructive category and structure.

Exposed Elements Categories	Building Structure	Resistance
Cemeteries and gas station	Mixed structure (timber masonry and/or similar)	Weak
Continuous urban fabric, churches	Masonry, mixed structure (masonry/reinforced concrete)	Medium
Discontinuous urban fabric, schools, other community buildings	Reinforced concrete	Strong
Towers	Steel, Reinforced concrete	Very strong

).

According to the authors [50], a vulnerability class from 1 to 4 (lower to higher vulnerability) was attributed to each building category and to roads for the four different intensity levels (

Table 8. Vulnerability classes from 1 to 4 (lower to higher vulnerability) attributed to exposed categories for the four different intensity levels.

Exposed Elements Categories	Low Intensity	Medium Intensity	High Intensity	Very High Intensity
Continuous urban fabric, churches	1	2	3	4
Discontinuous urban fabric and schools	1	1	2	3
Cemeteries	1	2	3	4
Cultivated areas	3	4	4	4
Wooded areas	1	2	2	3
Pastures and shrubs	1	1	1	1
Tracks	1	1	1	2
Roads	1	2	3	4
Motorway	1	1	2	3
Railway	2	3	4	4
Gas station	1	2	3	4
Transmission and communication towers	1	1	2	3

). The remaining non-built exposed elements were evaluated by referring to the vulnerability values of the buildings. The different vulnerability classes were attributed to the exposed elements based on the different expected responses of each category to the stresses induced by shallow landslide phenomena. Among others, the following considerations guided the attribution of the vulnerability values:

• buildings of public importance (i.e., schools) and/or strategic buildings and discontinuous urban fabric are supposed to have better quality building characteristics (e.g., due to a more recent construction age, according to specific Technical Codes and Regulations) than buildings in historical centres and churches [71], as well as in the more ancient villages that constitute a continuous urban fabric;
• viaducts/bridges and tunnels make motorways less vulnerable than ordinary roads, because those artefacts are generally not impacted by mobilised materials in cases of shallow landslides;
• the fragility of train tracks makes railways more vulnerable than ordinary roads.

3.5. Risk Assessment

The study area was discretised in terms of the risk value according to the classic risk equation:

R = H ∗ E ∗ V

(3)

where R is the risk value, and H, E, and V are, respectively, the hazard (Figure 8), exposure (Table 5), and vulnerability (Table 8) values. All these values are discretised into classes from 1 to 4. The consequential set of values (from 1 to 64) was further divided into 4 homogeneous classes by attributing the value 0 to the portions of territory not affected by landslides.

4. Results and Discussion

There are substantial differences in the observed landslide dimensions in terms of the geometric and kinematic parameters, sometimes exceeding even one or two orders of magnitude for area, volume, and energy.

4.1. Runout Assessment Result

Comparing observed and estimated runout values using site-specific parameters allowed us to evaluate the performance of the Legros equation (Figure 10). A good forecasting capacity is expressed by the statistical measurement (R² = 0.7369), considering, however, that the estimated runout distances can vary between 0.5 and 2 times the corresponding observed debris flows paths. The extension of the estimated prone areas to the propagation of shallow landslides constitutes a modest portion of the entire study area (1.83%); within this area, the hazard classes are distributed quite evenly (

Table 9. Absolute and relative extensions of the estimated prone areas to a propagation of shallow landslides. The percentage values of the four hazard classes refer only to the total runout hazard areas.

Hazard Classes	Low Intensity	Medium Intensity	High Intensity	Very High Intensity
Km2	0.105	0.091	0.087	0.095
%	28	24	23	25

).

It is important to note that the use of the Legros equation in the GIS tool was somewhat constrained. In fact, while the Rickenmann formula [72] yields excellent results in the back analysis, it cannot be applied predictively. In formulas by other researchers [73,74], which demonstrate predictive reliability equal to or even superior to that of Legros [68], the slope break point serves as a crucial geometric parameter. To date, automating the identification of the slope breakpoint has proven elusive. Consequently, the implementation of GIS tools proceeded using the Legros formula. Ongoing efforts may focus on resolving this challenge, aiming to automate the identification process and potentially incorporate alternative formulas into the GIS tools.

4.2. Hazard, Exposure, Vulnerability, and Risk Map Results

The results show how the areas with higher risk levels are concentrated in the continuous urban fabric zones, where more important and vulnerable exposed elements are situated. The risk level is heavily influenced by the vulnerability class attributed to each exposed element and depends on the intensity value of the expected events. This is clearly shown in the cases of the two debris flows toward Giampilieri Village (Figure 11, Figure 12, Figure 13 and Figure 14).

It is crucial to note that fragility curves for shallow landslides typically emerge from correlating damage to exposed elements with the thickness and/or velocity of the mobilised material. In contrast, the four intensity levels used in this study express exclusively the speed and the kinetic energy of the overall mobilised mass. Future improvements in GIS tools for the runout should prioritise estimating the thickness of the deposit along the path. This enhancement would enable attributing the intensity value in a more representative manner through a matrix considering both the velocity and deposit thickness. The availability of an exhaustive damage database would significantly strengthen the implementation of more meaningful fragility curves for the various exposed element typologies.

4.3. Comparison with the Current Hazard and Risk Maps

Comparing the risk map developed in this study with the hazard and risk maps from the Sicilian River Basin Master Plan (Piano di Assetto Idrogeologico—PAI), produced and issued by the Regional District Basin Authority [75], provides insights into the consistency and reliability of different methodologies used for landslide hazard assessment.

The seven Basin Authorities currently existing in Italy are the institutions charged for drawing up the land planning and intervention programming tools to face and mitigate hydraulic, landslides, and drought risks [76]. The Basin Authorities employ various methodologies, including qualitative matrix, geomorphological, quantitative statistical, and mixed methods, all based on parameters from the available landslide inventory, such as the type of movement or activity level.

The landslide hazard assessment of the Sicilian DBA currently relies solely on the landslide inventory that has occurred in the area, lacking a probabilistic analysis of the evolution of the slopes. For the intensity of the landslide phenomena, the speed and the kinetic energy are considered, which estimations are, respectively, based on the landslide typology and on the relationship between the areal dimensions of the landslide and its typology. Given the intrinsically neo-formative characteristics of shallow landslides, the current River Basin Master Plan of the Sicilian DBA provides only partial and not particularly detailed information regarding the landslide hazard of the region. The DBA hazard maps show where phenomena occurred in the past but not where new shallow landslides may be activated and in which areas they may propagate.

The landslide risk map of the regional DBA identifies four classes of elements at risk (from E1 to E4), with increasing vulnerability values. The assessment of the vulnerability of these elements takes into consideration the probability of them being affected by a landslide event, the presumed value of the element at risk lost in the potential event, and the possibility of endangering people’s lives. Only effects relating to the physical sphere are considered, i.e., the degree of loss of one or more elements due to the natural event. By conservatively assuming the vulnerability of each exposed element to be equal to 1, the DBA risk map is an expression of the formula R = H ∗ E.

A comparison between the DBA and current study hazard maps (Figure 15) is aimed at representing the innovative opportunities provided by the proposed approach, taking into account that

• the DBA hazard map is not properly a predictive instrument for this landslide typology, being only a photo of past occurred events;
• due the correspondence with the inventory of landslides already occurred, the detail level of the DBA hazard map is inversely proportional to its predictive capacity;
• the methodology proposed in this study is aimed at producing real hazard maps, i.e., predictive tools, starting from areas not (yet) affected by past phenomena, defined through the susceptibility analysis;
• the current study map is the result of a modelling process that inevitably provides an approximate output.

The comparison evidenced distinct hazard levels assigned to the two landslide areas of Giampilieri Village. Although characterised by a level of detail consistent with the modelling approach, the runout assessment used in this study allows for a more detailed delineation of varying hazard levels along the flow path. Similar evidence emerges by a comparison between the corresponding risk maps (Figure 16), with the additional consideration that the DBA has adopted a precautionary approach in attributing risk.

Both the compared map outputs are completely independent of any triggering factor, such as rainfall. Even the runout GIS tool does not consider rainfall or other parameter rainfall-dependent (e.g., rheologic characteristics) as an influencing factor, apparently reducing its usefulness and efficacy and missing an objective achieved by other methodological approaches [77,78,79]. Nevertheless, starting from the potential source areas, the runout GIS tool process may provide different hazard outcomes according to different susceptibility maps that are directly dependent by the intensity of the triggering factors; the greater the triggering factor intensity, the greater the initiation areas activity and, consequently, the greater the volume of the mobilised material of each shallow landslide. Consequently, hazard and risk maps produced applying the runout GIS tool may (indirectly) be dependent on the triggering factor intensity (e.g., rainfall) that acts on the susceptibility maps in terms of the extension and class of the prone areas. Therefore, the proposed approach well suits studies focused on different triggering factors, such as the evaluation of the cascading effect determined by the earthquake–landslide–debris flow disaster chain [80].

5. Conclusions

The shallow landslides risk analysis approach applied in this study is integrated with the second segment of a complete hazard analysis process at the basin scale. Even if, in this study, the runout assessment process used initiation areas from past events as starting data, in a typical back analysis approach, the proposed process is equally suitable for using potential initiation areas derived from susceptibility maps. The main findings of the study are

• the specific runout toolset, implemented with QGIS, significantly facilitated the runout analysis, involving the execution of numerous algorithms sequentially to identify the landslide path, the runout distance from several potential initiation areas, and the intensity assessment. The automatic execution of this chain of operations as a single process provides a time- and effort-saving advantage in GIS analysis. The quick execution time of the runout assessment process guarantees that the implemented tool is well suited for integration into land planning and landslide early warning systems;
• although the runout distances calculated with the Legros equation exhibit a certain level of approximation, linked to the extreme simplification of the process-conditioning factors, the estimate allows quantitative reasoning about the hazard level in the valley zones downstream of the initiation areas;
• through a qualitative risk evaluation matrix, the analysis allowed to discretise the territory into qualitative different classes. However, the introduced approximations in volume, speed, and, consequently, energy estimation and its distribution along the path make the process a preliminary approach;
• landslide hazard in the Messina area, which has experienced widespread fires in the past, often faces a generalised high risk due to specific urbanisation and the limited availability of hydro-geomorphologic risk reduction plans. The approach used in this study, considering the relative variability of the vulnerability level, provides results with greater detail than the more precautionary maps of the Sicily Region.

Further improvements will be addressed:

• to implement and integrate alternative formulas into the GIS tools;
• to improve the vulnerability assessment process;
• to collect multiple sets of landslide events triggered by different rainfall events of different intensities for a specific study area, with the aim of identifying different site-specific runout parameters for different rainfall scenarios.

Despite its apparent limitations, the process has proven to be particularly valuable in meeting the demand of public administrations for straightforward tools to manage and guide regional planning. Acknowledging an inherent level of uncertainty, the application of the methodology to other areas prone to shallow landslides will be instrumental in refining and strengthening the accuracy of the process. Further investigations and enhancements of the method are expected to contribute to fill the current gap in defining hazard maps for shallow landslides. Given the reliance on public datasets, the methodology can be successfully exported to other contexts and countries where decision makers are actively engaged in and committed to landslide risk assessment. This adaptability enhances the broader applicability and utility of the approach in diverse geographical and administrative settings.