**1. Introduction**

Landslide hazard assessment requires a multi-hazard approach, since the types of landslides that will occur usually have different characteristics with different spatial, temporal, and causal factors [1]. The first step towards the evaluation of landslide hazards on a regional scale (e.g., 1:25,000–1:250,000) is the assessment of the relevant susceptibility, which is defined as the likelihood of a landslide occurring in an area in relation to the local geomorphological conditions [2]. In addition, the landslide susceptibility map can be used as an end product in itself [1]. In order to develop a susceptibility map, it is mandatory to first compile an inventory map where the spatial distribution of existing slope failures is shown. It should additionally be pointed out that on a regional scale map is not feasible to discriminate in detail the type of landslide and delineate the runout per failure.

Having developed the landslide inventory map, the likelihood of slope failures i.e., susceptibility, can be assessed by both qualitative and quantitative methods. The former group of methods includes the knowledge-driven methods (direct and indirect mapping), and the latter group includes the data-driven and the physically-based ones [1]. Considering the regional and local scale maps, the knowledge and data-driven approaches are suggested to be applied; for the former approach a geoscientist i.e., geomorphologist, can directly determine the level of susceptibility based on his/her experience and information

**Citation:** Tavoularis, N.; Papathanassiou, G.; Ganas, A.; Argyrakis, P. Development of the Landslide Susceptibility Map of Attica Region, Greece, Based on the Method of Rock Engineering System.*Land* **2021**, *10*, 148. https://doi.org/ 10.3390/land10020148

Academic Editors: Enrico Miccadei, Giorgio Paglia and Cristiano Carabella Received: 29 December 2020 Accepted: 29 January 2021 Published: 3 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

related to terrain conditions, while the data-driven mapping statistical models are used in order to forecast likely to landslide areas, based on information obtained from the interrelation between the spatial distribution of landslide conditioning factors and the landslide zones [3]. The most widely applied data-driven approaches are [1]: bivariate statistical analysis, multivariate statistical models and data integration methods like Artificial Neural Network analysis. Bivariate statistical methods (e.g., fuzzy logic, Bayesian combination rules, weights of evidence modelling) are considered as an important tool that can be used in order to analyze which factors play a significant role in slope failure, without taking into account the interdependence of parameters. Multivariate statistical models evaluate the combined relationship between the slope failure and a series of landslide controlling factors. In this type of analysis, all relevant landslide parameters are sampled either on a grid basis or in a slope unit and the presence or absence of landslides is evaluated. These techniques have become standard in regional-scale landslide susceptibility assessment.

Nowadays, the majority of the studies considering landslide susceptibility mapping makes use of digital tools for handling spatial data such as Geographical Information Systems (GIS). Specifically, the GIS-based techniques are considered very suitable for the landslide susceptibility mapping, in which the predisposing factors (e.g., geology, topography) are entered into the GIS environment and combined with the spatial distribution of slope failures i.e., landslide inventory map [3–6]. For the purposes of this study, the semi-quantitative methodology of Rock Engineering System (RES) originally introduced by Hudson [7] was implemented in Greece, particularly in the Attica region for the assessment of landslide susceptibility. This region, which is a county with a size of approximately 3800 km2, was selected due to the following reasons:


Region. Furthermore, the generated landslide susceptibility map will serve for many authorities related to public works, as a dynamic map for the planning, design, and implementation of a long-term landslide reduction strategy as well as identifying the areas where more detailed investigations will be required for the planning of critical infrastructure.

(vi) taking into account that the next five to ten years, very important civil engineering projects are about to be constructed in Attica county (such as transports network elements: highways, railroads, metro-tunnels, hospitals, administrative buildings, security/emergency structures, residential buildings) the existence of a regional-scale landslide susceptibility map could be a very useful tool for supporting decisions in order to prevent the location of high-value constructions in unsuitable locations.

