Flood Susceptibility Assessment through Statistical Models and HEC-RAS Analysis for Sustainable Management in Essaouira Province, Morocco

Abstract

Floods are natural disasters that often impact communities living in low-lying areas in the northern and central parts of Morocco. In this study, our aim was to create a flood susceptibility map using three methods; the hierarchy process (AHP) frequency ratio model (FR) and the weights of evidence (WoE) model. We extensively examined the area identified by these approaches using a hydraulic analysis software called HEC-RAS (version 6.3.1). Our analysis focused on the Essaouira watersheds in Morocco, where we identified around 197 flood locations. Out of these, we randomly selected 70% for modeling purposes while the remaining 30% were used for validation. Ten factors that influence floods were considered, such as slope, elevation, proximity to rivers, drainage density, stream order, land use patterns, rainfall data, lithology (permeability level) index (TWI), and curvature. We obtained these factors from data sources. Finally, we generated a flood susceptibility map and evaluated its accuracy by calculating the area under the curve (AUC). The validation results confirmed that all three models were robust and effective with an AUC of 90. Moreover, the research uncovered a trend of vulnerability with the most susceptible area being in close proximity to the city of Essaouira along the Oued Ksob. A detailed analysis using HEC-RAS was conducted at this identified location, pinpointing the village of Diabat as highly exposed. These findings hold significance for flood management, empowering decision makers, scholars, and urban planners to make informed choices and implement strategies that can minimize the impact of floods in susceptible regions while minimizing potential damages.

1. Introduction

Floods are an occurrence that can have consequences for both people and property [1,2,3]. They happen frequently and impact various aspects of the environment, including buildings, possessions, and even human lives [4,5,6]. In fact, floods account for 40% of all losses caused by natural disasters [7]. Urban areas are especially vulnerable to the effects of floods due to their population density and size [8]. Floods can be triggered by rainfall or melting snow, which overwhelm the drainage system and cause rivers to overflow their banks. This results in a surge of water downstream and at the watershed outlet leading to flash floods [6,9,10,11].

Changes in land use, land cover, and urbanization activities combined with an increase in density in flood areas have the potential to contribute to flood damage [6,12,13]. Certain practices like deforestation for settlement construction and alterations to river courses can also play a role in the occurrence of floods [14,15,16]. In fact, severe flooding has impacted a number of individuals from 2000 to 2008 with a higher concentration observed in developing nations [17]. Floods can cause indirect harm, including damage, to transportation systems, cultural landmarks, natural ecosystems, agriculture, bridges, and the economy [18,19,20]. Additionally, floods lead to the loss of thousands of lives and the displacement of millions of people worldwide each year [21,22,23,24].

In recent decades, remote sensing and geographical information systems (GIS) have emerged as tools for managing hazards and creating geo-spatial maps using various statistical models [16,25]. Flood hazard studies have been conducted on various scales [16,26] in regions without hydrological stations. Remote sensing technologies are employed to extract information regarding floods and the extent of the resulting damages [27,28,29,30,31,32,33].

There are three methods for mapping flood-prone areas: physically based, empirical, and physical modeling [34,35,36]. Physical modeling involves conducting experiments on small-scale models to verify the accuracy of the model predictions [34]. It can be simulated in one, two, or three dimensions using models like Delft3D or HEC-RAS. These models help determine the flooding process. This can predict the flood depth and extent for return periods [37,38,39]. The physically based approach is highly accurate in predicting floods as it simulates the hydrodynamic processes of fluids like water. It requires input data including hydrological, topographic, morphological, and remote sensing data processed in GIS [6]. These models help determine the spatiotemporal flooding process and can predict the flood depth and extent for return periods [37,38,39].

The empirical approach can be categorized into three models: criteria analysis, statistical methods, and machine learning and artificial intelligence models [6,40,41,42,43,44]. Among the physically based, the physical, and the empirical models, the later approach was the one most employed between 2000 and 2021, with approximately 46.2% of studies utilizing it. The physical modeling approach accounted for 43.8% of studies, while the physically based approach represented 10% [34].

For flood analysis and mapping flood areas, one popular model within criteria analysis is the analytic hierarchy process (AHP) method [45]. AHP involves an evaluation process that includes pairwise comparisons of parameters, ranking them accordingly to find the solution to a given problem [6,34,46,47]. This approach has been successfully implemented using Geographic Information Systems (GIS). Remote sensing techniques [16,48,49] have proven to be an accurate method that can be easily applied in various regions across the globe [16,50,51]. The ratio (FR) method is an effective statistical technique with a simple yet comprehensible concept [52]. It enables base scenario analysis (BSA) to study the impact of different factors on flooding occurrences [13,24]. Furthermore, while the weights of evidence (WoE)-based method is commonly used for landslide mapping purposes, its application in flood modeling is new [53,54]. WoE utilizes the model for decision making under uncertain conditions. This makes it suitable for hazard mapping as it considers the uncertainties associated with hazard events and their relationships with landscapes [55]. This approach allows for the implementation of the BSA to the extraction of correlations between the classes of conditioning factors and flooding events [56].

Generally, techniques for mapping flood susceptibility rely on a range of conditioning factors that represent the characteristics of the area under investigation. The specific selection of these factors heavily depends on the scale at which flood susceptibility analysis is being conducted [57]. These factors typically include geology or lithology morphometric properties such as elevation and slope river network density, soil permeability, and land use/land cover, among others [57].

As hydrological methods require fieldwork and significant budget allocation for data collection purposes [58], Tith (1999) [59] conducted a study in Austin, USA, to determine the floodplains by combining HEC-RAS and GIS. The findings highlighted that this combination of simulation techniques offers capabilities. Similarly, Qafari (2004) [60] used HEC-RAS and GIS to simulate the behavior of Babolrood River. It was concluded that these tools possess great potential in simulating flood hazard zoning [61].

Morocco, a country situated in northwest Africa, is currently grappling with the impacts of climate change. One of the consequences is the occurrence of intense rainfall, which often leads to hazardous flooding situations [62,63,64]. Moreover, Morocco faces challenges in terms of unpredictable rainfall patterns [65,66,67]. Throughout history, several regions in Morocco have witnessed floods that have resulted in loss of human lives. Notable incidents include 100 deaths in 1950, 730 deaths in 1995, and 500 deaths in 2014. These flood events have also caused losses amounting to around USD 2.2 million in 2001 and USD 5.2 million in 2014 [63,67,68]. Fortunately, advancements in natural science studies specific to the arid areas of Morocco provide valuable insights that can aid authorities in effective spatial management strategies.

The Essaouira province in the western coastal region of Morocco has experienced severe flood occurrences in the past due to rapid urbanization. However, conducting an accurate assessment of flood susceptibility is a challenging process due to the lack of data and other obstacles. This paper aims to: (i) provide an inventory of the flood records and marks; (ii) assess the role of flood conditioning factors; (iii) demarcate the critical locations to flooding in the Essaouira province by combining two geometric simulation techniques—(1) geo-spatial analysis, AHP model, FR model, and WoE model; and (2) the HEC-RAS hydraulic model in the Oued Ksob near the Moroccan city of Essaouira; and (iv) validate the produced susceptibility maps and identify the most accurate one.

2. Study Area

The study area under consideration in this research encompasses the province of Essaouira, located on the western side of the High Atlas mountains in Morocco, and covers an area of approximately 6335 km². The administrative boundary of the province and the boundary of the watersheds within the region were taken into account during the delimitation of the examined area. The province was created in June 1975 and comprises 57 municipalities, including the city of Essaouira, which is the capital of the province. According to the 2014 census, the province’s population was 450,527, with five urban municipalities and a medina spanning about 30 hectares along the coast on a rocky peninsula in the ocean [69].

This area is particularly susceptible to drought due to its reliance on rainwater for recharge [70]. It has a diverse climate with variations in temperature and precipitation due to its location between the ocean and the mountains [71]. The province is part of Morocco’s semi-arid environments with an average annual rainfall of less than 300 mm [72]. In the summer, hot and dry winds from the east and southeast cause an increase in temperature, and the average temperature fluctuates around 20 °C with a thermal amplitude of up to 17 °C between the driest month (August) and the coldest month (January) [73].

The area has three important aquifers: Korimat, the Meskala syncline, and Bas Ksob. However, the lack of water resources and rugged terrain with steep slopes mean that the water is quickly transported to the sea without being used for irrigation [73,74,75,76,77]. Groundwater is the primary source of water supply for domestic needs and to a lesser extent for agriculture and industry. The bassin of Essaouira dominates the study area, created by the central Atlantic Rift during the Triassic period. The Bassin of Essaouira is an open synclinal area that extends between two anticlinals, the Jbel Amsittène to the south and the Jbel Hadid to the north [78]. Despite the presence of several synclinal cuvettes, it is affected by numerous tectonic accidents, including the Bouabout, Korimat, and Essaouira synclinal cuvettes, as well as the Akermod plain [79].

The study area’s watersheds, which encompass five major catchment areas, include Chaibat El hamra, Ouazzi, Ksob, Tidzi, and Igouzoulen, with the Ksob watershed being the largest (Figure 1). Additionally, a small portion of the Tensift watershed lies in the extreme north border of the study area.

