European Geologist Journal 44

Landslide susceptibility mapping using the Rock Engineering System approach and GIS technique: an example from southwest Arcadia (Greece)

by Nikolaos Tavoularis1, Ioannis Koumantakis1, Dimitrios Rozos1and Georgios Koukis2

 Department of Geological Sciences, School of Mining and Metallurgical Engineering, NTUA, 9 Heroon Polytechniou St., 157 80 Zografou Athens, Greece

2  Department of Geology, University of Patras, 265 04 Patras, Greece



The purpose of this study is to prepare a susceptibility map in a landslide-prone area in Greece using Rock Engineering System (RES) and a geoprocessing tool called Model Builder. The implementation of RES is achieved through an interaction matrix, where ten parameters were selected as controlling factors for the landslide occurence. The validation of the developed model was achieved by using field-verified data, showing excellent correlation between the expected and existing landslide susceptibility level. In conjunction with Model Builder, which can overlay different layers and produce landslide susceptibility maps, RES can act as a tool for calculating the instability index and getting a prognosis of a potential slope failure in relation to sustainable development planning processes in landslide susceptible areas.

1. Introduction

Over the last two decades, many studies of landslide susceptibility assessment have been made. It is believed that the accuracy of landslide susceptibility mapping increases when all determining parameters are included in the analytical process. Rock Engineering System (RES), which is a semi-quantitative rock engineering approach and the basic tool for representing the parameters and their interaction mechanisms, can thus be useful in decision making regarding land use and development planning processes in landslide susceptible areas by providing a tool for zoning landslide hazard. It is based on an interaction matrix which represents the key parameters as leading diagonal terms and their binary interaction mechanisms as off-diagonal terms (Figure 1).

Figure 1: Summation of coding values in the row and column through each parameter to establish the cause and effect co-ordinates (after Hudson, 1992).

RES was developed by Hudson (1992) to determine the interaction of a number of parameters in rock engineering design and calculate an instability index for rock slopes. In this study, an attempt is made to prove that RES can be implemented with the same efficiency in forecasting landslides, which are related to a variety of geomaterials [for instance soils, soft rocks, flysch formation (intercalation of different geological formations, etc.)].

Here, the RES approach was used for landslide susceptibility mapping in the wider area of Megalopolis, located in southwest Arcadia, which is part of the administrative region of Peloponnese, Greece. In the first stage, 21 landslide locations were identified in the study area from the literature field surveys and interpretation of satellite pictures. Secondly, ten data parameters (layers) were used as landslide conditioning and triggering factors for susceptibility mapping.

Next, the examined area was analysed using the RES methodology in conjunction with Geographical Information Systems (GIS), which facilitated the manipulation of these ten selected landslide data layers. To be more specific, a geoprocessing GIS function called “Model Builder” from ArcGIS (ESRI) was applied, providing automatic preparation of the landslide susceptibility map. The results of the RES analyses proved the instability of the field-verified landslide locations. Moreover, the verification results showed not only excellent correlation between the susceptibility map produced and the existing data of the 21 historical landslide locations but also indicated that many more potential slope failures could be taken place in the wider region of Megalopolis in the future. In conclusion, this contribution (the combination of RES with Model Builder resulting in landslide susceptibility mapping) provides originality to this study.

2. Materials and methods

2.1 Establishing the interaction matrix and matrix coding

The RES approach, which is based on an expert’s judgement, uses an interaction matrix in which the main parameters thought to govern a particular circumstance (e.g. slope failure) are selected and the interactions between them are considered (Hudson, 1992). This enables a comprehensive assessment of the factors and interactions, the advantage being that all potential influencing factors can be included initially. RES methodology reduces uncertainty because study of the interactions between the factors indicates the degree of influence of the factors in the system being considered – which are dominant and which have a much lesser or insignificant contribution – thus reducing the uncertainty.

The interaction matrix (Table 1) shows in its main diagonal cells the principal parameters considered responsible for controlling the potential instability of the examined slopes, while its off-diagonal cells contain the coded expressions of all possible binary interactions. For the purpose of the present work, a range of possible interactivity from 0 to 4 was adopted (Koukis and Ziourkas, 1991), where ‘none’ is coded 0 to indicate the most stable conditions, and other interactions are ranked ‘weak’ (coded 1), ‘medium’ (coded 2), ‘strong’ (coded 3) or ‘critical’ (coded 4 – the most favourable condition for slope failure).

