Ranking Sub-Watersheds for Flood Hazard Mapping: A Multi-Criteria Decision-Making Approach

: The aim of this paper is to assess the extent to which the Sad-Kalan watershed in Iran participates in ﬂ oods and rank the Sad-Kalan sub-watersheds in terms of ﬂ ooding potential by utilizing multi-criteria decision-making approaches. We employed the entropy of a drainage network, stream power index (SPI), slope, topographic control index (TCI), and compactness coe ﬃ cient (Cc) in this investigation. After forming a decision matrix with 25 possibilities (sub-watersheds) and 5 evaluation indices, we used four MCDM approaches, including the analytic hierarchy process (AHP), best–worst method (BWM), interval rough numbers AHP (IRNAHP), picture fuzzy with AHP (PF-AHP), and picture fuzzy with linear assignment model (PF-LAM, hereafter PICALAM) algorithms, to rank the sub-watersheds. The study results demonstrated that PICALAM exhibited superior performance compared to the other methods due to its consideration of both local and global weights for each criterion. Additionally, among the methods used (AHP, BWM, and IRNAHP) that showed similar performances in ranking the sub-watersheds, the BWM method proved to be more time-e ﬃ cient in the ranking process.


Introduction
Floods are the most severe and common environmental hazard worldwide, associated with an immeasurable impact on humans, properties, and infrastructure [1,2]. In 2019, approximately 317 significant natural disasters occurred worldwide, with floods accounting for 45% of these natural disasters. The continent of Asia has faced more flood disasters than anywhere else [3,4]. On the continent, Iran is located in hazardous flood covers approximately 10,358.9 km 2 ( Figure 1) The case study is one of the main headwaters of the Karun River, which is Iran's most affluent and only navigable river. The watershed contains 25 sub-watersheds. It is situated in a moderate climate in the summer (August) and a cold and snowy climate in the winter (February). The annual average temperature is +9.7 °C, and the annual average rainfall is 313 mm, which occurs mainly between February and April.

Data and Methodology
The main framework of the current study to assess flood hazards and rank sub-watersheds was adopted using a GIS-based multi-criteria decision analysis ( Figure 2). Due to the existence of GIS, scholars can handle large amounts of data and combine the intrinsic or value-based knowledge involved in MCDA. Based on an extensive literature review, our experience, availability of data, field observation, and duplication of factors, seven flood hazard-driven factors, i.e., entropy of the drainage network, slope, topographic control index (TCI), and compactness coefficient (Cc), were investigated to provide a flood hazard ranking map. Subsequently, a short description of the mentioned driving factors was carried out. Entropy of a drainage network: The drainage network and its features, such as density, length, and bifurcation, are the most important geomorphological features of sub-watersheds that influence soil erosion and flood peak through direct and indirect effects. The direct effects refer to the relationship between flood features (e.g., velocity and depth) and the drainage network, while indirect terms are used to describe the role of the drainage network as a geological index. The entropy of a drainage network is defined as the state of the art of methodology for describing drainage networks and their complexity, so a higher value of entropy tells us that the desired drainage network has more effect on the flooding degree than a drainage network with low entropy. To calculate the entropy of the drainage network, we used the box-counting method (for more descriptions, refer to Sepehri et al. [28]).
Slope: The slope criterion, as a hydrologic-related morphometric factor, can be used to identify areas susceptible to flooding in low-slope gradients. Low-slope areas behave as ponds that retain surface runoff as temporary storage [11,29,30].
Topographic control index (TCI): Topographic elements, such as slope, contributing areas, and the volumes of depression, are also key factors that significantly influence the flow direction, velocity of runoff, and the potential places where pluvial flooding can occur. For a single depression, the quicker it becomes inundated, the more prone it is to flooding. The time required to fill a depression is affected by the topography of its watershed, such as the watershed area, the watershed slope, and the ponding volume of the depression. A watershed includes the upslope contributing area and the depression itself. To evaluate the flooding risk for each depression, an integrated topography control index, namely TCI, was developed [31].
Compactness Coefficient (Cc): Cc refers to the shape of a sub-watershed, which is calculated based on the ratio of the perimeter of the basin to the circumference of a circular area, which equals the basin area. It is one of the most commonly used morphometric features of watersheds in flooding studies. When the values of Cc tend to be close to one, it means the shape of the watershed is closer to a circle, and, therefore, it has the lowest infiltration capacity and the highest sensitivity to flooding [32,33].
Stream power index (SPI): SPI is one of the most commonly used indicators in drainage network analyses, and it can be used to detect and quantify the power of a flow at a given point on a topographic surface [34][35][36][37].

