Determining Flood Zonation Maps, Using New Ensembles of Multi-Criteria Decision-Making, Bivariate Statistics, and Artificial Neural Network

Abstract

Golestan Province is one of the most vulnerable areas to catastrophic flood events in Iran. The flood severity in this region has grown dramatically during the last decades, demanding a major investigation. Accordingly, an authentic map providing detailed information on floods is required to reduce future flood disasters. Three ensemble models produced by the combination of Evaluation Based on Distance from Average Solution (EDAS) and Multilayer Perceptron Neural Network (MLP) with Frequency Ratio (FR), and Weights of Evidence (WOE) are used to quantify the map flood susceptibility in Golestan Province, in the north of Iran. Ten flood effective criteria, namely altitude, slope degree, slope aspect, plan curvature, distance from rivers, Topographic Wetness Index (TWI), rainfall, soil type, geology, and land use, are considered for the modeling process. The flood zonation maps are validated by the receiver operating curve (ROC). The results show that the most precise model is MLP-FR (AUROC = 0.912), followed by EDAS-FR-AHP (AUROC = 0.875), and EDAS-WOE-AHP (AUROC = 0.845). The high accuracies of all methods applied to illustrate their capability in predicting flood susceptibility in future studies.

1. Introduction

Among climate-related hazards, flood is considered to have the largest impact on humans and the environment, globally [1]. In the last decades, most cities have experienced a population explosion at an alarming rate, leading to a dramatic increase in flood potential [2]. Population growth and, therefore, urbanization results in the change of agricultural land and natural vegetation into mainly impervious built-up environments, leading to higher storm-water runoff and flooding [3,4]. Flood hazards are expected to increase in frequency and severity in many parts of the world through the impacts of global change on climate, severe weather in the form of heavy rains, and river discharge conditions [5]. For all of the above reasons, it is crucial to perform flood hazard, vulnerability, and risk assessments and incorporate the results in urban planning and management to mitigate future risks [6]. Producing flood zonation maps to diagnose at-risk areas is one of the primary steps for flood hazard mitigation [7].

Iran is one of the most flood-prone countries in the world, owing to its geological and geographical setting. Characteristics of the rivers and topographic conditions of Iran, especially in the north of the country, greatly contribute to the extent and intensity of floods [8]. Golestan Province in the north of Iran has recently experienced casualties in human lives and livelihoods due to severe floods [9,10]. Hence, it is of immense significance to identify flood-prone areas as the most vital preventive measure.

Geographic Information System (GIS) technologies have made significant contributions to the analysis of natural hazards [4,11]. Using a flood inventory map produced by flooding location data is a popular way of generating flood susceptibility maps. Flood inventory maps are recognized as the crucial tools for understanding the relationship between flood occurrences and flood influential factors [10]. In the past, various kinds of databases have been used to diagnose flood-affected locations [12] in the literature, such as governmental sources, field surveys [13], remote sensing, or aerial photographs [14].

Over the past few decades, researchers have developed various methods and models for flood hazard mapping, including bivariate statistics [12,13,14,15], multi-criteria decision-making (MCDM) [16,17,18], and neural networks [19,20,21,22]. MCDM methods are popular for flood vulnerability assessment in literature due to their relative simplicity and their capacity to optimize decision-making considering diverse types of data [16,23]. One of the most widespread MCDM methods is the Analytical Hierarchy Process (AHP), which is applied in a wide range of studies [17,24,25]. EDAS (Evaluation Based on Distance from Average Solution) is one of the latest methods of MCDM [26], in which its capability in flood prediction has not been investigated previously. Moreover, the accuracy of the MLP method has been reported in several research studies in different areas [27,28,29].

The objective of this study is to produce and validate flood susceptibility maps of Golestan Province based on the factors affecting river flooding. Since employing only one technique may not reveal the performance in flood susceptibility assessment and ensemble models have shown higher accuracy in previous research [22,30,31], two ensembles of EDAS, namely EDAS-FR-AHP and EDAS-WOE-AHP are evaluated in this study.