By coding the interaction matrix components and then summing the values in the row and column through each parameter, “cause” and “effect” co- ordinates are generated, indicating a parameter’s interaction intensity (Figure 2).

Table 1: Interaction matrix of the RES method.

Figure 2: Parameter interaction Intensity and Dominance (after Hudson, 1992).

The influential role of each parameter on slope failure (weighted coefficient influence) is revealed from a Cause [C] versus Effect [E] diagram (Figure 3), while the role of the system’s interactivity is expressed from the histogram of the interactive intensity (C+E) (Figure 4). The C+E values will be transformed into a percentage form acting as weighting coefficients, which express the proportional share of each parameter (as a failure-causing factor) in slope failure and are normalised by dividing by the maximum rating (4), giving the ai percentage.

Figure 3: Cause-Effect diagram for the 10 parameters of the Megalopolis study area.

Figure 4: Interactive intensity of parameters.

The next step is to compute the instability index (Ii) for the considered slope using the following equation:

Ii=Σai x Pij,

where i refers to parameters (from 1 to 10), j refers to the examined slope and ai is the weighting coefficient of each parameter given by the formula:

ai=1/4 * [(C+E)/(ΣiC+ ΣiE)]%,

scaled to the maximum rating of Pij (maximum value=4). Pij is the rating value assigned to the different category of each parameter’s separation, which also fits better to the conditions related to the parameter in question regarding the examined slope failure (Rozos et al., 2008). The instability index is an expression of the inherent potential instability of the slope, where the maximum value of the index is 100 and refers to the most unfavourable conditions.

2.2 Geological setting of the study area

The study area is located in Greece and specifically in the southwestern part of Arcadia in Peloponnese (Figure 5).

Figure 5: Modified geological map of Megalopolis area, scale 1:100,000 (Tavoularis, 2017).

In this particular region two alpine geotectonic units of the external Hellenides are present, namely (i) the Tripolis unit (shallow – water carbonates, Triassic – L. Eocene and flysch, L. Eocene – E. Miocene), and (ii) Pindos unit (pelagic limestones, radiolarites, the so-called “first flysch”, thin-bedded limestones, L. Cretaceous and flysch, Danian-Eocene). Pindos units overthrust Tripolis units, forming successive thrust sheets with movement direction from east to west. The neotectonic macrostructure of the broader area (SW Peloponnese) is characterised by the presence of large grabens and horsts bounded by wide fault zones, striking N-S and E-W (Ladas et al., 2007). In addition, this area was affected by the great Neogene phase of crustal extension.

The examined area is included in the 1:50,000 Geological Map of Megalopolis and has an extent of 614 km2. The main rivers are the Alfios and Elisson. The mean annual rainfall is around 1000 mm, while the maximum precipitation falls between November and March. Altitude values in the study area vary between 88 to 1340 m (Tavoularis, 2017).

2.3 Selection of the appropriate parameters controlling the slope failures

The majority of such studies follow five basic concepts for the chosen parameters, specified by Ayalew and Yamagishi (2005). Parameters should:

  • be representative of the entire study area,
  • vary spatially,
  • be measurable,
  • not account for double consequences in the final result,
  • have a certain degree of affinity with the dependent variable (the presence or absence of landslides).

In the above case study area, ten parameters were selected as independent controlling factors for the landslide occurrence and each factor was classified into 5 classes. The factors utilised for the RES methodology were:

P1 – distance from roads,

P2 – tectonic regime,

P3 – slope inclination,

P4 – slope orientation (aspect),

P5 – lithology,

P6 – hydrogeological conditions,

P7 – rainfall,

P8 – thickness of weathering mantle,

P9 – distance from rivers and

P10 – distance from tectonic elements.

The geodata were adjusted to the local conditions of Megalopolis area and rated for construction of the interaction matrix (Tavoularis, 2017). These parameters can be quantified more easily than those of time – and money – consuming ones (Table 2).

Table 2: The selected parameters and their rating.

2.4 Application of Model Builder in landslide susceptibility mapping

The susceptibility approach was designed by the USGS in the 1960s as a qualitative way to prepare landslide maps or to delineate zones affected by landslides, assessing the propensity of a given slope unit to generate a landslide based on spatial data (Brabb et al., 1972). Naturally, any landslide susceptibility prediction has a level of uncertainty. Sources of uncertainty include:

  • errors and incompleteness in the landslide and thematic information to complete the analysis,
  • an imperfect understanding of landslide processes and their geographical and temporal evolution;
  • limitations in the techniques used to determine the susceptibility,
  • the inherent natural variability of landslide phenomena.

