Flood susceptibility mapping using an improved analytic network process with statistical models

Flooding is a natural disaster that causes considerable damage to different sectors and severely affects economic and social activities. The city of Saqqez in Iran is susceptible to flooding due to its specific environmental characteristics. Therefore, susceptibility and vulnerability mapping are essential for comprehensive management to reduce the harmful effects of flooding. The primary purpose of this study is to combine the Analytic Network Process (ANP) decision-making method and the statistical models of Frequency Ratio (FR), Evidential Belief Function (EBF), and Ordered Weight Average (OWA) for flood susceptibility mapping in Saqqez City in Kurdistan Province, Iran. The frequency ratio method was used instead of expert opinions to weight the criteria in the ANP. The ten factors influencing flood susceptibility in the study area are slope, rainfall, slope length, topographic wetness index, slope aspect, altitude, curvature, distance from river, geology, and land use/land cover. We identified 42 flood points in the area, 70% of which was used for modelling, and the remaining 30% was used to validate the models. The Receiver Operating Characteristic (ROC) curve was used to evaluate the results. The area under the curve obtained from the ROC curve indicates a superior performance of the ANP and EBF hybrid model (ANP-EBF) with 95.1% effi- ciency compared to the combination of ANP and FR (ANP-FR) with 91% and ANP and OWA (ANP-OWA) with 89.6% efficiency.


