Geomatics, soft computing, and innovative simulator: prediction of susceptibility to landslide risk

1.   Introduction

In the international context, the detailed understanding and categorization of terrain based on their distinctive geomorphological features have long been regarded as critical aspects of geographical analyses. The creation of maps aimed at delineating landslide susceptibility zones represents a significant advancement in geospatial analyses, offering a comprehensive view of terrain vulnerabilities. These maps serve as invaluable tools, providing continuous assessments of slope conditions and enabling proactive measures to mitigate risks associated with landslides. The process of developing such maps involves the integration of various data sources, including satellite imagery, topographic surveys, and field observations, to accurately capture the intricate terrain characteristics. The resulting susceptibility maps not only highlight areas prone to landslides, but also offer insights into the magnitude of potential hazards posed to communities, infrastructures, and ecosystems. This information plays a crucial role in guiding land-use planning, infrastructure development, and disaster management efforts, thus allowing decision-makers to prioritize resources and implement targeted mitigation measures. The principal problem in the creation of such maps is the possibility to have upgraded and accurate systems based on near-real-time data. Besides the mapping traditionally conducted with geographic information systems technologies ^[1,2], artificial intelligence techniques are used to help better create these kinds of maps, and several authors in the literature have addressed this. For example, in various contributions ^[3,4,5], Nwazelibe et al. have demonstrated advanced mapping technologies for testing the performance of different algorithms in the landslide susceptibility mapping of Nigeria's Udi Province, which is a region known for incessant soil erosion and landslide events.

In regard to the Italian context, one of the most dangerous environmental risks is associated with the problem of landslides, which have occurred in Italy with increasing frequency in recent years, making the risk of landslides a real and widespread problem. For this reason, various bodies, and research institutes, such as the Research Institute for Protection (IRPI), have focused on studying and monitoring the evolution of landslide susceptibility. This institute provides in-depth information on landslides and contributes to understanding the evolution of risk in the country. Among its activities, the IRPI deals with the analysis of landslide phenomena, the development of prediction and monitoring models, and the promotion of risk mitigation and management measures. Additionally, the Higher Institute for Environmental Protection and Research (ISPRA) played a significant role in the analysis of this risk. In Italy, according to the ISPRA classification, the most affected regions include Valle d'Aosta, Liguria, Calabria, Campania, and Molise, with the latter presenting a significant risk, with 16% of its territory at risk of a landslide. These regions are characterized by geologically unstable terrains and climatic conditions that can increase the risk of landslides ^[6,7,8]. The report on hydrogeological instability in Italy 2021 published by ISPRA provided an overview of the danger linked to landslides, floods, and coastal erosions throughout the country. Furthermore, it presented risk indicators relating to the population, buildings, businesses, and cultural heritage. This report was compiled by ISPRA as a part of its institutional tasks of collecting, processing, and disseminating data regarding soil protection and hydrogeological instability in Italy. To obtain a complete map of landslide hazards across the entire national territory, ISPRA adopted a mosaic approach that combined information relating to the different hazard classes of the Hydrogeological Management Plans (PAI). This process involves the harmonization of the legends used, dividing the territory into five categories: very high danger (P4), high danger (P3), medium danger (P2), moderate danger (P1), and areas of attention (AA). According to the 2020–2021 mapping, many areas of the country, especially those classified as P3 and P4, presented a high risk of landslides. This is the result of a range of factors, including unstable terrains, coastal erosion, and climate change. The new national hazard maps were created using the PAI and hydraulic hazard maps according to the directives of Legislative Decree 49/2010. These maps considered the updates provided by the District Basin Authorities. Overall, 18.4% of the Italian territory was classified as having a high landslide hazard and either very high or medium hydraulic hazard (with a return period between 100 and 200 years). Compared to the previous edition of the report, a percentage increase of 3.8% of the surface classified as having a high and very high landslide danger and of 18.9% of the surface of medium hydraulic danger was observed. This increase is attributable to an improvement in the information provided by the District Basin Authorities, who have conducted detailed studies and mapped new landslides and recent flood events. Moreover, the report presented updated data on the dynamics of the Italian coasts in the period 2007–2019. It highlighted that 19.7% of the coasts were advancing while 17.9% were retreating. Despite the progressive increase in coastal defense structures, there was an increase in stable and advancing coastlines compared to the 2000–2007 period, with a 1% decrease in eroding coastlines. However, the situation varied regionally, with some coastal regions being more prone to erosion. At the level of Italian municipalities, 93.9% are at risk from landslides, floods, and coastal erosion. There are 1.3 million inhabitants at risk of landslides and 6.8 million inhabitants at risk of floods. The regions with the highest number of populations at risk are Emilia-Romagna, Tuscany, Campania, Veneto, Lombardy, and Liguria. Furthermore, almost 548,000 families are at risk from landslides and over 2.9 million from floods. Of the more than 14.5 million buildings in Italy, over 565,000 are located in areas with a high and very high landslide danger (3.9%) and over 1.5 million are in floodable areas in the medium scenario (10.7%). Moreover, the report presented data on the presence of industries and services in areas with a high and very high landslide danger, with over 84,000 companies and 220,000 employees at risk. Additionally, there are over 640,000 workers exposed to flood risks in the medium scenario (13.4% of the total). Finally, the report included an estimate of the cultural assets at risk of landslides and floods, with over 12,500 assets potentially subject to landslides and almost 34,000 assets at risk of flooding, ^[9]. It follows that monitoring landslide risk areas is essential to prevent disasters and adopt timely mitigation measures.