Determining the errors associated with the geomorphological, geological and other thematic information is of primary importance. Improving the understanding of the landslide processes is feasible, but requires time and resources often not available to landslide investigators (Guzzetti et al., 2006). Consequently, some parameters will have to be rated from the beginning. To overcome these difficulties, in this study Model Builder, an application of ArcGIS (ESRI), is used for the automatic preparation of a landslide susceptibility map by modifying each parameter easily and quickly at any time.

Model Builder is a visual programming language for building geoprocessing workflows that allows multiple processes to be combined. The model is represented as a diagram that links together sequences of processes and geoprocessing tools, using the output of one process as the input to another process. It enables the user to visualise work flow (in the form of flow chart diagrams) and author and automate geoprocessing tasks that would normally be executed in single steps in ArcMap. It also has the resultant advantage of allowing the user to document the steps involved in the development of a model. While the development of the initial version of a model might take a little more time than conducting the steps manually, it is extremely useful when conducting multiple runs of a model – the model can be run on different data or small changes in the model can be made and the model rerun to examine model alternatives and assumptions (National Land Service of Lithuania, 2008). By using this application, it becomes easier to test the susceptibility model (Figure 6).

2.5 Data analysis

The following step was the production of the landslide susceptibility map through the construction of different thematic maps associated with landslide-related variables. The data used for the preparation of these layers were obtained from the Hellenic Military Geographical Service topographical sheets (scale 1:50,000) and IGME geological map (Megalopolis, scale 1:50,000). All data layers were digitised either from the original thematic maps or derived from spatial GIS calculations and finally were converted into grids with cell size 30 × 30 m.

The next step was to assign weights and rank values to the raster layers (representing factors) and to the classes of each layer respectively. This was realised with the use of RES. Finally, the weighted raster thematic maps with the assigned ranking values for their classes were multiplied by the corresponding weights and added up (through the ArcGIS tool of weighted sum) to yield a simple map where each cell has a certain landslide susceptibility index value. After reclassification this map represents the final susceptibility map of the study area (Figure 7).

Figure 7: Landslide susceptibility map of Megalopolis area scaled in 1:100,000 (Tavoularis, 2017).

3. Results and discussion

3.1 Implementation of RES in Geological map of Megalopolis (scale 1:50,000)

RES was implemented in the area defined by the geological map of Megalopolis, taking into account the interactions of the examined principal parameters and the calculation of their weighting coefficients. This resulted in the determination of an instability index for each examined slope of the study area. The RES matrix shown in Table 1 provides interactions of the chosen parameters based on the ratings outlined in Table 2.

For example, regarding the effect of rainfall on the thickness of the weathering mantle, the runoff erodes the surface soil and weakens rock formations. Also, the infiltrated water increases the pore water pressure, alters the weathering of the existing surface (and subsurface) geomaterials and consequently increases the thickness of the weathering mantle (rating: 4). On the other hand, the thickness of the weathering mantle does not influence rainfall at all (rating: 0).

From Table 1, it can be seen that hydrogeological conditions and thickness of the weathering mantle are the most interactive parameters (C+E=37), while rainfall is the least interactive (C+E=19). This suggests that rainfall does not depend on the influence of the other parameters but is an independent agent, concerning the whole system.

Based on the above, the landslide susceptibility index (LSI) values in the resulting susceptibility map vary within the range of 0 and 100. LSI values were classified into seven categories, namely “Negligible”, “Low”, “Middle”, “High”, “Very high”, “Extremely high” and “Landslide”, according to the classification for landslide susceptibility shown in Table 3 (Brabb et al., 1972). The higher the LSI, the more susceptible the area is to landslides. In this study, the critical zones were those corresponding to an instability index higher than 49%, the “Very High” and “Extremely high” zones.

Table 3: Classification for relative landslide susceptibility proposed by Brabb et al. (1972).

3.2 Validation of the landslide susceptibility map

