Mass movements susceptibility mapping by using heuristic approach. Case study: province of Tétouan (North of Morocco)

Abstract

This research paper aims to model Mass Movements Susceptibility (MMS) in the province of Tétouan. First, we identified the characteristics and spatial mapping of the different types of MM (collapse, mudflows, and complex landslides) by means of the interpretation of satellite images and from fieldwork. Subsequently, we selected the predictive parameters controlling the occurrence of MM e.g. lithology, land use, fault density, hydrographic network density, slope degrees, slope aspects, and elevation. We used the heuristic method for Modeling Mass Movements Susceptibility (MMMS). The choice of this method compared to other methods (fractal, factorial, and neurons) is justified by the possibilities of intervention and the judgment of the expert who relies on the ground truth to select the parameters, to identify the classes, and to assign the weights to each one; unlike to other methods with steps that are done automatically and randomly. The results of the validation of the susceptibility map correspond to 70% compared to the field data and it includes five susceptibility classes (not susceptible, low, moderate, high, and very high). Indeed, the originality of this paper relies on the fact that the creation of our susceptibility map will eventually indicate the areas of roads, dwellings, the extension of urbanization, and dams, which are located in areas at risk of MM. Our map is also a powerful decision-making tool to conduct management plans and to guide the selection of sites to build new projects; which help mitigate the socio-economic impacts usually encountered when mass movements in Tétouan province are triggered.

Introduction

Mass Movement (MM) is a phenomenon of geological, geophysical, and/or seismic origin, where a mass of terrain descends on a slope (Cruden 1991). The term mass movement is typically used to describe a wide variety of processes that were translated by gravitational movements of soil, rock, and vegetation (USGS, 2004)¹. According to Varnes (1984) a mass movement is defined as: “The probability of a mass movement in a specific time and in a given region of a potentially harmful phenomenon”. Several studies have been conducted worldwide (Chalkias et al. 2014; Del Ventisette et al. 2014; Shahabi and Hashim 2015; Bozzano et al. 2020) and in the Northern Rif to identify, classify, and map mass movements (Maurer 1968; El Gharbaoui 1981; Millies-Lacroix 1968; El Kharim 2012; Mastere et al. 2015; Ait Brahim et al. 2018). Some of these mass movements can be caused by the nature (Okeke and Wang 2016) or induced by man-made activities. As most of them took place in populated areas, the results can be tragic, causing loss of life and economic damage. For this reason, and as Carrara et al. (1995) said: “The past and present mass movements are the key to the future”, it is mandatory to delineate areas prone to mass sliding in order to prevent and better manage these catastrophic events before they take place in the future.

The main goal of this research paper is to develop Mass Movements Susceptibility Mapping (MMSM) by using ground truth knowledge, remote sensing data, GIS-based techniques, and the application of heuristic approach. New technologies like GIS have made the impossible possible and facilitate mapping processes with a minimum of data (Brabb 1984). During the two last decades, many authors such as Van Westen et al. (2003); Guzzetti et al. (2005); Anbalagan et al. (2015); Roy and Saha (2019) have developed various statistical methods using GIS techniques. Recently, authors such as Lin et al. (2017) have grouped the most common methods used for MMS into three categories: statistical, physical, and deterministic. The statistical methods are processed in an automatic manner without a previous study of the relationship or correspondence with the ground truth and regarding the physical methods, provide precise results; however, they are expensive, especially for large scale. On the contrary, deterministic approach combines both the computing to run the model that are easily implemented, compared to the physical ones and require the expert definition of the weighting mass movements factors, which is not a dilemma in itself, since we are familiar with the study area. In fact, the application of an effective model that relies heavily on ground truth, expert knowledge, satellite imagery, and a qualitative approach to mapping MMS is a crucial task for further study of vulnerable zones in our area. Hence, this approach could eventually apply to other countries in the surrounding region of North Africa that are subject to mass movements and have limited resources for large investigations. The approach that is considered in this research paper is more robust than:

• 1 Qualitative methods where the creation of Mass Movements Susceptibility Map (MMSM) is carried out either solely by simple zoning and geomorphological interpretation in the field, from satellite images, or by interpretation of aerial photos (Guzzetti et al. 2006), without a mandatory inventory of mass movements (MM) from field surveys. Otherwise, by assigning relative weight to the parameters selected (Anbalagan 1992; Nagarajan et al. 2000; Saha et al. 2002) on the basis of in-depth knowledge of the field and the expert experience involved in this work (geologist and/or geomorphologist). The weight assigned to each parameter is related to its relevance in the genesis of Mass Movements (MM).

• 2 Quantitative methods that are based on the useability of probabilistic and statistical empirical models, which calculate a different weight as a function of the number of mass movements that we have mapped within the different classes identified for each controlling parameter (Dai et al. 2001; Lee and Min 2001; Lee et al. 2004).