Introduction
The rapid increase in urbanization and the effects of climate change have led to numerous environmental problems and disasters in recent years Malik et al. 2020b). Floods are one of the most destructive natural disasters (Doocy et al. 2013), causing significant loss of life worldwide and financial damage that equates to 40% of the total economic loss caused by all-natural disasters annually . Any type of land-use change resulting from urbanization leads to an increase in flooding due to the rise in the rainfall-runoff coefficient (Myronidis et al. 2016;Dammalage and Jayasinghe 2019). Preventing floods is complicated, but their occurrence can be predicted and partially controlled using appropriate methods and measures (Mansur et al. 2018;Malik and Pal 2020). Having access to accurate and current information is key to preventing floods or at least reducing their harmful effects. One of these key pieces of information is a flood potential map (Papaioannou et al. 2018). Morphological changes caused by human activities can have a significant impact on changing the conditions for floods (Yousefi et al. 2018).
So that, urbanization, deforestation, and the impact of climate change on rainfall characteristics in susceptible areas may all contribute to an increase in flooding (Kundzewicz et al. 2014). Many factors are important in determining flood risk zones, the most important of which are land use/land cover, geology, slope, drainage network and morph metric factors in general . Also, the river network in an area can be one of the most effective reasons for the morphological change of floods, and this morphological change will be caused by human interaction and communication with the rivers (Yousefi et al. 2019).
Present flood control measures have failed to mitigate flood problems. It is not possible to eliminate the flood problems, but with proper flood forecasting system and management, damages can be minimized (Das et al. , 2019. Flood zoning provides valuable information regarding the nature of the floods, their effects on the floodplain, and determining the riparian zone (Roopnarine et al. 2018). GIS is a useful tool for investigating multidimensional events, such as floods, at various spatial and temporal scales. Spatial analysis with GIS has been used to determine flood zones (Chowdhuri et al. 2020). Also, the use of remote sensing data, along with improved modelling techniques, such as data mining and machine learning, could help to determine flood-prone areas better.
Several variables affect the occurrence of floods in a watershed and determining the flood catchment areas is an important factor in investigating the floodplains. The most important factors and variables are physical and natural factors (e.g. elevation, slope, aspect, slope curvature, lithology, topographic position index and rainfall), hydrological factors (e.g. drainage density, river distance, topographic wetness index, and stream power index) and human disturbance factors (e.g. land use and road distance) (Al-Juaidi et al. 2018;Zhao et al. 2018;Darabi et al. 2019;Costache 2019a). Recently, several methods and models have been used to study flood hazards. The wildly applied methods and models for flood hazard mapping are Logistic Regression (LR) Lee et al. 2018;Rahmati et al. 2019;Malik et al. 2020a), Weights of Evidence (WofE) , Support Vector Machine (SVM) Tavakkoli Piralilou et al. 2019;Ghorbanzadeh et al. 2019c), Random Forest (RF) Chen et al. 2019;Darabi et al. 2019;Shahabi et al. 2019), Frequency Ratio (FR) (Lee et al. 2018;Rahman et al. 2019), Decision Tree (DT) (Chen et al. 2019) and Naïve Bayes Tree (NBT) (Khosravi et al. 2019). In addition to these models, multi-criteria decision-making techniques such as the analytic hierarchy process (AHP) (Satty 1980) and the Analytic Network Process (ANP), which is an extended version of the AHP and consider as a more appropriate approach for complex decision problems can also be used for flood hazard mapping (Ahmadlou et al. 2019;Kanani-Sadat et al. 2019). Although these models have criticized for their some inabilities and associated uncertainties, they have optimized by integrating with several statistical and mathematical techniques (Ghorbanzadeh et al. 2019b). The statistical techniques of the sensitivity and uncertainty analyses have been applied along with the AHP and ANP by using the Monte Carlo simulation (MCS) (Ghorbanzadeh et al. 2018b(Ghorbanzadeh et al. , 2019a. The frequency ratio (FR) also has combined with the AHP in the hybrid spatial multi-criteria evaluation (SMCE) method to enhance the overall accuracy of criteria weightings. Moreover, Ha et al. (2017) applied the AHP, which was integrated with fuzzy-TOPSIS for their study. The fuzzy logic has been integrated with both the AHP and ANP to evolve their results in a number of studies to minimize the uncertainty associated with decisionmaking procedure. Ghorbanzadeh et al. (2018a) applied the interval matrices for pairwise comparisons of AHP, which could be helpful for the cases that a lot of experts with a high variation of preferences are involved in complex decisions. The combination of GIS and multi-criteria decision-making has already been used in spatial modelling and disaster analysis of natural disasters such as floods (Paquette and Lowry 2012;Sol ın 2012;Kazakis et al. 2015). The primary purpose of this research is flood susceptibility mapping using improved ANP with statistical models. Moreover, we evaluated the efficiency of the FR, EBF and OWA methods in providing flood probability and susceptibility maps for Saqqez City in north-west Iran.

Description of the study area
The study area is Saqqez City, which is located in the Zagros mountains in the north of Kurdistan Province (51 45 46 51 east longitude and 35 46 36 36 north latitude) (Figure 1(b)). The city has an area of about 4514 km 2 and a population of 226.451(2016 Census). The average annual rainfall in the study area is 450 mm. Rainfall usually starts in the springtime, and during June, July, and August the city is at risk of flooding. There are many streams and also many springs which underlie the genesis of the rivers located across the study area. The notable rivers of the study area are Zarinehrood, Khorkhor, and Saqqez (Figure 1(c)). Due to a lack of flood control measures and improper flood management, the leading causes of flooding in the urban and rural areas are the rivers (Figure 2).

Proposed process
Flood susceptibility mapping (FSM) in this study was performed with a new approach, by combining the ANP multi-criteria decision model with FR, EBF, and OWA statistical models in five stages ( Figure 3). Collecting data and information from the study area to prepare a map related to Flood Conditional Factors (FCFs) and a flood inventory map. Multi-collinearity test was performed to ensure the independent effect of each FCFs. Determining the relationship between the locations of occurring floods and FCFs classes using Frequency Ratio (FR) and Evidential Belief Function (EBF) algorithms. Using a linear regression relationship and FR values for each class of factors, the initial weight values for FCFs were calculated. The importance (final weight) of each of the factors used in this study was determined using the initial weights in the Analytic Network Process (ANP) model. The calculated values of FR and belief (Bel) were applied to each of the sensitivity classes in each factor. For Flood susceptibility modelling, new combined models ANP-FR, ANP-EBF, and ANP-OWA were made by applying secondary weight to each of the factors. The performance of flood forecasting models was evaluated using Receiver Operating Characteristic (ROC) curve and area under the curve (AUC), standard error (SE) and confidence interval (CI) values.

Flood inventory map (FIM)
A flood distribution map is one of the most critical layers in flood susceptibility mapping, which is based on the analyzes of correlations between factors affecting the flood and the location of the flood (flood zones) in the region . Also, considering that in the study area, the conditions of floods (geographical, climatic, etc.) have not changed much compared to previous years; previous flood points can be used to predict future flood zones. In total, 42 flood events have been recorded in the Saqqez Catchment using the flood line and field investigation. The data regarding flood locations were extracted from the archive of the Water Resources Management Organization which provided the flood locations as points. It should be mentioned that this institution is the official authority in Iran in the domain of hydrology and water management. The use of flood locations as points is also recommended by the high-quality results and models performances higher than 90%, that were achieved in the previous studies related to the flood susceptibility assessment (Chapi et al. 2017;Choubin et al. 2019;Bui et al. 2019a;Ali et al. 2020;Yariyan et al. 2020c). Also, after obtaining data about the flood locations, using field research and regional documents, the occurrence of floods in each place was ensured. We used 70% of these flood location points to map potential flooding and the remaining 30% to validate the models (Costache 2019b).

Flood conditional factors (FCFs)
Determining the factors affecting flooding is one of the essential steps in mapping flood potential (Pham et al. 2016). In this study, the 10 factors influencing floods are slope, slope aspect, curvature, altitude, rainfall, distance from river, topographic wetness index, slope length, lithology and land use/land cover (LU/LC). By collecting information about 10 factors affecting the flood, all maps with a pixel size of 90 Â 90 and a scale of 1: 600,000 were converted to raster format. Each factor was divided into several classes ( Figure 4). The map classification related to each of the factors was done by the Natural break method to perform FR analysis and better display. ArcGIS and SAGA GIS software was used to provide the required layers by using a digital elevation model (DEM) with a spatial resolution of 12.5 meters. The DEM was downloaded from the ALOSPALSAR sensor from Landsat images. When high-precision topographic data are used, flood analysis has good predictive accuracy . Therefore, the elevation layer and its derivatives have a significant impact on the prediction of flood-prone areas. The altitude class map is based on the digital elevation model and comprises 5 classes (Figure 4(f)) (Bui et al. 2019b). The slope is the most critical factor for flooding due to its direct influence on surface runoff in watersheds. The effect of slope on the occurrence of devastating floods is related to runoff speed (Yariyan et al. 2020c). So that, on slopes with high slopes, the speed of water penetration is much lower, this will increase the speed and flow of water. The slope map of the study area is based on the digital elevation model and is classified into 5 classes  (Youssef et al. 2016). The slope aspect is a further factor influencing flooding. This map in 9 class (flat, northeast, east, southeast, south, southwest, west, and northwest) was created using the digital elevation model (Figure 4(e)) (Bathrellos et al. 2018).
The slope length layer produced by the combination of the catchment area and the slope layers affects the water flow. We used this layer to examine the influence of the topography on the flood flow. The topographic wetness index (TWI) is used to measure topographic control in hydrological studies (Chen and Yu 2011). The TWI factor indicates the amount of cumulative flow associated with a point in the watershed that tends to move to less sloping areas under the influence of gravity . Therefore, it cans gravity to increase water flow velocity and destructive power. Also, the amount of spatial distribution for runoff production will be identified by measuring TWI. Therefore, the TWI parameter related to water flow can be calculated through the following equation: where A s is the catchment area and b is the slope (in degrees). The TWI and slope length (SL) were prepared using SAGA GIS software version 3.2 ( Figure 4(c,d)). The rainfall layer was developed based on the collected data from seven rain gauge stations (Figure 1(b)) Baneh, Boein, Bastam, Kelahshin, Vezmaleh, Miradeh, and Saqqez for the period 1999-2018. The rainfall map in five classes from 290 to 550 mm have been classified and prepared by the data of the Meteorological Organization of Iran (Figure 4(b)). We used the inverse distance weighted (IDW) method in GIS to interpolate the precipitation stations in the study area. Distance from the river is an important and influential factor in reducing or intensifying the occurrence of floods (Khosravi et al. 2018) because the areas adjacent to the river are more prone to floods. In this study, the river network was extracted using the 'Arc hydro' integration tool in ArcGIS through DEM, and then a layer of distance from river was prepared using the 'Euclidean distance' tool ( Figure 4(h)). Land use/Land cover plays an essential role in the flooding process and directly or indirectly affects some hydrological processes such as permeability, evapotranspiration, and runoff (H€ olting and Coldewey 2019). When it rains in a basin, the amount of rain that enters the rivers depends on the conditions of the area, topography, and LU/LC. The land use/land cover map of the study area was prepared using the Landsat Eight image and OLI sensor on 2018.06.20. The generated LU/LC map has 7 classes, namely: garden, good rangelands, dry farming, moderate rangelands, poor rangelands, urban area, and water ( Figure 4(j)). The geology of the area plays a key role in mapping the potential of floods because many geological units are active in hydrological processes, also can have a great impact on the conduction and penetration of water flow. The geological layer was prepared using a 1:100000 map of Kurdistan Province and include 26 geological units (Figure 4(i); Table 1).
The information on formations in this area is shown in Table 1. The distance from the river map was prepared using GIS software and classified into 5 classes ( Figure 4(h)).

Analytic network process (ANP) model
The purpose of decision-making is to choose the best option by weighing up the influencing factors. Multi-criteria decision-making methods use multi-criteria rather than an optimal measurement criterion (Ghorbanzadeh et al. 2018a;Mitroulis and Kitsios 2019). All elements must be well understood and precisely defined to model and analyze a multi-criteria system (Yariyan et al. 2019). The analytic network process is one of the multi-criteria decision-making techniques that fall into the category of compensatory models (Alilou et al. 2019). This model is designed based on the hierarchical analysis process, and the network has replaced the hierarchy (Peng and Tzeng 2019).
One of the assumptions of the hierarchical analysis process is that the higher sections and branches of the hierarchy are independent of the lower sections and branches (Saaty 2004). In many decisions, however, decision elements cannot be modelled hierarchically and independently of each other, therefore, different elements are combined. Saaty (2004) suggests using the ANP technique to this end. In the ANP, the relationships between different levels of decision-making are considered one-way. The advantage of the ANP approach is that the various elements are measured based on their inter-relationships (Ghorbanzadeh et al. 2018b). Unlike the AHP model, which uses an entirely hierarchical structure for modelling, the ANP model uses a pairwise comparative measurement scale and performs modelling using the system perspective feedback (Alizadeh et al. 2018).
In this ANP model, the initial value for each factor is determined on a scale of 1 to 9. After the pairwise comparisons matrix is formed, the local priority vector is calculated using the following formula: Where A is the pairwise comparison matrix and kmax is the largest value in this matrix, which is calculated as follows: where a ij is the value of the factor and W j is its weight. The weights of any pairwise comparison matrix can be calculated separately. Moreover, the resulting weights are organized in an overall supermatrix for evaluating the final weights of criteria (Ghorbanzadeh et al. 2018b).
The resulting values from the limited matrix are considered as the final relative weight W ¼ ðw1, w2, ::, wnÞ that is calculated using the supermatrix. The whole process of the limit matrix is fully explained by Saaty (2007).
After carrying out the ANP, the value of the CR compatibility index is calculated (Equation (5)). If its value is CR 0.1, the matrix is compatible.
where n is the rank of the arbitration matrix, and RI is a random index type with the following values (Alizadeh et al. 2018) ( Table 2).

Frequency ratio (FR).
The frequency ratio is a bivariate statistical model that can be used as a simple spatial tool for calculating the probable relationship between independent and dependent variables and includes several classified maps. This method is used to map the flood susceptibility . The FR approach is based on the observed relationship between the distribution of flood point locations and flood-dependent criteria (Yariyan et al. 2020b). The FR of each class of each criterion is calculated according to the following equation.
In this equation, FR-the value of frequency ratio given to a class/category i of parameter j, Np ðLXiÞ is the sum of flood locations inside a class i of a predictor X, NpðXjÞ is the sum of locations inside a predictor Xj, m is the sum of classes contained by a flood predictor Xi, and n is the sum of flood predictors inside the study area. Flood susceptibility index (FSI) for a cell is equal to the sum of the frequency ratio of that cell in all variables. If there are M effective factors, the susceptibility index of flood occurrence is as follows: Finally, the rates obtained for each class were applied to this layer in the geographic information system for this class, and a raster map was used to predict the likelihood of flooding.
In this study, FR was used to obtain the true weight of each factor. First, by classifying the factors affecting the flood, the FR values for each class of factors were calculated (Table 3). Then, we determined the initial weight of each factor using a linear regression equation between the percentage of flood occurrence points and the  (8)) shows the effect of an independent variable on a dependent variable and measures the correlation between them .

Evidential belief function (EBF)
The Evidential Belief Function (EBF) data-driven model is based on Dempster-Shafer theory, which was first introduced by Dempster and developed by Schafer (Ghorbani Nejad et al. 2016). The most important advantage of the Dempster-Shafer theory is its ability to combine uncertainties from different sources of evidence and its relative flexibility to accept uncertainties. The Dempster-Shafer theory provides a framework for estimating EBF under Dempster law. The evidential belief function is a combination of disbelief, degrees of belief, uncertainty, and probability (Bui et al. 2012).
Higher uncertainty values indicate a stronger correlation between the factors affecting the flood and the likelihood of flood occurrence, and vice versa (Pourghasemi and Kerle 2016). The following Equations (9)-(14) are used to calculate these four indices: where Bel Cij (degree of belief) is the certainty value, Dis Cij (degree of disbelief) is the uncertainty value, N C ij \ D ð Þ are the density flood pixels in class D, NðC ij Þ is the total density of the floods that occurred in the study area, N D ð Þ is the pixel density in class D, Unc is degree of uncertainty, Pls is the degree of plausibility and N T ð Þ is the pixel density in the whole study area.

Ordered weight average (OWA)
There are different weight composition methods for calculating the final weight in multi-criteria decision-making methods. Here we used the ordered weight average (OWA), which is proposed by Yager (1988). We used the OWA method to weight the criteria and prioritize the criteria by combining the weighting with the prioritization of evaluation criteria. Prioritizing weights enable the direct control of criteria. The main feature of this method is the possibility to reclassify weights. The coefficients of weights are applied to the prioritized position of the criteria values for the chosen alternative instead of applying them to the criteria directly (Maldonado et al. 2018).
The OWA operator depicts n dimension space in one dimension (Liu et al. 2008). There are three basic operations (or integration functions) based on the bellman model for fuzzy set integration operations: fuzzy set sharing operators, fuzzy set community operators, average operators (Yager 1988). Figure 5 illustrates the decision process of OWA (Salman Mahini and Kamyab 2012).
In Figure 5, the x-axis is the maximum precaution, and the y-axis represents the trade-off between criteria. Weighted linear combination (WLC) is based on the weighted average, in which the criteria are combined in a standard common numerical interval using the weighted average . The ordered weighted average method can provide prediction maps and shapes for various evaluations by varying the order and benchmark parameters (Ferretti and Pomarico 2013).
OWA is defined as a hybrid technique with an ith location, a set of order weights v ¼ v 1 , v 2 , :::, v n ðv j 0, 1 ½ , j ¼ 1, 2, :::, n, and P n j¼1 wj ¼ 1) and a set of criteria weights w ¼ w1, w2, w n ðwj 0, 1 ½ , and P n j¼1 wj ¼ 1Þ (Motlagh and Sayadi 2015). In general, the OWA model with two parameters, AND ORness and TRADEOFF, is defined as follows: ORness ¼ 1 À ANDness TRADEOFF ¼ 1 À ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n P W order i À 1=nÞ 2 n À 1 s (17) where n is the number of factors, i is their order and the W orderi is the weight of the order factor.
The OWA operator is calculated with the following equation: where Zin ! ::: ! Zi1 is specified by the order of the criterion (xij), Vj is the weight of the order, and Wj is the weight criterion (Malczewski 2006). The OWA model requires standardized layers. We used a standardization method based on the minimum and maximum factor values (Equations (18) and (19)). Therefore, the following two equations are used to convert crude values to standard values: where S ij is the value of i for the criterion j, S min j is the lowest value for the criterion j and S max j is the highest raw value for the criterion j .  In this study, the weights of factors obtained using the ANP method ( Figure 6) and (Table 4) were calculated using Equation (21). Then, the order weights obtained from OWA were multiplied by the ordered quantities (Table 5).
where wj is the criterion weight, and a is the degree of optimism OWA obtained from Table 5. The OWA technique is an embedded type of multi-criteria decision-making and enables geographic information systems to perform data analysis and decision-making. Input for this model is raster maps of the study area and descriptive tables.

