Comparing multi-crit ria methods for landslide susceptibility mapping in Chania Prefectur , Crete Island , Greece

Introduction Conclusions References


Introduction
Landslides are the most commonly occurring geological hazards worldwide and the second most frequent natural catastrophic events, after hydro-meteorological events.Landslides are considered to be one of the most catastrophic natural hazards due to the fact that they cause extensive damages to constructions and infrastructures as well as thousand casualties annually (Saha et al., 2002;Akgun and Turk, 2007).
Landslide hazard assessment and risk reduction can be accomplished by providing the risk managers with easy accessible, continuous, and accurate information about the landslide susceptibility.Thus, an accurate susceptibility mapping can be the key information for a large variety of users from both private and public sectors, from governmental departments and the scientific community and from local to international levels.
A landslide susceptibility map depicts areas likely to have landslides in the future by correlating some of the principal factors that contribute to landslides with the past distribution of slope failures (Brabb, 1984;Fall et al., 2006).Crucial factors for the construction of reliable susceptibility maps are the quality and the amount of available data and the selection of the best methodology analysis.The process of creating these maps involves different approaches and methods distinguished as: qualitative, semi-quantitative and quantitative (Lee and Jones, 2004;Castellanos Abella and van Westen, 2008).Quantitative methods are based on mathematical expressions of the correlation between controlling factors and landslides.The two types of commonly used quantitative methods are the deterministic and the statistical (Ward et al., 1981;Nash, 1987;Terlien et al., 1995;Kramer, 1996;Atkinson and Massari, 1998;Aleotti and Chowdhury, 1999;Fall and Azzam, 2001a, b;Akgun and Bulut, 2007).The statistical approach uses statistical and probabilistic methods.According to these methodologies the previous landslides can be related to measurable elements of the landscape (Koukis and Ziourkas, 1991;Dai et al., 2002;Ohlmacher and Davis, 2003;Akgun and Bulut, 2007) and these elements can be used to predict possible future landslides.Typical multivariate statistical approaches used to map landslide susceptibility are the weight of evidence (WoE) modeling method, a quantitative "data-driven" method used to combine datasets (Lee et al., 2002;van Westen and Lulie, 2003;Lee et al., 2004a;Thiery et al., 2007;Mathew et al.,. 2007;Neuhauser and Terhorst, 2007;Poli and Sterlacchini, 2007;Rezaei Moghaddam et al., 2007;Dahal et al., 2008a, b;Sharma and Kumar, 2008;Barbieri and Cambuli, 2009;Ghosh et al., 2009;Regmi et al., 2010) and the logistic regression method (Wieczorek, 1996;Atkinson and Massari, 1998;Guzzetti et al., 1999;Gorsevski et al., 2000;Lee and Min, 2001;Dai et al., 2001;Dai and Lee, 2002;Ohlmacher and Davis, 2003;Ayalew and Yamagishi, 2005).
Furthermore, during the last decades, many studies (such as van Westen, 1994;Carrara et al., 1995;Aleotti and Chowdhury, 1999;Dai et al., 2002;Cevik and Topal, 2003;Ayalew and Yamagishi, 2005;Fall et al., 2006;Carrara and Guzzetti, 1995;Carrara et al., 1995;Dikau et al., 1996;Luzi et al., 2000;Miles and Ho, 1999;Miles and Keefer, 1999;Carrara et al., 1999;Fall, 2000;Refice and Capolongo, 2002;Fall et al., 2006;Kouli et al., 2010;Arnous, 2011) have been carried out to overview the use of GIS for landslide susceptibility assessment.In a geologically active country like Greece the intensive tectonic activity is not related only to earthquake events but also to other natural hazards.Active tectonics generates steep morphological structures, strained geological formations and intensive weathering procedures which are considered to be the main causal factors of landslide events.
Crete Island is located in the southern part of the Hellenic fore-arc, an area of important tectonic deformation and high seismic activity.This high seismic activity is caused by the subduction of the oceanic African lithosphere beneath the continental Anatolian-Aegean lithosphere (Papazachos and Comninakis, 1971;McKenzie, 1972;Angelier, 1979;Makropoulos and Burton, 1984).
Chania Prefecture, located at the western part of the Island (Fig. 1), suffers from severe landslide phenomena.The lack of a landslide susceptibility mapping in the approximately 2330 km 2 study area makes the risk management impossible for the local authorities.The landslide susceptibility analysis introduced in this study covers the entire Chania Prefecture and it is conducted by the application of the Weight of Evidence (WoE) statistical modeling method and the Weighted Linear Combination (WLC) semi quantitative approach.Besides the obvious benefit of providing landslide susceptibility maps, the study compares the efficiency of the applied methods.