We would like to overemphasis that our method is more efficient than the previous ones aforementioned since the formulas, which we used (Voogd 1982) combine the weight of the parameters determined by qualitative approach and the weights of the different classes of each controlling parameter by quantitative approach. In addition, the susceptibility map obtained is crossed-validated with half of the mass movements that we previously kept for validation purposes. Finally, the ROC curves were established and the areas under the AUC (= 0.79) were calculated to assess the degree of fit of the model to choose the best reclassification of Mass Movements Susceptibility Map (MMSM). We would like to draw the attention of the readers that the calculation has been performed, but not presented in this article.

Study area

The province of Tétouan is located at the far east of the Tangier peninsula (north of the Rif chain of Morocco). The area of our zone extends over approximately 2137 Km², between the parallels 35^∘ 29" 58.781’ and 35^∘ 55" 18.326’ North and the meridians 5^∘ 9" 44.071’ and 5^∘ 29" 55,667’ West (Fig. 1). On May 5-7, 2021, we conducted a field observation to verify our results and to note fresh/new mass movements in the study area (Fig. 2)

Fig. 1

Location of the study area and the distribution of mass movements

Fig. 2

a,b,c,d,e,f Examples of mass movements along the road networks of the province of Tétouan: a road N2, Tétouan (35^∘33’25.3"N 5^∘24’46.3"W), b road N16, Fnideq (35^∘50’08.7”N 5^∘21’19.0’W), c raod N16 (35^∘31’58.2”N 5^∘11’52.3”W), d mediterranean rocade, Wadi Laou (35^∘27’30.3”N 5^∘07’05.9”W), e road N2 (35^∘21’53.2”N 5^∘22’14.4”W),and f road P4702 (35^∘25’11.7”N 5^∘26’58.7”W). g mass movements in the middle of urban settlement from road N16 (35^∘34’42.2”N 5^∘18’13.6”W). h sliding and failures related to erosion along the highway A6 (35^∘46’24.3”N, 5^∘21’44.1”W)

Climate and precipitation

Climate

The study region benefits from two maritime influences: on one hand, a Mediterranean influence that is introduced into the chain by opening the limestone chain at the Tétouan Cluse; and on the other hand, an influence of the Atlantic Ocean, because it is exposed to the winds of the Western Atlantic and the Southwest. It is characterized by the succession, during the year, of two distinct seasons, a hot, dry summer that lasts several months, followed by a rainy, colder season corresponding to autumn, winter, and early spring. The rains are not permanent throughout this part of the year; but have alternating periods of rain and fine weather with a maximum precipitation (80%) of November and April. They are very violent with a maximum threshold in November, December, and January; consequently, they are causing a significant erosion on the slopes and sudden flooding.

Precipitation

In recent decades, torrential and irregular rainfall have been recorded in the Tétouan region, which have reached 150 mm/day, with an average frequency of 5 events per year (Smir Dam Stations, Torreta Bridge). Luino (2005) reported that in section 4. Extraordinary hydrological events of his paper that: “Since the end of 1960s, observation of the behavior of northwestern Italian basins during extraordinary hydrological events has shown that the number and type of instability processes triggered by rainfall are not only related to the morphological and geological characteristics of the area where the rain falls, but also to the distribution of intense rainfall during the meteorological event.” That is why we paraphrase it to this sentence: “...quantity and duration of instability during the rainy period”. He explained that when precipitation exceeds a certain threshold (which may vary depending on hydrological conditions), the instability activities begins to trigger reliably. In Tétouan province, between December 2009 and March 2010, the region received nearly 1,100 mm of precipitation (80% of annual precipitation), which was the highest cumulative total recorded in the past decade. The combination of that type of irregular but intense precipitation with low-strength geological formations (clays and marls) and vegetation cover heavily degraded locally by anthropogenic action (agriculture, livestock, urbanization, forest fires, etc.), induces after a much faster response time of mass movements (PATEAU 2014).

Geology

The province of Tétouan is part of the Tangier Peninsula, which forms the Western end of the Rif Chain, which are themselves part of the Enduring Alpine Ranges, grouped under the name of the Betico-Rifo-Tellian Arc (Delga et al. 1962). At the province of Tétouan level, the main structural areas are: internal domain, Flyschs nappes domain, and external domain (Suter 1980).

• The internal domain: it is mainly developed in the Tangier Peninsula. The Inner Domain is composed of three structural layers stacked on top of each other: the Sebtides, the Ghomarides, and the Dorsale Calcaire.

• Flyschs Nappes Domain: it is a set of units made up of allochton Flyschs series that encircle the Calcaire Dorsale and the Predorsalian zone outwards. The Flyschs nappes are mainly made up of lands of age ranging from the Lower Jurassic to the Oligocene (Delga et al. 1962). It includes: the Predorsalian nappes, the Mauretanian nappes, the Massylian nappes, the Numidian nappes.

• The external domain: it is a vast structural area largely carried by the Flyschs nappes and the internal domain. It is subdivided into three material sets ranging from the Triassic to the Miocene (Suter 1980), which are, from the inside to the outside of the chain: the IntraRif unit of Ketama-Tanger (IntraRif), the MezoRifaine zone of Senhaja-Aknoul (MesoRif) and the pre-Riman zone (PreRif).

Methodology

In the past two decades, many authors and scientists have used the heuristic approach to mapping Mass Movements Susceptibility (MMS) at large scale (Carabella et al. 2019). The heuristic approach is based on the judgment and knowledge of the expert considering all the various specificities of the case study. As has been stated by many researchers (Anbalagan and Singh 1996), the heuristic approach is a qualitative analysis of mass movement events that takes into account geomorphological and geological characteristics. This method is closely related to the expert’s experience and is the most useful and effective approach to mass movements caused by different factors. We chose this method for many reasons: 1) it is applicable on a large scale (Tétouan province), 2) it is suitable for the use of remote sensing data, 3) it depends heavily on the knowledge of the expert in terms of geological and morphological fields of the study area, and 4) it allows cross-validation with previous data on mass movements. The MMS depends mainly on the rigor with which environmental parameters are selected as input to run the model in order to generate MMS map. This method is divided into two processes: 1- direct mapping analysis and combination of qualitative maps; 2- assigning weights to all conditioning factors and their classes according to their field importance followed by cross-validation with half part of mass movements that were not used to train the model. In the first process, researchers determine susceptibility directly in the field that is based on the expert’s experience. In the second process, the expert uses his field knowledge to specify the different weight values of each class. For the assignment of weights to the 7 selected parameters and the 5 classes determined for each parameter, we combined the qualitative method and the quantitative one according to the following protocol in order to reduce the subjectivity of the method. The methodology adopted in this research paper is shown by the Fig. 3.

• 1 for each class (W_j) of the 7 controlling parameters selected, we used the quantitative method essentially based on the statistical study that considers the number of mass movements (MM) that we mapped by surface within each class (Dai et al. 2001; Lee and Min 2001; Lee et al. 2004). This number is recalculated based upon the frequency of MM existed in the study area. Afterwards, this latter is translated into the corresponding weight of each class. The values of each class range from 1 (Not susceptible) to 5 (highly susceptible).

• 2 for each parameter (W_j), we used the qualitative method (Anbalagan 1992; Nagarajan et al. 2000; Saha et al. 2002) by assigning a relative weight to the seven selected parameters (lithology, density of faults, density of the hydrographic network, landuse, slope degrees, slope aspects, and elevation) on the basis of our in-depth knowledge and extensive field experience. In addition to three subject matters expert analysts involved in this tremendous task (geologist and geomorphologist). The weight assigned to each parameter is related to its relative importance in the genesis of Mass Movements (MM). The values range from 1 (less important) to 5 (very important).

Fig. 3

Flowchart adopted to create the susceptibility map of MM by heuristic approach in Tétouan province

The main criteria used to assign the weights to the seven predisposition parameters selected come from the specific ground truth of our study area related to precipitation (heavy precipitation 150mm/day) which constitute the main dynamic factor triggering mass movement (MM). Hence, lithology (marlstone and mudstone are very indurated, and flysches moderately resistant), and land use (vegetation is highly degraded) are considered the major factors influencing the occurrence of MM in the province of Tétouan. They will be represented by the most important relative weight of value 5. The slope comes second with a weight of value 4 because the rugged and steep terrain indicates that any increase in the slope angle results in an increase in the tangential shear stress and consequently a decrease in the security rate that accelerates the MM. Afterwards, it comes the fault density by its role of water infiltration hence the attribution of the value 3. Then, the parameters of elevation and the hydrographic network density, which is impacted by the lateral undermining of torrents, or that of the drainage; and therefore, it serves to evacuate the earth debris. These two parameters are also important in the control of MM, but they come at second position compared to the parameters; and consequently, they receive a weight value of 2. To end, slope degrees and slope aspects are assigned a value of 1 as a relative weight (table). As a matter of fact, the spatial subdivision of each parameter into 5 classes aims, during the calculation process, to take into account only the weight assigned to each parameter class on the basis of real statistical data (approach quantitative) and the area (spatial distribution) of each parameter class. This homogenization of the number of classes, will make it possible to obtain a well-prioritized susceptibility level with very contrasting areas thanks to the cumulative effect of the weights assigned on one hand to the 5 classes of each parameter (qualitative approach) and on the other hand to the weights attributed to the 7 different predisposition parameters selected (qualitative approach).

• 1 Once the different weights are assigned, we apply either the mathematical formula proposed by Francipane et al. (2014), which takes into account the indexing of the weights of all the controlling parameters selected for the creation of the susceptibility map of the province of Tétouan :

  $$S = \left(I_{cop} * I_{prend}\right)*{\frac{\sum_{n}\left(I_{n} * P_{n}\right)}{\sum_{n}P_{n}}} $$

  Where S is the Susceptibility Index, I_cop and I_pend are the indexes of the discriminating parameters of coverage and slope, and respectively, I_n and P_n are the index and the weight of the n-th inducing parameter.

• 2 Or the method proposed by Voogd (1982), of which all the parameters and indexed classes are gathered in a single Index Map of the Mass Movements (IA_MT),and this using the formula that allows the calculation of the weighted linear summations:

  $${IA}_{MT} = \sum_{j=1}^{n}\left(W_{j}*w_{ij}\right) $$

  Where IA_MT: Mass movements susceptibility index; W_j: Value of the weight assigned to parameter j; w_ij: Value of the weight assigned to the class i parameter j; n: Number of parameters.

Data and materials

In this research article, a heuristic approach by developing a scoring scheme, which is based on the importance of each predisposing factor in triggering MM in the study area, was employed. The methodology requires selecting factors by using external sources such as satellites imagery, completeness and double checking with existing maps (Geological and Topographical maps), analysis, geo-processing, and interpretation by using Erdas IMAGINE 2014. The next step consists of numeric data ranking and weighting, integration of data into GIS environment by using ArcGiS 10.2 software and then generation of thematic maps of all the predisposing factors (lithology, landuse, fault density, hydrographic network density, slope degrees, slope aspects, and elevation). Next, the transformation into the grid format and the classification into five classes of all factors. Finally, the MMS is validated in regard to the existing distribution of field movements and a statistical significance test. Figure 2 shows the flowchart adopted for the elaboration of MMS by using heuristic approach in our study area.

Mass movements

The characteristics and spatial mapping of different types of MM in Tétouan province have been identified from the work of Millies-Lacroix (1968); Kirat (1993); Sitteri (1999); Ait Brahim et al. (2002); El Kharim (2012); Mansour and Ait Brahim (2005); Fonseca (2014); PROkOS et al. (2016); Elmoulat and Ait Brahim (2018); Ait Brahim et al. (2018), and completed by our analysis of high-resolution satellite images and our fieldworks. We gathered a catalog of 228 MM; these in particular were recently active, but the old ones were hard to detect. The different types of these movements encountered in our study are collapse, mud flows, and complex mass movements covering an area of 249.27 km², about 8.57% of the total area of the province (2137 km²). The map of the mass movements treated in our study area is represented by the Polygon-shaped Fig. 4-a to reflect the ground truth and by Fig. 4-b in form of centroids to be able to perform geo-processing in ArcMap 10.2. For further geoprocessing within the GIS environment, each single mass movement polygon was converted to its centroid to change its geometry into point. The collapses and complexes landslides were developed at the level of solid materials (limestone and dolomite) of the internal Dorsal domain in steep slopes. The mud flows are located mainly in the middle hills, with moderate slopes of Flyschs (predominantly mudstones). The MM appear at the low and open hills of the Senonian nappes of marls of the Tangier unit. The complex MM were resulted from lithological and structural heterogeneity and steady elevation (moderate to high slopes associated with Flysch nappes). Several of these MM were keenly studied by our predecessors in Rif areas such as Kirat (1993); Sitteri (1999); El Kharim (2012); Ait Brahim et al. (2002); Mansour and Ait Brahim (2005); Fonseca (2014); Elmoulat et al. (2015); PROkOS et al. (2016); Ait Brahim et al. (2018); Elmoulat and Ait Brahim (2018). Among these studies, the most important one in the province of Tétouan has been examined by El kharim (El Kharim 2012). The development of several MM on both slopes of Martil Wadi (G6, Fig. 4) were accentuated by anthropogenic factors in the city of Tétouan and its surroundings. The G5 mass instabilities are recorded in the Onsar region in the north of Tétouan (G5, Fig. 4). These movements of specific-character slopes present the concomitant of plastic-type processes with collapses. Thus, movements, flows and solifluxions, affecting the marls and Flyschs are associated with collapses by means of boulder falls or mass movements of above carbonate materials. The M’diq region had experienced complex landslides and collapses (G8, Fig. 4). As an example of these disastrous events we list : 1) mudslide and boulders on the eastern slopes of Zem Zem mountain (G9, Fig. 5), 2) mass movements in the North of Fnidaq (G10, Fig. 5) and south of Tétouan, 3) Douar Hayounech on the left bank of the oued el Kbir, 4) Tamezakht movements from 1990 and 1991 north of Douar Tamezakht located at the bottom of the Calcaire Dorsale on a slope with a slightly low slope (G1, Fig. 5), 5) the unstable slopes of the Oued Mhajrate (G3 and G4, Fig. 5), and 6) mass movements on the predominantly sandstone or limestone and mudstones of Koudiat Krikra straddle between the Nakhla and Martil dams (G2, Fig. 5).