In order to test the performance of the produced susceptibility map, it was compared with the distribution of the major landslide events that had occurred in the study area and the predicted map showed very satisfactory results. To be more specific, in the landslide susceptibility map of Megalopolis, 5% of the locations of actual landslides correspond to the “Very high” and 95% are associated with “Extremely high” landslide susceptibility (Figure 7, Table 4)”. The susceptibility map shows that slope failure incidents were located in areas where flysch formations, schist-cherts and Neogene sediments outcrop on slopes near major fault zones and thrust surfaces. Moreover, in order to examine the potential landslide risk in respect to settlements, villages and cities of the study area were overlaid on the susceptibility-hazard map. This correlation suggests that 45 settlements are entirely or partially located within “Very high” or “Extremely high landslide susceptibility” zones”.

Table 4: Calculation of Instability Index of Megalopolis area.

4. Conclusions

In this paper, landslide susceptibility was assessed by examining ten landslide parameters using RES and a GIS geoprocessing tool called Model Builder. Based on the selected parameters, all interactions that were revealed have been implemented through an interaction matrix. The outcome of this procedure produced the final susceptibility map for the southwestern part of Arcadia in Greece. The validity of this approach was tested using the slope failures that had occurred in this particular region. The instability index of all those slopes was found to be larger than 49 (out of 100), proving their instability (e.g. 21 out of 21 landslides were observed in either the “Very high susceptibility” or “Extremely high susceptibility” zone). In addition, it became clear that many more potential landslides could take place in the wider region of Megalopolis in the future.

Based on these positive results, we are confident in saying that the adaptability of RES to local conditions and to the given characteristics of existing geodata allow the use of an efficient tool in estimating landslide susceptibility hazard by adopting parameters that can be quantified more easily compared to other more expensive and time-consuming techniques. Moreover, experts (geologists, engineers) can use the RES approach before site investigations (geological–geotechnical) take place without knowing in advance if any slope failure has occurred in this area already.

As a consequence, RES could be an inexpensive and effective tool in ranking the instability in natural and man-made slopes and be useful in decision making regarding land use and development planning processes for zoning areas of potential landslide phenomena, such as those of southwest Arcadia.


Ayalew, L., Yamagishi, H. 2005. The application of GIS-based logistic regression for landslide susceptibility mapping in the Kakuda-Yahiko Mountains, Central Japan. Geomorphology. 65(1–2). 15–31. DOI: 10.1016/j.geomorph.2004.06.010

Brabb, E., Bonilla, M.G., Pampeyan, E. 1972. Landslide susceptibility in San Mateo County, California. US Geological Survey Miscellaneous Field Studies, Map MF-360, scale 1:62,500 (reprinted in 1978).

Guzzetti, F., Reichenbach, P., Ardizzone, F., Cardinali, M., Galli, M. 2006. Estimating the quality of landslide susceptibility models. Geomorphology. 81. 166–184. DOI: 10.1016/j.geomorph.2006.04.007

Hudson, J. 1992. Rock Engineering Systems: Theory and Practice. Ellis Horwood Limited: Chichester.

Koukis, G, Ziourkas, C. 1991. Slope stability phenomena in Greece: a statistical analysis. Bulletin of Engineering Geology and the Environment. 43(1). 47–60. DOI: 10.1007/BF02590170

Ladas, I., Fountoulis, I., Mariolakos, I. 2007. Large scale landslide susceptibility mapping using GIS-based weighted linear combination and multicriteria decision analysis – a case study in Northern Messinia (SW Peloponnese, Greece). 8th Greek Geographical Congress, 4-7 October 2007, Athens.

National Land Service (Lithuania). 2008. Spatial Analyis and Modeling. Training material GII-07, Civil Servants’ Learning Program: Distance Learning of Geographic Information Infrastructure.

Rozos, D., Pyrgiotis, L., Skias, S., Tsagaratos, P. 2008. An implementation of rock engineering system for ranking the instability potential of natural slopes in Greek territory. An application in Karditsa County. Landslides. 5. 261–270. DOI: 10.1007/s10346-008-0117-4

Rozos, D., Bathrellos, G., Skillodimou, H. 2011. Comparison of the implementation of rock engineering system and analytic hierarchy process methods, upon landslide susceptibility mapping, usin GIS: a case study from the Eastern Achaia County of Peloponnesus, Greece. Environmental Earth Sciences. 63(1). 49–63. DOI: 10.1007/s12665-010-0687-z

Tavoularis, N. 2017. The contribution of landslide parameters to the prognosis of slope failures. PhD thesis, School of Mining and Metallurgical Engineering, National Technical University of Athens (Greece) (under review).

This article has been published in European Geologist Journal 44 – Geology and a sustainable future.

Read here the full issue: