Impact of geo-environmental factors on landslide susceptibility using an AHP method: A case study of Fruška Gora Mt., Serbia

The paper considers the outcome of multi-criteria analysis of landslide susceptibility on the NW outskirts of Fruška Gora Mountain, Serbia. The area of the interest is known for landslide occurrences, and to focus on the most affected areas, it was necessary to consider some principal factors (lithology, slope inclination, rainfall, erosion, vegetation, altitude and slope aspect) and sort them by their importance to the phenomena. Prior to any criteria assessment, available data records had been assembled and refashioned as raster datasets. Thereafter, the criteria arising from an analytical hierarchy process (AHP) provided their weights of preference in the final model. In addition, the model was analysed for the information gain and classified in accordance to the optimal informativeness. Being tailored in the context of raster modelling, aided by the GIS spatial tools, our result gained substantial correlation to the control reference map (a digital photo-geological interpretation map of active and potential landslides).


Introduction
Over the past few decades, the geographic information system (GIS) has been applied to a variety of spatial-related problems.Throughout this period, the applications of its use have been proliferating, achieving more and more impressive results.Thus far, the im-provement has affected many fields of science and engineering, but the practice of GIS use has not been equally pervasive around the world (ALLEOTTTI & CHOWDURY 1999, CHACÓN et al. 2006).Namely, sufficiently developed GIS usage in developed countries contrasts with the situation in developing countries, which results in unbalanced insights regarding natural Impact of geo-environmental factors on landslide susceptibility using an AHP method: A case study of Fruška Gora Mt., Serbia phenomena.Engineering geology practice in Serbia is an example of the latter.This paper represents an attempt of the above-mentioned implementation of the GIS within an engineering geology scope.Environmental hazards, which are addressed within the field of engineering geology, affect both the social and the economic aspects of human lives.Hazards strike at different rates, with varying intervals and duration, leading to different outcomes.Hence, it becomes desirable, if not necessary, to predict their behaviour, number and severity prior to their potential triggering.It has been shown that GIS-based techniques are powerful tools for handling such challenges (CHACÓN et al. 2006;BONHAM-CARTER 1994).
Herein, one of the most widespread hazard phenomena (ALLEOTTTI & CHOWDURY 1999;SMITH 2001), more precisely their susceptibility, is to be considered.This addresses landslides and mass movements alike.