Fig. 4

Mass movements of Tétouan province in polygon-shaped (a) and in form of centroids (b)

Fig. 5

Some examples of mass movements (G) in Tétouan province: The Tamezakht Movement (G1). The movements of Koudiat Krikra near the Nakhla and Martil dams (G2; G3; G4). The Slope Movements in the El Onsar Pass (G5). The mass movements on the 2 slopes of Oued Martil (G6) with Tétouan and its surroundings (Kourrat Sbaa, Khandak Ezzerbouh Samsa, Bouanan, Kaa El Hafa etc)

Predisposing factors

Seven parameters were used in mapping MMS, i.e. lithology, landuse, fault density, hydrographic network density, slope degree, slope aspects and elevation.

According to many authors who have worked in the field of mass movements in the Rif and based on our own work for over than 30 years, we were able to identify geological, geomorphological, geotechnical, seismic, climatic, anthropogenic parameters in relation to mass movements. This research work has been limited only to the most relevant controlling factors, permanent and conditioning of the occurrence of mass movements namely: lithology, density of fault, density of hydrographic network, land use, slope degrees, slope aspects and elevation. Precipitation, earthquakes, anthropogenic actions have been ruled out because they are considered as causative.

These factors have been carefully selected for their importance (as predisposing factors) of mass movements in this region. All thematic maps have been converted to Grid format. This format is preferably used by many authors as Ruff and Czurda (2008) because, it has some advantages in the handling and processing of data compared to vectors and other types of formats. The resolution of the grid is 27m.

Lithology

The lithological nature of the soil is a major predisposition factor in the mapping of susceptibility to gravitational movements of the terrain, both by the nature of the materials (geotechnical quality of the rocks), and by the structure (slope and dip of geological layers). Mass movements in the province of Tétouan are related to the petrographic nature of the area (characterized by the presence of limestone, sandstone, clays, marls) and to the infiltration of water mainly at the level of clays and marls sensitive to changes in water content. The map of the lithology was made from the analysis and interpretation of the Landsat ETM-7 images that we performed on a textural basis. The geological map of the study area (Durand-Delga and Kornprost 1985), combined with our knowledge of the field, allowed us to validate the results from remote sensing techniques. All digitization activities were completed using ESRI ArcGIS 10.2. Thus, five lithological units could be distinguished (Fig. 6-a).

Fig. 6

Seven digital thematic maps of the main selected predisposition parameters of Tétouan Province

Land use

Remote sensing gives us the opportunity to access to mountainous areas of the province of Tétouan where extensive fieldwork is extremely difficult. Indeed, for these maps to be effective and efficient, it is necessary to have robust and reliable automatic methods capable of using the availability of data. Broadly speaking, automatic approaches are useful for producing land use maps from satellite images that are often based on image classification methods. For this reason, the methodology adopted in this study is the supervised classification. This classification based upon the use of samples for which land use is known as examples for learning process. Supervised classification usually provides better results, but it requires baseline data for learning that are expensive to obtain (field campaigns, photo-interpretation, ground truth, multi-source database, etc.). All supervised classification and geo-processing procedures were accomplished using ERDAS IMAGINE 2014. The results of the supervised classification were integrated into the GIS environment for grouping and vectorization (Fig. 6-b).

Fault density

According to previous works by many researchers and experts in the Western Rif such as Deffontaines et al. (1992); Ait Brahim (2003); Ait Brahim and Alaoui (2003); Mansour and Ait Brahim (2005) our knowledge of the terrain and the geological work carried out by our predecessors allowed us to retain only the faults present in the study zone that could play an instrumental role in the instability of the field (crushed land, water infiltration, emergence of springs, etc.). To highlight the fracture network at the provincial level, we used remote sensing techniques. The processing of Landsat-7 ETM+ satellite images generated a lineament map that was subsequently validated on the basis of existing geological maps and our knowledge of the field (Fig. 4-c). Each scene consists of seven 30 m resolution spectral bands and a 15m-resolution Panchromatic Band. Digital image processing was conducted using Erdas IMAGINE 2014 and ArcGIS 10.2 software. More than 421 types of tectonic movements were mapped and analyzed using Landsat ETM+7. The faults in our study area are represented on ETM+7 images by approximately straight structural lineaments whose texture, hues, and variations of lithological highlights are different. The created lineament map was subsequently validated based on the existing geological maps of Ceuta and Tétouan-Ras-Mazari on a scale of 1/50,000 and our knowledge of the field (Fig. 6-c). At the end, analysis of these structures allowed us to distinguish four broad categories of faults that are N-S, NE-SO, E-O and NO-SE.

Hydrographic network density