**Figure 1.** Simplified geological map of Attica region, based on the official Greek projection system (EGSA 87). Active faults were inserted in this map from the National Observatory of Athens (NOAFAULTs, https://zenodo.org/record/4304613# .YAmJbugza1Z).

Regarding the above-mentioned, the scope of this study that is part of the project "Landslide Risk Assessment of Attica Region (DIAS)", is twofold: (i) construct a uniform and updated geodatabase of slope failures induced the last sixty years in the whole territory of Attica Region, and (ii) compile a landslide susceptibility map, being the basic step to produce the upcoming landslide hazard and risk maps.

#### **2. Geology and Tectonic Setting of Attica**

Attica is located in the back-arc area of the Hellenic Arc. The geology of Attica comprises Alpine basement rocks, both metamorphic and non-metamorphic, and post-Alpine sediments (Figure 1). The Alpine rocks belong to the high-pressure metamorphic units of the Cyclades and Almyropotamos that extend from Penteli Mt, east Attica [10] to the southern Gulf of Evia and to the non-metamorphic units of Eastern Greece/Sub-Pelagonian units that outcrop in Parnitha Mt and in west Attica. The southern parts of Attica are also underlain by schists and marbles of the Cycladic Metamorphic Belt. An 8.2 Ma granodiorite outcrops in the Lavrion area of SE Attica. The post-alpine (syn-rift) formations consist of alternating beds of marls, lacustrine limestone marls and sandstones. Quaternary deposits are talus cones, sandy–clayey soils, scree, and unconsolidated clays. [11].

Rifting started in Middle-Upper Miocene and continues until the present day resulting in the formation of several basins. According to Freyberg [12], in the western part of the Athens Basin, the Pliocene formation (with a considerable thickness reaching locally more than 300–400 m) can be found, such as clays, sands and sandstones, and gravels in alternation with white limestone. The dating of the synrift ranges from Upper Miocene to Holocene times. There are also Quaternary volcanic formations consisting of loose volcanic extrusive rocks with tuff blocks, dacitic and andesite domes as well as alluvial fan deposits and steep talus cones covering parts of Aegina island, Poros island and almost the entire Methana peninsula.

The Athens basin is the main neotectonic feature in Attica, elongated in a NE-SW direction. An important tectonic structure is the NNE-SSW, west-dipping detachment fault that separates the metamorphic units to the east from the un-metamorphic units to the west [13,14]. The fault was active in Late Miocene-Early Pliocene and produced several hundred meters of debris-flow deposits. In addition, the active normal faults of Avlon-Malakasa, Afidnes, Milesi, Pendeli, Kakia Skala, Thriassion and Fili dominate the area [15,16]. These faults present characteristic features such as prominent scarp linearity, considerable scarp height, unweathered scarp appearance and fault-slip kinematics that are compatible with the regional stress–strain fields (N-S to NNE-SSW) [17,18].

Based on their morphotectonic features [16], all normal NW-SE trending major faults of Attica could be considered "active structures". Overall, the northern part of Attica is bounded by a series of north-dipping active fault segments, while the central part by south-dipping active faults, respectively [16,19–22]. The slip rates of active faults are less than 1 mm/year [15,21,22] and average earthquake recurrence intervals are expected in the order of a few thousands years.

An interesting part of the geological setting of the Attica region is Kithira and Antikithira islands, which are the southeastern islands of the Ionian Sea between Peloponnese and Crete and belong to the administration of Attica Regional Authority. The geological formations that are found there, comprise metamorphic rocks as well as carbonate rocks of Tripolis and Pindos geotectonic zone. Both islands are surrounded by N-S oriented active faults due to ongoing east-west extension in this area of the Hellenic Arc.

#### **3. Materials and Methods**

#### *3.1. Landslide Inventory of Attica*

The first step towards the compilation of a landslide susceptibility map is the development of a landslide inventory [23]. In this study, the generated inventory map, and the landslide geodatabase, cover a chronological period from 1961 up to the present.

The methods that were used for the generation of the inventory are classified into the following approaches:


**Figure 2.** *Cont*.

**Figure 2.** (**a**) Rockfalls at the coastal areas of Alepochori–Psatha (North-Western Attica), (**b**) earth fall on bank slopes subjected to undercutting by Chelidonous stream (North of Athens, Kifisia municipality). (**c**) Rockfalls occurred in Attica islands such as Spetses (e.g., Agriopetra) and (**d**) Kithira (e.g., Galani spring-Agia Pelagia). In Figure 2c, the blue circles around rocks emphasize the grea<sup>t</sup> possibility for rockfalls. (**e**) Complex slope failure in Salamina island (Porto Fino site), (**f**) A rock topple failure in Hydra (adjacent to Miaoulis statue), (**g**,**h**) An earth slide from Penteli area (Ntrafi site) at northeastern of Athens. The toponyms of each characteristic site are depicted in Figure 1.

> Following the terminology defined by the Working Party on World Landslide Inventory (1990) [25], the majority of the depicted slope failure sites hold information on location, dimensions-geometry, landslide-movement type, trigger mechanism, damage caused, slope and aspect, lithological composition, movement date, older activation, seismic risk zone, meteorological data, hydrogeological behavior, consequences, proposed remedial measures, the confidence of landslide identification, mass movement date–field survey date, bibliographic reference and characteristic photos for each slope failure. The developed landslide inventory map is shown in Figure 3, where slope failures are interpreted as points (220 sites), polygons (98 areas delineated based on the Oregon Protocol) and erosion lines based on data provided by the Hellenic Survey of Geology and Mineral Exploration (H.A.G.M.E.), assigning a unique identifier and a number of attributes to each landslide. Taking into account Varnes classification (1978) [26], the movement type of the 220 slope failures, shown as points, can be characterized as follows (Table 1):

**Figure 3.** Frequency of landslides of Attica for the period 1961–2020. Bin-size ten years. 34


**Table 1.** Movement type of the 220 inventoried slope failures in Attica Region (based on Varnes nomenclature). Each type is associated with a specific recorded number of failures. Each slope failure is depicted in Figure 4.

> The geodata within the DIAS database followed the EU Inspire Directive and is maintained in a digital format that can be adapted and updated for future use. Furthermore, from the DIAS geodatabase, some more extra remarks can be deduced about the frequency of slope failures per decade from 1961–2020 (Figure 3). It is noted that the number of recorded slope failures increased in the 2000–2010 and 2010–2020 decades in comparison to the pre-2000 data, and this can be explained due to intensive climate change and due to the execution of more detailed field and remote sensing surveys from public authorities, research institutes and consulting agencies.

> In the following Figure 4, the developed landslide inventory map is shown. In the legend of the map, the slope failures depicted with green circles correspond to landslides that have already manifested at Attica Region in the past. Slope failures in red polygon shapes are those that are delineated through the methodology described by the protocol of Special Paper 42, developed by the Oregon Department of Geology and Mineral Industries. Finally, erosion lines were provided by the Hellenic Survey of Geology and Mineral Exploration (H.S.G.M.E.) through a research project which proposed flooding mitigation measures in the Mandra area, west of Athens.

**Figure 4.** The developed by this study landslide inventory map of Attica Region for the period 1961–2020. Background hillshade image is derived from a high-resolution Digital Elevation Model.

#### *3.2. Assessment of Landslide Susceptibility at the Attica Region*

For the purposes of this study, a semi-quantitative heuristic methodology called Rock Engineering Systems (RES), originally introduced by J. Hudson (1992) [7], was applied to assess the landslide susceptibility. The Rock Engineering System approach has been used for a wide variety of rock engineering and other topics, such as surface blasting, natural slope instability, earthquake and rainfall-induced natural slope instability, road-cut induced slope instability, rockfall assessment, engineering geology zonation, coastal landslides, TBM performance, Metro tunnel stability, and many more applications in engineering modelling and design [27]. Furthermore, regarding recent findings of the implementation of RES generally in geotechnical engineering applications, it can be mentioned that:


In Greece, the RES methodology has been applied in different geological settings and scales. For example, Rozos et al. (2006) [32] have used RES for a study in Karditsa prefecture, Greece (scaled in 1:50,000), Rozos et al. (2011) [33] have compared RES and Analytical Hierarchy Process (AHP), Tavoularis et al. (2017) [34] tested RES on Malakasa (1995) and Tsakona (2003), Greece in site-specific scale (1:1,000 to 1:5,000), Tavoularis (2017) [35] implemented RES in a regional scale area (Geological Sheet of Megalopolis, Greece scaled in 1:50,000) in complex geological setting and tectonic regime environment.

In this study, an attempt is made to implement RES in a larger coverage area (scaled in 1:100,000) than those previously mentioned with many different geological settings (active faults, places adjacent to dormant volcanic eruptions, streams banks eroded by flash floods), densely populated and surrounded by many important infrastructure facilities.

#### 3.2.1. The RES Approach

A crucial problem of any engineering design is ensuring that all the necessary parameters are included and that the interactions among them are understood. John Hudson was the researcher that originally introduced the Rock Engineering Systems (RES) approach in 1992. The RES methodology is a synthetic approach which studies the problem (e.g., landslide), breaks it down into its constituent variables (e.g., predisposing parameters, estimation of landslide instability index), and assesses their significance (e.g., calculation of susceptibility analysis). In most slopes, that kind of analysis is complicated due to different interacting factors, complexity of geological formations, different scale of the instability events as well as a scarcity of detailed geodata. These problems can be solved through the use of RES, where its use can take into account the particular problems at any investigated site so as to identify critical sites in order to support decisions on land use and planning development [27].

For consideration of a specific engineering project–system (in our research the landslide susceptibility of the Attica region), some parameters are expected to show a greater effect on the project–system than others and some parameters will in their turn be significantly affected by the system. The RES methodology uses a table (i.e., interaction matrix) with xi rows and yj columns, in which the selected n parameters are selected as leading diagonal terms and the interactions between them are considered as off-diagonal terms. In Figure 5, the row passing through the parameter Pi represents the influence of Pi on all the other parameters in the system, whereas the column through Pi represents the influence of the other parameters on Pi. Afterward, we study this so-called influence by coding the off-diagonal components in order to express their importance. A semi-quantitative coding

method was used with values ranging from 0 to 4 corresponding to: 0-No interaction (most stable conditions); 1—Weak interaction; 2—Medium interaction; 3—Strong interaction and 4—Critical interaction (most favorable condition for slope failure), respectively. For eliminating the subjectivity, this coding method can be used by one or more experts familiar with the project being considered [7]. Next, the sum of each row (named as "cause-C") and each column (named as "effect-E") can be determined and designated as co-ordinates (C, E) in the diagram of Figure 5. The meaning behind this diagram is that C represents the way in which Pi affects the system; and E represents the effect that the system has on Pi, by indicating a parameter's interaction intensity (as the distance along the diagonal) and dominance (the perpendicular distance from this diagonal to the parameter point). By these two words, we quantify parameter significance inside the matrix system (i.e., landslide).

**Figure 5.** Interaction matrix. The dashed lines correspond to the terms "interaction intensity" and "dominance" respectively [7].

According to Hudson (1992) [7], there are many "constellations" that could occur, the two main ones being mainly along theC=E line or mainly along a line perpendicular to it. If the parameter points are scattered along the C = E line but close to it, then they can be ranked according to their parameter interaction intensity; in other words, they can be listed in order of interactive importance (Figure 6a). If, on the other hand, the parameter points are scattered on a line perpendicular to the C = E line, they will have similar interaction intensities but widely differing dominance values (Figure 6b). In the former case, it might be possible to use five or six parameters in such a scheme; in the latter case, all the parameters must be used.

The cause versus effect diagram reveals the influential role of each parameter on slope failure which is expressed by the term "weighted of coefficient influence". Respectively, the role of the system's interactivity is expressed from the histogram of the interactive intensity [cause (C) + effect (E)] against the parameters. This intensity is transformed into weighting coefficients, which express the proportional share of each factor in slope failure and normalized by dividing with the maximum rating (4), giving the ai%, as it is explained in the next paragraph.

**Figure 6.** Interaction intensity–dominance diagram, with different forms (**<sup>a</sup>**,**b**) of parameters constellations [34].

The next step is to compute the instability index (Ii) for each examined slope, by using the following equation:

Ii

$$\mathbf{k} = \boldsymbol{\Sigma}\mathbf{a}\mathbf{i} \times \mathbf{P}\mathbf{i}\mathbf{j} \tag{1}$$

where Pij is the rating value assigned to the different category of each parameter's separation, i refers to parameters (from 1 to 10 corresponding to this research, and generally from 1 to n, in other case studies where a different number of landslide parameters are selected), j refers to the examined slope and ai is the weighting coefficient of each parameter provided by the formula:

$$\text{Azi} = 1/4 \,\text{\*}\,\left[ (\text{C} + \text{E})/(\Sigma \text{iC} + \Sigma \text{iE}) \right] \% \tag{2}$$

normalized to the maximum rating of 4. It should be noted that the instability index is an expression of the potential instability of the slope, with values ranging between 0 (no slope failure at all) and 100 which refers to the most unfavorable conditions (i.e., landslide).

#### 3.2.2. Selection of the Parameters Controlling the Slope Failures

Ayalew and Yamagishi (2005) have reported that there are five basic concepts for the chosen parameters regarding the assessment of landslide susceptibility [36]. Parameters should: (i) vary spatially, (ii) be measurable, (iii) be related to the presence or absence of landslides, (iv) be representative of the entire study area, and (v) not account for double consequences in the final outcome. Ten parameters were selected as independent controlling factors for the landslide manifestation of the Attica Region, and classified into five classes. These factors which were utilized for the RES methodology are the (i) distance from roads, (ii) slope inclination, (iii) slope orientation (aspect), (iv) lithology, (v) hydrogeological conditions, (vi) rainfall, (vii) land use, (viii) distance from streams, (ix) distance from tectonic elements and (x) elevation.

In order to decide and consequently select the above-mentioned parameters, (a) we studied a huge amount of published and unpublished engineering geology reports, (b) we applied very interesting landslide research based on statistical analysis gained in Greek territory [37], and (c) we took into consideration the field observations conducted in Attica region in the frame of this study [38]. In the following paragraphs, the importance of each selected parameter for the initiation of a landslide and an analysis of what do we mean by the terms "dependence" and "independent" are provided.

The meanings of "dependance" variable and "independent" parameter are related to the role each one has inside the whole system we study. The system can be a slope failure or underground stability and support or the selection of the right type of tunnel boring machine or any other geotechnical engineering problem that can be addressed by using this semi-quantitative heuristic methodology of RES.

To be more specific, referring to "dependance" variable is meant the occurrence or not of a slope failure. For example, we study the interaction of ten landslide parameters and according to RES methodology, we calculate the weighted coefficient of each landslide parameter, estimating the instability index for each examined slope. If the calculated index is over a critical accepted threshold (as the one that we will present in the following section of this paper), this means that the selected parameters are crucial for the slope failure occurrence, and subsequently, measures must be taken in order to minimize their effect on slope instability. Otherwise, if the estimated index is under this critical threshold, then, no landslide is about to happen and immediately we conclude that those parameters that we selected are not crucial for landslide initiation.

On the other side, as "independent variables" are characterized the landslide controlling factors (such as geology, distance from roads, hydrogeological conditions, distance from tectonic elements), and each other is tested on how dominant or how interactive can be with the other selected landslide parameters.

RES studies the interaction of each parameter to the other and vice versa, by quantifying the different importance of these interactions. This is justified because some parameters will have a greater effect on the system (e.g., in our case the landslide susceptibility in Attica county,) than others and some parameters will in their turn be largely affected by the system. Thus, talking for example about the interaction of hydrogeological condition on lithology, it is meant how lithology can be affected by the permeability status that dictates the geological formation that constitutes the examined slope and vice versa how a specific type of rock or soil of the examined slope will affect the hydrogeological equilibrium of the slope. In another case, we examine how the distance from a road affects the amount of vegetation that exists around this. To be more specific, if a public authority plans to construct a new highway in a place where forest or a grassland area already exists in that particular zone, then it is proved that buffer zones of highway that are in a distance 50 or 200 m from the surrounded slopes affect the existence of vegetation dramatically [33]. Vice versa, the influence of vegetation on slopes that are in an x distance from roads is less important.

In the following paragraphs, a brief comment on the importance of each predisposing landslide parameter is presented.

#### (i) Distance from roads

During the construction of the road network, vegetation removal, and the application of external loads as well as extensive excavation are some of the most common human intervention actions which are taking place, and result in landslide triggering [39]. It should be mentioned that the digital record of Attica county roads for the generation of DIAS geodatabase was provided by the General Secretary of Civil Protection Agency of Greece. Buffer zones were created around the roads. According to many studies but mostly based on Rozos et al. [33], the slopes that are at a distance of 50 m from a road are more prone to failure.