On a hydrological level, the province of Essaouira faces significant disadvantages. Oueds such as Oued Ksob only have intermittent flow, solely occurring during the rainy season and in a torrential manner. The region is confronted with severe issues in terms of water resources. Regarding surface watercourses, it is worth noting that the province is intersected by the following wadis, each with an annual supply: Oued Tensift (250 Mm³), Oued Ksob (46 Mm³), Oued Igouzoulen (17 Mm³), Oued Igrounzar (46 Mm³), and Oued Zelten (17 Mm³) [80]. In order to mitigate the water scarcity, two dams were constructed to store the surface water; Moulay Abderrahmane dam on the intersection of Oued Zelten and Igrounzar (Figure 1), and Sidi Mohamed ben Aberhman El Jazouli dam (Figure 1) on Oued Igouzoulen.

Rainfall and flow data incorporating information were collected from seven meteorological stations inside and surrounding the study area, namely Abadla, Adamna, Chichaoua, Igrounzar–Zelten, Talmest, Igouzoulen, and Azrou (Figure 1). In order to perform hydraulic modeling, hydrological data are essential input parameters, including the maximum annual instantaneous peak discharge, which is required to estimate the magnitude of floods for different return periods (RPs). For this study, the Adamna hydrological station, which is close to Essaouira city (Figure 1), was utilized, and the available 52-year record (1970–2021) of instantaneous peak discharge data (

Table 1. Maximum annual instantaneous flow rate of Oued Ksob in Adamna station.

Year	Q Max (m3/s)	Year	Q Max (m3/s)	Year	Q Max (m3/s)	Year	Q Max (m3/s)
1970	325	1983	624	1996	1602	2009	2270
1971	250	1984	435	1997	619	2010	469.1
1972	73.7	1985	1330	1998	607	2011	142.5
1973	140	1986	141	1999	150.5	2012	391.1
1974	54.7	1987	583	2000	303.6	2013	92
1975	94.3	1988	1050	2001	1240	2014	287.5
1976	162	1989	69.2	2002	268.3	2015	22
1977	255	1990	80.5	2003	157	2016	51.2
1978	103	1991	568	2004	76.7	2017	6.7
1979	217	1992	10.9	2005	2550	2018	35
1980	613	1993	177	2006	421	2019	0.6
1981	182	1994	24.4	2007	175.4	2020	4.5
1982	198	1995	643	2008	1079	2021	62

) provided by the Hydraulic Basin Agency of Tensift (ABHT) was used. Previous studies, such as Namara et al. (2022) with 25 years [81], Prastica et al. (2018) with 16 years [82], and Al-Zahrani et al. (2017) with 15 years [83] demonstrated that this duration is sufficient for this type of study. Oued Ksob (Figure 1) flows in an east–west direction throughout its length and receives several tributaries on its left and right bank, and the branching of the network is equally pronounced between the upstream and downstream, resulting in five classes with Strahler classification.

The mean annual rainfall in the study area is approximately 280 mm, and it has two main seasons: (i) the wet season which includes winter and autumn, with the highest precipitation occurring between December and November (with marked peaks in autumn and winter that gradually decrease from February to May); and (ii) the dry season from April to September, which is characterized by low rainfall, with July and August being the driest months. The coastal zone has higher precipitation than inner areas, as well as higher humidity levels (above 75% of humidity). In Essaouira and other areas exposed to maritime influences, summer fog is particularly significant [84,85].

To evaluate the effects of morphologic details and geometric structure on flood inundation maps [86], it is crucial to establish the topographic predisposition of the watershed, since morphologic data affect flow conditions [87] and hydraulic simulations via geometric discretization utilized in the model [88].

The Oued Ksob, which borders the study area, is the most significant watercourse in the vicinity. Within the riverbed, there are small- and medium-sized pebbles, shaped stones, and substantial amounts of sand. The riverbanks of the oued can reach heights of up to ten meters above the riverbed. However, in certain areas, they experience erosion primarily due to the terrain’s nature and the strong flow of the river. Vegetation is absent within the riverbed, but beyond the riverbanks, it is typically present in the form of cereal crops. There are no specific obstructions impeding the water flow in the Oued Ksob, except for those created by ABHT and the mobile dunes at the outlet.

The Oued Ksob catchment has a Gravelius compactness index KG = 1.3 and is characterized by an amoeboid shape, abrupt slopes, and steep reliefs, with elevations ranging from 0 to 1697 m. Steep slopes are more prevalent in the upstream border, particularly in the southeast corner, which is closer to the High Atlas mountains. On the other hand, the extreme west (outlet) features gentle slopes and flat landscapes.

3. Methodology

The methodology used in this research involved several steps to assess flood susceptibility in the Essaouira province. Firstly, an inventory of flooded points was created, followed by the assessment of several conditioning factors that were responsible for the occurrence of floods. The next step involved combining remote sensing and geospatial analysis, as well as the use of the AHP, FR, and WoE models to produce flood hazard potential maps. The accuracy and predictive capability of various flood susceptibility models were evaluated using the receiver operating characteristic (ROC) curve and calculating the area under the curve (AUC), as demonstrated in studies by Linden (2006) [89] and Remondo et al. (2003) [90], and, the most accurate map was chosen later. Finally, the most susceptible area was identified and the hydraulic model of HEC-RAS was applied for further analysis. Field surveys and validation were used to support the preparation of the data sources used in the analysis (Figure 2).

3.1. Data Source

To develop flood susceptibility models for the study area, the lack of GIS data in digital format was a major challenge. To address this issue, several GIS raster layers were created using a freely available digital elevation model (ALOS PALSAR RTC DEM) with a 12.5 m resolution, land use data from 10 m resolution Sentinel satellite imagery, soil permeability data from 1/100,000 and 1/500,000 scale geological maps (

Table 2. Data sources.

Data Type	Data Denomination	Source	Scale/Resolution/Duration
Geological maps	Tamanar map	Ministry of Energy and Mines, Water and Sustainable Development	1/100,000
	Marrakech map		1/500,000
	Taghazout map		1/100,000
Meteorogical data	Adamna station	Hydraulic basin agency of Tensift (ABHT)	1977–2021
	Igrounzar–Zelten station		1968–2021
	Azrou station		2003–2021
	Talmest station		1984–2021
	Igouzoulen station		1997–2015
	Chichaoua station		1965–2021
	Abadla station		1969–2021
Satellite images	Sentinel	https://scihub.copernicus.eu/ (last access: May 2023) (Copernicus 2021)	10 m
	Digital elevation model	https://search.asf.alaska.edu/ (last access: May 2023) (JAXA/METI 2020)	12.5 m
	High-resolution ortho-imagery	https://earthexplorer.usgs.gov/ (last access: May 2023) (USGS-EROS 2018)	0.3 m

), and precipitation data from the Tensift Hydraulic Basin Agency. Flood risk analysis was conducted using a cell size of 12.5 × 12.5 m², while the HEC-RAS analysis utilized a cell size of 1 m × 1 m.

Given that high-resolution and accurate DEMs are crucial for precise hydrological analysis and the estimation of hydrographs, low-resolution data can significantly and indirectly impact hydrological modeling. To ensure the reliable and precise HEC-RAS analysis, the digital elevation model (DEM) used was extracted from elevation points obtained through surveying data using Leica differential GPS in RTK for the stream bed and surrounding area, contour lines, hydrographic network, and using a drone surveying Diabat village (which is the most susceptible area). These elevation points were then interpolated to generate a DEM with a resolution of 1 m, which provided high-accuracy elevation data, thereby enhancing the credibility and precision of the hydraulic model used in the analysis.

3.2. Flood Inventory

To accurately estimate the probability of the exceedance of flood events in a region, it is crucial to examine the historical records of past floods [56,91]. To accomplish this, it is essential to have an inventory map that documents flood occurrences in the area [92]. In our study, we created a flood inventory map for the Essaouira province by mapping flood locations using both documentary sources and field surveys. The historical flood inventory map does not impose any limitations on occurrence dates.

For our analysis, we divided the flood inventory map into two segments—70% for training and 30% for testing [56,93]. The training set consisted of selected flood locations and a separate layer representing floods was created as a variable.

3.3. Flood Conditioning Factors

To perform flood susceptibility mapping, it is crucial to identify the flood conditioning factors [4,56]. Therefore, in the first stage of this study, a flood-related spatial database was created. The conditioning factors were chosen based on the knowledge gained from the literature [94] and field investigations. Ten flood conditioning factors were selected for the susceptibility analysis, and the spatial database of these factors was compiled. The conditioning factors are as follows: slope (°); elevation (m); distance from river (m); drainage density (km/km²); stream order (which represented the flow accumulation); land use/cover (LULC); weighted rainfall; lithology (permeability level); topographic wetness index (TWI); and curvature (

Table 3. Classes, rates, and weights of conditioning factors.

Factor	Classes	Susceptibility Class Ratings	Weights (%)	Factor	Classes	Susceptibility Class Ratings	Weights (%)
Slope (°)	0–2	5	11.43%	Land use	Bare/spare vegetation	3	6.39%
	2–5	4			Built-up	4
	5–15	3			Dense vegetation	1
	15–35	2			Medium vegetation	2
	>35	1			Water	5
Elevation (m)	0–80	5	15.80%	Rainfall	269–290	1	10.88%
	80–350	4			290–330	2
	350–600	3			330–365	3
	600–950	2			365–380	4
	950–1697	1			380–407	5
Distance from river (m)	0–100	5	7.95%	Lithology (permeability level)	Very low	5	5.44%
	100–200	4			Low	4
	200–1000	3			Moderate	3
	1000–2000	2			High	2
	2000–8166	1			Very high	1
Drainage density (km/km2)	0–0.55	1	10.87%	TWI	<7	1	9.24%
	0.55–1	2			7–9	2
	1–1.47	3			9–12	3
	1.47–1.92	4			12–16	4
	1.92–2.38	5			>16	5
Stream order	Order 1	1	12.76%	Curvature	<−2.46	5	9.24%
	Order 2	2			−2.46–−1.13	4
	Order 3	3			−1.13–−0.14	3
	Order 4	4			−0.14–0.51	2
	Order 5	5			>0.51	1
	Order 6	5

). The choice of these factors was based on their theoretical relevance to flood risk.