Numerous works by several authors (Santangelo et al. 2013); have shown that the water system is one of the fundamental parameters that controls the triggering of gravitational movements of the land. Thanks to its mobility and erosive power, it generates instability on the slopes resulting in a long-term decrease in its resistance. This is how the removal of stop by erosion, causes the movements of the slopes or their reactivation. In our analysis, we integrate only surface water and not groundwater, since it relates to other types of untreated movements in this work. Indeed, in this study the river network parameter was automatically extracted from the Digital Elevation Model (DEM) using ArcHydro tool. This extension was developed by ESRI with a document that explains step-by-step how to analyze drainage. The DEM is used to generate data on flow direction, flow accumulation, and watershed delimitation. This data is then used to develop a vector representation of watersheds and drainage lines. Using this information, a geometric network is built and then converted into a grid. Using ArcGis 10.2 Spatial Analyst Tools module, the density of the hydrographic network shown in the map (Fig. 6-d), which takes into account the cumulative length of the drainage system and the water nodes (affluent and confluent), was calculated. It has five density classes, i.e.: very low, low, medium, high, and very high.

Slope degrees

The slope gradient is a permanent geomorphological factor crucial to understanding the dynamics of the slopes because it conditions the stability of the slopes. According to Pareta et al. (2012), when the slope gradient increases, it results in an increase in tangential stress; consequently, increasing the possibility of ground ruptures. Thus, a potential energy is available at the slope, to explain the initiation of a mass movement. The presence of the slope promotes runoff and gravity drainage, which catalyzes the appearance of phenomena such as debris flows and collapses (Capdeville et al. 2001). The calculation of the slope gradient takes into account the vertical component of the normal vector to the regression plane closest to the four points of the DEM delimiting the grid in question. This variable was calculated based on ArcGis 10.2 algorithm. Slope values are calculated using average maximum techniques (Burrough and McDonnell 1998) (Fig. 6-e). After the generation of the slope gradient parameter, several classes were obtained, which we reclassified into five classes, with a pitch of ten degrees areas with a gradient of slope between 0−10^∘ have relatively flat terrain that is located at the alluvial plains of main soft-material rivers. The strongest slope gradients are particularly located at the levels of the limestone of the dorsal of the internal Rif domain and the sandstone Flyschs nappes.

Slope aspects

This parameter expresses the orientation of the slopes in relationship to the North. It is determined by clockwise measurement of the angle between the North and the horizontal projection of the maximum slope line of a slope. The spatial distribution of MM is closely correlated with the exposure of slopes in a given region, since it is largely controlled the water and thermal regime (rain interception, evapotranspiration, surface and/or deep flow, type of vegetation and its rooting, etc.). The slope aspects map (Fig. 6-f) was generated automatically from radar Aster images previously used to build the slope aspects map. This map was made using the “Surface Analysis” of the Spatial Analyst, which is implemented in the ArcGis 10.2 software. The orientation of the slope aspects in Tétouan region mainly influences the degree of MM on Northeast facing for complex mud flows and landslide types. The south-oriented slopes with dry weather, SSE, SSW, show mainly runoff-dominated erosion. Finally, Block falls and complex slides are located on the west-facing slopes.

Elevation

This is an appropriate parameter to distinguish the different altitude levels in the terrain. We used it to identify the stage of relief development and quantify the relationship between altitudes and MM. It represents the volume of surface above a given altitude. In this study, the elevation map of our study area is generated from SRTM (Fig. 6-g). All of the above factors have been integrated into a GIS-based environment and stored in a raster grid format of 90 x 90m resolution. The grid of the study area is 496 lines by 496 columns. Then all the layers were divided into five classes.

The rating and weighting of each factor is shown in Table 1 above and the percentages of each factor’s area are given below and explained as follows:

Table 1 Relative weights of the each parameter (W_j) and relative weight of different classes (w_ij) influencing the MM genesis in the province of Tétouan

a- Mudstones and marlstones are weighted due to the high number of mass movements and their widespread outcrops. Mudstones (32% of occupied areas) and marlstones (32% of occupied areas) have a high value because of their geotechnical characteristics compared to other formations (Fig. 7-a and Table 1).

Fig. 7

Percentage of zones of seven predisposition factors in Tétouan Province

b- Vegetation is an important parameter influencing the events of mass movements. In addition, a diversity of vegetation cover is often changing the behavior of mass movements. Landuse is an indirect factor in slope stability due to the small number of mass movements and their limited occurrence in the study area. The water, builtup areas, and forest represent 2%, 3%, and 22% respectively. The Bare lands cover 27% of the study. These latter are vulnerable to severe weather and erosion, therefore, receive high in weight. At the end, agriculture (46%) also increases the susceptibility of mass movements and end up with the highest weight (Fig. 7-b and Table 1).

c- Fault density was used to generate the MMS map. The areas of major faults between the limestone and mudstones of the Tangier unit, as well as the layers of Flysch and the Tangier unit show a significant number of mass movements. Ait Brahim et al. (2018) mentioned that the study area contains significant faults with a high frequency of instability. The faults whose density is between 583m and 768m and between 416m and 583m represent the largest percentage (27%) for each class (Fig. 7-c and Table 1).