Analytic Hierarchy Process (AHP)
In 1980, Saaty [38] presented a multilevel and hierarchical structural technique that allows one to solve complex and subjective problems. In the technique, the weights of the criteria and sub-criteria are driven by pairwise comparisons. A matrix concerning the nth decision criteria or sub-criteria is created as follows: where shows the relative importance of the criteria/sub-criteria to the criteria/sub-criteria j. The level of importance is calculated based on the experts' knowledge through a 9point rating scale, which varies between 1/9 (lowest importance) and 9 (highest importance) ( Table 1). In the following stage, the elements of the A matrix must be synthesized and normalized to obtain the rank criteria/sub-criteria as matrix: Since the weights of the criteria/sub-criteria are calculated based on experts' knowledge, it is necessary to evaluate the consistency ratio (CR) of matrix A regarding its significant degree of uncertainty. If the CR is < 0.1, it can be assumed that the mentioned matrix has enough consistency. The rate of consistency can be expressed as follows: where RI (random index) is a random index that is determined based on randomly generated 500 matrixes [36], and the consistency index (CI) is calculated based on the maximum eigenvector ( ) of the A matrix and (number of criteria/sub-criteria), given as follows: In a pure AHP, each individual decision regarding the weights of criteria and subcriteria is expressed in rough values. In this state, decision-makers often face a dilemma in allocating relative importance to selected criteria and sub-criteria, and consequently, subjectivity, imprecision, and uncertainty rise in the process of decision-making. In recent years, a new concept in the IRN-based rough numbers theory has been introduced to deal with uncertainty and imprecision. In this regard, the combination of IRN and pure AHP, known as IRAAHP, means that scholars can better deal with the uncertainties of rough numbers. The IRAAHP algorithm is described as follows: The algorithm of IRNAHP

IRN Mathematical Model
Suppose that there is a set of decision-maker (DM) preferences as = ( , , … , ), and all its objects are defined in a universe and represent DM preferences.
Each object in an R set is defined by the interval = , , which is subject to ≤ (1 ≤ ≤ ). and denote the lower and upper limits of the class, respectively. If and satisfy < <, … , < and < <, … , < (1 ≤ ≤ ) , respectively, then = ( , , … , ) and = ( , , … , ) can be defined. In this state, the lower approximations ( ) of and are represented as follows: The upper approximations of and are expressed as follows: The lower limits of and and, consequently, their upper limits are expressed as follows: The upper limits are expressed as follows: where and represent the number of objects in the lower/upper approximations of and , respectively. The rough boundary (RB) for the lower/upper approximations of and can be defined as follows: Rough boundary : Rough boundary : Then, the rough boundary concatenation (RC) classes of and are expressed as follows: Finally, the interval rough number (IRN) of the object can be concluded as follows: Interval Rough Numbers AHP (IRNAHP) Mathematical Model In order to integrate the IRN and pure AHP, the elements of the matrix are transferred to the interval rough number , and the next stages of the pure AHP are repeated. Therefore, three to five stages are needed to obtain the IRNAHP mathematical model, as follows: Stage 1: Establishing the pairwise comparison matrixes by experts as follows: where and are the relative importance of the criteria/sub-criteria i to the criteria/sub-criteria j, which are selected based on Saaty's 9-point rating. If every expert has vagueness for selecting 2 values between the Saaty's 9-point rating, the ≠ is established. If there is no vagueness in selecting, the expert chooses one value, and, in this state, we have the following: = .
Stage 2: Calculating the consistency rate for every expert. This stage is similar to the AHP pure one, with the difference that there are two consistency rates, one of which is for upper approximations ( ) and the other one is for lower approximations ( ). The final CR is calculated based on ( + )/2.
Then, using and , the lower and upper limits of , = , , are obtained. The matrix can be shown as follows: Finally, the interval rough weighted coefficient, ( ) , can be calculated using Equation (26).
The best-worst method (BWM), developed by Rezaei [39], is one of the latest MCDA methods for dealing with multi-objective complex systems. Compared to AHP, BWM requires fewer pairwise comparisons to obtain the weights of criteria and sub-criteria, so it only requires pairwise comparisons, whereas in AHP, these comparisons are increased until one of the main sources of rising uncertainty is reduced.
The algorithm of the best-worst method (BWM) can be described as follows: Stage 1. Selection of a set of desired criteria/sub-criteria, and determination of the best/worst (most/least important) of them.
Stage 2. Determination of preferences of the best criterion over others (BO) and vice versa, i.e., preferences of all criteria over the worst criterion, using a 9-point rating scale (Table 1).
where indicates the preferences of the best criterion and shows the preferences of the worst criterion.
Stage 3. Finding the optimal weights ( * , * , … , * ) by solving the following model: Stage 4. Checking the consistency as follows: The consistency index is calculated and shown in Table 2. The consistency ratio varies between 0 and 1, where the lower values of the consistency ratio show a more consistent preference matrix. Gündodu et al. [27] introduced the picture fuzzy analytic hierarchy process and picture fuzzy linear assignment, a novel technique of hybrid multi-attribute decision-making using AHP and LAM under the picture fuzzy concept to reduce the hesitancy or uncertainty of decision-makers' judgments due to the lack of information or motivation on the examined problem. On Saty's 9-point rating system, the picture fuzzy is the most recent expansion of fuzzy logic and can be used as a valuable tool to deal with a decision maker's imprecision and uncertainty. (See the following for more information on the picture fuzzy analytic hierarchy process and picture fuzzy linear assignment algorithm.) Definition 1: The single-valued PFS of ÃU of the universe of discourse U can be defined as follows: where and the indeterminacy of to ̃ , respectively. Additionally, the degree of refusal of can be defined as follows:

Definition 2:
The basic operations of PFS can be defined as follows: • Multiplication by a scalar; λ > 0 Definition 3: The PFS-weighted geometric (PFWG) mean given that = , , … , , ∈ 0,1 , ∑ = 1 can be defined as one of the following equations: , … , Definition 4: To de-fuzzify, rank, and compare the PFS sets, the following score (SC) and accuracy (AC) functions can be used: After calculating the score and accuracy functions, the dominance rules can be written as follows: The stages listed below describe the picture blurry, which is related to LAM. The weights of criteria using picture fuzzy are determined in Stage 1, while the rank of alternatives is calculated in Stage 2.
Stage 1: Picture fuzzy analytic hierarchy process. 1.1. Using pairwise comparison matrices for the weights of the criteria The criteria are compared in this step by the decision-makers. They chose a value from Saaty's nine-point scale based on their job, and the value was then translated to the associated picture fuzzy digits (Table 3). Since the values chosen from Table 1 are based on experts' preferences, it is important to double-check the consistency ratio (CR) of each pairwise comparison, which is calculated using Equation (2), as well as the consistency index (CI) and random index (RI).

1.2.
To aggregate the decision-makers' assessments, we used a weighted geometric (PFWG) mean. There can be different comparison matrixes ̃ in decision-making situations because there are multiple decision-makers. It is necessary to employ geometric means (Equation (38)) to unify all comparison matrices ̃ in the next steps. Eventually, the final picture fuzzy weight ̃ must be calculated as follows:
It is required to export ̃ as a defuzzified value in this stage in order to use ̃ from the criteria as a weight for ranking the sub-watersheds (Stage 2).  (37), as shown in Table 5. 2.3. Defuzzification of the aggregated matrix using Equation (40) was performed to compare and rank options that are connected to each other.
2.4. Determination of the rank frequency matrix, which includes associated elements that show the number of times alternative m dominates on the nth criterion ( Table 6).
Determination of the weighted rank frequency matrix , which measures the contribution of the mth alternative to the overall ranking (Equation (45) and Table 7).
Construction of the linear assignment model based on and permutation matrix P (m*m) as follows: 2.7. Using Equation (46) to obtain the optimal permutation matrix ( * ).

