Soil Erosion and Landslide Susceptibility Mapping in Western Attica, Greece: A Rock Engineering System Approach

Abstract

Two of the eight main soil degradation processes with which soils worldwide are confronted are soil erosion and landslides. Specifically, landslides are a major threat in particular areas across Europe, often leading to serious impacts on population, property, and infrastructure. Regarding the abovementioned processes, the case study of the fatal Mandra flash flood (November 2017) in the Attica Region (Greece), which caused 24 deaths, and much infrastructure and building damage, is presented with the intention of assessing the relationship between soil erosion and landslide incidents. Investigations were executed from 2018 to 2022, and their outcomes were taken into consideration by the Technical Authority of the Attica Region. Soil erosion lines were delineated in a GIS and were validated using a previously generated regional Web-GIS landslide susceptibility map. The study presents soil erosion types from the Mandra fatal flash flood event and correlates them with already existing landslide susceptibility analyses for the Attica Region. The produced susceptibility map is a cartographic product on a regional scale (1:100,000) generated via a semiquantitative heuristic methodology named the Rock Engineering System (RES). The way in which both soil erodibility and landslide susceptibility maps were generated and validated could be the basis for proposing modeling approaches that can respond to new developments in European landslide policies.

1. Introduction

Soil erosion is a crucial triggering factor in land degradation worldwide [1,2] and specifically at the European level, with serious financial implications. To this end, the European Commission’s Thematic Strategy for Soil Protection recognizes soil erosion as a serious threat to the European Union’s (EU) soil resources [3]. Focusing more on soil erosion types, soil erosion by water and gully erosion are two typical causes of land degradation that subsequently lead to slope failures. For these reasons, different stakeholders need easy access to soil data and information of various types and scales to assess the state of soils [4]. To date, many researchers have used a variety of qualitative and quantitative techniques with erosion models integrating GIS applications to cope with soil erosion and land degradation issues [5,6,7,8]. To be more specific, soil erosion prediction models have been used to predict the hazard of soil erosion [9,10]. In addition, the most common erosion models are the Universal Soil Loss Equation (USLE) and its revised version (RUSLE), which estimate long-term average annual soil loss [11,12,13]. Considering spatial distribution, USLE/RUSLE models have limitations, which have been addressed using geospatial applications [14]. On the other hand, soil erosion assessment in large-scale field measurements may have some drawbacks, such as being time-consuming, expensive, or nearly impossible due to limited resources [15,16]. In addition to this, soil erosion assessmentmodels such as RUSLE/USLE have some drawbacks when predicting sediment pathways from hill slopes to water bodies, and in gully assessment [17].

Moreover, regarding gully erosion, which is a major land degradation process and probably the most severe type of water erosion [18] (and this type of erosion appears to occur very often in the broader study area of Mandra city), gully erosion susceptibility mapping is a valuable tool, as it is beneficial for identifying the spatial probability of gully incidents.

Traditional statistical models have been implemented to estimate gully erosion susceptibility mapping, such as analytical hierarchical process (AHP) [19], frequency ratio (FR) [20], logistic regression (LR) [21], and weight of evidence [22], but the prediction accuracy is relatively low when using these methods [23,24,25]. However, gully erosion susceptibility mapping prediction accuracy has been greatly improved using machine learning algorithms compared with traditional statistical models [26,27]. In particular, different machine learning algorithms have been used for gully erosion [28,29,30,31].

To be more specific, Eustace et al. [28] used random forest (RF) to model the relationship between gully erosion and some characteristic parameters, and they found that the the predictive capability of the generated model was modest.

Rahmati et al. [29] compared the behaviors of seven machine learning models (RF, BRT (boosted regression tree), BRANN (Bayesian regularized artificial neural network), L-SVM (Lagrangian support vector machine), RBF (Radial Basis Function)-SVM, P-SVM (parallelizing support vector machine), S-SVM (structured SVM)) in order to model the gully presence at the Kashkan–Poldokhtar watershed in Iran. They concluded that RF, RBF-SVM, BRT, and BP-ANN (backpropagation artificial neural network) can be implemented as machine learning models for gully erosion susceptibility mapping.

Arabameri et al. [30] used the random forest model, and found that the precision accuracy of the confusion matrix for this model was equal to 88.89%.

Finally, Pourghasemi et al. [31] used the combination of SVM and ANNs to model a gully pattern, and they concluded that these two algorithms could be useful tools for gully erosion modeling.

On the other hand, the application of machine learning algorithms to gully susceptibility mapping in large areas using very high-resolution datasets also has limitations regarding computing efficiency (e.g., pixel-by-pixel prediction is needed). In such a case, the spatial resolution must be reduced to improve the abovementioned computing efficiency. As a result, this procedure may reduce prediction accuracy [32]. In

Table 1. Some characteristic methodologies for the assessment of different types of soil erosion.

Methodologies for the Assessment of Different Types of Soil Erosion	Researchers	Number of References
Universal Soil Erosion Equation (USLE)	Wischmeier and Smith (1978), Panagos et al. (2015)	[11,13]
Revised USLE	Renard et al. (1997), Panagos et al. (2015)	[12,13]
Frequency ratio (FR)	Conforti et al., 2011	[20]
Logistic regression (LR)	Conoscenti et al., 2014	[21]
Analytical hierarchical process (AHP)	Arabameri et al., 2018	[19]
Weight of evidence (WoE)	Arabameri et al., 2019	[22]
Machine learning algorithms	Eustace et al., 2011; Rahmati et al., 2017; Arabameri et al., 2019; Pourghasemi et al., 2017	[22,28,29,31]

, a summary of the above-described methodologies is presented.

Thus, to predict soil erosion types (gullies included) with high accuracy, especially in areas where severe phenomena (e.g., extremely heavy rainfall, such as that which happened in Mandra on 14–15 November 2017) may occur, a more practical methodology still needs to be invented.

In this context, as land degradation (e.g., landslides) is associated with soil erosion, it is believed that finding tools and methodologies to cope with slope failures automatically succeeds in addressing issues caused by soil erosion. This statement is strengthened by the fact that researchers have obtained more reliable soil erosion susceptibility outcomes by using slope failure events, and soil erosion conditioning factors have been used in landslide susceptibility prediction [33].

Based on the high intensity of the Mandra fatal flash (with short duration) flood event of 14–15 November 2017, the erosivity of the rainfall and the runoff were (among other parameters) the main causes loosening the soil (and the weathering mantle of decomposed–fragmented rocks), resulting in weakening the slopes and as a result leading to the mass movements of solid and semisolid materials [34,35].

In addition, fires, which have recently manifested in the broader area of Mandra, can be a catalyst for erosion and landslides [36]. Thus, the consequences of fires, the fact that in the broader area of Mandra improper urban expansion has taken place in the last 40 years, as well as the triggering factor of heavy rainfall of 14–15 November 2017, have intensified the extension of gullies (soil erosion) and landslides. Another characteristic point that relates landslides to soil erosion and particularly gullies in the study area is slope modifications via human interference (e.g., cut and fill or artificial drains) [17,37].

Several researchers have noted the correlation between slope failure phenomena and soil erosion processes [8,38,39,40,41,42,43,44,45,46,47,48,49]. In addition, Lin et al. [50], Kou et al. [51] and Deljouei et al. [52] emphasized the role of vegetation in the interaction of slope stability and soil erosion.

Regarding the examined area, two characteristic types of slope failures have been recorded—(i) earth falls and (ii) rock falls, according to the Varnes classification [53]—whereas in the broader area of the Attica Region, there are also debris falls, earth slumps, earth slides and rockslide–debris avalanches [53]. Analyses of these slope failures have been developed by Tavoularis et al. [54].

Based on the abovementioned, a case study from the Mandra fatal flash flood (which took place on 14–15 November 2017) in the Attica Region (Greece) is presented with the intention of exploring the role of soil erosion in relation to land degradation (e.g., landslides). The combination of the high rapidity of the rainfall episode manifested on 14–15 November 2017 along with the industrialized character of Mandra city constitutes the basic reason why this flood and study area were chosen. Investigations performed by different stakeholders have been executed between 2018 and 2022, and the outcomes of these have been taken into consideration by the Technical Authority (Directorate of Technical Works) of the Attica Region in designing and implementing a priori mitigation measures (for debris flows and rainfall-induced soil erosion processes) against potential upcoming extreme rainfall episodes. Using a variety of tools, soil erosion types are defined and delineated in GIS maps and afterwards validated using the regional Web-GIS landslide susceptibility map of the Attica Region (DIAS project), which was completed in June 2021 by a research team [54], implementing a semiquantitative methodology named Rock Engineering System (RES). This map identifies specific zoning areas that are more susceptible to slope failure. The way this landslide susceptibility map is generated can offer a basis for modeling approaches that can respond to new developments in European policies (e.g., data, maps, technical reports), such as those of the European landslide susceptibility map version 2 ((ELSUS v2) and ELSUSv2_six_datasets and metadata) [55], or can contribute to the improvement of large-scale assessments, which can further generate landslide hazard and risk maps. Thus, the objective of this study is to explore the relation of soil erosion to landslides using methodologies that have been implemented in landslide susceptibility modeling. As per the author’s knowledge, no one else has predicted (at least at the Greek and European levels) soil erosion susceptibility using a semiquantitative methodology such as RES. Thus, implementing methodologies that have already been used in landslide susceptibility mapping can help us identify and estimate soil erosion hazards.

Soil movement is a sign of a soil erosion hazard. This includes riverbank erosion, gully erosion, debris-falls, rock-falls, and general slope failures that can damage the environment and human beings [5]. Annually, thousands of lives are lost due to mass soil movements worldwide [56]. However, Blaschke et al. [57] mentioned that the impacts of mass movements on soil erosion and land productivity are underestimated in the literature. Thus, less attention has been given in the research to soil erosion due to mass movement. Most soil erosion hazards prevail during a rainy season or after heavy rain [58,59]. Therefore, realizing susceptibility to soil erosion is very important for mitigation and hazard–risk minimization.

The current article is organized as follows: A brief description of the fatal flash flood that occured in Mandra (November of 2017) is first presented. Then, the soil erosion types that arose in that area are described. In addition, landslide susceptibility for the Attica Region is also briefly analyzed via a semiquantitative heuristic methodology named the Rock Engineering System (RES). The outcomes (e.g., soil erodibility and landslide susceptibility map) of this methodology are depicted to validate the correlation between landslide occurrences and soil erosion events manifested in the Mandra area. Moreover, some characteristic mitigation measures that have been designed are addressed against potential upcoming extreme rainfall episodes. The paper is concluded with suggestions for future research.

2. Materials and Methods

2.1. Location of the Study Area

The landscape characteristics of the study area are very diverse, due to its large land cover. The whole area is a mosaic of forested areas with dense vegetation cover consisting mainly of pine forest, lowland agricultural areas, areas of residential density (small settlements, settlement areas, etc.) and areas of industrial activity. The study area also presents a great diversity of topography, with hilly outcrops alternating with hollows and lowland sections, crossed by a number of small and larger streams that contribute to the main streams of the catchment areas into which the study area is divided [60].

In western Attica, in which the wider area of interest can be found, the prevailing winds are westerly. Geomorphological uplift prevents the movement of winds to the east, resulting in the condensation of water vapor and the occurrence of significant rainfall. According to the precipitation map of Greece, the area of flooding is located in the region reflected by the annual precipitation curves of approximately 400 mm and 600 mm rainfall height [60].