In a study by Costache et al. (2021) in small river catchments, an ensemble of FR and WOE methods with a fuzzy analytical hierarchy process (FAHP) was used. The authors reported that the FAHP-WOE and FAHP-FR models had better performance than WOE and FR models [31]. In another research study by Wang et al. (2020), the ANN, AHP, and FR methods were applied as both standalone and ensemble models for determining the flood potential. The researchers concluded that all ensemble models had higher accuracy than standalone models [29].

Considering that the Golestan Province experiences severe floods, the application of new MCDM methods would be rational and necessary to find a promising method for flood prediction in this region. An MLP-based model (MLP-FR) is also applied to compare the results. The application of novel ensemble models for flood susceptibility prediction is the main novelty of the present study.

2. Case Study

Golestan Province is located in the north of Iran with an area of approximately 20,438 ${\mathrm{km}}^2$ (Figure 1). It lies between the latitudes of 36°27′ to 38°14′ N and the longitudes of 53°40′ to 5°30′ E, and altitude of −40 to 3820 m. Golestan Province has a combination of mountainous and temperate Mediterranean climate with an average precipitation of 470 mm. This area largely consists of Cenozoic rocks, followed by Mesozoic, Paleozoic, and Proterozoic, in respective order. Furthermore, it is covered by mountainous regions, jungles, gardens, and coastal areas.

3. Materials and Methods

3.1. Research Overview

In the present study, 10 flood conditioning factors were selected, and the flood inventory map was achieved according to flood and non-flood locations. For the next stage, 70% of data were selected for the modeling process, and flood susceptibility maps of three ensemble models were produced. Besides, the models’ performance was determined by the remaining 30% of the data. As a result, the most accurate model was determined. Figure 2 illustrates the hierarchy diagram of the applied methodology, in which the details of each step are shown.

3.2. Flood Inventory Map

In order to predict the future flood occurrences in an area, it is necessary to have records of the past flood occurrences in that region [32]. In the current research, a total of 240 past flood locations were obtained from the Golestan Water organization from 2014 to 2019. Moreover, 240 non-flood locations were randomly selected in high-elevation regions with a low possibility of flooding. According to several research studies, these locations were randomly classified into two samples of training (70% of data) and testing (30% of data) [33,34,35,36] (Figure 1).

3.3. Flood Influential Criteria

In order to identify the flooding risk for any region, it is necessary to recognize the most probable hazard forming and vulnerability indicators. Besides, relative weightage should be given prudently to the chosen indicators by considering their contributing role in aggravating the flood risk. In the present study, 10 flood influential factors were selected (Figure 3) [5,19,37]. For simulation, all effective criteria are converted to a 30 × 30 m resolution raster. Flood criteria were classified based on the Natural Break technique [6,14]. Altitude and slope degree are significant factors, which have an inverse relationship with flooding. The slope aspect has an impact on hydrologic processes, rainfall direction, and evapotranspiration [38]. Curvature indicates the ground shape and is divided into three classes, including concave, flat, and convex. Distance from rivers influences flood magnitude and frequent lateral migration [39]. TWI shows the flow propensity for going downslope due to gravity (Equation (1)) [40] and is calculated as follows:

$$\mathrm{TWI}=ln\left(\frac{A_s}{\mathsf{\beta }}\right)$$

(1)

where A_s is the specific area ($\frac{m^2}{m}$) and β is the slope angle (degree).

Rainfall is the chief source of surface runoff. The intensity and volume of channel discharge are largely governed by rainfall. In this study, 10 years of meteorological data (2009–2019) from hydrometric stations were used. Types and sources of all applied factors are presented in

Table 1. Sources and types of the data used.

Dataset	Source	Data Type	The Scale of Source Data	Derived Factors
Digital elevation model (DEM)	United States Geological Survey (USGS) site	Raster	1:25,000	Altitude, slope, aspect, curvature, distance from rivers, TWI
Rainfall	Golestan meteorology organization	Vector	1:25,000	Rainfall map
Geological map	Golestan Regional Water Authority	Vector	1:100,000	Lithology, soil type
Land cover	Golestan Regional Water Authority	Vector	1:100,000	Land use

.

3.4. Feature Selection

Determination of the importance of influential factors in flood modeling is a necessary step. In this study, before the modeling process, Information Gain (IG) test was used to find the most effective influential factors in flood susceptibility determination [41]. Equation (2) shows the IG value formulation.

$$\mathrm{IG}\left(x·H\right)=\frac{Entropy\left(H\right)-{{\displaystyle \sum }}_{i=1}^n\frac{\left|H_i\right|}{\left|H\right|}Entropy\left(H_i\right)}{-{{\displaystyle \sum }}_{i=1}^n\frac{\left|H_i\right|}{\left|H\right|}log\frac{\left|H_i\right|}{\left|H\right|}}$$

(2)

where x is each attribute that belongs to dataset H with subsets H_i. This method helps to increase the prediction power and model’s efficiency [42]. The value of IG varies from 0 to 1, in which values close to 1 indicate a high prediction of flood intensity factors, while a value of 0 means that flood influential factors do not affect flood occurrence.

3.5. Flood Susceptibility Models

3.5.1. Frequency Ratio

Frequency ratio (FR), as a widely used bivariate statistical model, determines the probabilistic relationship between flood occurrences and different variables [14,43,44]. FR for each class of various factors is determined by Equation (3) [45].

$$FR=\frac{\raisebox{1ex}{$A$}\!\left/ \!\raisebox{-1ex}{$B$}\right.}{\raisebox{1ex}{$C$}\!\left/ \!\raisebox{-1ex}{$D$}\right.}$$

(3)

where A is the number of flood pixels per class; B is total flood pixels; C is the number of pixels per class, and D is the total pixels of the study area.

3.5.2. Weights of Evidence

Weights of Evidence (WOE) determines the relationship between various criteria and flood occurrences [11,46,47,48]. This method calculates the weights of each class as follows (Equations (4) and (5)):

$${\mathrm{W}}^{+}=\ln \frac{\mathrm{P}\left(\mathrm{B}|\mathrm{A}\right)}{\mathrm{P}\left(\mathrm{B}|\overline{\mathrm{A}}\right)}$$

(4)

$${\mathrm{W}}^{-}=\ln \frac{\mathrm{P}\left(\overline{\mathrm{B}}|\mathrm{A}\right)}{\mathrm{P}\left(\overline{\mathrm{B}}|\overline{\mathrm{A}}\right)}$$

(5)

where B and $\overline{\mathrm{B}}$ are indicative of the presence and absence of every criterion, in respective order. in addition, A and $\overline{\mathrm{A}}$ are reagents of the presence and absence of flood events. The standard deviation is determined by Equation (6).

$$\mathrm{S}\left(\mathrm{C}\right)=\sqrt{{\mathrm{S}}^2{\mathrm{W}}^{+}+{\mathrm{S}}^2{\mathrm{W}}^{-}}$$

(6)

where ${\mathrm{S}}^2{\mathrm{W}}^{+}$ and ${\mathrm{S}}^2{\mathrm{W}}^{-}$ are variances of positive and negative weights, computed using the Equations (7) and (8), as the following:

$${\mathrm{S}}^2{\mathrm{W}}^{+}=\frac{1}{\mathrm{N}\left(\mathrm{B}{\displaystyle \cap }\mathrm{A}\right)}+\frac{1}{\mathrm{B}{\displaystyle \cap }\overline{\mathrm{A}}}$$

(7)

$${\mathrm{S}}^2{\mathrm{W}}^{-}=\frac{1}{\mathrm{N}\left(\overline{\mathrm{B}}{\displaystyle \cap }\mathrm{A}\right)}+\frac{1}{\overline{\mathrm{B}}{\displaystyle \cap }\overline{\mathrm{A}}}$$

(8)

It is noted that N represents the number of unit cells. W_final is also determined by Equation (9):

$${\mathrm{W}}_{\mathrm{final}}=\mathrm{C}/\mathrm{S}\left(\mathrm{C}\right)$$

(9)

3.5.3. Evaluation Based on Distance from Average Solution

Evaluation Based on Distance from Average Solution (EDAS) is a novel MCDM technique that was developed by Keshavarz Ghorabaee et al. (2017) [26]. To evaluate different alternatives, the EDAS model employs positive and negative distances from an average solution. This technique includes six steps which are described below:

• Step 1. Determining the decision matrix:

$$\mathrm{X}=\left[\begin{array}{ccccc}{\mathrm{x}}_{11}& {\mathrm{x}}_{12}& {\mathrm{x}}_{13}& \dots & {\mathrm{x}}_{1\mathrm{m}}\\ {\mathrm{x}}_{21}& {\mathrm{x}}_{22}& {\mathrm{x}}_{23}& \dots & {\mathrm{x}}_{2\mathrm{m}}\\ \dots & \dots & \dots & \dots & \dots \\ {\mathrm{x}}_{\mathrm{n}1}& {\mathrm{x}}_{\mathrm{n}2}& {\mathrm{x}}_{\mathrm{n}3}& \dots & {\mathrm{x}}_{\mathrm{nm}}\end{array}\right]$$

(10)

where n shows alternatives and m is indicative of criteria. x_ij also shows the status of the i alternative in the j criterion.

• Step 2. Calculation of the average solution of criteria:

$$AV={\left[A{V}_j\right]}_{1\times m}$$

(11)

$A{V}_j$ is calculated using Equation (12).

$$A{V}_j=\frac{{{\displaystyle \sum }}_{i=1}^n{X}_{ij}}{n}$$

(12)

• Step 3. Computation of the positive and negative distance from the average value:

If the criteria have a positive impact on flooding, the distances from the average value are calculated using Equations (13) and (14):

$$PD{A}_{ij}=\frac{\max \left(0.\left(X_{ij}-A{V}_j\right)\right)}{A{V}_j}$$

(13)

$$ND{A}_{ij}=\frac{\max \left(0.\left(A{V}_j-X_{ij}\right)\right)}{A{V}_j}$$

(14)

where PDA is the positive distance from the average value and NDA is the negative distance from the average value.

On the other hand, in the case of the negative impact of criteria on flooding, positive and negative distances are computed using Equations (15) and (16).

$$PD{A}_{ij}=\frac{\max \left(0.\left(A{V}_j-X_{ij}\right)\right)}{A{V}_j}$$

(15)

$$ND{A}_{ij}=\frac{\max \left(0.\left(X_{ij}-A{V}_j\right)\right)}{A{V}_j}$$

(16)

• Step 4. Determination of the optimistic ($S{P}_i$) and pessimistic ($S{N}_i$) weighted sum of the positive and negative distances for alternatives:

For this step, Equations (17) and (18) are calculated as follows:

$$S{P}_i={\displaystyle \sum }_{j=1}^m{w}_{j}PD{A}_{ij}$$

(17)

$$S{N}_i={\displaystyle \sum }_{j=1}^m{w}_{j}ND{A}_{ij}$$

(18)

w_j is the weight of the j_th criterion, in which ${\displaystyle \sum }_{j=1}^m{w}_j=1$. It should be highlighted that the assigned weights are determined using the AHP technique [49]. The implementation of this method is explained in the literature [17,50,51].

• Step 5. Normalization of the SP and SN values:

Using Equations (19) and (20), the normalized values of the optimistic and pessimistic weighted sum of the positive and negative distances for all alternatives are computed, respectively.

$$NS{P}_i=\frac{S{P}_i}{{\text{max}}_i\left(S{P}_i\right)}$$

(19)

$$NS{N}_i=\frac{S{N}_i}{{\text{max}}_i\left(S{N}_i\right)}$$

(20)

• Step 6. Alternatives ranking:

The final score of each alternative ($A{S}_i$) is determined using Equation (21) and then they are ranked.

$$A{S}_i=\frac{1}{2}\left(NS{P}_i+NS{N}_i\right)$$

(21)

where $0\le A{S}_i\le 1$.

It should be noted that an alternative with the greatest $A{S}_i$ value is the best choice among all alternatives.

3.5.4. Multilayer Perceptron Layer

Multilayer perceptron layer (MLP) as an artificial neural networks model is a precise technique in natural hazard prediction [52,53]. It includes three stages as the following:

• Step 1. Building the architecture of the network:

In this study, the input layer contains ten neurons, corresponding to flood influential criteria. The output layer is composed of two neurons, indicating flood and non-flood pixels.

• Step 2. Training the network:

In this stage, in order to raise the accuracy of outputs, a backpropagation algorithm was employed [20]. In this procedure, weights of layers were computed via reverse calculation processes, and therefore, neuron numbers in the hidden layer were assigned by trial and error. This process continues until reaching the minimum RMSE [34].