There are various technologies and methodologies, geomatics in particular, to monitor these particularly susceptible areas. These methods, which mainly consist of on-site measurements, include traditional topographic monitoring, which involves the use of topographic instruments such as total stations and levels to measure changes in the terrain topography and detect even small changes in the ground surface, the Global Positioning System (GPS) method, which allows you to monitor ground movements precisely and continuously, making it particularly useful for monitoring slow and continuous deformations, and the most sophisticated and complex use of interferometric radar (InSAR, with all the its methods of use DInSAR, PSInSAR, ...), which uses satellite radar images to measure ground deformations, detecting millimetric changes in the ground surface. ISPRA uses different methods and approaches to evaluate the stability and danger of the slopes and calculate the susceptibility index to landslide phenomena by providing a Geographical Information System (GIS) with updated information on the susceptibility index derived from an information and statistical analysis of geological, geomorphological, and hydrogeological data, taking the various factors that can influence the stability of the slopes into account, such as the terrain morphology, geology, vegetation, rainfall, terrain slope, and other relevant parameters. Therefore, it is undeniable that in this context, monitoring meteorological conditions (with particular reference to rainfall) is essential, since intense rainfall can trigger landslides that are particularly dangerous for high-risk areas such as those adjacent to urban areas.

The constant vigilance of weather forecasts and the recording of current precipitation are crucial elements to recognize potentially dangerous situations. Mondini et al. ^[10] highlighted the importance of predicting and managing the risk of rainfall-induced landslides. Currently, both empirical models based on rainfall thresholds and models are used to predict the short-term occurrence of rainfall-driven shallow landslides. The article proposed the use of a strategy based on deep learning to link rainfall to the occurrence of landslides. The results of the study's proposed system indicated that it is possible to effectively predict the occurrence of rainfall-induced landslides over large areas and that their location and timing are mainly controlled by rainfall. This opens the possibility of developing operational landslide forecasts based on precipitation measurements and quantitative weather forecasts. However, these approaches may not be effective for operational forecasting over large areas.

Equally important is a geotechnical analysis: geotechnical deformation meters can be used to monitor geotechnical data, which are directly installed in the ground and can detect variations in the soil pressure, tension, and compression of the ground, providing valuable data on the state of the ground or seismometers and accelerometers to detect seismic tremors that could trigger landslides. Equally useful are the piezometry sensors that monitor the water level in the ground, which is a probable indicator of landslide risk if there is a sudden increase in the water level. It represents another type of approach and involves the measurement of parameters such as soil cohesion and porosity. This information is crucial to understand geotechnical conditions and to identify any significant changes. In parallel, in-depth geological analyses of the physical and geotechnical characteristics of a given area constitute an essential starting point to identify areas exposed to a potential risk of landslides. In this sense, some studies ^[11] have conducted a quantitative analysis of debris flow in a meteorized gneiss through a methodology aimed at identifying the key aspects related to risk management.

Therefore, to identify the most susceptible areas, the most common methodologies focus on how the deformation progresses over time ^{[12,13,14,15,16]}. While deformation progression is a crucial aspect of the problem, other factors must also be evaluated, including materials, failure geometry, human causes, deforestation, and weather events. As these factors change over time, research has focused on developing prediction models based on artificial intelligence techniques, in particular using various neural network configurations.

In this context, this article aims to present an experimental methodology capable of detecting susceptibility to landslide risk through the use of a WebGIS system with a Decision Support System (DSS) connotation; by exploiting cellular automata and soft computing techniques, this allows us to predict variations in the susceptibility index provided by ISPRA on large sample areas downstream of intense precipitation phenomena, even localized ones. In particular, starting from input data (such as the DEM Digital Terrain Model, degree of soil saturation, and ISPRA susceptibility index), using precise rainfall data (obtained from an innovative atmospheric simulator), discretizing the terrain (cellular automata), and using an appropriately optimized neural network, it will be possible to obtain the updated and predicted susceptibility index as an output, compared to the known one present in the susceptibility maps provided by ISPRA, all managed in a WebGIS environment.

To a certain extent, this system allows us to anticipate the possibility of imminent landslide events. The proposed system offers a number of key advantages in landslide risk management and land planning. First, it provides a long-term risk assessment, thus allowing the continuous monitoring of ground conditions and landslide susceptibility in a specific area. This long-term perspective is essential for prudent urban and regional planning, thus contributing to the sustainable management of resources and protection of the environment. Furthermore, the system supports informed decisions through accurate data and dynamic simulations. This information can influence urban planning, thus limiting development in high-risk areas and promoting responsible land management. Additionally, a significant advantage is the constant update on landslide risk, which allows authorities and experts to adapt mitigation strategies in real time, thus improving emergency response capacity.

2.   Materials and methods

The fundamental approach underlying the operation of the proposed method is based on the creation of a Digital Twin that reflects the surrounding territory, including all of its physical and topographical characteristics. This process employs the three-dimensional component of GIS by incorporating orographic data and a DEM. In particular, this approach involves the representation of the surfaces and volumes that constitute the terrain using three-dimensional cellular automata. Such cellular automata are equipped with specific state variables, derived from microphysics rules, which regulate their mutual interactions (Figure 1). Subsequently, this representation of the territory is subjected to simulations, with each iteration representing a specific moment in time. During this process, the system dynamically evolves, with each cellular automaton seeing its state variables changed with each iteration. Then, these variables are examined by a Pattern Detector, which is implemented using a neural network. This detector can identify any variations in the landslide risk susceptibility index, if present.

The proposed methodology is divided into four phases: STAGE Ⅰ, STAGE Ⅱ, STAGE Ⅲ, and STAGE Ⅳ (Figure 2).

• STAGE Ⅰ: This phase includes the development of the digital twin of the area in question by collecting data relating to the terrain, meteorological conditions, and all the variables relevant for the risk analysis. The WebGIS will process this data, allowing for the visualization of the information associated with the maps, creating real-time interactions, and managing the data with maximum precision and speed.

• STAGE Ⅱ: This phase involves the use of an innovative atmospheric simulator that emulates the behavior of the atmosphere by discretizing the particles in cubic cells, following the Smoothed Particles Hydrodynamics (SPH) model for their interaction. The output of this phase will consist of precise values of the climate parameters and the estimate of variations in space and time, thus allowing meteorological events to be anticipated.

• STAGE Ⅲ: This phase includes the development of the WebGIS forecasting system, which begins with the discretization of the terrain using cellular automata described by appropriate physical variables that follow a specific interaction law. The prediction of the diffusion of the properties of the cellular automata will be carried out through the optimization of a neural network which will detect the change in the characteristic susceptibility index, if present.

• STAGE Ⅳ: This phase is fundamental to evaluate the effectiveness of the proposed method by comparing the simulations carried out with the historical data collected. The simulations via WebGIS were compared with events that actually occurred in Calabria (Italy) in 2005.

2.1.   Data collection

In phase Ⅰ, all of the appropriate and necessary data for the subsequent phases of the proposed methodology were collected. First, an open-source American Standard Code for Information Interchange (ASCII) DEM from the Calabria region was used, which is fundamental to obtain terrain elevation data. As the name suggests, the DEM is stored in an ASCII text format. Moreover, the information from the DEM made it possible to attribute a fundamental parameter to the cellular automata, which is the slope (attributed only to surface automata).

Then, the data provided by ISPRA in relation to landslides in Italy were obtained, with a particular reference to the susceptibility index. As is known, the ISPRA Susceptibility Index is a parameter used to evaluate the susceptibility of an area to slope instabilities, including landslides. This index is an integral part of the ISPRA "Charter of Nature" and contributes to the understanding of the ecological vulnerability of a territory. ISPRA uses several factors to calculate the Susceptibility Index, including the geological, geomorphological, and climatic characteristics of the area. These factors are analyzed and weighted to determine the degree of susceptibility of a specific area to landslide phenomena. The overall evaluation of the index considers multiple variables that can influence the stability of the soil. The functional relationship (1) underlying the calculation of the Landslide Susceptibility Index may vary slightly based on specific local conditions; in general, it can be expressed as follows:

where:

A, B, C, D, E are the weights associated with geological, topographical, climatic variables, and so on. These weights may vary depending on the methodology used and the local conditions.

G represents the geology of the area and considers factors such as lithology, the presence of faults or geological fractures, etc.

R represents the slope of the land and takes the slope of the terrain in the area into account.

S represents vegetation and takes the vegetation cover in the area into account.

L represents rainfall and considers the rainfall levels in the area.

C represents land use change and takes changes in land use in the area over time into account.

The value in question presents variations both within the area taken into consideration and over time, which is in line with the updates made by the agency. This value is of a considerable importance in the process of the proposed methodology, since it contributes to modeling and simulations through the inference carried out by the neural network, thus incorporating the knowledge of the models used by ISPRA.

Soil moisture is one of the main key factors in the landslide triggering process. The infiltration of rainwater or melting snow causes an increase in pressure in porous soils, with a consequent change in their consistency due to a decrease in cohesion and internal friction. During periods of drought, deformation phenomena are less evident; however, in soils with a significant presence of clay, cracks can form more easily, facilitating the infiltration of subsequent rainwater and causing a loss of cohesion in the soil. Similar effects can occur along riverbanks following rapid changes in water levels in the surrounding basins due to the drag forces of fine sand and silt grains, which can lead to liquefaction of the soil. One of the main physical causes of landslides is intense and prolonged rainfall. However, having only the weather data provided by the control units of the Calabria Region (7 km distance) available and wanting the most accurate measurements possible in a limited area, the surface humidity value (provided as input) was calculated through downscaling operations carried out from an appropriate atmospheric simulator developed by the Geomatics Laboratory of the University of Reggio Calabria, starting from the rainfall value recorded in the area of interest. In fact, the simulator has provided precise precipitation values in the study area, thus allowing a limited area to be taken into consideration and proceeding with a more precise and accurate spatial analysis.

With regard to the functioning of the atmospheric simulator (phase Ⅱ) used to obtain highly detailed information on the quantity of rainfall in the specific area under study, please refer to the authors' publication ^[17,18,19].

2.2.   Cellular automata

For the discretization of the terrain and, therefore, phase Ⅲ of model analysis, it was decided to use the cellular automata model ^[20,21,22].