Case study area
The study area is located in the NW part of Serbia, on the mountain Fruška Gora, in the vicinity of Novi Sad.The site is contoured by the river terrace of the Danube on the north, the central mountain's ridge on the south and local ridges along east and west (Fig. 1).The area spreads over approximately 85 km 2 of hilly landscape, with intriguing geoenvironmental features.As such, it has been widely studied in many aspects, including slope stability, and this paper is another contribution in this regard.Prior to any modelling and for the sake of appreciating the phenomena thoroughly, it is advisable to consider all of the environmental factors that are important for slope stability.Thus, only natural factors have been regarded, despite the apparent influence of human activity on slope stability.Hence, considered factors comprise climatic features, lithology and other aspects of geological setting, geomorphological characteristics, hydrological, hydrogeological features and finally, engineering geological properties.A brief presentation of some essential properties follows.
Although it is not of significant altitude (with the summit being slightly over 500 m), this mountain exhibits some climatic variability, particularly in the distribution of rainfall regime and intensity, which varies drastically from the base to the high-ground.This implies that the mountain shape and disposition, rather than its altitude, influence the distribution of moist air masses and the overall precipitation.Namely, moist air abuts the northern slopes and condenses upward, as the temperature decreases (by 1°C/200 m).On the contrary, as it descends down the southern slopes, moist air abuts a warmer environment and accordingly provides less rainfall.This effect is magnified due to the asymmetry of Fruška Gora Mt. because northern slopes rise abruptly in alti-tude, and mild southern ones gradually subside to a plain.The studied area belongs to the northern, moister realm, with drastic changes in rainfall (gradient of approximately 36 mm/100 m).These changes have a great impact on slope stability and erosion.
The geological setting of the entire Fruška Gora Mt. implies a zonal lithological and structural setting because of the complex horst-anticline forming the core of the mountain.The study area encompasses the NW part of this anticline, with the typical succession (Fig. 5a), starting with low-grade the Paleozoic crystal schists in the anticline base.Scattered in the stripelike segments, these metamorphic associations occupy the higher ground.They are characterised as a green formation, composed of a mixture of altered magmatic and sedimentary rocks, with regional faults and folds of W-E trends within.Subsequently, a portion of the Triassic basal sediments (conglomerates and sandstones) gradually shifted towards limestone, implying localised subsidence of the paleo-relief at the time.This lasted until the early Jurassic, when another intense uplift occurred, followed by some minor volcanic activity.During the Jurassic period, the far more prominent movement was the one related to the closing of the oceanic basin on the south, which left peridotite (serpentinite) thrusts and diapirs as evidence of this event.This movement culminated in the early Cretaceous, followed by minor gulf formations of coral limestone sequences, known as the Bačko-Banatska zone.The Post-Mesozoic tectonics had re-established W-E trends of structures at regional scale and induced NW-SE oriented faults, traversing the former structures.The tertiary is chiefly represented by marine clastites, gaining more carbonate components as the basin turned more limnic during the late Neogene.This interval is characterised by diverse lithology, ranging from sands and clays to limestones, via marls and other transitional forms.The most significant and the most widespread Quaternary unit is loess.It covers the lower landscape, flattening it towards the Danube's alluvium and ending with steep cliffs facing the river.The most recent Quaternary unit includes the fluvial deposits of permanent and periodical flows, represented by gravels and sands or their loose aggregations (ČUPKOVIĆ 1997).
Predominant geomorphological entities are of the aeolian and fluvial origin.In conjunction with other processes, these processes sculpted the current landscape of the terrain.Fluvial and aeolian processes alternated in supremacy during the terrain evolution, meaning that they had different enrolments at different times.Apart from these, other morphological processes left some significant imprints, such as larger landslides, proluvial fans and cones, scree slopes and gullies.
Hydrological and hydrogeological regimes are not overlapping the meteorological one, even though such overlapping could be expected.This discrepancy is again influenced by complex geological features.The ultimate aquifers are the Paleozoic schists due to their super-sized voids (i.e., the outcomes of multiple tectonic actions).As such, they provide a more balanced regime (minute annual variations of water table levels) and more constant water temperatures during the year, implying minor direct influence of meteorological phenomena (PETKOVIĆ et al. 1976).The water balance analysis suggests that only 30% of the atmospheric precipitation runs out surficially, whereas the rest leaves in evapo-transpiration or as groundwater.Hence, it is more likely that springs and streams are governed by the groundwater regime.Generally, there could be three to four major groundwater horizons specified in the schistosity core, and many smaller localised accumulations in different aquifers.It could be inferred that the former are the most important for the overall water regime of the area, whereas the latter still could be significant by their local influence on the slope stability.
As for the engineering geological properties, it is necessary to stress that surveys have not yet been sufficiently detailed.Seismological features imply relatively stable ground, even though an active fault zone propagates throughout the midst of the mountain, creating possible seismic hazards, especially in loose rock masses.Apparently, the study area seems to be sufficient in size and distance from this zone to be regarded as uniformly affected by minute seismic inconveniences from time to time (PAVLOVIĆ et al. 2005)