In the current research, 70% of input data were chosen randomly for training and 30% for testing purposes to verify the network performance.

• Step 3. Testing:

The remaining 30% of data were selected to determine the model’s efficiency in this stage (Equation (22)). Besides, the two factors of sensitivity and specificity were used to evaluate the model’s predictive ability (Equations (23) and (24)).

$$Accuracy=\frac{TP+TN}{TP+FP+TN+FN}$$

(22)

$$specifity=1-\left[\frac{TN}{TN+FP}\right]$$

(23)

$$sensitivity=\left[\frac{TP}{TP+FN}\right]$$

(24)

where TP indicates correctly grouped flood pixels, TN shows correctly grouped non-flood pixels, and, FP and FN express numbers of incorrectly grouped pixels.

3.6. Validation of Models

Validation of the applied models was evaluated on both training and testing datasets using the ROC curve. Validating by training dataset indicates the goodness-of-fit of the models, while validating through testing datasets indicates the predictive capability of the models. ROC curve determines the model’s performance in diagnosing flood-prone areas and has been employed in several studies [24,54,55,56]. The area under the ROC curve (AUC) implies the precision of methods, varying between 0.5 and 1 [35,55].

4. Results

4.1. Importance of Influential Factors on Flood Occurrence

The values of IG for flood influential factors are presented in

Table 2. Importance of flood influential factors for flood modeling.

Flood Influential Factor	IG	Flood Influential Factor	IG
Altitude	0.69	TWI	0.42
Slope	0.78	Rainfall	0.35
Aspect	0.12	Soil type	0.23
Plan curvature	0.47	Lithology	0.51
Distance from rivers	0.73	Land use	0.57

. It can be observed that the slope factor had the highest importance on flooding due to the largest IG value (IG = 0.78), followed by distance from rivers (IG = 0.73), altitude (IG = 0.69), land use (IG = 0.57), lithology (IG = 0.51), plan curvature (IG = 0.47), TWI (IG = 0.42), rainfall (IG = 0.35), soil type (IG = 0.23), and slope aspect (IG = 0.12), respectively. As indicated, all the values are higher than 0, which implies the effectiveness of all factors on flood occurrences. Even though rainfall is an initial and key factor in flooding, it has one of the lowest IG values, which can be explained by the fact that morphological and topographical conditions may have more influence on flooding than rainfall intensity. The influential factors derived from DEM have a high impact on the flow accumulation process because the surface runoff moves to regions with low slope degrees and low altitudes where flow accumulates. Moreover, the distance from rivers factor had a high IG value because the rivers are the sources of flooding.

4.2. Coefficients Calculation

Results of AHP, FR, and WOE coefficient computations are displayed in

Table 3. Relationship between influential factors and flood points.

Criteria	Class	Number of Pixel	Number of Floods	C/Sc	FR	AHP	AHP × FR	AHP × C/Sc
Altitude	2027–3820	144,558	1	−2.147	0.121	0.32	0.039	−0.687
	1316–2027	233,726	1	−2.658	0.075		0.024	−0.851
	719–1316	365,764	12	−2.039	0.576		0.184	−0.652
	234–719	589,554	25	−1.651	0.744		0.238	−0.528
	−40–234	1,613,803	129	5.503	1.402		0.449	1.761
Slope	19.2–75.5	80,633	4	−0.282	0.870	0.09	0.078	−0.025
	12.4–19.2	229,861	10	−0.889	0.763		0.069	−0.080
	7.1–12.4	400,470	22	−0.186	0.964		0.087	−0.017
	2.3–7.1	589,786	32	−0.312	0.952		0.086	−0.028
	0–2.3	1,646,655	101	1.108	1.076		0.097	0.100
Aspect	Flat	12,050	6	5.290	8.736	0.02	0.175	0.106
	North	437,462	23	−0.420	0.922		0.018	−0.008
	Northeast	338,379	18	−0.311	0.933		0.019	−0.006
	East	316,204	20	0.492	1.110		0.022	0.010
	Souteast	366,613	16	−1.140	0.766		0.015	−0.023
	South	434,844	20	−1.038	0.807		0.016	−0.021
	Southwest	337,394	22	0.670	1.144		0.023	0.013
	West	323,159	21	0.637	1.140		0.023	0.013
	Northwest	381,300	22	0.061	1.012		0.020	0.001