Cellular automata are mathematical and computational models used to simulate the behavior of complex systems through the discretization of space and time. They were first introduced by mathematician John von Neumann and scientist Stanislaw Ulam in the 1940s. These models were later studied in detail by Stephen Wolfram and other scientists. They are composed of a regular grid of cells, each of which can be in a discrete state at a given instant of time. Each cell interacts with its neighboring cells according to a defined set of rules. These rules determine how the state of a cell changes over time based on the state of neighboring cells. The transition rules can be simple or complex, and the overall behavior of the system emerges from the combination of the rules of all the cells (Figure 3).

Cellular automata can be used in a variety of applications, including modeling natural phenomena and simulating complex processes. In the context of terrain discretization, cellular automata are used to represent terrains divided into small units or cells. Each cell in the case at hand has attributes such as elevation, slope, degree of saturation, and other relevant parameters. The interactions between these cells can be used to simulate the behavior of the ground in response to factors such as rain, erosion, landslides and more.

This method is particularly suitable for modeling complex spatial systems, where local interactions between units have a significant impact on the overall behavior of the system. Terrain discretization with cellular automata allows for a detailed representation of the landscape, making it possible to simulate geological and hydrological processes accurately. In the specific case, cellular automata are characterized by continuously evolving state variables during the simulation and by proximity interaction rules (local microphysics) between adjacent cellular automata. It should be noted that superficial cellular automata (i.e., those that intersect surfaces) have their state variables set by the inputs and not by proximity interactions.

The cellular automata used in this context to create a phenomenological Digital Twin of the area under study are three-dimensional, georeferenced, and have the following dimensions:

• Width: 1 meter.

• Thickness: 1 meter.

• Height: 1 meter.

State variables include:

• Slope, expressed in radians. Only relevant for shallow automata.

• Characteristic Susceptibility Index. Only relevant for shallow automata.

• Degree of Saturation Sr. $Sr$ = $Vw/Ve$, i.e., the volume of free water divided by the volume of voids (the spaces between waterproof materials).

• Sn, the number of hours out of a total of annual hours (8760) in which the soil had a saturation greater than 30%, considering the previous year as an instantaneous reference. The value varies from zero to one.

The selection of these specific variables is motivated by the fact that the output class, which classifies the susceptibility, is ultimately determined in our experimentation by artificial intelligence (AI) called to infer the law that establishes the correlation between the inputs (precipitation) and the output (change in the susceptibility index). Therefore, this choice is aimed at ensuring that the typically used neural networks work at their best. This approach does not aim to completely replace classical phenomenological prediction methods with the cellular automata simulator, but rather to serve as an intermediate layer that makes the most of the inferential capabilities of AI, especially neural networks, in the specific context.

Returning to the parameters previously defined and used in this work, we observe the following:

• The slope provides the neural network with information on the orography that conditions landslide events.

• The characteristic susceptibility index indirectly reports information obtained from classical methods, thus contributing to convergence during the learning phase.

• The degree of saturation plays the role of the geological parameter, which can vary abruptly (at each iteration), thus significantly influencing the susceptibility as it modifies the stability of the terrain and its intrinsic characteristics, such as roughness ^[23].

• Finally, the Sn index offers a simple and effective way to inform the neural network about the historical evolution of the phenomenon, avoiding the use of recurrent neural networks and the related architectural complexities, while including the cumulative effect of precipitation in the overall model in time.

Regarding the law used to model the local interaction law of cellular automata, the Green-Ampt Infiltration Model was selected, based on a simple law, which was focused on moisture diffusion.

The Green-Ampt method ^[24] is used to estimate water infiltration into the soil during a rainfall event. The basic mathematical relationship in the Green-Ampt method is used to calculate the initial infiltration rate (f) during the initial stages of a rainfall event. This initial infiltration rate is critical to understanding how water penetrates the soil. The main mathematical relationship (2) for calculating "f" is as follows:

where:

"f" is the initial infiltration rate (rate at which water penetrates the soil at the start of the rainfall event).

"${K}_{s}$" is the hydraulic conductivity of the soil in saturated conditions.

"${\theta }_{i}$" is the initial soil moisture content (before the onset of the rainfall event).

"${\theta }_{s}$" is the soil moisture content at saturation (maximum water holding capacity of the soil).