d- The distance closest to the water system is the highest probability for the triggering of mass movements. Broadly speaking, in relationship to shoreline undermining in meandering portions of streams. However, in the study area, there is a significant number of mass movements on either side of the banks of the straight-drew water system. The phenomenon is accentuated especially in the streams that are more in soft lands (mudstones and marlstones of Flych nappes) (Fig. 7-d and Table 1).

e- The slope is one of the factors that causes the increase in shear stresses at the level of the different rocks and soils present on the slopes. The steeper the slope, the more area becomes prone to mass movements. In our study region, slopes above 40 degrees receive the highest weight (4 and 5) (Fig. 7-e and Table 1).

f- The slope aspects map and our mass movements show that they are more abundant on north-facing slopes (23%, 5) and east-facing (21%, 4). For this reason, these directions get the higher weight (Fig. 7-f and Table 1).

g- At the study area level, a significant number of mass movements begin at high altitudes. Nevertheless, altitudes 706-1118 m and 1118-1920 m (4 and 5) have a high weight (Fig. 7-g and Table 1).

After determining all the predisposition factors of the ground movements, the susceptibility will be calculated for each pixel of the province of Tétouan.

Results

The calculation and application of the heuristic approach resulted in the susceptibility map divided into five classes (Very high, high, medium, low, and not susceptible). The final map classification was obtained using ArcGis 10.2. In order to highlight the easiness of the use of our map to our readers, we have created a map that contains only the dams, tourism infrastructures, roads, and certain urban areas the province of Tétouan (Fig. 8); in addition to this, we have drawn other layers, such as water, urban construction, cities, villages, airport, dams, roads on the final susceptibility map (Fig. 9).

Fig. 8

Map of dams, tourist infrastructure and roads in Tétouan Province

Fig. 9

Mass Movements Susceptibility Map of Tétouan Province

The mass movements susceptibility areas are explained by the combination of several factors: 1- Their geological setting (Mediterranean and mountains regions). 2- Their exceptional natural predispositions (climate is very varied and quite often extreme, mild in winter (150 mm/day), hot and dry in summer with changeable spring, and autumn either very wet or very dry; rugged terrain; dominance of argillaceous rocks). 3- Its anthropogenic activities (economic growth, demographic pressures, overgrazing, wood cuttings, forest clearing, forest fires, break up in slope following the installation of the road network as part of plans to develop the northern province and open up the Rif). Nevertheless, to minimize the mass movement vulnerability over the critical areas, an important number of urgent and mandatory measures could be conducted in order to slowdown the velocity of the mass movements over the study areas: 1- Reforestation of bare soils. 2- Forest protection from fire risks by raising public awareness and speeding up the response in the event of an incident or breakdown. 3- Maps creation of sustainable urbanization that prohibit the installation of anarchic housing and infrastructures in areas at high risk of mass movements, such as unstable steep slopes, banks of dried wadis, downstream of areas where rock blocks collapse. 4- Citizens especially investors awareness about areas at risk. 5- Landslides monitoring² particularly slow landslides and mudslides) by different methods of monitoring of their evolution such as the measurement of electrical resistivity (ER) (Crawford and Bryson 2018), geolocalization sensors with GPS (Benoit et al. 2015), high-resolution data of remote sensing images (Jaboyedoff et al. 2012) like LiDAR and InSAR to create a landslide of 3D representation and measure their evolution over time (Jaboyedoff et al. 2012; Travelletti et al. 2014). Recently, we have developed a monitoring method based on the use of the Internet of Things (IoT) (El Moulat et al. 2018). This model is ground-based remote monitoring system consist of more than just field sensors; they employ data acquisition units to record sensor measurements, automated data processing, and display of current conditions usually via Artificial Intelligence (AI) (Elmoulat et al. 2020; Elmoulata et al. 2021). This model is ground-based remote monitoring system consisting of more than just field sensors; they employ data acquisition units to record sensor measurements, automated data processing, and display of current conditions usually via Artificial Intelligence (AI). Our findings of this new monitoring approach are designed to detect “where” and more importantly “when” the slopes are ready to slide and can provide early indications of rapid and catastrophic movements. It reports also continuous information from up-to-the-minute or real-time monitoring, provides prompt notification of MM activities, advances our understanding of their behaviors, and enables more effective engineering and planning efforts.

Discussion

