The Prevention of Landslides Using the Analytic Hierarchy Process (AHP) in a Geographic Information System (GIS) Environment in the Province of Larache, Morocco

Landslides are one of the natural hazards that many countries around the world are facing. In Morocco, the Rif regions are the most affected by these phenomena. Each year they cause enormous damage to the road network and infrastructure, especially in our study region, the province of Larache. The study region is subject to several opening up and road construction projects, which is why it is necessary to predict and identify the most vulnerable areas beforehand, in order to propose measures and techniques which are adequate for protection and reinforcement. The main goal of this study is to develop a susceptibility map to ground movements using a multi-criteria spatial assessment approach, and in order to reduce subjectivity, we have used a method for analyzing such complex decisions, which is the analytic hierarchy process (AHP) implemented in the geographic information system (GIS). Seven factors have been considered as conditioning factors in the occurrence of landslides, which are: lithology, fracturing density, slope, aspect, land use, density of the hydrographic network, and altitude. To verify the results obtained, we performed a correlation analysis of ground movements, already inventoried and verified in the field, with the susceptibility classes that were calculated. This analysis is accompanied by a statistical study.


Introduction
Mapping the susceptibility of slopes to ground instabilities is a very important component in assessing the risks associated with mass landslides. This was because zoning according to the degree of occurrence of land movements can help managers and decision-makers during development work, and reduce the stakes for individuals, the environment and infrastructure.
Since the early 1970s, many scientists have attempted to assess the hazards associated with land movements and have produced risk zoning maps describing the spatial distribution of ground instability phenomena, applying different methods based on geographic information systems (GIS) [1][2][3]. Several models have been proposed in order to produce occurrence maps for landslides. Among these methods is the heuristic approach, a multi-criteria spatial method which is applicable on a regional scale. Being a semi-quantitative approach, the weight of each layer depends on the expert's judgment, so that the more precise this judgment is, the more compatible the produced map is with reality [4]. In order to reduce the subjectivity of the decision makers, we used the analytic hierarchy process (AHP) approach when evaluating the weight of each predictive factor in the occurrence of field instabilities. AHP is a theory of measurement and decision-taking which was developed by Saaty [5]. AHP thus makes it possible to solve a large number of decision-making problems quantitatively by developing a decision support model, represented in the form of a hierarchy which allows the comparison of quantifiable and intangible criteria [6].
Using this method, in the case of the mapping of slope susceptibility to landslides, each data layer used is divided into sub-factors (classes), which are weighted according to their importance before these layers are assembled to produce the final map. This approach is based on three principles: decomposition, comparative judgment, and synthesis of priorities [7].
The objectives of this study are: (i) the assembly of predictive factors and their structuring, (ii) inventory and mapping of landslides, (iii) evaluation of the relative weight to each predictive factor by applying the AHP model and finally, (iv) the superposition of different layers for the development of the map of susceptibility to landslides in a GIS environment.
The resulting map is subsequently validated by comparing the location of the mapped landslides with the susceptibility classes.

Geostructural Framework of the Study Region
The province of Larache is located in the northern Rif, between the area of the flysch in the North and the outer area. It consists mainly of (i) intrarifan lands represented by all the units of Tangier, Loukkos and L'Habt of lower to upper Cretaceous age; (ii) Meso-Rif lands formed mainly by the unit of Numidian sandstones; (iii) pre-Rif lands represented by part of the Ouazzane Nappe of lower Cretaceous to Middle Miocene Age. All of its units are separated by large overlapping accidents oriented NW-SE with vergence NE [8,9].
The north-eastern part of the region is characterized by reliefs in the form of ridges and often sandstone hills, which mainly belong to the Numidian tablecloth namely Jbel Soukna (1603 m), Jbel Smia (853 m), and Jbel El Kobba (801 m). The north-west part of the region is characterized by the presence of two front-pit basins located in the front of the anticlines (Boujediane anticline and Arbaa Ayacha anticline) forming part of a system of folds and of overlaps called "anticlinal ramp of the Arbaa of Ayacha" by [10]. This system of folds and overlaps shows with the Numidian massive of Jebel El Kobba a submeridian direction (NNW-SSE). These two basins of before the pit were the seat of a turbiditic syn-tectonic sedimentation, the age of which extends from the Oligocene to the Lower Miocene [11]. In the south-west part, known as Bas Loukkos, includes the argillaceous alluvial plain, the sandy plateau of Rehamna and the hills of Oulad Ogbane.
The structuring of the region is part of the Alpine tectonic evolution which marks the Rif range, the dominant Meso-Cenozoic and Mio-Plio-Quaternary lands record the different structures which mark the major Alpine phase: folds, thrust layers, overlaps and fracturing [12,13] (Fig. 1).

Fig. 1. Geological and structural map of the Rif (Morocco)
Source: the structures are traced according to [8,9,12]

Methodology and Approach
The zoning of the landslides risk is defined as the classification of an almost homogeneous surface area according to the potential of the risk of gravity movements. There are several methods of zoning landslide hazards. The selection and optimization of an appropriate method according to the characteristics of the site and in relation to a possible scale is a very important step [14,15].
The approach used in this work aims at evaluating and mapping the susceptibility of slopes to landslides on a regional scale, throughout the territory of the province of Larache. The results obtained will then be validated using an interpolation of the inventoried ground movements in relation to the calculated susceptibility classes. The main stages of our work can be summarized as follows: carrying out an inventory of landslides, using aerial images, Google Earth, previous studies, and site visits.
The identification and mapping of the various factors involved in the process of ground instability (the predictor variables). These factors concern four criteria: geological, topographic, hydrographic, and environmental which have been detailed in seven sub-criteria: lithology, fracturing density, slope, aspect, altitude, density of the hydrographic network, and the land use. This step is the most difficult and the longest, due to the scarcity of data and similar studies in the region. The stage was completed by the creation of a codified, structured and interactive database of raster maps of the factors involved in the process of ground instability, The estimation of the weight of each factor used a statistical method, the analytic hierarchy process (AHP), which is based on the comparison of the different conditioning factors, two by two, from the construction of a square matrix. The weight is assessed according to the relative importance of one factor over another, using an appropriate scale for this [16]. It is at this level that the expert must intervene, by referring to a deep knowledge of the ground truth, since he must assign a relative value (between 1 and 9) reflecting the importance of each factor in relation to the other. Once the comparison matrix is filled, on the "Set weights" function of the AHP application installed on ARCmap, a weight specific to each predictive factor is calculated by executing the "Compute" command.
Development of the susceptibility map to landslides: the corresponding posteriori probability map is generated by exercising the "Create map" instruction on the AHP tools. The resulting raster is only an expression of the hierarchy or the order of precedence calculated for the considered variables. The consistency of the decision process can be checked by the values of Consistency Index (CI) and the Consistency Ratio (RC), if RC is less than 10%, the matrix is acceptable [17]. This result can finally be judged and validated by the decision makers through a simple interpolation between the map obtained and the inventory of landslides which had occurred in the past.
The superposition of the factor maps by associating them with their relative weights, allowed us to establish a susceptibility map for landslides.
In our case, only the spatial component of the hazard was taken into account but it is possible to incorporate a temporal component in this analysis, in particular when a link is established between certain hydro-meteorological, seismological or other triggering parameters and the triggering of ground movements (the province of Larache is located in the northern Rif, between the area of the flysch in the North and the outer area) (Fig. 2).

Landslide Inventory
The observation of the ground supplemented by the analysis of the cartographic supports (geological and geomorphological maps, satellite images), allowed us to distinguish, according to the lithological nature of the materials put in motion, the shape of the shear surface also the magnitude of landslide, four major types of land movement: muddy flows, landslides, rock collapse, and complex landslides. The mapping of mass movements has systematically included the active forms where the risk is certain and the old forms stabilized or in a period of temporary remission. Figure 3 presents the centroids of mass movements that we used in the validation phase of the model.

The Conditioning Factors of Ground Instabilities
The causative factors of ground instabilities are classified into two categories, namely predictors and triggers. The predictor variables that we have taken into consideration are: lithology, slope angle, aspect, altitude, land cover, fracture density, and drainage density. The triggers or driving forces are: earthquakes, precipitation, and human intervention. The triggers ware only added if we are interested in the time component and since only the spatial component is considered in our study, only the predictor variables were mapped.
These variables were subsequently structured in a "Grid" format and then reclassified, after having developed a correlation between the different classes of each parameter and the presence of landslides in the present or in the past. Based on these statistical results, we assigned indices between 1 and 5 depending on the occurrence of each entity with respect to landslides (5 -very high occurrence, 4 -high occurrence, 3 -medium occurrence, 2 -low occurrence, 1 -very low occurrence). This facilitated the decision-making phase, during the binary comparison and the execution of the weight calculation matrix (Fig. 4).
Lithology. The lithology map was developed using existing geological maps (the Rif map at 1:500,000, Larache map at 1:50,000) coupled with a series of photoaerials, and photo-interpretation computer-assisted from optical satellite data (Landsat 7 ETM image), with field checks. In our work, the realization of the lithological map was carried out in two stages. First, we proceeded to identify the geologic units that have lithologic affinity. Subsequently, we carried out a requalification of these geological formations according to their resemblance in mechanical and geotechnical characteristics. We were able to distinguish four classes of lithological facies, namely marls, sandstones, shale, and alluvium. According to the occurrence of each lithological entity, and with respect to gravity movements, we attributed a possible index to each class, this being based on a statistical study of the interpolation of the points of instability mapped on the different lithological units. Then the map was reclassified into four classes to which we assigned an index (1,3,4,5). Fracturing. The inherited and recent tectonic activity can intervene in the activation of ground movements in weakened areas [18], either by pre-conditioning the event or as a trigger during a seismic crisis. The occurrence of mass movements is usually accentuated near resurgences. Indeed, water infiltrates at the level of the faults cause increases in pore pressure, which reduces the shear resistance of the Msoils. Once saturated, the mass movement is initiated with the formation of a deep shear plane.
The production of the index map linked to fracturing is based on a multisource study, according to the following steps: 1. Digitization of lineaments from analysis of Landsat 7 ETM+ satellite imagery (by combining directional filter techniques). 2. Extraction of the recent fracturing that affected the region based on the morpho-structural analysis of the hypsometric map coupled with those of the anomalies of the hydrographic network. 3. Digitization of all bibliographic flaws, in particular those coming from the geological map of the Rif range 1:500,000 [9,12,19]. 4. Correction of the lineament map, from morpho-structural analysis, previous work, and the ground knowledge. 5. Creation of the fracturing index map, in terms of density. To do this, we used the Spatial Analyst Tools module on ArcGis 10.x, using the "Line Density" function. It makes it possible to calculate the linear density (cumulative length) of any linear element per unit area (km²). 6. To each fracturing density class we assigned an index, depending on its occurrence with respect to landslides. The probability of the occurrence of landslides increases as the fracturing density increases, or also as the distance from the fault decreases.
Land use. The plant cover primarily exerts a phyto-stabilizing action on slopes with clay or rocky substrate developing a thin clay weathering cover. Anchoring by deep rooting can also limit the development of landslides with poor geotechnical characteristics (clayey, clayey-marly, etc.). However, in some cases, rooting can promote subsurface runoff (root mat), and causes the triggering of slope movements. As well as the absence of plant, or its scarcity on the slopes, accelerates the erosive effect by runoff, which is a motivating factor in this case.
The relationship between the spatial distribution of the different landslides, in the province of Larache, and the land use classes is supposed to be linked to human activity which manifests itself by destructive actions. The sum of these are degradation of the plant cover and the effects of water erosion, the loading of the head of the slopes by the construction works and the embankments, the abutment at the foot of the slopes by excavation and quarries, earthworks, filtering zones caused by the installation of septic tanks, and pipe leaks, modification of surface runoff by construction work and the installation of road layouts... etc., despite the planting trees actions carried out by the Department of Water and Forests and the Fight against Desertification, which are still humble.
However, integrating human activity into this modeling is difficult. Therefore, in this section we assess the role of plant cover in slope stability through land use mapping in the province of Larache.
The classification supervised by maximum likelihood of the Landsat 7 ETM+ optical data, allowed us to establish a land use map (Fig. 7) to which we attributed instability indices relating to the occurrence of occupancy classes of the ground in relation to the landslides. Four classes are distinguished: forests, matorrals, agricultural lands, bare lands and ponds and lakes.
Drainage density. In general, the potential for land instability increases near rivers, because watercourses can affect the stability of slopes, either through their erosive effect or by increasing the water level in the lower parts, which leads to saturation of the materials and an exaggeration of the pore pressure [20]. In our study region, the proximity of watercourses is one of the factors that influences the occurrence of landslides. A density map of drainage was created by interpolation using the "Line Density" function, in Spatial Analyst Tools in ArcGis 10.x.
Slope angle. Plays an important role in the susceptibility of slopes to landslides [21], since gravity movements are directly related to the degree of slope. This parameter is frequently used in the preparation of landslide susceptibility maps [17]. According to the correlation between the landslides, the slope map was reclassified into five classes with indices (1,2,3,4,5). However, the occurrence of the "slope angle" with respect to landslides, at the level of the studied area, is less important than the lithology because the geological units are of mediocre geotechnical quality which can generate landslides even on low slopes.
Aspect. Parameters related to factors such as the exposure of slopes to the sun, drying winds, precipitation (humidity or degree of saturation) and discontinuities can control the frequency of landslides, although some authors [22] have mentioned that exposure does not have a significant influence on slope instability. However, several researchers have shown the link between the orientation of slopes and the frequency of ground movements. In this study, the exposure was subdivided into five intervals (NNE, E, S, SW, NW). Classes with a strong occurrence are NW and NNE.
Altitude. The relationship between landslide activity and altitude is very debatable. Some researchers use elevation as a control parameter for ground instabilities [23]; other researchers find that the activity of slope movements, in a specific basin, is greater at certain altitudes [24]. While in most parts of the Northern Hemisphere the relationship between altitude and precipitation are direct, on the other hand changes in altitude have a close relationship with air temperature. As a result, changes in these two parameters (rainfall and temperature) influence the quantity of water storage, the type and density of plant cover as well as evapo-transpiration, etc. According to the correlation study with the mapped instabilities, the high occurrence of mass movements corresponds to the 150-500 m class.
All the data layers are classified into five classes according to their occurrence in landslides. Table 1 summarizes the obtained results.

Foundation of the Matrix of Binary Comparisons and Calculation of the Weight
The analytic hierarchy process (AHP), implemented on ARCmap, calculates the weighting for each criterion (W i ) by taking the eigenvector corresponding to the largest eigenvalue of the matrix, then by normalizing the sum of the components to the unit as follows [25]: where n -the number of criteria.
On the basis of the studies carried out on the conditioning factors, we carried out a comparison, two by two, of the occurrence of the conditioning parameters with respect to the ground instabilities. Each factor is evaluated with respect to any other factor by assigning a relative value between 1 and 9. When a factor i on the column is greater than the factor j on the row, this value varies between 1 and 9. Conversely, the value varies between 1/2 and 1/9.
We performed several trials of binary comparisons for the seven highlighted predictors. For the seven considered parameters we have a 7 × 7 comparison matrix, similar to the pair wise comparison matrices of a symmetrical nature. In this regard, only 25 values are needed to fill in, which corresponds to the upper diagonal and Table 1. cont triangular half of the matrix. The diagonal squares of a pair wise comparison matrix always take a certain value from 1 to 9. The squares of the upper and lower halves are symmetrical. Therefore, the corresponding values are reciprocal. Once the matrix is built, the AHP program generates the weights of the thematic layers of all the causal factors classified on the basis of the weights of the entered classes and on the basis of the binary made comparisons.
The results consist of the weights factor and a calculated consistency ratio (CR). Tests with a CR greater than 0.1 were automatically rejected, while a CR less than 0.1 was often acceptable.
In this study, we kept the result with the lowest CR, which corresponds to 0.04. This ratio indicates a reasonable level of consistency in the pair wise comparison for our model, which is good enough to recognize the weights factor. Therefore, the weight corresponding to the lithology is the highest, while the altitude is the lowest.

Creation of the Susceptibility Map to Landslides
The general principle of the multi-criteria heuristic approach is based on the combination of the relative weighting of several sets of conditioning variables.
For the mapping of the susceptibility to landslides, two weighting levels were taken into consideration: 1. In the first, a relative weight was attributed to the classes of variables. It corresponds to their estimated importance in the process of landslides. 2. On the basis of the same principle, the predictors were weighted by adopting the pair wise comparison.
Once the criteria scores were normalized and the weights were calculated, we were able to generate the susceptibility index map, using the WLC (Weighted Linear Combination) method [26]. This is the simplest method which combines the criteria to form a single evaluation score. In the WLC method, each criterion is multiplied by its weight from the pair wise comparison. The calculated weights are integrated into a linear sum, according to the following equation [26,27]: where: S -the final susceptibility index map, W j -the weight assigned to the variable j, X ij -the weight assigned to the class i of the variable j.
Weights can have a considerable influence on the result. Since the weights of the factors are added one by one, the final scores of the resulting map are expressed on the same scale. In addition, the weights assigned to each criterion determine the level of compromise with respect to the other explanatory factors.
The execution of the command "Create map" involves the application of the formula (2) [15]. The index maps of the posteriori probability derived from the AHP approach comprise five classes of susceptibility: very strong, strong, medium, low, very low (Fig. 5). The resulting landslide hazard zoning map indicates that high to very high-risk areas cover approximately 40.99% (1118.6 km 2 ) of the total area while approximately 21.83% (583.22 km 2 ) were classified as being at moderate risk and 37.64% of the study area (1027.32 km 2 ) is classified as barely susceptible or not.

Validation of Results
The interpolation, points of instabilities in the field on the susceptibility classes, shows that the classes with "high susceptibility" and "very high susceptibility" coincide with the greatest number of movements already inventoried in the region, especially since these two classes extend over very small areas compared to other classes. The histograms in Figure 6 present the statistical results of the validation of the resulting maps. The obtained results are considered satisfactory for the validation of the susceptibility maps of developed mass movements and the AHP model used. These landslides susceptibility maps provide logical and relevant information on current instabilities and also on possible mass movements that may arise in the future, so these results appear to be in perfect harmony with the morphological evolution of the region.

Conclusion
In the present study, we used the AHP analysis method in a GIS environment to produce a landslide susceptibility map at the level of the province of Larache (regional scale). The explanatory parameters used in the production of this map were carefully selected by analyzing the relationship between the parameters and the genesis of landslides in the region. A database is produced and structured by digitizing ground data, satellite images, DEMs, topographic and geological maps, taking advantage of the various functionalities that ArcGis10.X offers us. The weighting of the affecting factors and landslides were obtained by means of the AHP method. The effect of the classes of the data layer, as well as the weight associated with each factor, are quantitatively determined by the AHP method.
The obtained results in this modeling have demonstrated that the use of the AHP method can be considered among the practical and quite realistic approaches to define the weights of the factors in the case of evaluation of the susceptibility to ground instabilities on a scale of regional order.
These results revealed that landslides in the province of Larache are strongly influenced by lithology, degree of slope and land use. The most susceptible areas to landslide are located in the parts belonging to the Tangier and Loukkos units. The produced map will serve as a basic document in order to plan appropriate risk mitigation measures and future development strategies based on the identification of areas prone to slope movements in the region.
The final landslide susceptibility map can be used in possible land use planning and the prevention of landslides.