Material and Methods
Initially, our approach addresses the optimal selection of scale versus complexity of the problem to be modelled.At this point, the chosen mid-scale (1:50000) tolerates not only the extent of approximation or even exclusion of some factors but also certain subjectivity in the selection of class intervals (SÜZEN 2004).The former corresponds to the combination of the heuristic and semi-quantitative approach, which is believed to provide high quality insight in regard to susceptibility, hazard or risk assessment ( VAN WESTEN et al. 2006).Thus, the landslide susceptibility has been accessed in a somewhat subjective manner, but the inevitable subjectivity adhered to herein appears to be quite desirable.
Modelling was tailored in a sense of multi-criteria analysis, whereas natural environmental factors have been chosen as individual criterions.Refashioned by the scope and the requirements of this research, the modelling procedure generally fits the usual patterns of similar problems (Fig. 2).Initially, it addressed the digital elevation model (DEM), essential for the other implemented models.Additionally, it encompasses the following: rainfall distribution, altitude, aspect and slope models, linear erosion pattern, vegetation distribution and finally, the lithology model.
• DEM was derived from digitised topographic map at 30-m resolution (the same as in all of the following data sets).To increase precision and flatten the outliers such as hill peaks or minor depressions, the DEM was adjusted with tools at hand (standard tools Impact of geo-environmental factors on landslide susceptibility using an AHP method 93 Fig. 2. Schematic display of the research procedure and resources (light grey pieces stand for source data, medium grey for the analyses and calculations, whereas dark-grey ones mark the modelled outcomes).
of Watershed analysis, ArcGIS package).Its original digital number (DN), extending from 0 to 255, was reclassified to a range from 0 to 100.Resembling the percentages, reclassified DEM eased the merging with other models, including the models derived from DEM itself.
• The slope inclination model is considered to be of great importance because of its direct physical relation to the slope instability.Moreover, there are empiric proportions on how the inclination relates slope stability in homogeneous masses (GIRAUD 2007).Practically, the slope model is another morphometric feature derived directly from DEM.Because thorough classification would require extensive field work, general assumptions on slope inclination intervals were adopted as imposed by a host of authors (ABOLMASOV 1997;VOŽENÍLEK 2000;GIRAUD 2007,).They suggest distinction between several classes, particularly 0-5°, 5-10°, 10-15° and >15°, where the last two classes include very susceptible slopes, responsible for approximately 95% of landslide occurrences in similar area.All of the classes are assigned DN values according to the proneness they induce.Greater angles make slopes more prone to sliding, and lesser angles make slopes less prone to sliding.It should be stressed that this proportion couples all types of slope instabilities, not only landslides (which, in particular, rarely exceed 15-20°) but also landfalls, creeps, debris flows, etc.Moreover, it is apparent that the critical slope angles entirely depend on geological setting (lithological content and structural features).This is exactly why the criterion has been conferred with references of authors that have addressed this problematic more thoroughly (Fig. 3a).
• Slope aspect is usually a factor with minor influence and tends to be excluded in cases where several other factors contribute in greater proportion.Here it will be considered, but will be given appropriately inferior weight in the final model.Generally, this model substitutes environmental phenomena such as changes in the moisture content or changes in depth of the eluvial crust, both of which are caused by seasonal and diurnal variations of the solar path.At the longitudes in question, it could be expected that northern slopes would retain higher moisture content, whereas southern ones would have more profound eluvium, which leads to a dilemma regarding how to classify this model.Nevertheless, the study area is rich in forestall vegetation, which reduces the crust depth and provides rigidity, but amplifies the moisture absorption.Therefore, the slope stability is primarily influenced by the moisture content as far as the aspect model is concerned (the variation in depth of eluvium could be justifiably excluded from the model).The slope aspect parameter discriminates between these phenomena as follows: NW slopes indicate the most adverse conditions (highest moisture content in the soil) and yield the highest DN values.Conversely, SW slopes provide the most favourable conditions (lower moisture, lower DN values), while moderate conditions characterise NE and SE slopes (Fig. 3b).
• The altitude model follows from the DEM as another morphometric feature.Terrain is classified into four altitude entities (Fig. 3c).Governed by the natural break divisions, a certain DN value has been assigned to each altitude interval.Because it is more likely that lower slopes offer better conditions and seem less susceptible to sliding, they have been assigned lower DN values than the higher ground.
• The rainfall distribution model depicts the data of the Hydrometeorological Survey of Serbia.Furthermore, this model yields abrupt differences while using records of average rainfall per month, particularly in July.This is when the rainfall varies most substantially, both over single month and single day periods.These particular cases have been processed (Fig. 3d) to simulate the most disadvantageous conditions.The matter of classification is quite delicate, considering that the distribution significantly differs from Gaussian.For this instance, a natural break method (Jenks' method) was used to depict rainfall classes (WEBSTER & OLIVER 2001).Four classes have been accepted, and the appropriate DN values were allocated to each class.
• Geomorphological processes implying flowing water as the medium govern the linear erosion pattern, particularly the drainage pattern of fluvial or proluvial origin.Because lateral erosion propels slope instabilities, geomorphological processes are coupled with the occurrences of landslides.The basic principle suggests that the slope segments in the vicinity of the drainage should suffer a greater impact on their stability than remote, higher segments (ridges).The modelling involved the calculation of the pixel distance (Fig. 3e) from the drainage pattern vector, as well as the classification of the later outcome.The initial distance calculation has been corrected to assign respectively lower DNs to higher slope segments and vice versa because lateral erosion dominates in the area below 200 m (Fig. 4).
• Vegetation cover was modelled due to its beneficial contribution to slope stability, especially because the forestall flora dominates the terrain.Some simple approaches and techniques were engaged to display the model of vegetation distribution.For instance, the normalised distribution vegetation index (NDVI) applies well.It required satellite imagery data, particularly the 3rd and 4th channel of the Landsat TM imagery (RAVI 2002).It uses the differences in reflective responses of vegetation versus bare soil or rock (VINSENT 1997).With regard to raster processing, a simple function combines the imagery, creating the vegetation distribution model.Naturally, the presence of vegetation lowers the slope susceptibility to sliding versus bare soil, so DN values had been specified for two classes: with and without significant vegetation cover (Fig. 3f).
• Lithology predominantly determines the susceptibility pattern.However, the geological setting is complex and diverse (Neogene formations especially).To simplify this model, improvised classifications were used, and some different entities were merged into a single class.For example, the common DN value has been allocated to a solid rock masses (Paleozoic) and alluvion (Quaternary) because of their approximately equal unlikeness to host landslides.Accordingly, those are the lowest DN values in this model.In contrast, loose and clayey grounds had the highest DN values (Fig. 5a, b).

Results and Discussion
The final model gathered all of the raster data sets, yielding a susceptibility pattern based on the superposition of their weights (their relative influences).As a convenient procedure that addresses multi-criteria hierarchical structures (GENEST 1994, SAATY 2003, the analytical hierarchy process (AHP) ERCANOGLU et al. 2008) has been applied to perform the weighting of the influence for each raster model via a pairwise ratio scale.The foregoing procedure has been fashioned as a visual basic application (VBA) macro for the ArcGIS package (MARINONI 2004).Prior to obtaining weights, the procedure applied the gross estimation of the factor's preference, based on experts and their experiences (VOŽENÍLEK 2000;ESMALI 2003;KOMAC 2005;ÇAÐÝL et al. 2006;ERCA-NOGLU et al. 2008).For this instance, a nine-point scale (the range from 1/9 to 9) has been chosen to reflect the pairwise relations between input raster sets, yielding a two-dimensional reciprocal and inconsistent matrix -the comparison matrix (Table 1).
Normalisation of the matrix and averaging by rows generated the priority vector (SAATY 2003, GENEST 1994), which represents the distribution of the weights (Table 2, shaded columns).Thereafter, it was necessary to establish the procedure for shifting from an inconsistent to a near-consistent matrix.Versatile solutions proposed by different authors had been considered (GENEST 1994, LAININEN 2003, SAATY 2003, PO et al. 2007).However, they all prove to be slightly different from the outcome of the original technique (SAATY 1977).It was feasible to control the matrix consistency on the simplest basis, i.e., by Saaty's consistency parameters CI and CR (consistency index and consistency ratio, respectively) and criterion (CR<0.1).In this way, the initial subjectivity of the weights distribution (Table 1) has been unbiased up to a certain level, leaving the refined weights depicting the final pattern (Table 2, shaded columns).
Finally, the priority vector or, more appropriately, the linear distribution function of the weights, ap- where M i corresponds to the influence factor's DN values respective to their appearance in the matrices in tables (M 1 = lithology, M 2 = slope,…,M 7 = aspect).As the indices relate to the corresponding factor, the calculus of their weighted DN values generates the final raster -the model of susceptibility (Fig. 5c).
This model depicts the spatial distribution of the susceptible zones, separated into four classes: low, mild, moderate and high susceptibility.According to the presented approach, the first class (black in Fig. 5c) corresponds to the areas where the input raster had the smallest contribution, and the fourth class (dark grey in the Fig. 5c) stands for the zone with the highest overall contribution.
In addition, we need to become conversant with the classification criteria used for the final raster.Namely, it has been speculated what number of classes would be optimal for landslide susceptibility maps (CHACÓN et al. 2006).However, there is no decisive formulation on this topic, so we employed the entropy approach (PÁSZTO et al. 2009).
We have tested the final raster for the entropy function behaviour over the entire DN span (from 1 to 256 classes) with natural breaks intervals to define the highest information gains (Fig. 6).The trend curve approximates the behaviour of the entropy function, proving the criteria for information analysis.We considered only the span from 2-20 because the trend is fluctuating at higher values.Information gain reaches its maximum when raster is displayed in nine classes, with sub-maximum at four classes, leaving us to choose between these two cases.The former would be quite uncommon in relation to the usual practice, so we have chosen a four-class model (low, mild, moderate, and high susceptibility).
The slope stability map (PA-VLOVIĆ et al. 2005) which has been processed through the remote sensing techniques and the fieldwork was used for comparison as a control reference.
The susceptibility map (Fig. 5c) suggests that the area could be qualified as predominantly stable ground, with the second class (mild susceptibility to sliding) widespread.Severely endangered areas are the northern slopes (by the riverbank of the Danube), as well as the internal areas in the NE.
The basic statistical comparison between the model and the landslide inventory (Fig. 5b) resembles the proposed visual evaluation.Figure 7-a reveals that active landslides mostly occupy pixels with values of the fourth class (high susceptibility).The same applies to the category of temporarily inactive landslides; overall, these landslides contain approximately 20% of the fourth class pixels and 14% of the third.However, the terrain, with its portion of the second class exceeding 50% (Fig. 7b), predominantly remains stable but locally inclines to instabili-Impact of geo-environmental factors on landslide susceptibility using an AHP method 97   ties because by area it is nearly 30% third class (moderate susceptibility).
The numbers of the pixels within the first and the second class are surprising due to the imperfections of the model.Despite these unpredicted inaccuracies, the model yields fairly acceptable predictions of susceptible areas, which precedes possibilities for more detailed assessments, whether concerning susceptibility, hazard or risk.

Conclusion
The final model shows substantial correlation in the most critical areas, as detailed in Fig. 5c.Particularly, the landslide occurrences in valleys in the NE domain parallel the remote sensing data.Every single landslide polygon, either characterised as active or temporarily inactive, encompasses fair quantities of pixel values within the classes of moderate and high susceptibility.Nevertheless, it is obvious that criterions in use were too rigorous in some regions, most likely due to the emphasised influence of the lithology model.In particular, the entire northern outskirts are qualified to be moderately to extremely prone to instabilities, which might not be the case.Priority vector weights have been combined in such a manner to create a dilemma in analysing the results: at some points, our model fits the crucial remote sensing evidence but contradicts it in other regions.Seemingly, this could be avoided by more detailed lithology modelling or filtering the existing model, which would lead to an even more subjective approach, stressing the relevance of heuristics in the analysis.
To improve the final model, further directives must also address the perpetuation of other input data sets.This concerns the DEM, which could be processed in higher resolution, acquiring additional precision for the directly dependable models, such as models of altitude, slope, aspect and even linear erosion pattern.Furthermore, the inputs could be more temporally correlative (this is especially relevant to the rainfall distribution record).Finally, some superior methods in lithological quantification could be used and would minimise the perplexity and subjectivity in the classification.However, the majority of the capitalised adjustments and directives require extensive data, which were unavailable because they demand more detailed surveys.
In essence, the result errs in favour of safety, indicating that large portions of the total area could be potentially triggered.However, this remains arguable and should be regarded cautiously, even though it should not differ from the actual state of susceptibility over the terrain (according to authors such as KOMAC 2005or ERCANOGLU 2008, where the similar approach matched very rigorous criteria).

Fig. 4 .
Fig. 4. The cross section of the riverbed (bold line) in comparison to the theoretical curve of erosion basis (dashed line).Note that the interception point borders vertical and lateral erosion preference at approximately 200 m, wherein vertical erosion dominates higher grounds and lateral lower grounds.

Fig. 6 .
Fig. 6.Information gain of landslide susceptibility model based on the entropy function (the highest gain appears in classes 4 and 9) and four-high susceptibility); d, susceptibility model in comparison to the landslide inventory.

Fig. 7 .
Fig. 7. a, histogram of pixel distribution within the landslide polygon area per each susceptibility class (continuous grey line stands for the active landslides, dashed for temporarily inactive ones and the bold black line is overall); b, log-normal distribution of data from previous plot with standard deviation, sd, and arithmetic mean, m.