The conditioning factors were classified as nominal, ordinal, and scale, and the scale factors were classified using the popular method of natural breaks (Jenks). Once the dataset was prepared, each conditioning factor was transformed into a grid spatial database with a size of 12.5 × 12.5 m, and the grid of the province of Essaouira was constructed using 1722 pixels. All the information was organized in a database using ArcGIS-Pro software.

-

Slope: The slope degree can serve as a surface indicator to identify the likelihood of flooding [21]. The slope of an area affects the movement of surface runoff and water infiltration, where regions with flatter slopes may experience prolonged flooding, whereas areas with modest to higher slopes facilitate the faster drainage of floodwater [67,95]. The slope map of the study area, was generated using the ALSA Satellite DEM and ArcGIS-Pro software.

-

Elevation: The study area’s flooding pattern is significantly affected by the elevation gradient. It plays a crucial role in controlling the direction of overflow and determining the flood depth [48,96]. Furthermore, varying elevations alter climatic conditions, which consequently impact the vegetation and soil characteristics [97]. Figure 6 presents the altitude map of the area, generated using ArcGIS-Pro from the ALSA Satellite DEM.

-

Distance from river: The distance from rivers is a crucial factor in identifying the flood-prone areas. Previous research has demonstrated that, as the distance from the river system decreases, the probability of flooding increases [33,98,99]. The Euclidean distance tool in ArcGIS-Pro was utilized to generate this map.

-

Drainage density: The time of concentration of runoff and the probability of flooding are linked to the drainage density in an area. Previous studies [67,100,101,102] have established that the volume of flow accumulation from upstream to downstream in the river basin is affected by drainage density. A higher drainage density indicates a greater surface runoff, which in turn increases the likelihood of flooding [103].

-

Stream order (representing the flow accumulation): The accumulation of water flow in a cell is a measure of the amount of water drained from upstream cells. The volume of water flowing downhill is directly proportional to the risk of flooding as it indicates higher discharges [67,99]. Initially, we tried using the flow accumulation map as an indicator but found it unsatisfactory due to its linear nature. Consequently, we replaced it with the buffered stream order map, which provided a similar variation. The buffer was created proportionally to the stream orders and the major riverbed of the drainage system.

-

Land use/cover (LULC): Identifying zones with high susceptibility to flooding can be accomplished through the consideration of land use/land cover, as pointed out by Norman et al. (2010) [104]. The type of land use plays a key role in determining the infiltration of rainwater into the soil, which then affects the resulting runoff [99]. Forests generally promote infiltration through the root system of trees and plants, while impervious surfaces such as roads and buildings hinder the infiltration of water, resulting in increased surface runoff [56]. To produce the land use map, Sentinel images were utilized and the supervised maximum likelihood classification was applied for ranking purposes.

-

Weighted rainfall: Rainfall plays a role in determining the likelihood of flooding in an area. It directly influences river flow. When heavy rain occurs within a period, it can lead to flash floods especially in semi-arid regions [105,106]. Flooding is primarily caused by precipitation. This study focused on two factors related to rainfall; the amount of rain within 72 h each year and the annual frequency of rainstorms. To calculate the number of rainstorms, historical flood data were analyzed to determine the amount of precipitation required to trigger such an event at existing monitoring stations. The count of rainstorm events exceeding this threshold within a year was considered as the number of rainstorms [107,108,109]. The total count of rainstorm events between 1984 and 2021 was used to determine the rainfall within 72 h each year. Rainfall data from stations including Abadla, Adamna, Chichaoua, Igrounzar–Zelten, Talmest, Igouzoulen, and Azrou were collected for this analysis. The final map depicting these findings was created using inverse distance weighted (IDW) interpolation in ArcGIS Pro.

-

Lithology (permeability level): The hydrological processes in watersheds are influenced by lithology, which is determined by rock permeability, the thickness of the layers, and their outcrop size [110]. The lithology map was generated by vectorizing geological maps 1/100,000 and 1/500,000, and reclassifying the rock formations into zones with varying permeability levels ranging from very low to very high. For instance, schistose limestone, crusted limestone, and crystalline formations are major flood hazards as they increase the destructive potential of floods, thus assigning the least permeability. On the other hand, alluvial and continental deposits are highly permeable, with a greater capacity for infiltration, promoting groundwater flow and reducing flood intensity. Therefore, these formations were assigned a high permeability level [32].

-

Topographic wetness index (TWI): The topographic wetness index (TWI) is widely used to analyze the distribution of conditions [111,112,113]. This helps understand how topography affects the saturation levels and runoff, making it a valuable tool for assessing flood potential in watersheds. The TWI map is calculated using an equation outlined by Moore et al. 1991 [114]:

$$TWI=Ln\left(\frac{As}{tan\beta }\right)$$

(1)

where As represents, at each cell, the specific upslope contributing area of the watershed and β denotes the local slope (Rahmati et al. 2016) [61].

-

Curvature: The curvature is a significant factor in assessing the flood potential of watersheds [113]. The map depicting the curvature was generated using a digital elevation model and analyzing the slope using ArcGIS Pro.

3.4. Susceptibility Modeling and Validation

The input data for the Essaouira basin from diverse sources were analyzed, processed, and merged in a GIS setting to create three models for evaluating the flood susceptibility. These models include the analytic hierarchy process (AHP) for multi-criteria analysis, frequency ratio (FR), and weights of evidence (WoE).

3.4.1. Analytic Hierarchy Process (AHP)

The input data for the Essaouira province, obtained from various sources, were processed, analyzed, and integrated within a GIS environment. A multi-criteria analysis was conducted using the analytic hierarchy process (AHP) model. The AHP model was developed based on the conditioning factors and aimed to determine the suitable parameters and weight values for assessing flood hazards [26]. Its purpose was to identify areas with high susceptibility to flooding and enable comparative analysis across different basins [99,115,116]. The AHP model primarily considered natural factors, which were selected based on their relevance to flood hazards [117].

AHP, a semi-quantitative technique [118], is a widely used multi-objective and multi-criteria decision making (MCDM) approach that enables the user to achieve scale advantages drawn from a set of alternatives [119]. It has been applied in various fields of earth sciences and natural disasters, such as flood hazard zonation [120], landslide susceptibility assessment [118,121], and groundwater resource potential determination [61]. The AHP model employs two methods, questionnaire and non-questionnaire, to obtain the weight of each factor. In the questionnaire method, ten questionnaires were given to ten hydrology specialists to assign a proper weight to each conditioning factor based on Saaty’s standard Table [119] (Table 3). AHP is based on pairwise comparison, which makes judgements and calculations easy [122]. In the AHP model, the effect of each conditioning factor on flood susceptibility is determined relative to each other based on dual evaluation. The rate determines the weight of each parameter to flood susceptibility, and ranks give the weight of each class of flood-conditioning factors on the flooding (1–5). For some parameters, numbers 1–5 were reversely assigned, meaning a higher weight was given to the first slope class if it had the highest effect [20]. Each of the 10 factors was reclassified and assigned a rating based on its impact on flood risk. The ratings ranged from 1 to 5, as indicated in Table 3. These values reflected the intensity of the flood hazard, ranging from very low (1) to low (2), medium (3), high (4), and very high (5) (Table 3).

However, the AHP method has a drawback: the selection of indicator weight values based on expert opinion, which may introduce the subjectivity and cognitive limitations [123,124]. Nonetheless, this weakness can be reduced by assessing the consistency of the ratio. Saaty (1980) [119] determined that the threshold for the consistency of the ratio should be less than 10% to obtain a coherent value between the weighted variables [6].

The weights assigned to the factors in the Essaouira province were established through the application of the analytic hierarchy process (AHP) methodology, as depicted in Table 3. AHP is a mathematical technique employed to analyze intricate problems involving multiple factors. The factors were arranged in a hierarchical structure, and a pairwise comparison matrix was constructed to facilitate the meaningful comparisons between their assessments. The importance of each factor was determined using a set of five scales, which can be found in

Table 4. Numerical expression of the relative importance of factors.

Importance	Scale
Very important	1
Moderate	1/2
Less important	1/3
Moderately less important	1/4
Much less important	1/5

.

The values in each row of the comparison matrix indicate the correlation between two factors under consideration (Figure 3). To find the average of each column, we need to add up all the values in the column and then divide that sum by the number of elements in that column. This process allows for the normalization of the matrix, as described by Saaty (1990) [46].

The average of the rows in the normalized matrix represents the weight (w) assigned to each factor. In the case of the Essaouira province, the AHP model identified elevation and stream order (flow accumulation) as the most influential factors, followed by slope, rainfall, and drainage density.