Plan Curvature	Concave	1,393,068	66	−2.062	0.831	0.1	0.083	−0.206
	Convex	1,360,360	76	−0.238	0.980		0.098	−0.024
	Flat	193,977	26	4.477	2.352		0.235	0.448
Distance from river	>3000	53,994	2	−0.615	0.650	0.17	0.110	−0.105
	2000–3000	342,862	10	−2.246	0.512		0.087	−0.382
	1000–2000	1,045,667	22	−5.660	0.369		0.063	−0.962
	500–1000	1,131,426	71	1.032	1.101		0.187	0.175
	<500	373,457	63	8.907	2.960		0.503	1.514
TWI	1.4–8.1	133,167	3	−1.642	0.395	0.05	0.020	−0.082
	8.1–9.8	544,028	17	−2.730	0.548		0.027	−0.137
	9.8–11.1	1,016,028	49	−1.443	0.846		0.042	−0.072
	11.1–12.6	875,180	55	0.863	1.103		0.055	0.043
	12.6–19.02	379,002	44	5.000	2.037		0.102	0.250
Rainfall	54–258	309,529	14	−0.914	0.794	0.07	0.056	−0.064
	258–417	689,864	46	1.215	1.170		0.082	0.085
	417–595	446,783	12	−2.813	0.471		0.033	−0.197
	595–751	864,585	24	−4.139	0.487		0.034	−0.290
	751–1000	636,644	72	6.423	1.984		0.139	0.450
Soil type	Entisols-Rock Outcrops/Aridisols/Inceptisols/Playa	1,345,296	76	−0.105	0.991	0.05	0.050	−0.005
	Mollisols	1,225,408	88	2.819	1.260		0.063	0.141
	Alfisols	327,844	2	−3.290	0.107		0.005	−0.165
	Salt flats	46,135	1	−0.974	0.380		0.019	−0.049
	Silty loam	2723	1	1.862	6.443		0.322	0.093
Litologhy	Paleozoic	226,1742	2	−7.890	0.016	0.04	0.001	−0.316
	Mesozoic	331,254	10	−2.127	0.530		0.021	−0.085
	Proterozoic	1166	1	2.708	15.052		0.602	0.108
	Cenozoic	353,244	155	15.488	7.698		0.308	0.620
Land use	Dense forest—Mountainous areas	702,647	7	−5.113	0.175	0.07	0.012	−0.358
	Forest lands—Agriculture	28,730	3	1.053	1.832		0.128	0.074
	Fruit trees—Agricultural lands	1,610,365	118	3.986	1.286		0.090	0.279
	Herbaceous plants—Groves	487,559	21	−1.404	0.756		0.053	−0.098
	Urban–Coastal areas	118,104	19	4.584	2.822		0.198	0.321

for each class. As indicated, according to the AHP technique, altitude, distance from rivers, plan curvature, slope, rainfall, and land use had the highest influence on flood occurrence, in respective order. On the other hand, soil type, TWI, lithology, and slope aspect had the lowest impact on flooding. The CR value for pairwise comparison was equal to 0.019 (<0.1), which expresses the reasonable compatibility of weighting.

Higher FR values imply a greater relationship between flood criteria and occurrences [14,43]. The results showed that the last class of altitude factors had the greatest effect on flooding due to the highest FR value. On the other hand, the other four classes had lower FR values. It can be concluded that the flood probability decreases in the case of altitude increases. Likewise, the FR weight for the slope factor was the highest in the last class and it had a lower value in other classes. It is due to the fact that the lower slopes have more possibilities of flooding and vice versa. Regarding the slope aspect factor, the greatest FR value was obtained by the flat category. However, the lowest value was determined by the southeast class. This finding is in accordance with the findings of Termeh et al. (2018) [19] and Rahmati et al. (2015) [44], as it is found that the flood potential in the flat aspect is higher than in other categories. Results of plan curvature indicated that the flat shape had the greatest FR value, followed by convex and concave. Flat shapes had more possibility of flooding, compared to other categories, due to keeping surface runoff for more time. For the distance from rivers factor, the last two classes had the greatest influence, owing to higher FR weights. In the case of the TWI factor, the last class had the highest weight, while the first class achieved the lowest value. As is shown for the rainfall criterion, the second and last categories positively influence flood occurrence. The growth in rainfall intensity had no impact on flooding, owing to the fact that by increasing altitude, precipitation would increase as well. The most susceptible category of soil factor was silty loam, followed by Mollisols. Also, other categories negatively affect flooding. In the case of the lithology factor, Proterozoic and Cenozoic classes had a great effect on flooding, followed by Mesozoic and Paleozoic, respectively. For the land-use factor, the highest weight was observed in urban-coastal areas. In general, these residential zones are more prone to flood risk, due to their proximity to rivers and having large populations [57]. Furthermore, the dense forests-mountainous areas had the least FR values and less influence on flood events.