Morphometric Parameters
Morphometric parameters of drainage basins have been successfully applied to simulate Earth's surface processes and landforms, incorporating hydrological, geological, and geo-morphological setups at different scales [23,[40][41][42][43]. It can be observed that the entropy of the drainage network ranges between 0.03 (sub-watershed #6) and 0.145 (sub-watershed #13). According to the results of a TCI, the highest value of this criterion is related to sub-watershed #12 (−0.81) and vice versa, and sub-watershed #15 acquired the lowest value (−1.29). The highest and lowest values of the slope criterion are related to sub-watershed #15 (40.65) and sub-watershed #3 (12.22), respectively. The values of SPI indicate that sub-watershed #15 (−1.47) and sub-watershed #4 (−3.73) have been ranked in the first and last positions, respectively. According to the results of the Cc, it can be concluded that sub-watershed #5 (2.78) gained the highest values and sub-watershed #18 (1.58) received the lowest values, respectively (Table 8). Ranking of sub-watersheds using the AHP technique During the analysis, the weights of the morphometric parameters are calculated based on their importance in the case study of area floods. The weight assignment for each criterion was calculated based on the local characteristics of each criterion, and the opinions of three experts in the field of hydrology science are shown in Table 9. Based on the criteria weights, the most important criterion regarding the occurrence of floods in the study area is related to the entropy of the drainage network. The Cc weight was defined as the least important criterion. The remaining criteria, i.e., SPI, TCI, and slope, are next in order of importance. After the weight assignment to each criterion, the total weight for each sub-watershed was computed using a simple weighted sum, which is as follows: where is the weight of the criterion jth and is the normalized value of the criterion jth. The weights of each criterion are given in Table 9. Additionally, for each pairwise comparison matrix, the calculated consistency ratio (CR) showed that the mentioned matrices have enough consistency.
In the BWM, based on the judgments of the decision-makers, the entropy of the drainage network was chosen as the best criterion (the most important flood-related criterion), and, at the opposite point, the Cc was selected as the worst criterion (the least important flood-related criterion). Then, using the best/worst criterion, the BO and OW vectors were determined, and, finally, the average weight of the criteria was calculated as follows: entropy of drainage network (0.462), SPI (0.203), TCI (0.133), slope (0.136), and Cc (0.067) ( Table 9). Finally, using Equation (1), the ranking of the sub-watersheds was acquired. On the other hand, in the IRNAHP, the entropy of the drainage network and the Cc were known as the most and least important flood-related criteria (Table 9).
Based on Table 10, the local and global weights of each criterion were calculated (Table 10), and the weights of the alternatives (sub-watersheds) are shown in Table 11. Consequently, by using the mentioned weights, the weighted frequency rank matrix was determined (Table 12). Finally, according to the weighted frequency rank matrix, the permutation matrix in Table 13 was acquired to rank the sub-watersheds based on the flood risk degree. Finally, by using Equation (42), the ranking of the sub-watersheds was acquired.   1th 2th 3th 4th 5th 6th 7th  17th 18th 19th 20th 21th 22th 23th 24th

