A comparative assessment of multi-criteria decision analysis for flood susceptibility modelling

Abstract The purpose of this study is to develop a reliable model for identification of flood susceptible areas. Three multi-criteria decision-making techniques, namely Analytical Hierarchy Process (AHP), Technique for Order of Preference by Similarity to Ideal Solution (TOPSIS), and Attributive Border Approximation Area Comparison (MABAC) methods combined with weight of evidence (WOE) were used in Babolroud Watershed, Iran. First, 50 flood locations were identified, of which 35 (70%) locations were selected randomly for modelling, and 15 (30%) locations were used for validation. Using GIS with eight conditioning factors including rainfall, distance from rivers, slope, soil, geology, elevation, drainage density and land use, the flood susceptibility maps were prepared. Area under receiver operating characteristic curve (AUROC) for test data of AHP-WOE, TOPSIS-WOE-AHP and MABAC-WOE-AHP methods were 72.8%, 90.5% and 81.6%, respectively, which indicate the reasonable accuracy of models. High accuracy of the proposed new model (MABAC) clarifies its applicability for preventive measures.


Introduction
Among various natural disasters, floods are one of the most devastating cases, causing many damages and human losses in the world (Youssef et al. 2016;Sandink et al. 2016). Floods can be defined as intense rainfall and snowmelt that overflow into rivers and flood plains, consequently covering the surrounding areas (Kron 2002). A flood can damage environmental ecosystems, transportation networks, structures, economy and human life (Shirzadi et al. 2012;Yu et al. 2013). Heavy and prolonged rainfall, landslides, climate and land use changes, and weak dam construction are some of the main reasons for flood occurrences (Barredo 2009;Ahmadlou et al. 2018).
Severe flood occurrences and subsequently huge environmental losses have an increasing trend in Iran, especially in the northern cities and southern margin of Caspian Sea (Boudaghpour et al. 2015). Among these areas, the Babolroud Watershed in Mazandaran Province is always vulnerable to flood hazards due to extensive rainfall, and geographical conditions (proximity to Caspian Sea, presence of low-landing plains, and presence of permanent rivers). The increasing urbanization and conversion of forest and agricultural lands into residential areas results in a hydrological process imbalance (Sadeghi-Pouya et al. 2017;Khosravi et al. 2019). Also, due to climate change, heavy rainfalls occur in a short time and cause flood occurrence, if the intensity is more than capacity of soil permeability. Severe flood occurrence in March 2019 in Mazandaran Province leads to huge damages according to Iranian Red Crescent (https://reliefweb.int/report/iran-islamicrepublic/information-bulletin-1-golestan-mazandaran-flood-iran-timeframe-16-26).
Iran Country government estimated US $80 million damage caused by this flood (https://www. irna.ir/news/83277780/). Difficulties in data collection, lack of information about floodprone areas and not taking the necessary preventive actions lead to these damages.
One of the leading solutions for flood hazard preventive measures is preparing flood susceptibility mapping (Dang et al. 2011;Bubeck et al. 2012). Effective management by applying appropriate approaches causes to minimize damages to people and environment. The risk level of flood in different parts of the area is determined using a flood susceptibility map. By matching these maps with land use maps, it is possible to identify endangered areas, including cities, villages, bridges, factories, etc., to take the necessary measures to protect them.
Generally, study of natural disasters requires multi-temporal datasets to diagnose the spatial changes and the processes of hazardous events (Martinez and Le Toan 2007). GIS and remote sensing (RS) techniques provide a vast application in the natural disaster risks analysis, due to their efficiency and cost effectiveness (Haq et al. 2012;Jaafari et al. 2014). Therefore, in several studies (Dewan 2013;Pradhan et al. 2011;Kazakis et al. 2015), these methods have been suggested to evaluate flood hazard zoning.
Multi-criteria decision analysis (MCDA) has been recognized as an essential tool for analyzing complex decision problems that often involve incomparable data or criteria (Hwang and Yoon 1981;Malczewski 2006). AHP (Saaty 1977) is a popular method of MCDA , based on experts' knowledge in assigning weights. Fernandez and Lutz (2010) have assessed the efficiency of MCDA and GIS for mapping the flood-susceptible areas in Tucuman Province, Argentina. Their study showed that the AHP method in the GIS environment is a powerful technique to generate natural disaster hazard maps with reasonable accuracy. Zou et al. (2013) recognized the AHP technique as an understandable, convenient, and cost-effective method for flood hazard evaluation. Subramanian and Ramanathan (2012) indicated that the AHP technique is suitable for regional researches. Other techniques of MCDA for evaluating natural disaster risks are Technique for Order of Preference by Similarity to Ideal Solution (TOPSIS) and Multi-Attributive Border Approximation Area Comparison (MABAC). TOPSIS has been increasingly applied in flood hazard management because of its ability to deal with multiple attributes (Brito and Evers 2016). Mojtahedi and Oo (2016) combined a non-parameter resampling bootstrap technique and TOPSIS to rank the flood risk of Australia's states. Zhu et al. (2018) developed a TOPSIS-based model for reservoir flood control operation. In another study, Luu et al. (2019) applied the new multiple linear regression-TOPSIS for Vietnam's flood risk assessment. The fuzzy-TOPSIS technique has been used to evaluate flood risk in several studies (Jung et al. 2011;Lee et al. 2014). On the other hand, MABAC is one of the latest MCDA methods, developed at the University of Defence in Belgrade research centre (Pamu car and Cirovi c 2015). A few studies were carried out in different fields using the MABAC method (Gigovic et al. 2017;Liang et al. 2019). To the best of our knowledge, this method has yet to be used to assess flood hazards and is new in the modelling of flood susceptibility. The main goal of present research is preparing flood susceptibility maps in a flood affected watershed of Babolroud, for the first time by MCDM methods. In general, MCDM methods are one of the most popular approaches in flood hazard determination, in which applied by many researchers in recent studies (Rahmati et al. 2016;Rahman et al. 2019;Ogato et al. 2020). Improving the accuracy of these models has always been one of the main concerns of researchers and decision makers. While, all papers have shown no universal guideline in finding the best model in a region due to the models limitation, some researches (Shafapour Tehrany et al. 2017;Hong et al. 2018;Costache et al. 2020;Rahman et al. 2019) proved that combining several methods may lead to higher accuracy in flood susceptibility assessment. On the other hand, applying just one method is not adequate to predict the flood potential in a specific area. Therefore, this study aims to apply several hybrid models to increase the accuracy of flood-prone prediction. In previous studies, MABAC has not yet been used for the flood susceptibility assessment, which means that this methodology lacks in literature. Thus, the purpose of present study is to evaluate the flood susceptibility through a novel hybrid model which was generated by the combination of MABAC, WOE and AHP. This method is a novel technique in flood assessment, and this gap in research was filled by determining the predictive power and comparing it to other hybrid methods (TOPSIS-WOE-AHP and AHP-WOE). The main reason for using these methods is to develop a high-performance model for estimating the flood potential in Babolroud Watershed. Therefore, the study's contributions are applying three kind of models to determine flood susceptible areas; and generate a new and unique flood hazard map for Babolroud Watershed, using hybrid models. For this purpose, the ROC Curve was used to validate and determine the accuracy of each model. Modelling and data processing were done in ArcGIS 10.5 environment.

Study area
The study area is located in the Babolroud Watershed ( Figure 1), with an area of approximately 2344 km 2 that occupied about 9.88% of Mazandaran Province in the north of Iran. The study area lies between the latitude of 35 56 0 to 36 44 0 N and the longitudes of 52 11 0 to 53 3 0 E. The area's minimum elevation is À40 m, and the maximum is 3991 m above sea level. High differences in elevation indicate the morphological imbalance of the watershed. The average annual rainfall is about 627 mm. Land use patterns include rangeland (mostly in the height of more than 2000 m above sea), forests, garden lands, bare lands and residential areas. About 46.92% of the watershed is mountainous and the rest is coastal plains. Babolroud Watershed has two different types of climate, due to its plain and mountainous areas. North part of the watershed has temperate Caspian climate because of presence of low lands. Central and southern parts of watershed have mild mountainous climate due to their high elevation. Some features of this type of climate are atmospheric precipitation and a little difference in the annual temperature. Some details about flood inventory map, flood conditioning factors, different flood susceptibility models, models performance evaluation and finally feature selection method are explained in this section.

Flood inventory map
The future flood event in a region can be predicted by analyzing the records of its occurrence (Manandhar 2010). A flood inventory map is considered the most crucial factor in studying the relationship between flood events and different effective factors. In this study, a flood inventory map was prepared using historical data of 50 flood occurrence between 2012 and 2017. These locations were validated in the field by Mazandaran Water Resources Department. Flood locations were divided randomly into two datasets of training (almost 70%) and predicting (nearly 30%) data (Ohlmacher and Davis 2003;Pradhan 2010;Pourtaghi and Pourghasemi 2014;Khosravi et al. 2019). For this aim, 50 flood locations were identified in the study area, of which 35 (70%) locations were selected  randomly for modelling, and the remaining 15 (30%) were used for validation purposes ( Figure 1). Some photographs of the recent flood in Babolroud Watershed that causes a lot of damages are shown in Figure 3.

Flood conditioning factor
For preparing a flood susceptibility map, determining the flood conditioning factors is crucial (Kia et al. 2012). The flood conditioning factors were determined by literature review, personal observation, experiences and experts' knowledge (Kourgialas and Karatzas 2011;Ouma and Tateishi 2014;Ogato et al. 2020;Suthirat et al. 2020). Accordingly, rainfall, distance from streams, slope, soil, geology, elevation, drainage density and land use were considered as important flood hazard factors in the Babolroud Watershed. All the mentioned factors were converted to a raster grid with 30 m Â 30m grid-cells for application models.
Among these factors, the rainfall factor is the leading cause of flood in an area. Because of complex spatiotemporal dynamics of rainfall, it is difficult to describe in modelling (Pradhan 2010). Thus, a long term of average annual rainfall (e.g. 10 years) was applied in flood hazard modelling. For this aim, rainfall data for period of 2008-2018 was collected at eight stations. For interpolation of rainfall, different methods such as simple kriging, ordinary kriging, Inverse Distance Weighted (IDW) with the power of 1-5 and spline were used. According to Table 1, the lowest root mean square error (RMSE) and mean absolute error (MAE) were calculated and the map of rainfall has been produced with IDW in ArcGIS software (Jahangir et al. 2019;Ogato et al. 2020). Then to create a continuous raster rainfall data, it was converted to a raster layer. To produce a spatial map of rainfall intensity, Fournier's index had been used by meteorological information of eight hydrometric stations of Babolroud Watershed from 2008 to 2018 as shown in Equation (1) (Morgan 2005).
MFI is the modified Fournier index, P i is the average monthly rainfall and P is the average annual rainfall. This index shows the sum of the average monthly rainfall intensity at a station (Kourgialas and Karatzas 2011;Hernando and Romana 2015). In the present study, the IDW with the power of 2 interpolation method which is the most appropriate method was used to construct a raster spatial rainfall map ( Figure 4a).
Streams to pass the rainwater are vital factors that identify their density and distance from structures effective in predicting and reducing flood damage. Due to most villages in the study watershed in the river area, the amount of financial and human losses during floods reaches its maximum. The distance map in the ArcGIS environment was prepared using Euclidean Distance function from Spatial Analyst tools ( Figure 4b).
The slope factor has a significant role in the passing and stillness of surface waters. This factor is inversely related to the flood phenomenon. In ArcGIS software, the slope raster map was prepared using the digital elevation model (DEM) with a spatial resolution of 30 m, which was downloaded from the United States Geological Survey (USGS) site and slope generation tools ( Figure 4c).
Soil attributes such as permeability, thickness and infiltration rate directly affects the rainfall-runoff process (Rimba et al. 2017). The flood hazard increases with a decrease in soil infiltration capacity, which causes an increase in surface runoff. When water is supplied at a rate more than the soil's infiltration capacity, it moves as runoff on sloping land and can cause flooding. There are some types of soil in the study area, namely Chromic Vertisols, Inceptisols, Coastal Sands, Salt Flats, Water Body, Alfisols, Mollisols were, Rock Outcrops/Entisols and Rock Outcrops/Inceptisols ( Figure 4d).
Geological formations have an essential role in conducting or absorbing surface water due to their permeability. geology layer was classified into five groups: areas with impermeable rock layers, areas with clay soilsbasaltgranite, areas with fine sandmixed soils, areas with conglomerate-travertine-coarse-grained texture and areas with gypsumlime structure ( Figure 4e).
The land use factor indicates the type of land can be a factor in accelerating or decelerating flood intensity. Land use change is one of the significant contributors to flooding, as urban expansion increases, impervious cover increases and forest cover decreases contributing to increased run-off (Hall et al. 2013). The land-use classes of the area were classified into five groups: urban-coastal areas, herbaceous plantsgroves, fruit treesagricultural lands, forest landsagriculture and finally, dense forestmountainous areas ( Figure 4f).
The elevation controls the movement of the overflow direction. The elevation raster map was produced using DEM with a spatial resolution of 30 m and slope generation tools in the ArcGIS environment ( Figure 4g).
Drainage density is defined as the ratio of the length of drainage per area. Drainage density is an inverse function of infiltration. In other words, greater drainage density shows high runoff for the area and less prone to flood occurrences. The drainage map is overlaid on the basin map to determine the total length of streams to the basin's total area. A line density module was used to compute drainage density of the basin using the Spatial Analyst extension in ArcGIS.

AHP model
The Analytic hierarchy process established by Saaty (1977) is a technique for solving multiple criteria decision problems by setting their priorities and enables the user to achieve scale advantages drawn from a set of alternatives (Saaty 1977). This method is implemented using the following steps: At first, the pairwise comparison judgments were determined, using the nine-point scale (Saaty 1977). For constructing the matrix, each flood hazard factor was rated considering with the other factors, which were given values from 1/9 to 9 ( Table 2). After that, all values from the pairwise comparison were normalized by dividing each entry in a column by adding all entries in that column to sum it to 1. Finally, to avoid the occurrence of any incidental judgement in the pairwise comparison matrix, the consistency ratio (CR) was generated using the consistency index (CI) value (Equations (2) and (3)).
Where k max is the largest eigenvalue of the pairwise comparison matrix and n is the number of factors.
where CI is the consistency index and RI is a random consistency index. The random consistency index of 1.45 is for the nine flood hazard factors (Saaty 1977). When the   (Saaty 1977).
Description Definition Degree of importance The two elements are equally important Same importance 1 One element is relatively preferable to the other Relative preference 3 One element is preferred over another High preference 5 One element is much preferred over another Too much preference 7 One element has a great advantage over the other.
Extremely high preference 9 Intermediate values in judgments 2, 4, 6, 8 calculated CR was equal to (or less than) 10%, the estimated weighting coefficients were acceptable.

WOE model
Weights of evidence method are based on log-linear Bayesian theory, which was first used in the field of mineral exploration in 1988 (Bonham-Carter et al. 1988). This bivariate statistical model is based on probability theory and can estimate the relationship between a flood occurrence and different effective factors. It has been applied for producing natural hazards susceptibility maps (Jaafari et al. 2014;Youssef et al. 2015;Khosravi et al. 2018). The model determines the weights of classes for each effective factor as follows (Bonham-Carter et al. 1988): where P is a probability, B is the presence, B is the absence of each effective factor that has an impact on flooding. Also, A is the presence and A is the absence of flood, respectively (Xu et al. 2012). Weight contrast is the difference between W þ and W À (Dahal et al. 2008). The magnitude of this contrast indicates the spatial correlation between the effective factors in flood occurrence. The value of C is positive and negative, for positive and negative relations, respectively (Pourghasemi et al. 2013). The standard deviation for W is determined as Equation (6): S 2 W þ and S 2 W À are the variance of positive and negative weights, respectively and calculate as follows: N denotes the number of unit cells. Finally, the final weight (W final ) is calculated as follows: Where C is the contrast and S ðCÞ is standard deviation.

MABAC model
MABAC model was first applied by Pamu car and Cirovi c in 2015 and developed at the University of Defence in Belgrade for weight criteria and evaluation options. The casual relationship between the criteria is estimated using the modified fuzzy DEMATEL method. A DEMATEL method is a pragmatic tool for visualizing the structure of complex casual relationships (Buyukozkan and Çifçi 2012). Based on this method, criteria are divided into two categories of cause and effect. The extent of this effect is used to estimate the criteria' weight (Dalalah et al. 2011). MABAC model consists of six important steps, as described below.
Step 1. Calculation of the initial decision matrix (X) In each problem, there are m alternatives and n criteria. Each of the options is displayed in the form of vectorization as A i ¼ ðx i1 :x i2 :x i3 : . . . :x in Þ where x ij specifies the status of the i option in the j criterion. Based on this, the initial decision matrix (X) is shown below.
X ¼ x 11 x 12 x 13 . . . (10) Step 2. Normalize the initial decision matrix elements (N) In this step, the decision matrix is normalized to neutralize the effect of the different scales of the criteria. Equations (11) and (12) are used to normalize the positive and negative criteria, respectively.
Where x ij , x þ i and x À i are the elements of the initial decision matrix (X). x þ i and x À i are the maximum and minimum value of the observed criterion according to the alternatives, respectively, and defines as: Step 3. Calculation of a balanced normal decision matrix (V) In this step, the elements of balanced normal decision matrix are determined by Equation (15).
Where v ij are the elements of the balanced matrix (V), n ij are the elements of the normalized matrix (N) and w i are the weight coefficients of the criteria, which are determined by the AHP method in this study.
Step 4. Determining the border approximation area matrix (G) The border approximation area for each criterion is calculated by Equation (16).
Where v ij are the elements of the weighted matrix (V) and m is the total number of alternatives. After determining the value g i for each criterion, a border approximation area matrix (G) is estimated according to Equation (17).
Where n is the total number of criteria according to a selection of alternatives.
Step 5. Calculation of the distance of the alternatives from the border approximation area This distance is equal to the difference between the weighted matrix elements (V) and the value of the border approximation area matrix (G), as indicated in Equation (18).
After calculating the value of matrix Q, the belonging of each alternative can be specified.
The upper approximation area (G þ ) is the region which consists of the ideal alternative (A þ ), and the lower approximation area (G À ) is the region which consists of the antiideal alternative (A À ). The belonging of alternatives to the approximation area is determined based on Equation (19).
Based on the MABAC method, for the alternative to be selected as the best in the set, it means to have as many criteria as possible belonging to the upper approximation area (G þ ).
Step 6. Ranking the alternatives In the final step of the MABAC method, the value of criterion functions for each alternative is determined based on the sum of the distance of the alternatives from the border approximation area (q i ) as shown in Equation (20).
Where n is the number of criteria, and m is the number of alternatives. 3.3.4. Topsis model TOPSIS method was first proposed by Hwang and Yoon (1981). The basic principle is to calculate the distance between the evaluation object and the ideal solution and rank the evaluation object according to the distance (Izadikhah et al. 2014). The ideal object was defined as the object closest to the positive ideal solution and also farthest from the negative ideal solution (Lin et al. 2008).
TOPSIS is a practical method that compares options according to the value of data in each weight of the criteria. The fuzzy TOPSIS method is more compatible than the nonfuzzy TOPSIS method due to fuzzy sets. The use of fuzzy function due to flexible composition and other integrated methods will lead to better results in flood zoning and provide results closer to reality. To determine a flood zoning map based on the AHP-TOPSIS method, the following steps are performed.
Step 1. Formation of initial decision matrix (X) This matrix is the same as the initial decision matrix in the MABAC model.

Step 2. Normalize the decision matrix (N)
To normalize the decision matrix, the division of each value by the root of the sum of squares of the values of the column element is used (Equation (21)).
Where X ij is the value of the decision matrix in column number i and row number j, and r ij is the dimensionless value of this entry.
Step 3. Formation of a normal weighted matrix In this step, the weight of the criteria obtained from the AHP method is multiplied by the normal matrix's values to obtain the weighted matrix.
Step 4. Finding the ideal and non-ideal solution In this section, the criteria are divided into two groups of positive and negative.
According to the study's purpose, the factors that increase the potential of flooding occurrence in the region are positive criteria including rainfall, drainage density, soil type, geology and land use. In contrast, the factors that lead to a decrease in flood potential are negative criteria, including the distance from the streams, slope and height.
Step 5. Calculating the distance from the ideal and non-ideal solution In this step, each pixel's distance from its positive and negative ideal is calculated (Equations (22) and (23)).
Where d þ i and d À i are distance from the ideal and non-ideal solution for pixel i, respectively, v ij is the value of pixel i corresponding to the criterion j, v þ j and v À j are positive and negative ideal values of criterion j, respectively.
Step 6. Pixels ranking The similarity index is determined using Equation (16), and the score of each pixel is determined based on the AHP-TOPSIS method. The closer this index is to the number 1, indicates higher pixel potential in terms of flood risk.

Performance assessment
Validation of the prepared susceptibility maps is crucial in the identification and development flood-prone areas . In this study, the accuracy of the produced flood susceptibility maps by different models was verified using the receiver operating characteristics (ROC). The ROC curve is a standard method to evaluate the result's accuracy and has been used by many researchers (  The ROC curve indicates the model's sensitivity to the percentage of pixels that were correctly predicted by the model versus the percentage of pixels predicted to total pixels. High sensitivity indicates many true predictions (true positives), while the high specificity shows a low number of false positives. In this curve, the false positive rate (1-specificity) is displayed on the X axis (Equation (25)) and the true positive rate (sensitivity) is displayed on the Y axis (Equation (26)) ). The area under the curve (AUC) of the prediction rate describes how well the model predicts the flood (Pourtaghi and Pourghasemi 2014). The qualitative relationship between AUC and prediction accuracy of model can be classified into five classes as follows: 0.5-0.6 (poor); 0.6-0.7 (average); 0.7-0.8 (good); 0.8-0.9 (very good); and 0.9-1 (excellent) (Yesilnacar 2005).

Feature selection method
Determination of importance of conditioning factors in flood modelling is a necessary step. In this study, Information Gain (IG) test was used to find the most effective conditioning factors in flood susceptibility mapping (Tien Bui et al. 2012): 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 (Chapi et al. 2017). The value of gain ration varies between 0 and 1, in which the value close to 1 indicates high prediction power of flood factors, while a value of 0 means that flood conditioning factors do not effect on flood occurrence. Table 3 shows the results of the IG test to determine the importance of each flood conditioning factor. As it indicated, elevation (IG ¼ 0.612) with higher value had the most impact on flood modelling, followed by distance from rivers (IG ¼ 0.597), drainage density (IG ¼ 0.544), slope (IG ¼ 0.475), rainfall (IG ¼ 0.436), soil type (IG ¼ 0.311), geology (IG ¼ 0.127) and land use (IG ¼ 0.081), respectively.

Application of the AHP and WOE methods
The results of prioritizing the effective factors in flooding using the AHP method and the relationship between flood occurrences and each class of effective factors using the WOE model are presented in Table 4. As indicated, all of the quantitative and qualitative factors were classified into five classes based on natural breaks in the ArcGIS environment (Lin et al. 2020;Ogato et al. 2020). In the next columns, the number of flood pixels (training floods) and total pixels for each class were specified. Using WOE model, the final weight (Equation (9)) of each class was obtained. These weights imply the level of sensitivity of the classes to flooding.
According to the AHP analysis, among the eight effective factors in determining the flood potential, the rainfall, distance from rivers and slope with the weight of 0.349, 0.208 and 0.128 respectively, have the most impact on flood hazard (Fernandez and Lutz 2010;Youssef et al. 2015). On the other hand, geology, elevation, land use, drainage density and soil type factors ranked next. The calculated consistency ratio of pairwise comparison was 0.084 which was less than 0.1, so the estimated weighting coefficients were acceptable.
Floods occur mainly in rainy conditions. The results show that with increasing rainfall, the sensitivity of the floors to flooding increases. For rainfall factor, the last class including 74-79.85 mm had the highest C=Sc (7.695) and the maximum impact on flooding. For the factor of distance from the rivers, the fourth and fifth classes with the lowest distance from rivers had the greatest influence on the occurrence of floods. The first two classes, due to more distance from rivers, had a negative impact on flooding. For the third effective factor or slope, the last (0-4) class had the highest influence on the floods, with a maximum C=Sc of 9.213. It is clear that by increasing the slope, runoff flows rapidly and the probability of flooding decrease. In the case of soil type factor, the last class or sandy soils with more permeability had the most weight (6.868) and soils with lower permeability had a lower influence on flooding. For the geology factor, the impermeable rock layers related to the Early Cretaceous period had the highest sensitivity to floods due to the low speed of water transfer to the groundwater with the highest weight (6.842). Another influential factor is elevation. The last class had the maximum influence on the occurrence of floods. It is indicated that with increasing the elevation, the sensitivity of the flooding decreases. The results of land use factor showed that the herbaceous plants and groves have the highest weight (7.450). As it indicated, although urban and rural areas have the highest level of flood risk, a small number of floods occurred in these areas due to their great distance from the rivers. Also, dense forests and mountainous areas had a positive impact on flooding due to their lower distance from rivers. The results of the drainage density parameter showed that areas with higher drainage density are more sensitive to the occurrence of flooding, so that the last class has the highest final weight (12.539).

Flood susceptibility mapping using three hybrid models
After determining the weight of the effective factors in the occurrence of a flood using the AHP method and multiplying it by the weight of the mentioned factors, which was obtained using the WOE method, according to Equation (24), weight maps were collected and the final flood susceptibility map (FSM) of the study area was prepared in the GIS environment. If the final weight of each factor class (last column of Table 4) is multiplied by the value of each pixel of MABAC and TOPSIS methods (Equations (20) and (24)), the final flood susceptibility maps based on the MABAC-WOE-AHP and the TOPSIS-WOE-AHP methods were prepared.
The risk level of each factor class was divided into five categories (Very Low, Low, Moderate, High and Very High) using Quantile method which is the most popular method for classification of susceptibility maps (Osaragi 2002;Tehrany et al. 2015;Khosravi et al. 2019). The final susceptibility maps are shown in Figure 5.
As indicated in Figure 5, a large number of floods occurred near the rivers. A common feature of flood zoning maps in all three methods is the very high risk of waterways. Waterways are the main route for passing floods and are the first place of flooding risks. Researches had shown that a large part of flood damages is due to the development of construction near the rivers and the lack of attention to the structures locating programs in the Mazandaran province (Boudaghpour et al. 2015;Sadeghi-Pouya et al. 2017). According to flood susceptibility maps, northern part of the watershed has very high and high flood risk due to presence of low lands and higher rainfall and more moisture, because of low distance from Caspian Sea. Also, in these parts, underground water level is very high and most of soils are saturated. Therefore, a heavy rainfall leads to a large amount of runoff flows in the rivers and finally flood has been occurred. While, low and very low classes are mostly located in the southern part of watershed, due to presence of mountainous areas, steep slopes, and also less moisture. One of the positive points of the map (b) and (c) is the distinction of waterway lines in the flood zoning mapping, while the map (a) is continuous and waterway lines are not visible in it. Map (a) is affected by the rainfall, slope and elevation maps and its risky areas (high and very high) are compatible with the high and very high levels of flood potential in mentioned factors. On the other hand, combined MABAC and TOPSIS methods had affected less from these three factors, due to the fuzzy normalization.
According to Table 5, based on the AHP-WOE method, 62.05% of the watershed was at high and very high flood risk levels. The TOPSIS-WOE-AHP and the MABAC-WOE-AHP methods estimated 55.23% and 59.99% of the Babolroud Watershed at high and very high risk levels, respectively. Also, the estimated value of low and very low classes for the AHP, TOPSIS and MABAC based methods were 19.32%, 24.69% and 16.85%, respectively. As a result, all methods predicted more than half of the watershed in high and very high risk classes and less than 25% of region in low and very low classes. The maximum difference in estimating flood-prone areas between three methods was 10.9%, which shows similarity in flood risk prediction.
Accurate verification can reassure the readers about the more accuracy of one method than the others. In the next section, the ROC method was used to validate the results.

Validation of the flood susceptibility maps
For validation of three flood susceptibility maps, the prediction rate curve was produced by 30% of the data (test data, which were not used in modelling) which contains 15 flood points. The results of the prediction rate curve are indicated in Figure 6. The AUC value for the AHP-WOE, the TOPSIS-WOE-AHP and the MABAC-WOE-AHP models was 0.728 (72.8% of accuracy), 0.905 (90.5% of accuracy) and 0.816 (81.6% of accuracy), respectively. The highest degree of accuracy was for the TOPSIS-WOE-AHP model due to the highest area under the curve. It can be highlighted that all the models have a reasonable performance of flood prediction in the study area. Also, a novel MABAC-WOE-AHP model estimated the floodprone areas with very good accuracy, which indicates the high efficiency of this model.
As can be seen in Figure 6, when X-axis value is equal to 10, it means that 90% of total non-flood points classified correctly in the non-flood class, and only 10% of flood points classified incorrectly as non-flood points. Also, the sensitivity is equal to 85%, 40%, 35%, and for the TOPSIS, the MABAC and the AHP based models, respectively. By increasing the sensitivity, the number of flood points that classified correctly in flood classes increase, so the accuracy of TOPSIS model is more than other two models.

Discussion
In Iran especially in the North, enormous floods have been happening every year. Due to its devastating consequences, the first step leads to less causalities is anticipating flood prone areas. One of the most important preventive measures is preparing flood potential mapping. Thus, flood prone areas can be identified, and the proper measures can be launched to decrease flood damages. In the present study, the capability of three ensemble models based on MCDA was tested to predict the flood susceptibility in the study area.
However, flood hazards are affected by different influential factors. Moreover, not all factors have the same impact of flooding; thus it is essential to identify appropriate flood conditioning factors. According to feature selection method, all factors were found effective on flooding, due to the IG values more than 0. These findings are compatible with results of Rahman et al. (2019) and Costache et al. (2020). Elevation was selected as the most effective factor, followed by distance from rivers, drainage density, and slope. In contrast, the least effective factor in flood modelling was land use. These results are Figure 6. The ROC curve for the flood susceptibility map of the study area using different methods consistence with other studies (Tehrany et al. 2015;Khosravi et al. 2018Khosravi et al. , 2019. While, Termeh et al. (2018) reported that land use, slope and distance from rivers have the greatest impact on flooding, respectively. Besides, Ahmadlou et al. (2018) concluded that elevation was among the most important flood factors, which is also concurred with our result.
Analysis of the results show that the TOPSIS-based model has the highest capability in predicting flood potential in the study area, followed by the MABAC and the AHP based models. The main reasons for having less accuracy in the AHP method are: (1) in the TOPSIS and the MABAC methods, the exact value of each pixel for all criterion was calculated, while in the AHP method, all pixels of a class in each criterion have the same values.
(2) According to Shafapour Tehrany et al. (2017) and Rahmani et al. (2019) ensembles methods predict the flood potential more accurately. So, combining the TOPSIS and the MABAC methods with the WOE and the AHP might lead to higher precision than the AHP-WOE method and output might be more reliable.
The results also show that the TOPSIS method outperforms the MABAC method. These methods have some similarities and differences. In both models, pixels values were normalized and the normal weighted matrix was formed, by using the AHP method. The main difference between two models is the way of ranking the pixels and their level of risks. In the TOPSIS method, the pixels were ranked by calculating the distance from positive and negative ideal solutions. Impacting factor for this distance, is the maximum and the minimum value of each criterion. Meanwhile, in the MABAC method, border of ideal positive and negative border was determined by Equation (16). This border causes discontinuities in the flood zoning map. However, map based on TOPSIS method is continuous, which leads to risk level of each region determined more clearly and more reliable. The reason of the higher value of the AUROC in the TOPSIS method can be considered as more coherence in the flood zoning map.
In other works, different kinds of MCDA and bivariate statics methods were employed and compared with each other. Khosravi et al. (2019) used several methods (FR, WOE and AHP) to determine the flood susceptibility at Haraz Watershed, which located in Mazandaran Province. They indicated all methods have a reasonable accuracy due to the high AUC value. In another study in Ethiopia, Ogato et al. (2020) shows the AHP method predicts the flood prone areas with a high accuracy. Also, Suthirat et al. (2020) indicates that AHP-GIS flood hazard map was in a good agreement with both the historical flood events and the annual maximum rainfall in Ayutthaya city in Thailand. For analyzing the flood risk in Vietnam, Luu et al. (2019) exerted a TOPSIS-based method. Results revealed that the proposed model can be used for sustainable management of flood prone areas with high accuracy. Using the WOE model by Rahmati et al. (2015) leads to evaluate the flood risk with reasonable result (AUC ¼ 74.74%) in Golestan Province in Iran. Application of the MABAC method for selecting the optimal locations for wind farms for the first time was validated with some other MCDA techniques by Gigovic et al. (2017). The results showed high capability of the MABAC method in its selection.
Our results are in line with previous studies, which combined the AHP and the TOPSIS with the WOE methods had high efficiency in our study area. Also, the novel proposed method (MABAC-WOE-AHP) can be useful for mapping flood potential with similar settings, due to high performance.
Considering that a large part of the basin is in the middle to high flood hazard, it is essential to manage floods. In general, the flood management cycle consists of four components: prevention, preparedness, response and rehabilitation, and therefore, in order to apply proper management. In flood hazard management, all aspects of the issue should be considered, from policies in development applications affecting flood intensification to coping strategies such as construction of control structures, flood warning, evacuation, rescue, and rehabilitation measures.

Conclusion
One of the crucial steps to reduce the harmful effects of floods is to diagnose flood prone areas and prepare flood susceptibility maps. The need to have a valid and accurate technique to determine the flood-prone areas prompted the authors to assess several methods to understand the methods' efficiency. The study focussed on analyzing flooding hazards using GIS-based and the MCDA methods in Babolroud Watershed, Iran. Therefore, combined AHP, TOPSIS and novel MABAC methods with weight of evidence model were used. For this, a flood inventory map including 50 flood locations was constructed. These locations were divided randomly into two groups of training data (35 locations) and test data (15 locations). Then, eight data layers containing rainfall, distance from rivers, slope, soil, geology, elevation, drainage density and land use were derived from the spatial database. The results of the AHP-WOE, the TOPSIS-WOE-AHP and the MABAC-WOE-AHP methods indicated respectively 62.05%, 55.23% and 59.99% of Babolroud Watershed was located in high and very high risk of flooding area. In order to validate the results of these methods, the ROC curve was employed for test data. The validation of the results showed that the AUC value for the AHP-WOE, the TOPSIS-WOE-AHP and the MABAC-WOE-AHP methods was 0.728 (good accuracy), 0.905% (excellent accuracy) and 0.8162% (very good accuracy), respectively. The map of this analysis was indicated the distribution of high-risk areas, which is essential for flood management strategies. Once the spatial distribution of flood hazard zoning is determined, the responsible organizations, planners, and engineers would be more cautious with preventive operations such as strengthening structures in floodplains against floods, developing flood safety guidelines, providing general education, providing flood insurance, and flood mitigation plan and implementation. By comparisons in the present study, for the future work, applying some machine learning models such as artificial neural network, decision trees, and support vector machine and integrated with AHP can help to scientists and governments to achieve more precise data.