The study concerns the area upstream of the residential area of Mandra and extends from the ridge of the adjacent mountain range (Pateras Mountain), with an altitude of 841 m, to the settlement of Mandra (Figure 1). The study area covers the catchment areas of the two streams Agia Aikaterini and Soures (covered area equal to 45 km²) located upstream of the settlement of Mandra (latitude: 38°4′30.1548″ N, longitude: 23°30′3.0636″ E, elevation above sea level: 79 m).

The altitude of the Agia Aikaterini basin ranges from 171 m to 838 m. The lowest altitudes are found in the eastern–southeastern part of the basin and the highest altitudes can be found in the western and southwestern parts. The gradients found vary, with the steepest gradients occurring to the west, east and north, and the gentlest gradients occurring to the northwest, northeast and south [60].

The altitude of the Soures basin ranges from 100 m to 674 m. The lowest altitudes cn be found in the eastern part of the basin and the highest altitudes in the western and northwestern parts. The gradients found vary, with the steepest gradients found in the westnorthwest, southeast and northeast and the gentlest gradients in the central and eastern parts of the basin [60].

According to the Koppen classification, the climate of the study area belongs to the Csa category, i.e., Mediterranean climate type with warm summers [60]. In addition, it is noted that: (a) The average air temperature ranges from 5.89 °C (January) to 26.36 °C (July). The minimum monthly temperature is 1.52 °C (January) while the maximum monthly air temperature is 32.04 °C (July). (b) The annual rainfall is 538.89 mm. The monthly rainfall ranges from 10.11 mm (July) to 85.88 mm (December).

A characteristic of the climate of the study area and the region of Western Attica (where the study area is located) is the high rapidity of rainfall. The rainfall duration–intensity–return period curves, constructed with data from the study area (Laboratory of Hydrology and Water Resources Utilization, NTUA), show a rainfall intensity of 33.2 mm/h for 2 h of rainfall (rainfall height 66 mm) with a return period of 100 years [60].