According to the WOE model, the last category of altitude had the greatest weight and positively influenced flooding. At the same time, the other categories had negative weights. As indicated for the slope factor, by increasing the degree, final weights decrease. Regarding the slope aspect, the highest weight was obtained by the flat category. In contrast, the southeast category had the least WOE weight. For distance from rivers, classes of <500 m and 500–1000 m had the highest effect. Although the last two classes of the TWI factor had positive weights, the other three classes had negative impacts on flood occurrences. Similar to the FR method, the second and last categories of rainfall factors had positive effects on flooding, but other categories had adverse impacts. The WOE analysis for the soil type displayed the more susceptible category related to Mollisols, followed by silty loam. In contrast, other categories negatively affect flood events. The result of lithology revealed that the Cenozoic class had the greatest effect on floods, but the Paleozoic category had a negligible impact. In the case of land use factor, the highest and lowest WOE values were obtained for urban-coastal and dense forests-mountainous areas, in respective order.

4.3. EDAS-FR-AHP and EDAS-WOE-AHP Methods

In order to produce flood zonation maps of EDAS-FR-AHP and EDAS-WOE-AHP methods, values of FR-AHP and WOE-AHP were multiplied by pixel values of the EDAS technique. Ultimately, flood zonation maps were categorized into five classes using the Quantile algorithm in GIS [13,23,58]. The following maps indicate the flood-prone areas by the EDAS-FR-AHP and EDAS-WOE-AHP models. As is shown, very low and low categories covered 29.82% and 37.62% of the region. The moderate category was also spread over 17.3% and 20.50% of the Province, and finally, very high and high categories covered 52.88% and 41.88% of the Golestan Province as shown in Figure 4 and Figure 5.

4.4. MLP-FR Method

After several trials and errors, the minimum RMSE value stood at 0.015 and remained stable. Resultantly, the final MLP-based structure was obtained as the following (Figure 6):

In order to evaluate this model’s performance, statistical measures were calculated for both training and testing datasets, which are displayed in

Table 4. Values of statistical measures.

	Training	Testing
Sensitivity	0.864	0.851
Specificity	0.912	0.893
Accuracy	0.892	0.876

. Regarding the results, the MLP-FR method had high precision in both training and testing samples. Hence, this method is suitable for predicting the flood-prone zones in this study area.

The result of the MLP-FR model is presented in Figure 7, which is also grouped into five classes. Concerning this map, very low and low categories accounted for 27.17%, while the moderate category covered around 19.02% of the area. Additionally, the largest part of the study location (53.81%) was covered by high and very high classes of flood potential.

4.5. Validation

In order to represent the validation results, both success and prediction rate, with the help of training and validating datasets were investigated (Figure 8). Regarding the success rate, MLP-FR was the most accurate model among the three ensemble methods (AUC = 0.934), followed by EDAS-FR-AHP (AUC = 0.901) and EDAS-WOE-AHP (AUC = 0.883). Besides, MLP-FR had the highest precision (AUC = 0.912) based on the prediction rate, followed by EDAS-FR-AHP (AUC = 0.875) and EDAS-WOE-AHP (AUC = 0.844), respectively.

5. Discussions