"t" is the time elapsed since the start of the rain event.

"${t}_{0}$" is a corrective term for the time that takes the initial soil conditions into account.

It is important to note that the Green-Ampt method simplifies some of the complexities of the soil behavior and that several variations of the method can be used.

Once the phenomenological Digital Twin of the area of interest had been created, we opted for the implementation of a pattern detector using a neural network. The choice of the neural network focused on the model known as Self-normalizing Neural Networks (SNN), which was enriched with a layer that makes use of the Scaled Exponential Linear Units (SELUs) activation functions ^[25,26,27].

2.3.   Self-Normalizing Neural Networks with SELU

The SELU activation function represents an important development in the field of artificial neural networks, as it was designed to overcome some of the limitations associated with the Rectified Linear Unit (ReLU) activation function and its variants.

A distinctive aspect of SELU is its ability to self-normalize. During the process of training a neural network, weights and biases are adjusted such that the network's outputs have a mean close to zero and a standard deviation close to one. This self-normalization process is critical to deal with the vanishing gradient problem, which can significantly manifest in deep neural networks during training. In other words, SELU helps keep gradients stable, thus allowing for the creation of deeper and more efficient neural networks.

Comparing SELU to the ReLU function, which is one of the most common activation functions, we notice some key differences. The ReLU is defined as $f\left(x\right) = max(0, x)$, which means it returns zero for negative values and the input value for positive values. SELU, on the other hand, is defined more complexly as $f\left(x\right) = scale* \left(max\right(0, x)$$+min(0, alpha* \left(exp\right(x) - 1)\left)\right)$, where "scale" and "alpha" are positive constants (Figure 4). This additional complexity comes from the need for self-normalization. Among the advantages of SELU, we find its ability to improve convergence and stability in deep neural networks.

The structure of the SNN is as follows:

• One input layer (DensData) with 16,800 inputs.

• One layer with linear activation function (16,800 neurons).

• A DropOut layer with a 30% activation rate, which is used only during the training phase to prevent overfitting.

• One layer with SELUs activation function.

• A linear output layer with a single output.

For the creation of the training set, a strategy was adopted that integrates with the overall process. Initially, areas such as the main study area were identified, considering variables such as lithology, volumetry, and relief morphology, together with the annual distribution of rainfall. For each of these areas, different years were selected, with the aim of including at least one landslide event, even if it was minor. Subsequently, simulations were run with a number of iterations corresponding to one year, using the precipitation history as the input. This process allowed us to populate the variable called "Sn", which represents the historical evolution of soil saturation. The periods of maximum rainfall were recorded, and if no landslides were reported during these periods, the values of the cellular automata simulations (in the 24 hours following the rainfall) were used as the input for the training set, associating the index of ISPRA characteristic susceptibility as the output. It is important to highlight that the values of the state variables coming from the shallow and deep cellular automata contribute differently to the neural network due to their specific geographic locations. This aspect further enriches the learning process of the neural network ^[28,29].

Finally, in phase Ⅳ, to conduct a complete performance analysis of the results obtained (both with the atmospheric simulator and with the proposed method), various analyzes were performed. With regard to the atmospheric simulator presented in this study, it was necessary to compare the average parameters obtained from 100 tests performed via the atmospheric simulator with those obtained from the control units of the Calabria Region. The results indicate an average discrepancy varying between 5% and 10%, which can be considered non-significant for the objectives of the proposed application. This process of comparing data between simulations and real measurements is fundamental to evaluate the accuracy and reliability of our simulator and to guarantee the validity of the results obtained in the simulation of the phenomena examined.

2.4.   Index validation of percentage difference

With regard to the validation of the results obtained, it was decided to conduct the simulation using data referring to the previous year in relation to which the landslide event actually affected the area under study, in such a way as to verify the correctness of the possible change in the predicted susceptibility index compared to the known ISPRA one. Furthermore, a mathematical analysis was conducted using the Index Validation of Percentage Difference (IVDP) method, which is used to compare the calculated susceptibility indices (${S}_{calculated}$) with the observed or historical ones (${S}_{observed}$) and takes the percentage discrepancy between the two indices into account.