MMS map shows that the urban centers and major rural settlements of the Province of Tétouan are in hazard zones between moderate and high (3 and 4), or even very high for a number of rural communes in the mountainous area on the 5 classes. Hence, the problem of their vulnerability needs to be addressed. In the field, there are MM activated by mostly anthropogenic activities (excavations, overweight, disturbance of pipelines) in several urban centers in the various urban municipalities of Tétouan Province. In addition, urban extensions are carried out in areas of stability, precarious and not previously developed by the urban agency (El Kharim 2012; Ait Brahim et al. 2002; Fonseca 2014; PATEAU 2014; PROkOS et al. 2016; Ait Brahim et al. 2018). 1- The susceptibility map of mass movements shows between the cities of Tétouan and Fnidek in the North, the presence mainly of a hazard zone between high and moderate (4 and 3) that encompasses the Tétouan Cluse and along the Mediterranean coast (Fig. 9). The percentage of susceptible areas are given by the f(Fig. 10). This area is characterized by the highest population density, residential buildings, tourist complexes (kabila, ksar Rimal, Marine Asmir etc.), ports (Mdik, marine Asmir), national roads (N16, N2, N13) and secondary roads (P4701, P4703), motorways (A6), and dams (Asmir, Hassan bel Mahdi, Ajrass). The hazard zone between high and moderate (4 and 3) is justified concretely at the level of the city of Tétouan and its surroundings (on the 2 slopes of the Oued Martil) by the impact of several mass movements(G6, Fig. 4) on the dwellings of the districts of Kurrat Sbaa, Khandak Ezzerbouh Samsa, Kaa El Hafa and roads N13, N16, N2, P4701 (Fig. 9). Impact of mass slides and collapses (G8, Fig. 11) and mudslide and boulders (G9, Fig. 11) and unstable areas (G7,G10, Fig. 11) on construction, tourist infrastructure, the A6 auto-road and the RN 13 national road in the Gabo Negro, M’diq, Jbel Zem and Fnidek region (Fig. 9) (El Kharim 2012; Fonseca 2014; PATEAU 2014; PROkOS et al. 2016; Ait Brahim and Elmoulat 2018). 2- The MMS map also shows, south of Tétouan (on 70% of the province’s surface) the presence of the largest hazard zone between very high and high (5 and 4). This mainly mountainous area that is characterized by a low population density and urban centers compared to the area located at the North of Tétouan (see map). However, the N2 national road between Tétouan and Chefchaouen (G3, Fig. 11) and the secondary roads N4702 (G3 and G4, Fig. 4) and P4704 (G12, Fig. 11) are experiencing disruptions due to complex landslides. The unstable slopes of the WadiHajera, Nakhla, El Kbir and Laou accentuated by the silting of banks during periods of flooding that contribute to the siltation of the Nakhla and Martil dams (G2,G3 and G4, Fig. 11), Oued Laou and My Bouchta dams (G13 and G12, Fig. 11). Villages that continue to develop are located within this area with very high and high level of risks (Zinat, beni Moussa, Ignezaouen, Ahoumdan, Ardifen etc.).

Fig. 10

Percentage of susceptible areas

Fig. 11

Examples of Vulnerability to Mass Movements (G) of Urban and Rural Centres, Tourist Complexes (Kabila, Ksar Rimal, Marine Asmir), Ports (Mdik, Marine Asmir) and Important Infrastructure (Dams, ports, highways, national roads) - G7: Cabo Negro Tourist Complex, G8: Mdiq (A6 motorway, and urban centre), G9: Jbel Zem-Zem (auto-road A6), G10: Fnideq (Urban extension on unstable areas,G11: Semsa dam, G12: My Bouchta dam, RN13, P4704; G13: Oued Laou Dam, Route P4105; G14: Idrene Slide

Conclusion

The novelty of this study is the use of the heuristic indexed method for modeling MMS. The selection of this method in regard to other methods (fractal, factorial, neurons, etc.) is justified by the possibilities of expert’s intervention who relies on the ground truth, and therefore he/she has the freedom in selecting the controlling factors, identifying the selected classes, and assigning the corresponding weights. This approach made this study different form other automatic methods that run the calculation without previous analysis and/or correspondence to the ground truth. Our findings confirm that 70% of our results correspond to field data used for validation purposes. The originality of this research remains in the fact that the created susceptibility map will indicate, by cross validation with ground truth, the areas of roads, dwellings, extension of urbanization, dams etc, that are located in risky zones. This latter is raising the possibility of increased awareness and providing more information and/or intervention in order to overcome the issue of their vulnerability for better socioeconomic developments; either by stabilizing the slopes, or comforting, or by choosing new bypass routes etc.

In closing, The MMS map developed in this research study can serve as an important document to be considered for the construction project conducted by the Ministry of Housing and Urban Affairs (mass movements, foods, and earthquakes) with a budget of one million Euros for the main urban and rural centers in order to point out: 1-Safe Zones, which could be opened up for regular urbanization; 2 Areas at risk, requiring general planning arrangements (stabilization of slopes, treatment of soils, etc.); 3- Zones requiring special arrangements to be made by constructors (strengthening soil structures etc.); 4- Strictly prohibited zones for any type of construction, because they are inserted in areas presented with non-repairable risks.