To maintain the consistency of the eigenvector matrix, in analytic hierarchy process (AHP), it is crucial to evaluate the consistency ratio (CR) using Equation (2). The CR compares the consistency index (CI) (Equation (3)) of the matrix with the index (RI) of a generated matrix [101].

CR = CI/RI

(2)

The CR represents the measure of consistency while CI represents the consistency index. If the CR value is less than or equal to 0.10., this indicates a level of consistency in the matrix. However, if the CR value exceeds 0.10, it suggests inconsistency. This necessitates a reassessment of the judgements.

$$CI=\frac{{\mathsf{\lambda }}_{max}-n}{n-1}$$

(3)

The calculation of CI relies on obtaining ${\mathsf{\lambda }}_{max}$ from the comparison matrix and considering the number of factors involved ($n$). In this case, with ${\mathsf{\lambda }}_{max}$= 10.81, $n$ = 10 and a random index (RI) value of 0.091, we computed CI accordingly. By determining ${\mathsf{\lambda }}_{max}$ and examining the consistency across all judgements, we ensure that our AHP method adheres to Saaty’s suggestion that CR should be equal to or less than 0.1 as per his work in 1990 [46].

After analyzing the matrix, we found that the CR value is 0.061, which is below the threshold of 0.1. This suggests that the factors were consistently weighted. We then combined the 10 chosen factors with their weights using an equation. Subsequently, we calculated the flood risk index using Equation (4):

$$AHP=\sum _{i=1}^{10}(wi\ast Xi)$$

(4)

where $Xi$ represents the classification of each factor at a given location, $wi$ denotes the weight assigned to each factor, and (n) represents the total number of factors.

Utilizing Equation (5), the AHP analysis was conducted in the Essaouira province:

AHP = [11.43 × (Slope)] + [15.80 × (Elevation)] + [7.95 × (Distance From Rivers)] + [10.87 × (Drainage density)] + [12.76 × (Stream order)] + [6.39 × (Landuse)] + [10.88 × (Weighted Rainfall)] + [5.44 × (Lithology)] + [9.24 × (TWI)] + [9.24 × (Curvature)]

(5)

Assessing flood susceptibility allows us to examine and evaluate how likely flooding is to occur and determine how accurately our chosen model predicts it [125]. In studies, researchers have used the proportion of impacted features compared to the number of vulnerable features as a method to evaluate the susceptibility of a region to flooding [126,127].

3.4.2. Frequency Ratio Model (FR)

Various fields of natural hazards, such as groundwater spring and qanat potential mapping, landslide susceptibility mapping, and flood susceptibility mapping have utilized the FR method, which is one of the most commonly used bivariate statistical methods [20,56,120,128,129]. The advantages of this model are its simplicity in implementation and the ease of understanding its results [130,131].

The FR method is utilized for evaluating the impact of each influencing factor on flood occurrences [132,133]. The flood susceptibility index (FSI) is calculated by summing up the frequency ratios of all the given parameters or factors, as expressed in Equation (7):

$${FR}_{ij}=[N_{pix}\left({SX}_i\right)/\sum _{i=1}^m{SX}_i]/[N_{pix}(X_j)/\sum _{j=1}^n{N_{pix}(X}_j)]$$

(6)

$${FS}_i=\sum _{j=1}^n{FR}_{ij}$$

(7)

In Equation (6), $N_{pix}\left({SX}_i\right)$ represents the number of pixels with floods within class $i$ of the factor variable X, while $N_{pix}\left(X_j\right)$ represents the number of pixels within factor variable $X_j$. The parameter variable $X_i$ has m number of classes, and there are $n$ number of factors in the study area [20,134].

3.4.3. Weights of Evidence (WoE)

The method of weights of evidence (WoE) is an approach to produce hazard maps. It determines the importance of each factor by examining the relationships between flood locations and conditioning factors [54]. WoE offers advantages over statistical methods as it relies on data-driven techniques and utilizes the Bayesian probability model [135]. Although WoE has been commonly used for landslide mapping, its application in flood susceptibility mapping is relatively new [54,136].

To apply the WoE method, it is necessary to calculate two parameters, namely positive weight ($W_i^{+}$) and negative weight ($W_i^{-}$). The weight of each conditioning factor (B) is determined based on whether or not the disaster event (A) occurs within an area. This can be represented using the following equations;

$$W_i^{+}=\ln \frac{P\left\{B|A\right\}}{P\left\{B|\overline{\mathit{A}}\right\}}$$

(8)

$$W_i^{-}=\ln \frac{P\left\{\overline{B}|A\right\}}{P\left\{\overline{B}|\overline{\mathit{A}}\right\}}$$

(9)

In this context, the symbol P represents probability while ln stands for natural logarithm. The presence or absence of a conditioning factor is denoted by B and $\overline{B}$. Similarly A and $\overline{A}$ indicate the presence or absence of a flood event according to Xu et al. (2012), respectively [137]. A positive weight ($W_i^{+}$) signifies that the conditioning factor exists in disaster locations and its magnitude reflects a relationship with disaster occurrence. Conversely, a negative weight ($W_i^{-}$) indicates the absence of the conditioning factor, which shows the level of the relationship. By calculating the difference between $W_i^{+}$ and $W_i^{-}$, we can determine the association between the conditioning factor and disaster occurrence—a measure known as weight contrast. To assess confidence in this contrast being “real”, we calculate the contrast as the ratio of contrast to deviation. Regmi et al. (2014) [54] explained that this studentized contrast serves as a test to check whether it significantly deviates from zero. For formulation regarding the WoE method, refer to Bonham Carter (1994), Pradhan et al. (2010), Regmi et al. (2013), and Tehrany et al. (2014) [54,56,135,138].

In order to verify the accuracy of the flood susceptibility models, we established flood partitions as validation subsets and compared them in terms of their spatial distribution. To evaluate the performance of each model, we used the area under the curve (AUC) and ROC curves [89,90] to calculate the AUC values.

3.5. HEC-RAS Analysis

To model free surface flow in rivers with floodplains, numerical computation are essential [139]. Hydrodynamic modeling requires an accurate representation of the topography of the river channel and neighboring floodplains to simulate the streamflow as realistically as possible [140]. Various numerical tools, such as one-dimensional (1D), two-dimensional (2D), or three-dimensional (3D) approaches, are available for river and floodplain modeling, each with distinct capabilities and data requirements [141,142,143,144]. HEC-RAS model was selected for the study area based on meteorological and spatial data availability. HEC-RAS is a Windows-based hydraulic model developed by the U.S. Army Corps of Engineers that can perform one- or two-dimensional hydraulic computations for a full network of natural and constructed channels [145]. It can calculate water surface heights at all points of interest for given input values. The software (version 6.3.1) can simulate unsteady- and steady-state water flow in a single river or a complex network of interconnected reaches that includes hundreds of channels [146]. Additionally, the model can simulate different flow types, including supercritical, subcritical, and mixed flows. For a subcritical flow at each 1D cross-section, the Bernoulli Equation (10) is employed [147]:

$$Z_2+Y_2+\frac{{\mathsf{\alpha }}_{2V_2^2}}{2_g}=Z_1+Y_1+\frac{{\mathsf{\alpha }}_{1V_1^2}}{2_g}+h_e$$

(10)

where the elevation of the principal channel inverts is represented by $Z_1$ and $Z_2$, while $Y_1$ and $Y_2$ represent the depth of water at the cross-sections, and V¹ and V² represent the mean flow speed (total discharges/total surface flow), with α₁ and α₂ as velocity coefficients. The gravitational acceleration is denoted by $g$, and he represents the energy head loss.

HEC-RAS utilizes a set of parameters such as Manning roughness coefficients, cross-sections, and boundary conditions (upstream and downstream) to perform hydraulic simulations [148]. However, to determine the width of the floodplain area for different years of return periods, HEC-RAS alone may not be sufficient to generate a flood inundation map. Therefore, an extension in ArcGIS Pro was used for this purpose to generate flood hazard maps.

To estimate flooding, a range of spatial databases, including the digital elevation model (DEM), land use, and flow data of Oued Ksob, were collected from various sources such as the Hydraulic Basin Agency of Tensift, the US Geological Survey Earth Explorer, Alaska Satellite DEM, high-resolution ortho-imagery, and measured elevation points. To create a flood simulation in the downstream part of the Ksob basin surrounding Essaouira City, four important basic steps need to be followed.