Models validation
The Receiver Operating Characteristic (ROC) curve was used to evaluate the produced models. The ROC curve is a graphical representation of the balance between the negative and positive error rates for each possible value and is based on the set of pair points (X and Y) calculated using Equations (22) and (23). The X and Y were obtained from the comparison matrix and varied between 0 and 1 (Yariyan et al. 2020a). The area under the ROC curve is called the area under the curve (AUC). The AUC gives the value of the system's prediction by describing its ability to accurately predict events (floods) and non-events (non-floods). The values are a percentage of true reality and the percentage of false reality of the graph is calculated according to the following equations.
where TN, FP, TP, and FN are true negative, false positive, true positive, and false negative (Leman et al. 2020).

Multi-collinearity test of geo-environmental factors
Multi-linearity (or collinearity) is a type of statistical phenomenon in multiple linear regression analyses in which two (or more) predictor variables are related or interrelated (Kroll and Song 2013). It is necessary to choose a method that can predict the impact of variables independently in such research (Gudiyangada Nachappa et al. 2019;Yousefi et al. 2020). In this study, the linearity of flood conditioning factors was investigated using the Variance Inflation Factor (VIF) and Tolerance (TOL). If the TOL <0.1 or the VIF ! 5, there is a linear relationship between the conditioning factors (Equations (24) and (25)) where R 2 j is the coefficient of determination of the regression model .