Geological framework of Chania prefecture
Crete Island is characterized by an extremely complicated geological structure with intensive tectonic fragmentation.Although many researchers through the years have studied the geological evolution of the island (Tataris and Christodoulou, 1965;Bonneau, 1973;Creutzburg, 1977;Fitrolakis, 1980;Tsiampaos, 1989), there are still many unsolved issues.Based on current theories, the island constitutes of repeated tectonic covers consisting of variant geotectonic zones' geological series.
Considering the surficial geology of Chania prefecture (Fig. 2), Quaternary deposits forming depositional plains dominate especially at the northern part of the study area, as well as along the coast line.Miocene to Pliocene sediments crop out in the central and the northwestern part of the study area and carbonates of the Tripolis nappe in the northeastern part.Dissected hills of phyllites and quartzites, a late Carboniferous to late Triassic package of sedimentary rocks composed mostly of quartz-rich siliciclastic sediments, with minor limestone, gypsum, and volcanic rocks (Krahl et al., 1983)  Geological, topographical, land use and precipitation data were collected and introduced in a Geographic Information System (GIS).By processing the primary data layers the following landslide casual factors were selected and prepared as secondary data layers: unified lithology, faults proximity, elevation, slope gradient, aspect and curvature, roads and rivers proximity, precipitation, and land use types.
The geological maps of the Institute of Geology and Mineral Exploration (IGME) presenting lithological and structural units, at a scale of 1 : 50 000 were digitized.The Digital Elevation Model, with a cell size of 20 m, was generated from the topographic maps of the Hellenic Military Geographical Service, at a scale of 1 : 50 000.Many terrain attributes such as slope angle, aspect, curvature and stream network were computed using the digital elevation model.Furthermore, the Corine Land Cover map 2000 (CLC2000 100 m, version 1) of the European Environment Agency (EEA, Copenhagen, 2000; http://www.eea.europa.eu) was adopted for the assignment of the land use classes.Precipitation data covering a time period of 25 yr were also introduced in the GIS.The landslide inventory map of 108 landslide locations was created using data coming up from field investigation, previously documented events and Google Earth satellite images interpretation (Fig. 3).At the inventory map all the landslides are introduced as points because of the small scale enforced by the great extent of the study area.All thematic maps were generated using ArcGIS Desktop 9.1 software.All the raster layers were resampled using a cell size of 20 × 20 m.
The formations of the first four classes compose of alternating layers with grate variations in their mechanical and permeability properties setting the preparatory causal factors (Terzaghi, 1950;Popescu, 1994) for the manifestation of translational landslides or creep movements.On the contrary, at the areas occupied by limestones and marbles only rock falls, determined by the fragmentation degree of the rock mass, can be expected (Fig. 3).

Faults proximity
Tectonic structures associate with extensive fractured zones and steep relief anomalies, especially in rock formations.As a result these zones present favorable conditions for the manifestation of landslides.Therefore, major tectonic structures, digitized from the IGME geological maps, were included as a secondary data layer in this study.
For the introduction of the tectonic structures influence, multiple buffer zones were generated around the faults lineament, at distances of 100 and 200 m (Fig. 4).The three classes generated by means of the above described procedure were rated accordingly, as presented in Table 1.

Land use types
Land use is also related with the triggering and preparatory causal factors of the landslides.For instance, the strong root system of the woody vegetation provides both hydrological and mechanical effects that generally stabilize slopes.As a result, the landslide events recorded in woodlands are much less than those of the un-vegetated or irrigated cultivated areas.The role of vegetation in the slope stability was introduced by rating the main land use classes of Chania prefecture, as mentioned at the Corine Land Use map 2000.The eight main land use classes were ranked as shown in Table 1.

Rivers proximity
The density of the drainage network is considered as an important factor in characterizing landslide susceptible areas.Fluvial erosion is one of the most common triggering causal factors of the landslides (WP/WLI, 1994), and it usually affects the slopes toe.
In this case study the drainage network was automatically extracted from the DEM and the drainage tributaries were classified according to Strahler's system (Strahler, 1957(Strahler, , 1964)).
A drainage network buffer zones data layer was produced considering the streams order (Fig. 4).Generally the streams of each order were multiply buffered at two to three distances giving rise to three to four classes, respectively.The multiple buffer zones applied in each steam order as well as the rating values are analytically presented in Table 1.

Slope angle and curvature
Slope geometry is a vital controlling factors in slope stability (Huang and Li, 1992;Wu et al., 2001), and a digital slope image is therefore a fundamental part of a hazard assessment model (Clerici et al., 2002;Saha et al., 2002;Cevik and Topal, 2003;Ercanoglu and Gokceoglu, 2004;Lee et al., 2004a, b;Yalcin, 2008).Slope angle was extracted from the DEM, and the derived slope angle map saw values ranging from 0 • to 84 • .The produced slope image was reclassified into six slope angle classes (0 • -5 • , 6 • -15  (Kayastha et al., 2012).Slope curvature was classified into three classes namely convex, flat and concave.

Slope aspect
The slope aspect is defined as the downslope direction of the maximum rate of elevation change.Aspect can influence indirectly the landslide initiation because it controls the exposition to several climate conditions (duration of sunlight exposition, precipitation intensity, moisture retention, etc.) and as a result the vegetation cover (Wieczorek et al., 1997;Dai et al., 2002;Cevik and Topal, 2003;Suzen and Doyuran, 2004;Komac, 2006).
Considering that in Greece, because of the major mountain chains orientation, the NW oriented slopes are more violently affected by rainfalls, the slope aspect was classified accordingly.The nine slope aspect classes presented in Table 1 allow the application of rating values covering all possible climate conditions exposition.

Elevation
The elevation is a parameter secondarily related to the landslide occurrence.For instance, the high relief areas are usually occupied by the most cohesive formations (e.g., limestones), they are exposed in unfavorable climate conditions (e.g.intensive rainfalls), they are uninhabited and landslide events rarely take place and even more relay are recorded there.The low elevation areas occupied by Quaternary and Neogene sediments, depending on the roughness of the morphology, can present low to very high landslide frequency.So there is not any general rule relating the landslide occurrence with the elevation.The reclassification of the elevation data layer and the rating values of the classes are highly depended on the expert's opinion.The expert has to correlate as many geological and morphological factors as possible with the elevation in order to introduce correctly the elevation data layer to the Geographic Information System (GIS).In this study the higher rating values were given to the classes covering the elevations from 700 to 1000 m (Table 1), presenting steep morphology and mainly occupied by flysch formations and the phyllites -quartzites unit.

Precipitation
High precipitation is characterized as the physical processes constituting the main triggering causal factors of landslides (WP/WLI 1994; Koukis et al., 1997;Polemio and Sdao, 1999;Sdao and Simeone, 2007).For this work, precipitation data for a timeperiod of about 25 yr were collected from three different agencies; the Department of Hydrology and Water Resources -Crete Region (DHWR-CR), the National Agricultural Research Foundation located in Chania in Crete Island (NAGREF) and the Branch of the Institute of Geology and Mineral Exploration situated in Rethymnon in Crete Island (IGME) (Kouli et al., 2009).The average annual precipitation of the study area ranges from 800 to 1400 mm.Three classes were defined, using an equal interval of 200 mm (Table 1).

Roads proximity
Road construction is related to extensive excavations, application of static and dynamic loads, vegetation removal etc. along natural and engineered slopes.These landslide triggering actions (WP/WLI 1994) were considered in the design of the landslide susceptibility maps by introducing a road network buffer zones data layer.
Multiple buffer zones were applied, within distances of 20, 50 and 100 m, to the road network producing four classes (Table 1).

Methodologies
At a first step, an expert based, relative weighting -rating approach (WLC) was used for the landslide susceptibility mapping.According to this method, certain weights (Table 1) were assigned to each factor by taking into account the specific landslide occurrence parameters of the study area.The lithology-based weighting-rating system introduced in Kouli et al. ( 2010) was adopted as it was proved to be satisfactory accurate.After that, the rating values assigned to the different classes of the factors thematic layers were introduced as attribute information in the GIS and an "attribute map" was generated for each data layer representing a causative factor.The vector data layers were reclassified using the assigned rates and the related raster data layers were produced.Finally, the reclassified raster layers were used as input parameters for the raster calculator function.By means of the raster calculation the lithology-based weighted factor maps were overlaid (summed) to provide the landslide susceptibility maps.
Moreover, the weight of evidence (WoE) method was applied using the 80 % of the inventory map landslides as the training set and the rest 20 % for accuracy assessment purposes.
The weight of evidence (WoE) modeling method was firstly developed for medicine purposes (Spiegelhater and Knill-Jones, 1984), but it was also applied to the identification of mineral deposits potential (Bonham-Carter et al., 1988, 1989), to landslide susceptibility mapping (Lee et al., 2002;van Westen et al., 2003;Lee and Choi, 2004;Thiery et al., 2004Thiery et al., , 2007;;Mathew et al., 2007;Neuhauser and Terhorst, 2007;Poli and Sterlacchini, 2007;Rezaei Moghaddam et al., 2007;Dahal et al., 2008a, b;Sharma and Kumar, 2008;Barbieri and Cambuli, 2009;Ghosh et al., 2009;Regmi et al., 2010;Kayastha et al., 2012;Xu et al., 2012;Armas, 2012) and most recently to groundwater mapping (Lee et al., 2012).Bonham-Carter et al. (1988, 1989) gave a detailed description of the mathematical formulation of the method.Weights for each landslide predictive factor are calculated based on the presence or absence of landslides within the different classes of a causative factor (Kayastha et al., 2012) (1) where W + and W − are the weights for the presence or absence of landslides within a certain class of a causative factor map, P {A|B} is the conditional probability of A occurring given the presence of B, F signifies the presence of a landslide, L is a class of a causative factor, and the bar above a symbol signifies the negation (Kayastha et al., 2012).Applied to the landslide susceptibility estimation, the technique of log-likelihood ratios aims to identify the degree of influence, expressed as "weight" that each variable has on the development of a landslide event (Barbieri and Cambuli, 2009).Weights are calculated based on the spatial development of landslides in the thematic maps used as evidence.
The statistical significance of the weights can be verified based on their variances, S 2 (W + -W − ), which can be estimated as (Bonham-Carter, 1994): The contrast, C = W + -W − for a given class provides a measure of the correlation between the class of a causative factor and the occurrence of landslides.It measures the degree of correlation between point pattern and binary map, and can be tested for statistical significance (Agterberg and Cheng, 2002).For a spatial relation, the value of C is positive, and when a spatial relationship is lacking, the value is negative.When the landslides are randomly distributed within the study area, then W + = W − and C = 0.The variance of the contrast, S 2 (C), is given by the sum of S 2 (W + ) and S 2 (W − ), and the studentised contrast, C/S(C), gives a measure of confidence (Neuh äuser and Terhorst, 2007).The higher the studentised contrast the stronger the relationship between the causative factor class and the occurrences of landslides.If the studentised contrast is high, there is likely a positive spatial association between the occurrence of landslides and the causative factor class while if the studentised contrast is small, the causative factor class disfavors the occurrence of landslides.

Landslide susceptibility maps -results
The landslide susceptibility maps with both WLC and WoE methods were produced within a raster/grid GIS and were reclassified into low, moderate, high and very high susceptible zones, using the natural breaks classification method.The Landslide Hazard Index (LHI) for each grid cell is given by the summation of the raster thematic maps after their multiplication by the corresponding weights.In the case of WLC method, the Landslide Hazard Index (LHI) for each grid cell was extracted by the summation of the raster thematic maps after their multiplication by the corresponding lithology-based weights (Kouli et al., 2010) as shown in Table 1.The LHI is presented by the expression as given below: where n is the total number of causative factors.
In the case of WoE application, the hazard index maps were obtained by combining the contrasts of each causative factor according to the following equation: where C i j is the contrast value for class j for the causative factor i and n is the total number of causative factors.
The resulting landslide susceptibility maps are shown in Figs. 6 and 7, respectively.The weights, variances, contrasts and studentised contrasts for all parameter classes after the application of the WoE method are presented in Table 2.According to this table and having in mind that high weights indicate high probability of landslide occurrence, the distance to roads and geology casual factors mainly govern the landslide phenomena in the study area while distance to rivers, slope curvature and elevation play an important role.As it seems the most significant causative factors classes with a positive impact on landslides are: (a) slope angle of 16 The above mentioned impact of the causative factor on the landslides occurrence appears to be reasonable.But a detailed examination of the WoE method results reviles some irregularities.For example, the existence of forests appears to have positive impact on the landslide occurrence, fact totally unacceptable.Obviously numerous roads intersect areas characterized as forests at the Corine Land Cover map confusing the WoE method results.Another obscurity of the WoE method is that the C i j values are strongly related to the quality of the landslide sample.Usually the landslides are recorded when they damage the road network, villages and, in general, when they interrupt the human activities.So usually the available sample of landslides is not representative.These irregularities appear to introduce a random error at the Landslide Hazard Index values.
According to Rossi et al. (2010), a combination of the two landslide susceptibility maps was produced using the logistic regression approach.For this purpose, the presence or absence of the landslides in the study area was the dependent variable, while the two landslide susceptibility maps were the independent explanatory variables.The resulted combined map is shown in Fig. 8.

Validation of the applied methods -success rate curves
In order to proceed to the validation of the landslide susceptibility maps the success rate curves (Chung and Fabbri, 1999;van Westen et al., 2003;Kayastha et al., 2012) were adopted.The "areas under the curves" constitutes one of the most commonly used accuracy statistics for the prediction models in natural hazard assessments (Begueria, 2006).The rate explains how well the model predicts the landslide (Chung and Fabbri, 1999) and the area under the curve can be used to assess the prediction accuracy qualitatively (Fig. 9).The cumulative percentage of observed landslide occurrence was plotted against the real cumulative percentage in decreasing LHI values to obtain the success rate curve for the study area.The rate curves were created and the "areas under the curves" were calculated for the three cases of hazard maps using the existing landslide location data.For example, in the case of the WoE, the 10 % of the study area can explains 67 % of all the landslides in the success rate and was classified as a zone of very highly susceptibility while at the same time the same 10 % of the study area can explain 53 % of all the landslides in te case of WLC method.Therefore, it is obvious that the WoE method resulted in better predictions (success rate of 87.4 %) compared to the WLC method (success rate of 84.7 %).The combined model shows 85 % area under curve reducing the uncertainties introduced to the WLC method by subjective knowledge of experts.

Conclusions
For the Chania Prefecture area two landslide susceptibility maps were produced adopting the weighted linear combination (WLC) and the weights of evidence (WoE) methods.A third map was produced from the combination of the two methods using the logistic regression approach.The produced landslide susceptibility maps are substantial for the land degradation management and planning of the study area and they clearly present the areas prone to landslides.The validation procedure confirmed the index setup procedure of both WLC and WoE methods and as a result the accuracy of the indicated hazardous areas.So both methods provided accurate susceptibility maps which they can be used safely from the local authorities for slope management and land-use planning.It must be noted that the WoE method provided more accurate results despite the limited inventory data.The current study also proved that, despite the subjectiveness introduced to the WLC method by utilizing the knowledge of experts, the produced susceptibility maps can be proved to be satisfactory accurate.The combined model resulted to the smoothing and balancing of the output susceptibility maps but it did not overcome the prediction accuracy of WoE.Nevertheless, in cases where the quality of the available inventory data is considered satisfactory and the landslide sample is statistically adequate, the combination of both methods can be proved significant for the crosschecking of the results and the tracing of any possible mistakes.Spec. Rep.-Transp. Res. Board, Nat. Acad. of Sciences, 247, 76-90, 1996. Wieczorek, G. F., Mandrone, G., and DeCola, L.: The influence of hillslope shape on debrisflowinitiation, in: Debrisflow Hazards Mitigation: Mechanics, Prediction, and Assessment, edited by: Chen, C. L., American Society of Civil Engineers, New York, 21-31, 1997.WP/WLI (International Geotechnical Societies' UNESCO Working Party on World Landslide Inventory): A suggested method for reporting landslides causes.Bulletin of the International Association of Engineering Geology, 50, 71-74, 1994. Wu, S., Shi, L., Wang, R., Tan, C., Hu, D., Mei, Y., and Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ; Guzzetti 75 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | cover the central part of the study area.Carbonates of the Trypalion nappe are exposed in the central-eastern part of the Chania Prefecture, extending from the north to the south, with a NE-SW direction.Limestones of the Plattenkalk zone are mainly exposed to the most southeastern part of the study area.The tectonic regime of the study area is characterized by faults in NW-SE and E-W directions.77 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | 3 Data layers-landslides casual factors

79
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

81
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper |

85
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | slopes of concave shape, (d) areas with elevation 200-400 m, (e) distance from rivers 20-100 m, (f) land uses of forest and heterogeneous agricultural areas, (g) the geological formations of phyllites-quartzites and flysch, (h) distance from faults 200-400 m, (i) precipitation of 800-1000 mm yr −1 and (j) areas proximal to roads (distance <100 m).Respectively, the most significant parameter classes with a negative impact on landslides are: (a) slope angle 0 slope aspect E , (c) the geological formation of Plattenkalk limestones, (d) elevation less than 200 m, (e) distance from drainage more than 100 m, (f) land use of permanent crops and scrub and /or herbaceous vegetation associations, (g) distance from faults greater than 400 m, (h) precipitation more than 1200 mm yr −1 and (i) distance from roads > 500 m.

Table 1 .
Xu, R.: Zonation of the landslide hazard in the forereservoir region of the Three Gorges Project on the Yangtze River, Eng.Geol.The causative factors and the corresponding weighting and rating values adopted in this study for the extraction of the Landslide Hazard Index (LHI) with the WLC method.

Table 2 .
The parameters obtained after the application of WoE method for the extraction of the Landslide Hazard Index (LHI).Italics are used for maximum and minimum C values i.e. for the classes of maximum positive and negative relation with landslide phenomena.Bold type show the major landslide causal factors based on weight values.