• The first step involves preparing a triangular irregular network (TIN) in ArcGis Pro with 1 m resolution, which is created from the measured points obtained using professional GPS Leica and a drone in Diabat village, measuring a total of 29,205 points (Figure 4A).
• The initial stage of hydraulic modeling in HEC-RAS involves inputting the geometric data that characterize the study area. This includes creating a stream centerline, main channel banks, and material zones known as “channel geometry”, as well as choosing cross-section lines that run perpendicular to the river and are used to extract the elevation values from the terrain model (AQUAVEO, 2008). In order to obtain detailed flood maps, a total of 82 sections were considered across the oued with 50 m spacing and more than 4 km long, as shown in Figure 4B. It is important to note that reducing the number of sections may result in lower quality inundation maps [88]. Cross-section water surface elevation can be observed in Figure 4B. The roughness of the various cross-sections was determined using Manning’s n-values, which were assigned based on the land use characteristics of the study area in addition to an intensive field survey [149].
• To evaluate the extent of extreme events such as floods or low flows, the frequency analysis of the maximum annual peak flow is a popular statistical approach utilized in hydrology. This analysis determines the probability of exceedance for a given flood amplitude based on the design of flood sluicing buildings, bridges, culverts, and other flood control projects [150]. In this study, seven methods of allotment types were selected to determine the best distributions for hydrologic frequency analyses, including Weibull, Gamma, Lognormal, Pearson type 3, GEV, Gumbel, and Normal (

  Table 5. Discharges for different scenarios (2-, 5-, 10-, 20-, 50-, 100-year return periods) of seven different frequency distributions for the Adamna hydrological station.

  T(year)	GEV	Gumbel	Weibull	Normal	Lognormal	Pearson III	Gamma
  2	291	324	220	413	175	224	211
  5	730	805	657	871	675	708	681
  10	1060	1120	1050	1110	1370	1100	1080
  20	1410	1430	1480	1310	2450	1510	1510
  50	1920	1820	2090	1530	4730	2060	2090
  100	2350	2120	2580	1680	7330	2470	2540

  ). El Adlouni et al. (2007) [151] proposed two selection criteria, the Akaike information criterion (AIC) [152] and Bayesian information criterion (BIC) [153], which are used to decide which distribution is the best match based on the differences between the appropriate distribution and the empirical probability. The distribution associated with the lowest BIC and AIC values is considered to be the best match [154] as defined in Equations (11) and (12):

  AIC = −2 log(L) + 2k

  (11)

  BIC = −2 log(L) + 2k log(N)

  (12)

  with likelihood function (L), number of parameters (k), and sample size (N).

  Figure 4. (A) The constructed DEM from measured points; (B) cross-section lines, river, banks, and flowpaths.

  Figure 4. (A) The constructed DEM from measured points; (B) cross-section lines, river, banks, and flowpaths.

The roughness coefficients are estimated based on field observations with an adequate margin of safety in order to account for any changes that may affect the morphology of the riverbed, as well as the surface geology map provided in the report. The Manning’s “n” coefficient ranges between 0.024 and 0.044.

4.

For this study, hydraulic modeling was conducted using the HEC-RAS software (version 6.3.1) under unsteady conditions. The software solves the full 2D St. Venant equations for an unsteady open channel flow. The normal depth upstream and downstream boundary conditions were set at a friction slope of approximately 0.00375 for simulation in HEC-RAS. The flood hydrograph considered in this study corresponds to the adjusted flow rate for a 100-year return period, which was used to define the upstream boundary condition (inflow).

Following the input data in the HEC-RAS model, the output was exported to GIS in the form of a RAS GIS Export file. Flood depth mapping for the 100-year return period was generated using GIS.

4. Results

4.1. Flood Inventory

Similarly to other basins in the kingdom, the Tensift–Ksob–Igouzoulen basins experienced widespread rainy days with localized and unevenly distributed thunderstorm episodes. These rainfall events generated varying degrees of runoff in the region’s river systems, leading to an improvement in the water reserves in dam reservoirs.

These rains resulted in flooding in certain areas. These floods are primarily due to the contributions from pre-urban basins and the presence of water intakes on the riverbanks, causing the riverbeds to disappear under the accumulation of sediment.

The province of Essaouira experienced flooding in the city of Essaouira due to the malfunctioning of the stormwater drainage system. Other floods, caused by water inputs from peri-urban basins surrounding the centers (particularly Smimou, Essaouira, Igouzoulene, Ksob, Mejji, Tamanar, and others), led to runoff that posed a threat to the safety of residents and their properties.

A Geographic Information System (GIS) database was created to inventoried flood events in the Essaouira province. This database was constructed by mapping flood locations using both documentary sources and field surveys. A total of 197 flood locations were identified and used to generate a flood inventory map.

The flood inventory map was divided into two subsets: 70% for training and 30% for testing, following the methodology suggested by Ohlmacher and Davis (2003) [93]. The training set consisted of 138 randomly selected flood locations, which were used to create the flood layer as the dependent variable. The remaining 59 flood events were reserved for testing purposes.

Over the years, 11 localities in the Essaouira province have experienced recurrent flood events, leading to the loss of human lives, property damage, crop destruction, livestock loss, and the deterioration of health conditions due to waterborne diseases. These localities include Essaouira, Smimou, Sidi Jazouli, Tafedna, Tamanar, El Hanchane, Lagdadra, Had Dra, Tidzi, Zaouia Ben Hmida, Imintlit, Oulad Mrabet, and Mejji. Figure 5 displays images of the floods that occurred in these areas, along with their corresponding dates.

4.2. Flood Conditioning Factors

After preparing the dataset, each conditioning factor was transformed into a grid spatial database with a size of 12.5 × 12.5 m. The grid covered the entire Essaouira Province and consisted of 39,833,555 pixels, equivalent to 6224 square kilometers.

Initially, all the derived conditioning factors were continuous variables. However, the Kernel method was employed to classify each factor into different categories. The following ranges were used for classification: slope (0–75.21), elevation (0–1697), distance from the river (0–8166), drainage density (0–2.38), stream order (order 1–order 6), land use (sparse vegetation, built-up areas, dense vegetation, medium vegetation, and water), rainfall (269–407), lithology-permeability (very low to very high), TWI (1.19–19.79), and curvature (−42.88–41.6).

All the factors’ maps were classified into five classes, except for stream order, which had six classes. This was due to the large Tensift basin, which reached the sixth class in the Oued Tensift.

• The slope factor map (Figure 6) reveals that the downstream and coastal areas of the basin have a very gentle slope compared to the upstream region. This characteristic makes these areas more prone to flooding as runoff takes longer to drain, resulting in increased inundation.
• The elevation map, which is a crucial factor in flood susceptibility assessment, was categorized into five classes: 0–80, 80–350, 350–600, 600–950, and 950–1697 m, representing classes 1–5, respectively. The highest terrains are found in the southeastern borders, around Jbel Hadid in the north and Jbel Amseten in the central part. In contrast, the rest of the study area has relatively lower elevation.
• According to the map that measures the distance from the river, it is evident that areas within 200 m of the river are more susceptible to flooding. Conversely, areas located at a distance of 400 m or more from the river face no flood risk. The regions with vulnerability to flooding are primarily limited to the river networks (Table 3; Figure 6).
• The drainage density factor map provides a detailed depiction of the drainage density throughout the entire study area. This indicates that the intensity of hazards is particularly high around Oued Tensift and Oued Ksob, where the density ranges between 1.5 and 2.38 m per square kilometer (Figure 6). Lower density classes are concentrated in other parts of the study area.
• The stream order map factor, representing flow accumulation, demonstrates a range from order 1 (for small catchment areas and drains) to order 6. The highest values coincide with the water flow of the main Oued Tensift in the large Tensift watershed. The Oued Ksob is assigned the fifth order (Table 3; Figure 6). To identify areas with higher flood hazards, a value of 5 was assigned to the high order of this factor for AHP analysis.
• Land use in the Essaouira province exhibits strong contrast in terms of class size, properties (density, area, location, etc.), and temporal variability. Generally, it comprises sparse vegetation and bare soil. The land use data in the Essaouira province area were classified into five classes, as displayed in Table 1. The city of Essaouira, situated in the center of the province, with various infrastructures and impermeable materials, amplifies the occurrence of downstream floods. Agricultural activities are prevalent along the riverbeds in the form of terraces and bare soil. Residential areas dominate this class. The type of land use influences the rate of precipitation infiltration into the soil. Residential areas are more exposed to flood hazards compared to other land types. Therefore, a rating of 5 was assigned for water, 4 for built-up areas, 3 for sparse vegetation, 2 for medium vegetation, and 1 for dense vegetation (Figure 6).
• Weighted rainfall: Using the previously outlined methodology and examining the historical flood events from 1984 to 2021, it was found that the minimum precipitation value resulting in a flood and causing one fatality was 25.3 mm in 2007 in the Sidi Jazouli commune. Based on this threshold, Table 6 below presents the number of days when precipitation exceeded 25.3 mm, which may or may not lead to a flood. Table 1, along with the factor map (Figure 6), indicates that the Adamna station, which is the closest to Essaouira city, recorded the highest number of such occurrences, followed by Igrounzar, which is located upstream. Consequently, the value of weighted rainfall also increases. The average annual rainfall ranges from 161 to 319 mm, while the weighted rainfall varies from 199 to 406 mm with decreasing values from west to east. The areas around Essaouira city exhibit the highest values of weighted rainfall due to both the elevated number of occurrences and the influence of the oceanic effect on the precipitation distribution.

Figure 6. Factors used in the susceptibility analysis for Essaouira province.

Figure 6. Factors used in the susceptibility analysis for Essaouira province.

Table 6. Mean rainfall and weighted rainfall values for different meteorological stations.

Station	Mean Rainfall	Number Values >25.3 (1984–2021)	Weighted Rainfall
Abadla	161.91	20	199.29
Adamna	319.58	101	406.96
Chichaoua	176.25	32	213.66
Igrounzar	286.18	87	373.07
Talmest	254.04	71	360.39

Table 6. Mean rainfall and weighted rainfall values for different meteorological stations.