Multi-collinearity test of geo-environmental factors
The results of the multi-collinearity analysis showed it is not multi-collinearity between the conditioning factors and the effective factors in the flood risk mapping ( Table 6). The highest VIF and the lowest TOL were 1.762 and 0.567, respectively.

Result of the analytic network process (ANP) model
The results of the ANP indicated that the distance from river (H) has the highest weight of all the factors (Figure 7). Following the distance from river (0.147), the order of weights of the other factors is as follows: aspect (0.143), altitude (0.139), slope (0.13), rainfall (0.114), curvature (0.098), LS (0.082), TWI (0.065), geology (0.049), and LU/LC (0.033) ( Figure 6). The value of the compatibility index (CR) is 0.00041 (less than 0.1), which indicates that there is no incompatibility between the factors and the results are thus valid. Based on the results of the ANP and the importance of different conditioning factors, we could suggest the following general equation for a flood susceptibility map: where A, B, C, D, E, F, G, H, I and G are slope, rainfall, slope length, TWI, aspect, altitude, curvature, distance from river, geology, and LU/LC, respectively. The values of A-G are obtained based on the FR, EBF, and OWA. This new approach aims to determine the weights of factors with high accuracy because the initial weight of the factors is determined objectively.

Frequency ratios for different classes of flood condition factors
The factors with the highest value of frequency ratio are shown in Table 6 in bold.
In Figure 6, the x-axis represents the percentage of flood occurrence points in each class and the y-axis represents the calculated FR values in each class. By calculating the coefficient of determination for each of the regression relationships, the initial weight of each factor was determined according to the R coefficient on a scale of 1 to 9. So, the higher the R 2 , the higher the weighting factor. Then, Super decisions software was used to implement the ANP model and to form a pairwise comparisons matrix.