1. For each area in which there are calculated susceptibility indices (${S}_{calculated}$) and observed or historical ones (${S}_{observed}$), the percentage difference (DP) between the two indices was calculated (3):

The percentage difference represents how much the calculated index differs (in percentage) from the observed one.

2. The average of the percentage differences on all the points or areas considered (4) was calculated:

where ${\mathit{\Sigma}} \left(DP\right)$ represents the sum of the percentage differences for all points or areas and n is the total number of points or areas considered.

3. A low value of IVDP indicates a good agreement between the calculated and observed susceptibility indices, while a high value indicates a significant discrepancy.

3.   Case study and results

The study area was located in Favazzina (Figure 5), which is a fraction of the municipality of Scilla, in the region of Reggio Calabria, Italy. The aim was to make a comparison between the susceptibility index calculated by ISPRA with the susceptibility index provided as the output from the entire methodology. In this area, two recent landslide events, classified as debris flows, were recorded. These events originated as translational landslides involving the surface and altered material of the metamorphic substrate. The landslides occurred at positions higher than the channels and generated debris flows as the soil eroded and additional material was added to the movement.

The most significant landslide events in this area date back to May 2001 and March 2005, both with grave consequences involving several vital infrastructures.

On May 12, 2001, two surface translational landslides occurred at the head of the Favagreca River. These landslides originated at heights of 567 meters and 558 meters above sea level, corresponding to two incisions. The unstable masses merged at approximately 300 meters above sea level, forming a primary channel before reaching the SNAM gas pipeline station, the SS 18 main road, and the railway. This event led to the derailment of the Turin-Reggio Calabria intercity train.

On March 31, 2005, a similar incident took place in the valley near Favazzina. Three shallow translational landslides occurred at approximately 370 meters, 242 meters, and 170 meters above sea level, transforming into debris flows. This debris flow caused considerable damage to the transportation infrastructure, including the SS18 state road and the railway, thus resulting in the derailment of the ICN Reggio Calabria-Milan intercity train.

Geometrically, the triggering areas of the May 12, 2001 event had a prismatic shape, with a sliding surface located at a depth of about 1.5 meters. The study area is situated in a region with a Paleozoic crystalline substrate, exhibiting intense alteration conditions. The lower-middle areas of the slope have highly and moderately altered rocks, while completely altered rocks predominate above 300 meters above sea level. About 60% of the area is covered by susceptible debris flow materials, known as class Ⅵ gneisses (Figure 6).

Loose alluvial deposits with gravel and sand are present along the main waterways, and beach deposits composed of sand and gravel were observed between the sea and the base of the slope. The map indicates that the trigger areas of the 2001 debris flow mainly involved class Ⅵ gneiss, which was still visible on-site. However, in the 2005 trigger areas, rocks from classes Ⅵ, Ⅴ-Ⅳ, and Ⅲ emerged, possibly due to the 2005 debris flow removing much of the previously visible class Ⅵ rocks.

Tectonically, the study area is intersected by fault segments primarily oriented in NE-SW and WNW-ESE directions. The main NE-SW fault system gradually aligns northwestward, influencing the Favazzina slope's morphology. The older WNW-ESE oriented fault system, which is morphologically less evident, contributed to the formation of the hydrographic network, including the Favagreca canal flowing towards the coastal plain of Favazzina. Figures 6 and 7 depict the geological map and widespread landslide areas of the studied region, sourced from the National Geoportal.

Input data were initially acquired at a given instant, including the following: the surface moisture value calculated by scaling from spatial scaling operations performed on the precipitation value recorded in the area of interest; land use (forests, crops, burned areas...) to be encoded as an appropriate state variable for surface cellular automata; and the characteristic susceptibility index provided by ISPRA (Figure 8).

Subsequently, the domain of interest was represented using 3D cellular automata, as shown in Figure 9, where it is possible to visualize the three-dimensional view of the DEM and observe how the cellular automata intersected with the surface.