Producing accurate flood zonation maps is considered a crucial way for preventive measures and flood mitigation across the world. The latest destructive flooding in Golestan Province led authors to apply several new models to identify flooding regions. Hence, three ensemble models, namely EDAS-FR-AHP, EDAS-WOE-AHP, and MLP-FR were employed to generate flood zonation maps. Flood occurrences are affected by different effective factors, demanding to recognize them. According to the feature selection method, the slope factor was the most effective factor, followed by distance from rivers, altitude, and land use. Among different factors, the least effective one was the slope aspect. This result is in agreement with some previous studies [9,19,23,43]. The validation findings showed that the proposed MCDM-based methods (EDAS) predict the flood potential of the study location with high precision. The main advantage of the techniques is to present optimistic and pessimistic solutions that provide flexibility in the final evaluation for decision-makers. Also, the EDAS approach is performant in solving stochastic problems due to the computation of the average solution by the arithmetic mean [26].

The upper crescent part of Golestan province had low altitude and low slope degrees. As a result, areas with low altitudes and low slopes receive more runoff during and after rainfall. These are the main reasons behind the large proportion of high and very high susceptible areas to flood. According to the maps generated by three models, the areas with very low to moderate flood susceptibility are mainly located in the southern and southwestern parts of the region, where the Alborz Mountain range acts as a barrier and prevents the entry of humidity derived from the Caspian Sea into these areas. Thus, these regions have low rainfall and a dry climate. As a result, the possibility of flood occurrence in this part of Golestan Province is low. Besides, the areas with high flood susceptibility are located in the northern and northwestern parts. Evaporation of the Caspian Sea increases the humidity in these regions giving rise to heavy rainfall that can lead to floods. The proximity of the water table to the ground surface, as well as the saturated soil, increases the flood intensity.

A comparison of bivariate statistical methods such as FR and WOE was performed by several researchers. In a study in Golestan Province, the flood susceptibility was assessed using FR (AUC = 76.47%) and WOE (AUC = 74.74%), in which both models had reasonable accuracy and almost similar results [44]. It should be highlighted that the ensemble models of this study had higher accuracy, which proved their capability in the study area. Likewise, for flood zonation determination across Mazandaran Province in Iran, the FR model showed higher prediction performance than the WOE method, which is in line with the reported studies [14]. In another study in Golestan Province, the use of a novel MCDM model combined with bivariate statistics was investigated and both MAIRCA-FR and MAIRCA-WOE indicated reasonable precision [10]. Research in Bangladesh found that the precision of ensemble methods was more than standalone methods [59]. In another research study in the China basin, it was revealed that multilayer-based methods had better efficiency in comparison to the MCDM methods, such as TOPSIS and VIKOR [23], which concurred with the findings of the current study. This is because of the independence of flood factor weights based on the experts’ judgment [37,60]. The evaluation of the validity of standalone and ensemble of the MLP method in a study in Golestan Province indicated that the accuracy of the standalone method was not nearly as much as that of the ensemble model. It was also shown that the northern parts of Golestan Province were prone to severe floods rather than other areas [61].

Considering that in the current study, FR-AHP and WOE-AHP were combined with a new MCDM model (EDAS). The FR method was also applied as an ensemble with the MLP technique, and the precision of the flood susceptibility was improved compared to the previously used methods.

6. Conclusions

This research focused on the determination of flood risks across the Golestan Province in Iran, using three ensemble models, including EDAS-FR-AHP, EDAS-WOE-AHP, and MLP-FR. A flood inventory map was generated by 240 flooding and non-flooding locations, and then they were divided randomly into two categories of training (70% of data) and testing (30% of data). Ten flood effective criteria were chosen for the modeling process. Results showed that 52.88, 41.88, and 53.81% of the study region was covered by high and very high flood hazard risk for EDAS-FR-AHP, EDAS-WOE-AHP, and the MLP-FR, respectively. Accordingly, the MLP-FR model was the most efficient, owing to the highest AUROC value in training (0.934) and testing (0.912) datasets, followed by EDAS-FR-AHP and EDAS-WOE-AHP, with the success precision of 0.901 and 0.883 and predictive precision of 0.875 and 0.844, respectively. Newly employed combined methods (EDAS-based techniques) show reasonable precision in flood zonation prediction, as well as the MLP ensemble model. Resultantly, the MLP-based model was more precise than the MCDM-based model in this research owing to over-fitting prevention and noise problems, and also no dependence on the experts’ judgment for weights assignments, causing uncertainty in MCDM models, which is the most important limitation of this work. It is recommended to use EDAS-based models in other study areas in further research to ensure its validity in predicting flood potential areas.