Taking into account the climatic Atlas of Greece in the period 1971–2000 constructed by the National Meteorological Service (EMY, http://climatlas.hnms.gr/sdi/ (accessed on 6 October 2023)) for the wider area where flood phenomena occured (Soures catchment area), the annual rainfall ranges from 396 mm in the east (Mandra) to 635 mm in the west (Mount Pateras) [60].

2.2. The Fatal Flash Flood Event of Mandra (11/2017)

On Tuesday 14 November 2017, heavy rainfall occurred in the wider area of western Attica (in Greece), which continued at high intensity into Wednesday 15 November 2017. The total amount of rainfall during the core of the event (Pateras Mountain) exceeded 200 mm over a period of 6 h, with the highest intensity mainly at 5:00 to 8:00 on 15/11, representing a very heavy and relatively short period of rainfall. A remarkable fact of the flooding episode was the great rapidity that was recorded. This rainfall corresponds to approximately 1/3 of the total annual rainfall of the western part of the basin (Mount Pateras) of the two streams and 1/2 of the annual rainfall of the eastern part of the basin in the area of the city of Mandra [62]. As a result, increased water flow and the consequent erosion and solid discharges entered the settlement of Mandra, causing extensive damage to houses and roads and the deaths of twenty-four (24) people. The study area (a part of it is depicted in Figure 2) covers the catchment area of the two streams Agia Aikaterini and Soures located upstream of the settlement of Mandra (e.g., the Soures catchment area is equal to 23.0 km² and the Agia Aikaterini catchment area is equal to 22.0 km²) [63].

2.3. Geological Setting

The main lithostratigraphic characteristics of the formations found in the study area, listed from the most recent to the oldest, are as follows: (a) alluvial deposits; (b) modern lateral deposits; (c) old alluvial deposits; (d) marls, clays, sandstones, travertines, marly limestones and conglomerates; (e) limestones; (f) dolomitic limestones [61]. In addition, the study area contains some faults and a well-developed hydrographic network. The geological mapping of the catchments of the two streams that run through and north of the settlement of Mandra have provided useful information about the geological formations, their specific physical characteristics and their tectonic stress. This information is decisive for determining the extent of the contribution of each formation to the creation of flood conditions and the increased supply of sediments of various sizes in the area. More specifically, the largest part of the catchments, which is the mountainous part, is made up of dolomites and dolomitic limestones, which are highly fragmented, resulting in crumbling and mantles of alteration. The fractures and foliations of the fragmented rocks decrease their mechanical strength and increase their capacity to generate loose geomaterials, which lead to the formation of soil through the process of erosion. In addition, the erosivity of rainfall and runoff were the main causes for the loosening of the weathered mantle of decomposed–fragmented rocks, resulting in the weakening of the slopes and, as a result, leading to mass movements of solid and semisolid materials. All these materials can be easily eroded, and in the event of heavy rainfall they can provide materials for transport by water through streams.

2.4. On Site Geological Findings

The smoothest parts of the catchment areas (lowland) consist of deposits, i.e., materials that have been deposited in these areas in an earlier geological time (Upper Pleistocene) (Figure 3), and their deposition process is similar to the process that deposited the sedimentary materials during the flood event of 15 November 2017 [61]. That is, it was found that there were several successive flood events in the area in the recent geologic past that created the deposits in the lowland sections. These materials—by water movement and by mechanical means—have been welded together, forming very cohesive horizons, which are impervious to water, but also very resistant to erosion. These cohesive horizons are found in the lowlands at depths ranging from a few centimeters (20–50), i.e., below the soil mantle, to a maximum of 1 to 3–4 m. The presence of these cohesive horizons prevents streams from digging their beds deeper, so that in the event of increased flow, their cross-section is insufficient to drain the water. Consequently, they overflow, and the water flows uncontrollably over the surface of the adjacent areas, resulting in the creation of many furrows or small depressions that erode the soil mantle, adding new material to the solids [61]. Based on the above observations, it is considered appropriate to provide interventions at critical points in the hydrographic network that will limit the potential for erosion caused by streams and the downstream transport of erosion materials. In addition, in the lowland parts of the basins, where streams are unable to create a bed of sufficient cross-section due to the cohesive underlying materials, or where the bed is almost nonexistent, it is necessary to define the stream and to regulate it by creating a bed of sufficient cross-section. This can be done either by excavation or by creating artificial raised walls on the sides of the streams. The geological engineering characteristics of the formations that form the basin of the Mandra area, as well as their susceptibility to erosion, were studied to estimate the solid transport. The assessment of solid transport is very important for the planning and management of protection measures. The approximation of the total solid transport of the considered catchments is based on the calculation of the loss of soil mass by means of the measured erosion forms, and on the other hand on the calculation of the total volume of the flood deposits, the extent of which has been mapped by means of high-resolution satellite imagery and in situ measurements [61].

2.5. Causes of the Mandra Flood

The main causes of the catastrophic flood were [61]:

(1)

Urbanization, which significantly increases flood risk as it leads to the impairment of critical stream cross-sections;

(2)

The complete disappearance of small watercourses or their conversion into roads or parking places;

(3)

The construction of sub-dimensioning projects, which cannot take the hydraulic load in cases of heavy rainfall, as in the case of Mandra;

(4)

The dramatic reduction in the capacity of the soil to absorb part of the rainwater;

(5)

The very high concentration of atmospheric precipitation in a small area (such as that of Mandra) and the intensity of the rainfall.

This localized rainfall, which occurred in the mountainous parts of the region, caused a flash flood in the catchments of the area, and specifically the streams that flow into the lowland area of Mandra and the industrialized area of Mandra. These streams activated alluvial riffles formed by the supply of sediments within the Quaternary, resulting in flooding with a large amount of sediment. The areas inundated by water are located on top of the geological formations of the ripples. The development of settlements perpendicular to the water flow and on these alluvial ridges cut off the smooth drainage of water to the sea, resulting in water flooding the settlements. The development of roads and railways perpendicular to the flow is also estimated to have cut off the smooth drainage, resulting in the floodwaters being amplified in width.

2.6. Soil Erosion—Erodibility

Soil erosion or soil loss is the amount of soil or rock material that becomes detached under the influence of rainfall and subsequent surface runoff over a given period [14]. In our case, the soil erosion process occured due to extremely heavy rainfall, runoff, and gravity [61]. Research claims that soil erosion is extremely sensitive to rainfall [65,66]. The process of water erosion leads to the decomposition, transportation and deposition of sediment in a separate place. Erosion depends on factors such as [61]: (i) the erosivity of the water (rainfall intensity, volume and kinetic energy of the flow); (ii) morphological gradients; (iii) the presence or absence of vegetation cover; (iv) the resistance of the soil to erosion (erodibility); (v) climatic conditions and (vi) human interference—natural soil erosion is locally accelerated with anthropogenic interventions (road construction in stream beds, embankments, etc.), which disturb the soil structure and produce loose materials. The total land loss during the flood event of November 2017 amounted to more than 150,000 m³ [61]. The heavy rainfall that led to the catastrophic flood of 15 November 2017 in the urban area of Mandra was also accompanied by extensive soil erosion in large parts of the examined catchments. Soil erosion by water includes two main phases: the detachment of soil particles and their subsequent transport by water flow. When the transport energy of the water is reduced (reduction in flow rate, reduction in flow surface slope, etc.), the deposition phase of the transported materials follows. The typical forms of soil erosion that have occurred in the Mandra catchments and have led to extensive solid transport are as follows [61]: (a) surface (sheet) erosion, (b) rill erosion, (c) gully erosion, and (d) slope and subslope erosion.

2.6.1. Surface (Sheet) Erosion

Surface erosion is caused by the detachment of aggregates or soil components by raindrops [61]. Some of the detached particles move through the hydrological network, while others are deposited lower down. The erosive capacity of raindrops depends on their kinetic energy and momentum, which is a function of the size of the droplets combined with the direction and intensity of the rainfall, and its final velocity. The resistance of soils to erosion (erodibility) also plays an important role in the action of raindrops. Non-weathered rocks and cohesive soils have a low susceptibility not only to surface erosion, but also to all forms of erosion. Soils covered by vegetation are largely protected from erosion. Soil loss also increases significantly with the increase in the slopes of the morphological relief. Surface erosion, although less visible, is an important factor contributing to soil loss due to the destruction of the soil structure, which is then easily washed away by surface runoff. Although the calculation of the total soil loss due to surface erosion in the Mandra catchments requires systematic measurements over a super annual period, it is estimated that its contribution to the production of solids during the phase of the intense rainfall episode (15 November 2017) was much lower compared to the following forms of erosion. The reason is that the topographically higher parts of the catchments, where the rainfall intensity is highest, are structured by carbonate rocks, which generally exhibit low erodibility. In addition, most of the basins are covered by vegetation, which acts as a protection against surface erosion.

2.6.2. Rill Erosion

Rill erosion is the development of surface (sheet) erosion and is caused by the concentration of rainwater, which cannot percolate down to the subsoil, in small soil cavities, which progressively overflow and form irregular longitudinal surface runoff [61]. As the volume of runoff gradually increases, both the carrying capacity of the already extracted material and the potential for new material to be extracted along the water path increase. The erosive power of water depends on the volume, velocity, and rate of supply of surface runoff, slope gradient and soil erodibility. This creates many small and shallow (<1 m) rills that follow the direction of the maximum slope of the soil surface, forming water flow axes of usually triangular crosssections. The resulting sedimentary material is coarser than that produced by surface erosion and consists mainly of medium-sized gravels and coarse sand. The development of rill erosion in the catchment area of Mandra is particularly pronounced in areas structured by the engineering geological unit of Holocene deposits, where vegetation is also absent (Figure 4).

2.6.3. Gully Erosion

This type of erosion is developed by the flow of water along the axes of the maximum slope gradient of the slopes, where gullies of more than one meter in depth are formed [61]. The mechanism of their formation is similar to that of rill erosion, but on a larger scale, due to the greater concentration of water and the increase in its erosive and transport energy. The cross-section of the trenches depends on the erodibility of the soil with depth. In easily corroded materials, there is usually continuous dredging, and the cross-section is typically triangular, but deep erosion slows down considerably when the bottom of the trench reaches a corrosion-resistant formation. In this case, erosion progresses mainly laterally, and the cross-section of the gully is U-shaped. In the catchments of the Mandra rivers (e.g., Agia Aikaterini and Soures), intense trench erosion developed after the heavy rainfall that occurred on 15 November 2017, which contributed significantly to the production of fine and medium-sized sediments. Gully erosion developed in the upper and mainly middle elevations of the two basins, especially along valleys with a strong gradient, where the increased water volume significantly increased its erosive energy. Typical forms of gully erosion in the Mandra catchments are shown, which were developed mainly within the Holocene deposits (Figure 5).

2.6.4. Slope and Subslope Erosion (Riverbank Erosion)

In the central beds of the streams of the two catchments of the Mandra, as the elevations and gradients decrease, gully erosion develops into slope and subslope erosion, taking the form of a riverbank erosion feature (Figure 6). The strong flow of stream water detaches and carries materials from the foot of the slopes downstream, widening or deepening the bed [61]. The amount of erosion depends on the force of the water, the slopes of the bed, the meanders and bends it creates, and the erodibility of the geological formations through which the stream flows. The presence of vegetation or the local accumulation of coarse material in the streambed often diverts the main flow of the stream towards the slope, resulting in severe erosion. Slope and subslope erosion can be localized, as they usually develop perpendicular to the flow direction of the stream. Subslope erosion is also created in the beds of streams, in places of abrupt changes in elevation, where small waterfalls are formed. These forms of erosion can be observed in almost all the main and some secondary branches of the hydrological network of the Mandra catchment areas. The intense runoff that caused the major flooding in Mandra in November 2017 led to extensive slope and subslope erosion, resulting in the large widening of stream valleys and the removal of large volumes of material that were transported downstream. The slope and subslope erosions in the Mandra basin primarily affect the Holocene deposits of the stream slopes, but also affect the Pleistocene deposits, with differential erosion of the less cohesive horizons.

2.7. Modelling of Soil Erosion

2.7.1. Susceptibility to Soil Erosion

To investigate the intensity of the erosion that can be seen in the examined area of Mandra, a particular methodology was implemented by the Regional Authority of Attica/Directorate of Technical Works, and the Hellenic Survey of Geological and Mining Exploration (H.S.G.M.E) [61], as follows.

First, the erosivity of surface runoff (such as that derived from the flash flood on 14–15 November 2017), as well as the erodibility of the geological formations in the broader area of Mandra, determine how susceptible the geoenvironment is to the creation of gullies [20,21], a type of soil erosion that appears in the examined area of Mandra.

A soil erosion map was derived by a research project that was accomplished by the Hellenic Survey for Geology and Mineral Exploration (H.S.G.M.E.) in 2018, under the auspices of the Technical Authority of the Attica Region, which supervised this project (e.g., checking, reviewing, and approving the deliverables (e.g., technical reports and cartographic products)). The process of how those soil erosion lines were delineated is described more analytically in the following lines. The soil erosion map was produced based on the synthesis of three individual maps: (a) the engineering geological map, (b) the vegetation map, and (c) the morphological (slope) map.

Concerning the first map, the engineering geological map examined and classified the formations that structure the area, according to their lithological composition and hydrogeological behavior. The influence of tectonics was considered using the existing hydrographic network. Geological data were taken from the literature and are supplemented by field surveys. For the soil erodibility map, the roles of lithological composition, infiltration capacity and permeability were considered.

Infiltration is a property of soils. Continuous infiltration can result in the saturation of the soil cover or the weathering of the mantle. In the area under consideration, this property is very important, because Neogene formations and the weathering mantle can be saturated by prolonged rainfall and then become loose and easily eroded.

There should be no confusion between the concepts of permeability (which is a measure of the ability of the formation to allow water to pass through its mass) and infiltration, i.e., the ability of water to percolate through the soil mass without necessarily moving to the depth required to reach the water table. Infiltration refers to soils and not to rocks, and usually, if infiltration is present, it is referred to as deep percolation. In soils, continuous infiltration leads to saturation. The role of the soil cover on the rocks, which is usually a decomposed, clayey mantle with permeability, is therefore very important. The existence and thickness of the ground cover were estimated from field observations, particularly in areas recently affected by extensive forest fires, and from the 1:50,000 scale soil maps of the Institute of Forestry Research [61,68].

Regarding the second (vegetation) map, land use is based on mechanical and hydrological characteristics, as well as vegetation to control the slope’s stability [69,70,71]. During rainfall, land use influences the behavior of the soil as well as the magnitude of potential geomaterial that is about to move. In this case study, the land use classification from the research of Tavoularis et al. was used [54], wherein five classes are identified: (i) barren areas, (ii) urban areas, (iii) forest areas, (iv) shrubby areas—natural grassland and (v) cultivated areas. For the requirements of this work, the barren area class was implemented to estimate the soil erosion susceptibility.

Finally, considering the third (morphological/slope) map, the value of the slope was set as 5%, because this value is mentioned in the literature as the limit of the processes of deposition of fractional materials under conditions of diffuse flow (mudflow) and sedimentation on the surface of the earth (creation of cones of craters). The slope map distinguishes between areas with a slope of more than 5%, where erosion processes are more intensive, and areas with a slope of less than 5%, where erosion processes are less severe [61,68].

Thus, the susceptibility to erosion of the geological formations of Mandra was assessed based on the observations derived from technical geological mapping of the erosion phenomena observed and recorded in the wider catchment area (from technical site investigations in 2017 [67], 2018 [61], 2022 [62]), combined with the physical and mechanical characteristics of the geological formations. The grading of the units was carried out based on their physical and mechanical characteristics and estimated permeability.

In the present case study, the most interesting unit (associated with the highest soil erodibility) was the Holocene deposits, which cover an area of 2.05 km² [61]. Those geomaterials derived from river deposits are brown to reddish-brown in color, and generally loose to slightly cohesive, covering the smooth parts of the catchment. They originate from the erosion of older formations and deposition via surface runoff and gravity. They comprise mixed phases of materials, with a corresponding average proportion of coarse to fine-grained materials and rapid changes in grain size during both horizontal and vertical development. In places, they may represent lateral or alluvial deposits. Their thickness varies locally from a few centimeters to two meters, but is usually on the order of 50 cm. The coarse phase consists mainly of mudstones, gravels, and coarse-grained sands of mainly dolomitic and calcareous composition. The fine phase, which is the binder, consists mainly of sandy to silty sands (CL, SC) with a variable but generally small proportion of clay. They are characterized as soft to medium cohesive soils and are rarely very soft or stiff. The plasticity of the soils is generally low. They are characterized by low to medium water permeability and very high erodibility.

The above (three) parameters are the criteria used for assessing the susceptibility to soil erosion, and assuming equivalent weighting and linear correlation between them, the thematic maps were superimposed and corresponding susceptibility maps of the Agia Aikaterini and Soures catchments were obtained. Moreover, the previously mentioned findings are in line with those derived in the Hydroment Associate’s site work (2021) undertaken on those two examined watersheds (Agia Aikaterini and Soures), which was executed under the auspices of the Region of Attica’s Authority, which is responsible for designing and implementing mitigation measures against flooding [62]. The combination of the three final maps produced an erosivity map, part of which is depicted as soil erosion lines in the following map (Figure 7).

The described soil erosion lines were validated not only on-site by technical visits executed by the Technical Authority of the Attica Region (2017, 2022), but also using a pre-existent web-GIS landslide susceptibility map, the creation of which was accomplished via a research project briefly entitled “DIAS” [54]. Specifically, the generated soil erosion lines overlaid on the landslide susceptibility map show that those lines achieve a consistent coincidence with the zones of “Extremely high (slope failure) susceptibility” and “Landslide”.

The abovementioned soil erosion types recorded during the in-field site investigation were validated through a previously generated landslide susceptibility map [54], which is briefly described in the next section. In this map, a strong correlation was found between the defined soil erosion lines and the landslide susceptibility zones that were produced via a semiquantitative methodology named the Rock Engineering System (RES).

2.7.2. Rock Engineering System (RES) Methodology

This approach is mainly based on the correlation of mechanisms between landslide parameters through a matrix table and uses slope failure factors that can potentially be identified during the preparation phase of a civil engineering project. The purpose of using RES is to generate a landslide susceptibility map after the estimation of a landslide instability index [54]. RES is a semiquantitative approach, and the principal tool for representing the selected parameters and their interaction mechanisms (Figure 8). It is based on a matrix that depicts the key parameters as leading diagonal terms, and their binary interaction mechanisms as off-diagonal terms. The RES was developed by Hudson in 1992 [72] to determine the interaction of different rock engineering parameters and calculate the instability index for rock slopes. Since then, it has been applied to different aspects of geotechnical engineering such as rock stability problems, landslide susceptibility analysis and rock engineering [73]. Thus, by implementing the RES methodology, experts select the most important causative and triggering factors responsible for slope failures, rate their interactions, estimate their weighted coefficients, and calculate the instability index of the examined slope.

The calculation of the landslide susceptibility index is based on the following actions: (a) selection and rating of parameters, (b) generation of interaction matrix, (c) rating of the interactions among the selected parameters, (d) construction of a cause–effect diagram (the use of a cause–effect diagram can help the expert to understand the role of each independent landside factor in slope stability), (e) calculation of the weighting coefficient of each parameter and (f) estimation of the landslide instability index. In Tavoularis et al. [54], ten parameters were used as independent agents (e.g., thematic layers) to model the landslide susceptibility, implementing the RES approach. These parameters were: (i) distance from roads, (ii) land use, (iii) slope inclination, (iv) slope orientation (aspect), (v) lithology, (vi) hydrogeological conditions, (vii) rainfall, (viii) elevation, (ix) distance from streams and (x) distance from tectonic elements.

The landslide susceptibility map was generated in ArcGIS 10.2.2, using the abovementioned thematic layers, based on geodata provided by the Hellenic Survey of Geological and Mineral Exploration and the Digital Elevation Model from Hellenic Cadastre S.A. [54].

Afterwards, every data layer was digitized and converted into grids (cell size: 20 × 20 m). Furthermore, taking into consideration the RES methodology, weights and rank values were assigned to the reclassified raster layers and to the classes of each layer. Finally, the weights of the raster thematic maps were multiplied by the corresponding weights and added up (through “weighted sum”, which is a tool of the GIS environment) to generate a landslide susceptibility map, wherein each cell has a particular index value. After the reclassification of the above map, the final susceptibility map of Attica Region was generated and divided into five susceptibility zones according to Brabb et al.’s classification [76]: “Low–Middle”, with Instability index (Ii) < 25, “High” with 25 < Ii < 42, “Very High” with 42 < Ii < 53, “Extremely High” with 53 < Ii < 70”, and “Landslide” with Ii > 70%. From this classification, the main conclusion is that the higher the LSI is, the more susceptible the area is to landslides (instability index higher than 70%). The overall methodology is depicted in Figure 9.

The generated landslide susceptibility maps can be combined with the zoning maps of potential high-flood-hazard areas already produced for the Attica Region, authorized by the Greek Ministry of Environment and Energy, and with the flooded area maps delivered by the Copernicus Emergency Management Service Mapping [63]. The outcome of the DIAS project is accessible to the public, through open-access web-based platforms (https://gis.attica.gov.gr/en/node/1216 (accessed on 1 May 2023), https://maps.nccs.nasa.gov/arcgis/apps/experiencebuilder/experience/?id=29bd25e78fff45f0a6dbfd0328b4d03e (accessed on 6 October 2023)), to aid awareness of landslides among different stakeholders (e.g., landslide experts, government authorities, planners, decision-makers, citizens). Moreover, the DIAS project can facilitate the work of Civil Protection Authorities by providing inputs for prevention and preparedness. For further details on the methodology briefly described here, readers are kindly encouraged to read the following relevant research [54,73,74,75].

3. Results

The landslide susceptibility map generated for the entire county of the Attica Region is shown in Figure 10. The map has been validated using a confusion matrix and 220 slope failures [54].

Mapping Performance Evaluation

The aforementioned soil erosion types were subjected to expert-based cross-checks and validated through on-screen visual interpretation of the generated landslide susceptibility map of the Attica Region. The soil erosion types are depicted as lines and have been overlaid on the landslide susceptibility map (on landslide susceptibility zones) to identify the correlation of potential slope failure zones with soil erosion lines (Figure 11). It was found that there is a consistent coincidence between all of the derived soil erosion lines and the areas that are characterized, regarding landslide susceptibility, as “very high susceptibility”, “extremely high susceptibility” and “landslides”, based on the generated landslide susceptibility map of the Attica Region.

The verification was executed by implementing a frequency ratio statistical analysis, where the relationship between the spatial distribution of landslide susceptibility zones and soil erosion lines was studied. Specifically, the ratio is that of the area where soil erosion lines manifested to a particular landslide susceptibility zone. According to Pradhan et al. [39], a value of 1 is an average value. A value lower than 1 means a lower correlation of soil erosion lines, whereas a value greater than 1 means a higher correlation.

All the derived soil erosion lines were found to be in the “very high susceptibility”, “extremely high susceptibility” and “landslide” zones. To estimate the frequency ratio, the frequency distribution of soil erosion lines was calculated for each landslide susceptibility zone of the Mandra area. Moreover, the area ratio for each landslide susceptibility zone (measured in pixels) as well as the meter ratio for each soil erosion line were computed. Finally, the frequency ratios for soil erosion lines associated with the “very high susceptibility”, “extremely high susceptibility” and “landslide” zones were calculated by dividing the frequency of soil erosion lines by the landslide susceptibility zone area (measured in pixels).

The frequency ratio generated for the Mandra area is shown in

Table 2. Frequency ratio values of soil erosion lines to landslide susceptibility zones.

Landslide Susceptibility Zones	Pixel in Domain	Pixels (%) (a)	Soil Erosion Lines (m)	Pixels (%) (b)	Frequency Ratio (b/a)
42–53%	777.203	42.32	917.35	2.3	0.05
53.01–70%	933.057	50.81	21.920	55.11	1.08
70.01–100%	126.032	6.87	16.935	42.59	6.2
Total	1836.292	100.00	39.772	100.00	1

. For a lower probability of landslide occurrences (e.g., 42 < Ii (instability index) < 53), the frequency ratio is equal to 0.05, which indicates a poor relationship between the soil erosion lines and the generated landslide susceptibility zones. For extremely high landslide susceptibility (e.g., 53 < Ii < 70) and beyond (e.g., 70 < Ii < 100), the frequency ratio is found to be greater than 1, which indicates strong relationship between landslide susceptibility zones and soil erosion lines.

Taking into consideration the abovementioned factors, this kind of landslide susceptibility model (e.g., DIAS project) could be used as a soil erosion prediction model. Since soil erodibility reflects the soil’s susceptibility to erosion, the accurate mapping of susceptibility to erosion hazards is crucial to avoid economic losses and life losses. Thus, implementing the outcomes of landslide susceptibility mapping in validating soil erosion areas can help in identifying and consequently estimating soil erosion hazards and vulnerability.

4. Discussion

Although Kue et al. [51] reported that quantitative characterizations of the interaction between soil erosion and landslides are rare, in the present study, an attempt was made to quantitatively characterize the interaction between soil erosion and landslides using the RES methodology. Kue et al. [51] calculated the area and volume of 5420 shallow loess landslides and compared these results against the Chinese Soil Loss Equation (CSLE)-derived soil erosion rates of 15 subcatchments. Their analysis revealed satisfactory linear fitting between the areas of landslides and soil erosion. In this study, based on a heuristic approach, RES can be used to combine different parameters and study the interactions of those parameters to quantitatively estimate the instability index of an examined slope that can be associated with landslide susceptibility zones. In our case study, some of those zones coincided with soil erosion lines. These lines were derived not only from site investigations (in three different chronological periods: 2017, 2018, 2022), but also from a particular methodology (e.g., soil erosion susceptibility, which is described in detail in the manuscript).

In the studies of Lee et al. [38] and Pradhan et al. [39], the soil erosion map was verified using landslide locations via the integration of USLE with GIS to model the potential for soil erosion. In Rozos et al. [8], the Revised Universal Soil Loss Equation (RUSLE) was implemented to predict sites susceptible to slope failures caused by soil erosion. In those three papers, the verification that followed was the same as in the present study, which means that the generated soil erosion map was verified by comparing the soil erosion hazard zones with the spatial distribution of slope failures. Furthermore, one additional common characteristic is that the assessment of the dynamic soil erosion process can be correlated with another equally important and related threat of landslides [47]. Finally, from the abovementioned studies, considering the present study, it can clearly be said that, by adjusting the factors that are responsible for soil erosion such as lithology and land use, the rate of soil erosion can be minimized. As a result, gully and landslide inventory maps are very useful in identifying the distribution of sediment sources and the landscape characteristics that are connected to their existence [47]. By defining that distribution, this study could contribute to monetary estimates regarding the cost of removing sediments, while trying to implement flood mitigation measures for new potential heavy rainfall episodes.

In Tavoularis et al. [54], different types of slope failures [77] were found to manifest in the entire region of Attica, such as falls, rockfalls, slides, and debris. Among them, regarding the Mandra area, shallow landslides are the most prominent type that affects that region. Shallow landslides are very much related to gradual soil erosion, since they easily affect soil materials transposed by erosion [78]. Thus, the interaction of soil erosion and slope failure processes contributes to the loss of fertility, and the alteration of the landscape, not to mention damage to infrastructure and human facilities. Moreover, soil erosion and shallow slope failures are important geological issues, so specific studies are crucial to generate models to cope with hazards and vulnerability to these phenomena [79]. Taking into consideration the abovementioned factors, a well-structured landslide susceptibility map, such as that previously described, may be helpful in identifying and characterizing areas (in this study on a regional scale) of potentially increased land degradation. Therefore, it is important to gather and study spatial information on soil loss interacting with slope failures.

4.1. Prevention and Control Actions

It is scientifically documented, legally established and empirically proven that, to be effective in dealing with the flooding action of a stream, the study of flood control works should focus on the entire catchment area and include the necessary forestry works of stream regulation, which have the following positive results [61,62]:

• A reduction in the amount of solid material transported, with a corresponding reduction in the erosive capacity of the flood waters and the volume of the flood wave;
• The velocity of the flood wave is reduced, resulting in a delay in its occurrence downstream and a reduction in its destructive momentum;
• The effects of erosion on unprotected soils are reduced;
• The natural environment is protected and enhanced, especially through planting and soil protection projects;

The types of projects proposed include the following:

• Construction of small dams to grade the bed and retain the slopes;
• Construction of dams for the retention of debris;
• Construction of culverts in places where the existing road network is eroded by streams in the study area;
• Settlement of part of the hydrographic network of the study area by constructing an artificial bed with a dike;
• Implementation of horticultural works;
• Forestry measures for the management of the overall forest complex in the study area;
• Opening of forest roads to reach the sites of the proposed projects.

Particularly, given that the purpose of the interventions is to reduce flood risk and prevent catastrophic events, it was considered appropriate to limit the proposed measures to small dams, which are mainly of a preventive nature and aim to manage runoff effectively. These measures were considered to be the most appropriate proposal for the area, at least because they meet the requirements of being as cheap as possible and having a very short implementation time. To be more specific, the types of recommended dams are as follows (Figure 12):

(a)

The sediment barriers are intended to be between 3 m and 8 m high, with a reinforced concrete or unreinforced concrete construction material. The purpose of these dams is to counteract the axial erosion of the bottom of the streambed, by reducing the drag force of the water and retaining the sediment;

(b)

Construction of graduation dams–slope stabilization—these dams are proposed to be made of concrete (reinforced or unreinforced) or of reinforced wire mesh (e.g., sarsenet). The construction of the dams will be carried out either in places where there is evidence of gradual erosion, or in places where there is axial erosion of the bottom of the bed, in combination with the abovementioned sediment barriers. The height of these dams according to the theoretical assessment carried out is proposed to be between 1 m and 2 m.

4.2. Research Implications, Limitations, and Future Directions

To the author’s knowledge, the relationship between soil erosion and landslide events has not yet been approached using an expert semiquantitative methodology, named the Rock Engineering System (RES). Apart from the already known models that have been implemented before today (e.g., USLE, RUSLE, etc.), and which models have limitations (as was mentioned in the Section 1 (Introduction)), soil erosion identification can also be approached via a tool such as that of RES in conjunction with landslide susceptibility mapping, which, in our case study, has been generated for a region of approximately 3800 km². This map has already been validated by 220 slope failures recorded in the last 60 years (1960–2020) in the entire region of Attica. The validation via too many slope failure incidents, along with the possibility that different experts can rate the selected parameters and code their interactions through their expertise and experience to reduce subjectivity, make the RES approach a simple and reliable tool in investigating quickly and in a scientifically appropriate manner the susceptibility of land degradation. In this direction, the implementation of geographical information systems technology enables the continuous monitoring and evaluation of soil erosion susceptibility to hazard reduction.

In addition, this study could assist the local authorities of the Mandra municipality in defining areas susceptible to soil erosion, and as a result it can contribute to the development of resilience for future flooding protection, and minimize further damage or prevent the occurrence of phenomena such as slope failures.

This study could also contribute to some quantified monetary estimates about the costs of removing sediments during the process of implementing mitigation measures against upcoming potential flood episodes.

In conclusion, the implementation of the landslide susceptibility model in this study can contribute to the online repository of scientific information in the EU Soil Observatory of the European Soil Data Centre (e.g., datasets, maps).

In this direction, the implementation of artificial intelligence and machine learning methodologies using (free) open-access Web-GIS platforms (such as those of the DIAS project), accompanied by a variety of geo (including soil) data, could lead to the further validation of the European Landslide Susceptibility Map (ELSUS), offering more accurate regional susceptibility maps through the evaluation of information downloaded from the European Soil Data Centre (ESDAC), further permitting the identification, correlation and quantification between land degradation and soil erosion.

5. Conclusions

This study investigated the correlation between soil erosion and land degradation in the Mandra area (Attica, Greece). In particular, soil erosion types were depicted on a soil erodibility map over an area that suffered from a disastrous flash flood (November 2017). The reliability of this map was assessed by comparing the slope failure predictions generated in a research project (DIAS) and validated via several field works focused on designing flood mitigation measures under the auspices of the Directorate of Technical Works of Attica Region. The main conclusions of this work can be summarized as follows: (a) There is a very strong correlation between soil erosion and potential landslide susceptibility zones. (b) The landslide susceptibility model (DIAS) could be used as a (preliminary) guide for investigating and identifying soil erosion issues. Based on the abovementioned, we can offer an alternative way to identify soil erosion. Thus, implementing the RES approach could be of great use to different stakeholders in designing appropriate mitigation measures against floods and landslides.