Discussion
Based on the results of the AHP technique, the sub-watersheds #22, #1, and #6 have very high flood probabilities, while the sub-watersheds #23, #11, and #14 have very low flood probabilities (Figure 3). Based on Figure 3, it can be observed that the sub-watersheds #22, #1, and #6 receive the lowest value of entropy of the drainage network, and the SPI criteria in these sub-watersheds are nearly high. These results are consistent with the findings of Sepehri et al. [28] and Zhang et al. [44] in that they emphasize the drainage network as the most important geomorphology feature of a watershed that has the main effect on flooding and sediment yield, so that a watershed that has a drainage network with more complex features (such as more entropy/length/stream power) has more susceptibility to flooding and sediment yield than a watershed with a low-complexity drainage network. The flood observation points, which were prepared by the General Department of Natural Resources of Hamadan Province and are shown as red circles in Figure  3, represent that the accuracy of the AHP output regarding flood ranking is acceptable, so that the high number of flood observation points are located in sub-watersheds #22, #1, and #6. In the BWM method, the results show the same result as with the AHP technique, so that sub-watersheds #22, #1, and #6 are positioned at the first rank of flooding degree, while in contrast, sub-watersheds #23, #14, and #11 have the last rank of flooding degree. However, these same results show that the BWM method is less time-consuming than the AHP technique because, in this study, to reach the weights of the used criteria, the BWM method only needed 21 pairwise comparison matrices (for three experts), whereas in the AHP technique, this number of pairwise comparisons is significantly increased, i.e., 60 pairwise comparisons. It is obvious that in other studies with a higher number of criteria, the difference between the pairwise comparisons will be dramatically increased, and, in this state, the uncertainty of the considerations will be increased [19,25,26]. In the IRNAHP method, the sub-watersheds #22, #1, and #6 with the highest elevation are located in the first three ranks, respectively, and, in contrast, sub-watersheds #20, #15, and #11 are specified in the last rank and are the most susceptible to flooding, respectively. Even though these same results through the IRNAHP can be acquired by solving complex algorithms rather than the AHP technique, it must be noted that in MCDA studies that involve a high number of criteria and are associated with a large amount of uncertainty and subjectivity, the experts face a dilemma in choosing a crisp value as the initial weight for each criterion. Therefore, in these studies, IRN can be used as an effective tool to exploit uncertainty [26,45]. The output of PICALAM shows that the sub-watersheds #6, #22, and #1 are posited in ranks 1, 2, and 3, whereas the sub-watersheds #4, #15, and #13 are located in the last ranks. The first three ranks are similar to the previous methods, but the last three ranks are quite different. According to Figure 3, it can be observed that in sub-watersheds #4, #15, and #13, despite having the highest value of entropy of the drainage network and the lowest value of SPI, other flood-related criteria have the lowest relationship with the flooding degree. Additionally, in the PICALAM, the rank of the middle sub-watersheds is not the same as the three remaining methods. For example, the sub-watersheds #16 and #9 in PICALAM have been ranked as the most susceptible areas to flooding (ranks 4 and 5). The flood observation point in sub-watersheds #16 and #9 proves that these sub-watersheds have been ranked correctly. Therefore, it can be concluded that the PICALAM method has better performance than other methods. In PICALAM, there are two main advantages over other methods that can enhance its performance. (1) In the PICLAM, the values of the nine-point rating scale will be transferred to their corresponding fuzzy numbers so that these numbers are based on human reasoning and reduce uncertainty [27,46]. (2) One of the main characteristics of watersheds is their spatial variability, which means that the role and importance of each criterion subjected to the desired objective (here, flooding degree) will differ from one sub-watershed to the next. On the other hand, each criterion has an interrelationship with other criteria that cannot be considered due to a lack of data or an increase in the number of used criteria. For example, one of the main criteria used in this study is TCI, which represents the capacity of each sub-watershed to retain water and is calculated based on surface features (i.e., slope, upslope contribution area, and water volume). It is obvious that these parameters also depend on soil conductivity and land use, which were not considered in this study. In PICALAM, experts can use global and local weights to assign weights to commonly used criteria. For global weights, on which the AHP, IRNAHP, and BWM methods are based, the weight of each criterion is calculated as general for all sub-watersheds without taking into account the properties of watershed spatial variability. For local weights, the weight of each criterion can be calculated separately without considering other sub-watersheds. Therefore, PICALAM allows for greater flexibility and accuracy in decision-making due to the uncertainty and imprecision of complex decision analyses.

Conclusions
In this study, a range of MCDM approaches, such as the AHP, BWM, IRNAHP, and the PICAHP and LAM, were used to rank Sad-Kalan sub-watersheds in terms of flooding degree. To accomplish this, five essential criteria related to flooding, namely entropy of the drainage network, SPI, slope, TCI, and Cc, were chosen. The mentioned MCDM approaches were employed to assign a specific weight to each criterion based on its significance in determining the degree of flooding. Subsequently, by combining the weighted criteria, a flood hazard map was generated. The results indicated that among the MCDM approaches utilized, the PICAHP and LAM method demonstrated superior performance due to its consideration of the spatial variability of watersheds. Additionally, the results revealed that the AHP, BWM, and IRNAHP methods exhibited comparable performance in ranking the sub-watersheds. However, it should be noted that the BWM method stands out for its time-efficient operation in comparison to the other two approaches. Overall, it is important to acknowledge that this study was conducted to assess how human linguistic terms and their quantity influence the uncertainty of flood degree ranking in a case study with only five flood-related criteria. It is worth noting that the performance of the MCDM approaches used may vary significantly in other studies that involve a higher number of criteria and greater complexity. Furthermore, it is essential to recognize that flood hazards, as the initial component of flood risk studies, primarily focus on the physical and climatic aspects of floods. For future studies, it is advisable to incorporate socialeconomic criteria to facilitate comprehensive flood management.

Conflicts of Interest:
This manuscript has not been published or presented elsewhere in part or entirety and is not under consideration by another journal. There are no conflicts of interest to declare.