Station	Mean Rainfall	Number Values >25.3 (1984–2021)	Weighted Rainfall
Abadla	161.91	20	199.29
Adamna	319.58	101	406.96
Chichaoua	176.25	32	213.66
Igrounzar	286.18	87	373.07
Talmest	254.04	71	360.39

• Lithology (permeability level): The analysis of the geological map of the Essaouira basin reveals the predominance of Quaternary and Cretaceous formations, along with some Triassic magmatic formations. Based on the attributed permeability classes and their respective percentages, the following categories were identified: very low permeability (5%) for “Clays” and “Marls”, low permeability (23%) for “Conglomerates and Consolidated Dunes” and “Doleritic Basalt”, moderate permeability (58%) for “Dolomitic Limestones”, “Limestone, Yellow Marl, and Sandstone”, and “Lumachelic Clayey Limestones”, high permeability (12%) for “Colluviums/Medium Terrace”, “Sandstone and Red Clay”, and very high permeability (2%) for “Barkhanes” and “Mobile Dunes”. As illustrated in Figure 6, the areas with high and very high permeability classes (14%) are located in the river valleys (oueds) and the coastal dunes.

In the AHP analysis, a rating of 5 was assigned to very low permeability, as impermeable formations tend to exacerbate flooding. On the other hand, alluvial and coastal deposits received a lower rating due to their higher capacity for infiltration, facilitating groundwater flow and reducing the severity of the floods.

• Topographic wetness index (TWI): The topographic wetness index (TWI) serves as an indicator of areas with potential for flooding and is utilized in predicting the occurrence of saturation excess overland flow. For Essaouira province, the TWI was classified into five classes, displaying similar variations to the elevation and slope factors since they were derived from the same data source.
• Curvature: the curvature factor was classified into five classes, as depicted in Figure 6. Classes 4 and 5 (ranging from −0.14 to >0.51) corresponding to flat terrain (with a value of around 0) tend to promote the occurrence of floods. This leads us to consider urban floods, which occur when city or village landscapes are unable to absorb excess water following prolonged periods of intense rainfall, river overflow, or storm surge. Such urban floods have occurred numerous times over the years in Essaouira city and in villages like Tafedna, Tamanar, Smimou, El Hanchane, Lagdadra, and Ounagha.

4.3. Susceptibility Modeling and Validation

Once the conditioning factors were chosen and standardized with consistent cell sizes, an analysis was conducted in ArcGIS Pro software to generate the ultimate flood susceptibility maps for the three models. To create a flood susceptibility rating map, the map needed to be divided into distinct categories. In this research, the flood risk was classified into five categories: very low, low, medium, high, and very high, employing the natural breaks method.

4.3.1. Analytic Hierarchy Process (AHP)

After the application of weighting to all conditioning factors, the final flood susceptibility map was generated. The map was classified using the natural breaks (Jenks) grading method, as illustrated in Figure 7. In the context of the AHP method, this classification approach is considered the most suitable for delineating the flood susceptibility zones [123]. The respective areas of each flood susceptibility zone and the proportion of flood events associated with each zone are presented in

Table 7. The percentages and covered areas of the flood susceptibility model classes.

	AHP		FR		WoE
Susceptibility Classes	Area (km2)	Percentage (%)	Area (km2)	Percentage (%)	Area (km2)	Percentage (%)
Very low	0.45	0.01	5148.87	82.77	1608.53	25.86
Low	1980.53	31.84	837.18	13.46	2921.66	46.97
Moderate	4147.67	66.67	202.69	3.26	1124.74	18.08
High	92.12	1.48	24.36	0.39	490.95	7.89
Very high	0.00	0.00	7.68	0.12	74.90	1.20

.

Regarding the distribution of flood susceptibility classes in terms of area percentage, the lowest share was observed in the very high (0%) and very low (0.01%) classes. The high flood susceptibility class covers an area of 1.48%, while the low susceptibility class represents 31.84% of the study area. The largest share (66.67%) is attributed to the moderate flood susceptibility class. Consequently, based on the incorporated conditioning factors and their corresponding weights, 1.48% of the Essaouira province is identified as having very high to high flood susceptibility according to the AHP analysis.

When analyzing the factors, it is evident that elevation has the significant impact, on the occurrence and spread of flooding in the study area with a weight of 15.80. Elevation plays a role in determining how floods move and their depth. It also affects how the extent of floods is controlled. On the one hand, the slope factor carries a weight of 11.43. The slope of an area influences the surface runoff and water infiltration. Areas with low slopes tend to experience periods of flooding, while those with moderate-to-steep slopes facilitate the faster drainage of floodwater.

Regarding the spatial distribution, the AHP analysis results indicate that there are no areas in the Essaouira province classified as very highly susceptible to flood occurrence. However, there are areas with high susceptibility to floods. These include Oued Ksob near Diabat village and Essaouira city, Guazoua village located in the south of Essaouira, Oued Tamanar, Oued Smimou, Had Dra village, and Souira Guedima city situated near the outlet of the major Tensift river.

4.3.2. Frequency Ratio Model (FR)

Based on the results of the frequency ratio (FR) model, it was determined that 0.51% of the study area in the Essaouira province is categorized as having high to very high susceptibility to flood occurrence, as shown in

Table 8. Spatial relationship between each conditioning factor and flooding extracted by FR and WoE method.

Factor	Class		Frequency Ratio Index (FR)	Weight of Evidence (WoE)
Slope (°)	1	0–2	3.140	1.144
	2	2–5	1.563	0.446
	3	5–15	0.258	−1.353
	4	15–35	0.072	−2.632
	5	>35	0.070	−2.634
Elevation (m)	1	0–80	7.298	1.988
	2	80–350	0.701	−0.356
	3	350–600	0.296	−1.216
	4	600–950	0.070	−2.634
	5	950–1697	0.070	−2.634
Distance from river (m)	1	0–100	6.240	1.831
	2	100–200	1.624	0.485
	3	200–1000	0.475	−0.744
	4	1000–2000	0.488	−0.717
	5	2000–8166	0.070	−2.634
Drainage density (km/km2)	1	0–0.55	1.662	0.508
	2	0.55–1	1.135	0.127
	3	1–1.47	0.950	−0.051
	4	1.47–1.92	1.074	0.071
	5	1.92–2.38	0.070	−2.634
Stream order (flow accumulation)	1	Order 1	0.634	−0.456
	2	Order 2	4.283	1.455
	3	Order 3	16.292	2.791
	4	Order 4	14.191	2.653
	5	Order 5	36.101	3.586
	6	Order 6	10.425	2.344
Land use	1	Bare/spare vegetation	0.398	−0.921
	2	Built-up	49.693	3.906
	3	Dense vegetation	0.972	−0.028
	4	Medium vegetation	0.133	−2.019
	5	Water	1.176	0.162
Weighted rainfall	1	269–290	0.070	−2.634
	2	290–330	0.254	−1.371
	3	330–365	0.373	−0.985
	4	365–380	1.047	0.046
	5	380–407	3.666	1.299
Lithology (permeability level)	1	Very low	0.500	−0.693
	2	Low	0.174	−1.748
	3	Moderate	0.257	−1.358
	4	High	3.177	1.156
	5	Very high	18.888	2.939
TWI	1	<7	0.665	−0.407
	2	7–9	1.297	0.260
	3	9–12	2.604	0.957
	4	12–16	3.143	1.145
	5	>16	0.070	−2.634
Curvature	1	<−2.46	0.410	−0.891
	2	−2.46–−1.13	1.171	0.157
	3	−1.13–−0.14	1.207	0.188
	4	−0.14–0.51	1.148	0.138
	5	>0.51	0.722	−0.325

. Conversely, over 96% of the province exhibits very low to low susceptibility to floods.

Analyzing the factors (Table 8), the land use factor, specifically the built-up class, was found to be the most influential factor in the occurrence of floods. The presence of built-up areas hinders water infiltration, leading to increased surface runoff and higher flood susceptibility. However, the fourth class of slope (15–35°) was identified as the least influential class for flood occurrence, as water cannot easily stagnate in this slope. As a result, a value of 0.07 was assigned to classes with no intersection of the inventoried points, as indicated in the gray-colored table.

For the slope factor, the first class (0–2°) had the highest weight of 3.140, while the class of 15–35° had the lowest weight of 0.072. In terms of elevation, the range between 0 and 80 m exhibited the highest weight of 7.298, indicating high flood susceptibility within this elevation range. Distance from the river, which can range from 0 to 100 m, displayed a relationship with the likelihood of flooding. On the other hand, areas located beyond 1000 m from the river showed a tendency to have an effect on flood incidents. These findings suggest that floods primarily happen near the edges of rivers and occasionally extend into regions.

Regarding drainage density, the range of 0–0.55 had the highest weight (1.662), while the range of 1–1.47 had the lowest weight (0.950). For stream order, represented by flow accumulation, the fifth class (order 5) had the highest value of 36.101, while the first order (order 1) had the lowest weight of 0.634.

In terms of the land use factor, the built-up class, encompassing roads and constructions, exhibited the highest value of 49.693, indicating maximum flood susceptibility. Dense and medium vegetation acquired weights of 0.972 and 0.133, respectively, suggesting a negative correlation with flood occurrence. Vegetated areas can reduce runoff and mitigate flooding. In the case of rainfall, the highest value of 3.666 was observed in the fifth class (380–407 mm), which is expected since it represents the highest weighted precipitation value. Regarding permeability, the classes of high and very high permeability associated with sand and oueds’ alluviums exhibited the highest values of 3.177 and 18.88, respectively. This lithology indicates a higher flood susceptibility in the study area, as it is typically found in areas where oueds (seasonal rivers) are present.