Next, the simulation was started using the cellular automata, whose simulated states represented the inputs for the described neural network. This network returns as the new value (if changed) of the susceptibility index as the output. In fact, the output of the pipeline is provided by the pattern detector. The values of the state variables of the individual automata come to be the input of the latter module of the susceptibility classification process. As said in the methodology section, the pattern detector was an implemented SNN equipped with a layer with an activation function SELU. The only existing output produces the simulated susceptibility index.

The classification values used in this study were drawn from the ISPRA dataset. Specifically, a numerical scale was employed to designate the severity of a landslide risk within various areas. When a warning area was identified, it was attributed a value of 1.0 on this scale. For regions characterized by a moderate level of warning, the assigned value was escalated to 2.0; similarly, for areas posing a medium warning, the value increased to 3.0. Higher degrees of risk were denoted by assigning values of 4.0 and 5.0 for areas with high and very high warnings, respectively. Moreover, to enhance the versatility of this classification system, two additional categories were incorporated. First, for areas deemed not to require immediate attention, a value of 0.0 was assigned. Second, an index value of 6.0 was designated when evidence of recurring patterns leading to landslides under similar circumstances existed. The implementation of the proposed methodology yielded insightful outcomes for the area under investigation. Specifically, the analysis conducted resulted in a computed value of 5.8, as illustrated in Figure 10, which serves to underscore the severity of landslide risk within this region.

3.1.   Discussions and validation

To validate the results obtained, it was observed that the iterations corresponding to the date of March 29, 2004 returned an index of 5.8 in the platform (Figure 7). This value indicates a significant deviation from the indices provided by ISPRA (4.9), thus signaling, a priori, a marked deterioration in the state of stability of the slopes in 2005, which actually occurred in the same year with the consequent landslide recorded. We decided to call the new obtained index "Fast reclassification index", as it allows us to obtain a quick reclassification of the area under investigation through the data input provided. The representation is typical of the GIS environment with a polygon representing the cellular automata having inherited that value. Finally, the IVDP showed that for the results obtained in the last months of simulations before the landslide event, the percentage value was on average around 3.2% and 4.7%, showing a good adherence between the known values and the calculated and simulated values. The adherence of simulated results to known data suggests that the model used is reliable.

4.   Conclusion

The results obtained suggested the need to further explore the prospects offered by the presented technique. Although the study was conducted on a specific location, the validation period of the results was extended to one year, in accordance with expectations. A signal of worsening of the stability conditions emerged precisely around the time of the landslide events, thus confirming the temporal relevance. Future investigations should concern both the understanding of the micro-interactions between atmospheric conditions in a more detailed way and the consideration of a complex and non-uniform distribution of lithology in the study areas. The main objective of our work was to introduce an innovative methodology to predict susceptibility variations, which has the advantage of allowing the use of mathematical models with emergent properties and especially of artificial intelligence with neural networks, which is a field often considered complex according to literature. Further refinements of the method will be dedicated to the integration of contextual information derived from the geological characteristics of the terrain, seismic data, and the impact of human activities. In regard to the construction of the training set, we will exploit monitoring tools for suspicious slopes and the SBAS-InSAR and PS-InSAR techniques ^[31,32] to expand the number of areas suitable for providing useful information for identifying any ongoing deformations, thus associating an estimated value of susceptibility to avoid training the network only with boundary values. In addition to the initial parameters mentioned, it would be equally significant to consider the land use and seismicity of the area. However, in an initial testing phase, we decided to exclusively focus on a few variables handled by a simple diffusion relationship. The aim of this work was to present the method, with the intention of carrying out more detailed and accurate future implementations. These future implementations will not neglect factors such as geotechnical mechanical properties, geomorphology, the seismicity of the terrain, and all possible natural and anthropogenic causes that can contribute to landslides. Finally, it is important to underline that the choice to have a single classification output is guided by the desire to test the effectiveness of the method in general terms and to confer, in this work, a greater reliability on the temporal domain compared to the spatial one, considering the complexity of the model and the richness of the training set.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare no conflict of interest in this paper.