Weight of different classes of flood condition factors based on EBF
In Table 3, N ðC ji Þ is the number of pixel units in each class, N ðC ji \ DÞ is the number of flood points in each class, and Wi represents the weight of the ANP. After applying the Bel weights to each class of flood factors, the ANP-EBF hybrid model was implemented based on Equation (24).
Then, to implement the EBF method to integrate different factors to reduce uncertainty in flood risk mapping, the values (Bel) were calculated for each class of effective flood factors (Table 4). For the combined modelling, all flood conditioning factors were classified based on the EBF weight obtained from Bel.

Flood susceptibility modelling (FSM)
By using the relationship between the flood points and the conditioning factors associated with floods, we created flood susceptibility maps using the hybrid method of the ANP (Figure 8 and Equation (26)) with FR (Table 3) and EBF (Table 4) models. Therefore, the weights of FR and EBF of each susceptibility class of each factor (Tables 3 and 4) were applied in Equation (25). The result of executing the above command in ArcGIS software is the flood risk map in the interval 0-1 (Figure 8(a,b)).
We implemented the OWA model using ANP weights, the weights of the factors, and the order weights. This second set of weights allows for direct control over the risk level. The order weight is a set of weights assigned to a given location (pixels). Applying the weight of the order based on the grading factors affects the minimum to the maximum value for each location. After standardizing and entering layers into the IDRISI software, we implemented the OWA analytical model and the result was a flood risk map with values between 0 and 1 (Figure 8(c)).