In terms of the topographic wetness index (TWI), the range between 12 and 16 had the highest weight of 3.143, indicating a high flood susceptibility within this TWI range. Additionally, a flat curvature had the highest value of 1.207, suggesting a high flood susceptibility in flat areas. This aligns with the natural characteristics of floods, which tend to occur more frequently in flat regions.

Analyzing the spatial distribution (as shown in Figure 8), the very high class, covering an area of 7 km², is distributed across several locations. This includes the city of Essaouira, where urban floods have been reported; Diabat village, located near the left limit of Oued Ksob Oued; Smimou village; and the outlet of Oued Ouazzi. The high class, encompassing an area of 24 km², is primarily visible in Tafedna village along the Igouzoulen outlet; the Tidzi and Had Dra villages; as well as the communes of El Hanchane and Lagdadra, which are situated in relatively flat areas. The remaining area, predominantly at higher altitudes and upstream of the watersheds, is occupied by the very low to low flood susceptibility classes.

4.3.3. Weights of Evidence (WoE)

In the previous section, all the parameters were calculated for each conditioning fac-tor, representing the relationship between the classes of each factor and flood occurrence. Analyzing the final results of the weights of evidence (WoE) in Table 8, it is evident that the slope with the built-up range of the land use factor had the highest weight among all other factor classes. This indicates that areas with built-up land use exhibit maximum susceptibility to flooding in the catchment. Conversely, the lowest value among all factors was −2.632. Therefore, a value of −2.634 was assigned to classes with no intersections with the inventoried points (colored in gray in the table).

For the slope factor, the first class (0–2°) had the highest weight of 1.14, while classes 15–35 acquired the lowest weight of −2.632. In terms of elevation, the range between 0 and 80 m had the highest weight, indicating high flood susceptibility within this elevation range. The distance from the river, ranging from 0 to 100 m, showed a positive correlation with flood occurrence, while areas that were more than 1000 m away from the river had a negative influence on flooding. This suggests that flooding mostly occurs near the riverbank and occasionally at a distance from the river. Regarding drainage density, the highest weight (0.508) was assigned to the range of 0–0.55, while the lowest weight (−0.051) was for the range of 1–1.47. In the case of the stream order representing the flow accumulation, the fifth class (order 5) had the highest value of 3.586, while order 1 acquired the lowest weight of −0.456.

In terms of the land use factor, the built-up class, which includes roads and constructions, had the highest value of 3.906, indicating maximum flood susceptibility. Dense and medium vegetation acquired weights of −0.028 and −2.019, respectively, showing a negative correlation with flood occurrence. This suggests that vegetated areas can reduce runoff and consequently decrease the risk of flooding. For rainfall, the highest value was observed in class 5, ranging from 380 to 407 mm, which is expected as it represents the highest weighted precipitation.

Among the different degrees of permeability, the class of high and very high permeability associated with sand and oueds’ alluviums had the highest values of 1.156 and 2.939, respectively. This lithology represents the highest flood susceptibility in the study area, as it is typically found in areas where oueds are present. In the case of the topographic wetness index (TWI), the range between 12 and 16 had the highest weight, indicating high flood susceptibility within this TWI range. Similarly, the flat curvature had the highest value of 0.188, indicating high flood susceptibility in flat areas. This aligns with the natural characteristics of floods, which predominantly occur in flat regions.

Based on the visual interpretation of flood-susceptible locations (refer to Figure 9), the area of 74 km2 classified as “very high” susceptibility is distributed across several regions. These include the big Oued Tensift, its outlet, and its affluent of Lahrichate. The Essaouira city, where urban floods were observed, which also falls within this category. Additionally, Diabat village, located near the left limit of Oued Ksob, as well as Smimou village and the outlet of Oued Tafedna are part of this high susceptibility area.

On the other hand, the area of 490 km2 classified as “high” susceptibility is predominantly visible upstream of Tafedna village along the Oued Igouzoulen. It also encompasses Oued Tidzi, Had Dra village, and the drains near the communes of El Hanchane and Lagdadra communes, along with the outlet of Oued Ouazzi.

The remaining classes, ranging from “very low” to “low” susceptibility, cover the remaining area, particularly at higher altitudes and in the upstream regions of the watersheds.

4.3.4. Validation (AUC-ROC)

However, visual assessment alone is not sufficient for comparison and judgment, necessitating the use of a proper validation method. To validate the resulting probability maps, the area under the curve (AUC) method was employed, and both the success rate and prediction rate curves were measured (refer to the Figure 10). AUC is widely recognized as a popular method for evaluating the effectiveness of the generated method, as it provides both success and prediction rates [155].

The validation process involved assessing the presence of flood locations (both training and testing) within each interval and measuring the resulting success and prediction rates. Spatial confrontation was carried out to validate all the flood susceptibility models, using independent flood partitions as validating subsets. The receiver operating characteristic (ROC) curves, as proposed by Linden (2006) and Remondo et al. (2003) [89,90], were computed for the predictive models, and the corresponding AUC values were calculated.

Figure 10 illustrates the AUC values obtained during the validation process for all models. It is worth noting that all landslide susceptibility models with AUC values exceeding 0.9 were considered to be outstanding, particularly the AHP model.

4.4. HEC-RAS Analysis

According to the results of the AHP, FR, and WoE models, one of the most susceptible areas was the vicinity of the Oued Ksob outlet, which also happened to be the most vulnerable area due to its high population density. The selection of the Oued Ksob reach was specifically chosen for the downstream section of Oued Ksob, extending from the Adamna gauging station to the outlet in the extreme west near Essaouira city. This particular sector exemplifies both the highest likelihood of hazard occurrence and the presence of valuable assets (Figure 7, Figure 8 and Figure 9).

The flood hazard mapping was conducted using HEC-RAS and the GIS environment. The analysis of the flood frequency utilized the maximum instantaneous annual discharge data for the Oued Ksob, as presented in Table 1.

Table 9. Result of the best-fitted distributions using the AIC and BIC criteria of the Ksob watershed (P(Mi) = a priori probability; P(Mi|x) = a posteriori probability (method of Schwarz); AIC = Akaike information criterion; BIC = Bayesian information criterion).

	NB	XT	P(Mi)	P(Mi|x)	BIC	AIC
Weibull	2	2582.145	14.29	72.04	732.778	728.875
Gamma	2	2537.99	14.29	26.61	734.77	730.867
Lognormal	2	7327.21	14.29	1.35	740.728	736.826
Pearson type 3	3	2474.973	14.29	0	758.398	752.544
GEV	3	2350.298	14.29	0	774.425	768.571
Gumbel	2	2120.29	14.29	0	783.641	779.738
Normal	2	1679.527	14.29	0	809.607	805.704

displays the results of the Akaike information criterion (AIC) and Bayesian information criterion (BIC) for each selected distribution. The quality assessment of the frequency distribution relies on the values of these criteria. By comparing the values provided by AIC and BIC, the most reliable distribution for the adjustment can be determined among the seven distributions. The results indicate that the Weibull distribution is the most suitable statistical distribution for the studied variables (Figure 11). It exhibits the lowest values for both criteria (AIC = 728.875; BIC = 732.778), signifying the best fit.

Furthermore, Table 9 presents the frequency analysis of the extreme flows calculated using the Weibull statistical model (Figure 11), resulting in a value of 2582 m³/s for a 100-year return period. These large discharges of water are considered highly destructive, particularly in a vulnerable environment like the Oued Ksob region.

It is important to note that the flow corresponding to 2582 m³/s for the 52-year flood is relatively high considering that the basin size is only 1764.30 km². This is a characteristic of hydrology in arid zones, where upstream basins generally contribute more water than downstream basins, although there are specific conditions in the Ksob basin where the downstream area receives more precipitation. As the slopes decrease, floodwater spreads and infiltrates the alluvial terraces and mobile dunes of Essaouira, resulting in a decrease in floodwater levels at the downstream hydrological stations.

The selection of Manning’s coefficient value used for floodplain mapping was based on a sensitivity analysis of this parameter in relation to the resulting floodplains. A simulation was conducted using the HEC-RAS software (version 6.3.1) for a 100-year return period.

The flood maps shown in Figure 12a,b provide valuable information regarding floodwater levels, flood amplitude, and floodplain areas for the 100-year return period, with and without considering obstructions (bridges and some constructed anti-flood dikes or longitudinal levees for defenses reasons). These maps are essential for visualizing and monitoring flood events, allowing real-time flood alerts, and ultimately reducing the potential loss of human lives. In both Figures, the blue-colored areas indicate locations within the study area that would be submerged at a discharge of 2583 m³/s, highlighting areas that are susceptible to damage. Figure 12 further demonstrates that the extent of endangered areas in the high-hazard class increases from upstream to downstream areas, primarily due to the decreasing slope of the land in the downstream direction, which leads to the greater accumulation of floodwater and consequently higher levels of risk.

The flood mapping in the study area shows that the actual hazard levels range from 0 to 10.2 m, with flooded areas primarily located along the Oued Ksob. The floodwater level is higher near the riverbanks, reaching depths of up to 7 m and causing severe damage. Even floods with depths of less than 0.8 m can still result in significant damage. Near the outlet, higher values are observed due to the geomorphological characteristics of the terrain, characterized by lower depressions that facilitate the extensive lateral expansion of floods. Furthermore, the digital elevation model (DEM) map reveals that, in almost all the right-side areas before bridge 2 of the Oued Ksob, the elevation is higher than the left side. Consequently, the low-hazard areas (in terms of severity) on the right side of the oued are fewer than on the left side, and this distinction is accurately reflected in the flood hazard map. Figure 12A,B also show that the floodwater extends further into the river sides, particularly in the downstream direction, as a result of the slope and flow velocity.

Without obstructions (Figure 12A), these floods could reach the city of Essaouira. The dune located next to the Oued Ksob outlet acts as a natural obstruction to the surface water runoff and flow into the Atlantic. This natural barrier leads to water retention during heavy rainfall events and contributes to the occurrence of flash flood hazards. Furthermore, the city of Essaouira experiences strong winds that transport sand grains to the outlet, which can accumulate in the water and alter the flow direction, as observed in previous years. Examining the historical course of the Oued Ksob reveals that, prior to the 1950s, it followed three channels: one passing by the lighthouse and two heading west. Currently, the Oued Ksob follows a single bed due to the construction of a gabion wall and the stabilization of the dunes.

The simulated floodplain depths for the 100-year return period indicate significant overflows from both riverbanks, as observed on the southern side near Diabat village. This area is particularly vulnerable to flood inundation, even with the presence of obstructions. It is evident from Figure 12B that many houses and agricultural fields are highly exposed to flooding in that part, and the water flows behind the obstruction, as seen in Figure 13 of station 1424, where the constructed obstruction is visible on the left side.

The agricultural area, residences, and roads near the Oued Ksob in Diabat village are the most extensively impacted by floods. The river continues to pose a hazardous threat during severe flooding events, as seen in the case of the 2009 flood. Understanding these patterns can facilitate better stream restoration planning to minimize losses and damages during future flood events.

To validate the accuracy of the simulation results, Figure 14 presents the areas within the study area that were affected by previous floods. For instance, Figure 14a illustrates the flooded area resulting from the 2009 flood event with a discharge of 2270 m³/s. This figure demonstrates that our modeling results align with the observed extent of the flood, with the exception that there were no obstructions present at that time.

5. Discussion

The utilization of AHP, FR, and WoE models allows for the classification of flood susceptibility levels without indicating the actual water depth. As a result, the province of Essaouira has been categorized into five groups representing varying degrees of susceptibility, ranging from very low to very high hazard (Figure 15). Among all the models, the areas classified as having very high and high susceptibility to flooding occupy a smaller proportion compared to the other classes. Conversely, the moderate, low, and very low susceptibility classes cover larger percentages of the study area.

In general, the severity of the flood susceptibility is particularly significant in certain localities, including Essaouira, Smimou, Sidi Jazouli, Tafedna, Tamanar, El Hanchane, Lagdadra, Had Dra, Tidzi, Zaouia Ben Hmida, Imintlit, Oulad Mrabet, and Mejji. These areas are characterized by their low-altitude plains and gentle slopes. The way hazards are spread out shows that there are agricultural and infrastructure areas like roads and bridges that face significant risks and could be seriously damaged. Proximity to the river significantly increases the exposure to hazards, with the magnitude of the risk ranging from high to very high.

Upon analyzing the results obtained from comparing the selected factors, it is evident that the three most influential factors in determining flood risk and contributing to the occurrence of floods, as identified by all the models, are the slope, elevation, and stream order (flow accumulation).

The area around the Ksob outlet near the province of Essaouira is identified as the most vulnerable region based on the results obtained from the three models and considering the population density. In order to further assess this area, we conducted an analysis using HEC-RAS. Remarkably, the outcomes of the three models closely align with the findings of the HEC-RAS analysis, which indicates the models’ effectiveness in mapping flood-prone areas in different regions, even in the absence of specific hydrological data from measurement stations.

The flood map (Figure 12) illustrates the extent of water resulting from the centennial flood, indicating the presence of protective measures. The embankments effectively prevent overflow onto the major bed of the Oued Ksob and subsequently towards the city of Essaouira. These embankments are situated on both the right and left banks between the footbridge and bridge 2, and only on the right bank downstream.

According to the Tensift Hydraulic Basin Agency and the local residents we have contacted, the Oued Ksob no longer causes significant flooding issues in the city of Essaouira. However, there are occasional problems of bank erosion in certain areas. The two Chaabats_minor water tributaries that conduct water to the wadi or oued_that pass through the boundaries of the urban development plan converge and create an obstacle to mobility during periods of heavy rainfall. The other watercourses within the city are part of the internal drainage system. The flooding problems experienced in Essaouira are primarily due to the combination of urban fringe basins’ runoff and the inadequacy of the urban drainage system (in terms of limited capacity, high tide, swell, low altimetric calibration, etc.).

Nevertheless, overflow from the right bank of the Oued Ksob also contributes to the flooding. The sandy nature of the riverbed and its banks, combined with the velocity of water during floods, cause varying degrees of bank erosion and significant accumulation of vegetation, which aids in sand fixation.

It is worth noting that new constructions have been established along the banks of the Oued Ksob, particularly at the entrance of the city of Essaouira. Moreover, some of these new constructions encroach upon the major bed of Oued Ksob (Figure 16) in the village of Diabat. This development occurred after 2007, as aerial photos and the historical imagery from Google Earth indicate the absence of these constructions at that time. The problem is that these constructions are the most exposed to flooding according to the results of our models and simulations.

Fieldwork observations of the water levels during various flood events, particularly the 2009 flood, along with surveys conducted among the local population, were carried out to spatially validate the extent of flood areas in the province of Essaouira and at the Ksob outlet. The results obtained from these surveys aligned with the findings of the models with some minor differences. These results can be utilized for future land-use planning in the Essaouira basin. However, it is worth mentioning that the models did not offer any details regarding the depth of the water in those regions. This limitation restricted our ability to evaluate the level of danger and how much infrastructure might be submerged.

Particular attention should be given to the high flood risk in Diabat village, especially in residential, agricultural, and infrastructure areas situated within the main riverbed in order to prevent disasters resulting from severe flood events. It is recommended that hydrological data measurements be considered by a regional intergovernmental committee, which could issue warnings and facilitate the evacuation of people, livestock, and property in the event of high floods, thereby minimizing the associated consequences.

Installing an early warning system to provide alerts during heavy rainfall is crucial to mitigate the damage caused by rapid river section filling. Additionally, in the develop-mental plans for Essaouira province, it would be prudent to construct a dike at the upper limit of Diabat village, complementing the existing one downstream to mitigate the impact of floods originating from the upstream basin. This infrastructure would help reduce the velocity of floodwater and prevent water entering Diabat village.

Furthermore, it is recommended that slope development measures are implemented, such as the construction of benches, protective sills, and dikes, along both tributaries and the main river. The purpose of these measures is two-fold: to control water erosion and to reduce surface runoff by promoting infiltration.

To address the geomorphological characteristics of the region, it is essential to establish the flood risk management strategies (FRMSs) formulated by local authorities. These strategies should define the prohibited and regulated zones for construction purposes. The implementation of the FRMSs will require action to reduce the vulnerability of existing infrastructures. The main challenges ahead are primarily financial, as current funding is primarily directed toward large basins and agricultural plains in other regions of Morocco.

6. Conclusions

In conclusion, the utilization of the AHP, FR, and WoE models has enabled the classification of flood susceptibility levels in the Essaouira province without explicitly indicating water depth. The province has been categorized into five susceptibility groups, representing varying degrees of hazard. While the areas with high and very high susceptibility occupy smaller proportions, the moderate, low, and very low susceptibility classes cover larger percentages of the study area.

The overlay between flood susceptibility with highly exposed areas is particularly significant in certain settlements located in low-altitude plains and gentle slopes. These areas, including Essaouira, Smimou, Sidi Jazouli, Tafedna, and others, are exposed to high levels of susceptibility_n some areas, aggravated by the high level of vulnerability of social and buildings’_putting residential, agricultural, and infrastructure areas such as roads and bridges in danger of significant damage. Proximity to the river greatly increases the exposure to flooding as the historical record of flood losses demonstrate.

The analysis of the selected factors across all models consistently highlighted the importance of slope, elevation, and stream order (flow accumulation) in determining flood risk and contributing to flood occurrence.

Furthermore, the area around the Ksob outlet near the province of Essaouira emerged as the most vulnerable region based on the results of the three models and population density. The outcomes of the models were closely aligned with the findings from the HEC-RAS analysis, demonstrating the effectiveness of these models in mapping flood-prone areas even without specific hydrological data.

To address these findings, it is crucial to pay special attention to the high flood risk in Diabat village, particularly in residential, agricultural, and infrastructure areas within the main riverbed. Measures such as establishing a regional intergovernmental committee to consider hydrological data, installing an early warning system, and constructing additional barriers in Diabat village could help mitigate the impact of severe flood events.

Additionally, implementing slope development measures, defining prohibited and regulated zones, and establishing flood risk management strategies (FRMSs) formulated by local authorities are recommended to reduce vulnerability and effectively manage flood risks. However, financial challenges may arise in implementing these measures, as funding is currently focused on other regions in Morocco.