Validation
In the ROC curve from the area under the curve (AUC), standard error (SE), and confidence interval (CI) were used to evaluate the efficiency of the models (Figure 9). Evaluation of the models showed that the hybrid models ANP-FR, ANP-EBF, and ANP-OWA yielded performance coefficients of 91%, 95.1%, and 89.6%, respectively. These values indicate the good performance of the hybrid models for flood risk assessment in the study area. However, the ANP-EBF model with the highest accuracy (AUC ¼ 0.951, SE ¼ 0.0494) compared to the ANP-FR (AUC ¼ 0.91, SE ¼ 0.0833), and ANP-OWA (AUC ¼ 0.896, SE ¼ 0.0632) models is best suited to determine the risk of flooding in the region (Table 7). Evaluation statistics of the applied methods are shown in Table 4. According to the results shown in Table 4, the ANP-FR has the lowest standard error (0.49) compared to the ANP-EBF (0.083) and ANP-WOA (0.063) models.

Discussion
Evaluation of the performance of the hybrid statistical models and multi-criteria decision-making showed that the ANP-EBF, ANP-FR, and ANP-OWA models perform well in flood hazard zoning in the Saqqez watershed. We evaluated the performance of the hybrid models based on the AUC, SE and CI criteria, which showed that the ANP-EBF hybrid model was able to map the flood susceptibility with an AUC ¼ 95.1%, SE ¼ 0.0494 and CI ¼ 0.779 to 0.998 efficiency coefficient, compared to the slightly lower values of 91.0% and 89.6% achieved by the ANP-FR and ANP-OWA models, respectively. Because the ANP model is a comprehensive technique that allows all direct and indirect factors to influence decision-making, its combination with other models increases their accuracy. ANP method allows for more complex relationships between decision levels and features to be evaluated because levels do not require a strict hierarchical structure. Because of the interdependence characteristics that exist in natural hazards such as floods, it is also important to consider the interrelationships between the criteria in the decision-making (Alilou et al. 2019;Kanani-Sadat et al. 2019). The ANP method allows us to examine the interdependencies between the criteria levels and is, therefore, an appropriate tool for multi-criteria decision-making (Sun et al. 2016 Also, based on the results of the study (Arabameri et al. 2019), the EBF algorithm has a better performance than FR in predicting flood sensitivity. In addition to the above, according to research (Chowdhuri et al. 2020), the EBF-LR combination for predicting flood sensitivity also has the highest accuracy. These can confirm the findings of our study. Therefore, the use of an EBF in flood risk assessment is useful and reliable. Since most of the floods are caused by the flooding of the main river channel, the distance from the river is one of the most important factors in the flood sensitivity analysis. Areas close to the river are more vulnerable to flooding due to the proximity to the mainstream and rapid response to floods (Chapi et al. 2017;Moradi et al. 2019). Our results showed that the 0-600 meters class in the factor of distance from river has the highest flood risk. Shafizadeh-Moghadam et al. (2018) showed in their study, that the distance from the river factor is one of the most important variables to determine the risk of flooding and that the areas close to the river are most sensitive to flooding. In the aspect factor with its high importance in flood occurrence (0.143), the northwest class has the highest sensitivity (FR ¼ 1.86, Bel ¼ 0.2529) for the flood. By comparing the maps obtained from the models and Figure 4(e), the accuracy of these results can be realized. The flood sensitivity map of the study area shows the northern part of Saqqez City is more sensitive to flooding. Observing Figure 4(f), we find that the northern parts of the study area with the lowest altitude are more prone to floods. This is completely consistent with the analytical calculations of the FR algorithm for the altitude factor. So that the highest sensitivity in the altitude class is less than 1600 meters. One of the primary reasons for this is the possibility of more concentrated water flow at low altitudes. Our result is consistent with the results of Fern andez and Lutz (2010); Tehrany et al. (2015); Termeh et al. (2018); Costache et al. (2020) who also stated in their study that with increasing altitude, the sensitivity of the classes to flooding decreases. Examination of the slope factor shows that the highest effective weight factor in flood sensitivity the class <15%. The higher concentration of water at low slopes increases the sensitivity and the risk of flooding. Vojtek and Vojtekov a (2019) found that one of the effective factors in flood sensitivity is the slope angle and that low slopes are very prone to flooding. As our calculations show, high FR and Bel values are in low slope classes. Therefore, it can be found that by establishing a linear relationship to reach the final weight through the FR algorithm and slope factor sensitivity classes, the high importance of the slope factor (0.13) has been calculated correctly. Therefore, the findings of the study (Khosravi et al. 2018) also emphasize the high importance of the factors of distance from the river, slope, and altitude in the occurrence of floods, which can confirm the weights calculated for these factors in our study.
Rainfall is one of the most important influencing factors in flood sensitivity. According to the rainfall map of Saqqez City, areas with precipitation of more than 450 mm in the northern part of Saqqez City are more prone to flooding. Also, according to the results of the ANP model, the rainfall factor with a weight of 0.114 is very important in the occurrence of floods. In addition, the obtained FR and Bel values for rainfall sensitivity classes indicate the presence of high sensitivity in the class> 450 mm. Therefore, the correspondence of the results in Tables 3 and 4, Figure 7, and the output maps (Figure 8), indicates the high importance of the rainfall factor and the complete agreement of the calculations. After Slope, Altitude, Distance from river, Rainfall, and Aspect factors, the curvature factor is the most important in flood occurrence. Relationships between flood locations and sensitivity classes indicate high sensitivity in the Flat class. Therefore, the prediction of floods on slopes and low altitudes and flat places in this study indicates the high accuracy of calculations and results. The relationship between FR and Bel values in the slope length factor sensitivity classes indicates high sensitivity in class <11. These calculations are fully correlated with the values obtained for the slope factor and the sensitivity at low slopes. So that on a slope with a lower degree and length, the degree of flood sensitivity is higher.
Investigation of the effect of TWI on flood sensitivity showed that areas with high TWI had more of an effect on flood sensitivity in the study area. As the TWI shows the distribution of moisture distribution in different parts, increasing this index increases soil moisture content and increases the probability of flooding. TWI is one of the most important factors affecting flood and that areas with a high TWI are highly sensitive to flooding. The results of Tables 3 and 4 for the TWI factor also show in the classes with the highest value (>16), high sensitivity. The relationship between flood and geological factors in this study shows the high sensitivity of Moderate-grade, regional metamorphic rocks (pCmt1) to flood. In general, the PC group was more effective. So that, after pCmt1, pCmb (Marble) and pCgr (Precambrian granite to granodiorite) groups had the highest flood sensitivity, respectively.
Investigation of the effect of land use/land cover on flood sensitivity in Saqqez City showed that orchards and agricultural lands have a more significant effect on flooding than other land-use types. As agricultural and orchard land-use types are extensive near the main river and downstream of the basin, these areas are most affected by flooding and inundation. So that floods near agricultural lands, despite having nutrients, can be one of the factors of quality and growth of agricultural products.

Conclusion
In this study, flood risk zoning was determined using a combination of statistical methods (FR, EBF, and OWA) and a multi-criteria decision-making model (ANP) in the Saqqez watershed in Kurdistan Province, which is one of the most flood-prone areas in this province. The environmental variables we used to model flood susceptibility were altitude, slope, aspect, distance from river, rainfall, TWI, LS, plan curvature, LU/LC, and lithology. The results of the combination of the statistical method of FR and ANP (ANP-FR) showed that this hybrid method was more efficient than the ANP-EBF and ANP-OWA methods. The analysis of the importance of factors using ANP showed that the variables of distance to river, aspect, and altitude were the more important in determining the flood zone in this area. The results of flood hazard zoning showed that the northern areas of the watershed are highly susceptible to floods, and these are generally are low-lying areas close to the river. Today, financial bottlenecks in watershed management require that a watershed be properly budgeted, and that flood protection and management measures are prioritized to maximize effectiveness. Therefore, surveying, analyzing, and mapping flood susceptibility zones can play an important role in better management of watersheds, budget allocation, and organization of structural and non-structural operations. Moreover, given the fact that, during the last decades, climate change creates the premises for severe flood events across many regions of the globe, we consider that flood susceptibility estimation through highly accurate models is mandatory. Therefore, since the accuracy of the three hybrid models, used in the present study, exceed 89% we have a higher degree of confidence that the present methodology can be successfully used also in other areas, across the Earth, which is exposed to flood occurrence. Thus, the novel hybrid models and the high accuracies of the methods make from this manuscript valuable research which will be of great interest for the potential readers across the entire globe.