Flood Hazard Mapping Using Fuzzy Logic, Analytical Hierarchy Process, and Multi-Source Geospatial Datasets

: Iran is among the driest countries in the world, where many natural hazards, such as ﬂoods, frequently occur. This study introduces a straightforward ﬂood hazard assessment approach using remote sensing datasets and Geographic Information Systems (GIS) environment in an area located in the western part of Iran. Multiple GIS and remote sensing datasets, including Digital Elevation Model (DEM), slope, rainfall, distance from the main rivers, Topographic Wetness Index (TWI), Land Use/Land Cover (LULC) maps, soil type map, Normalized Difference Vegetation Index (NDVI), and erosion rate were initially produced. Then, all datasets were converted into fuzzy values using a linear fuzzy membership function. Subsequently, the Analytical Hierarchy Process (AHP) technique was applied to determine the weight of each dataset, and the relevant weight values were then multiplied to fuzzy values. Finally, all the processed parameters were integrated using a fuzzy analysis to produce the ﬂood hazard map with ﬁve classes of susceptible zones. The bi-temporal Sentinel-1 Synthetic Aperture Radar (SAR) images, acquired before and on the day of the ﬂood event, were used to evaluate the accuracy of the produced ﬂood hazard map. The results indicated that 95.16% of the actual ﬂooded areas were classiﬁed as very high and high ﬂood hazard classes, demonstrating the high potential of this approach for ﬂood hazard mapping.


Introduction
Natural hazards frequently cause a variety of physical (e.g., injuries, casualties, and property damages) and non-physical (e.g., psychological and mental) disturbances worldwide [1]. The number of disastrous events has dramatically increased in recent decades, and consequently, the number of casualties/injuries, economic losses, and intense environmental disturbances [2,3]. Flooding, as the most frequent natural hazard, is not an exception, and over the last few decades, it has caused considerable negative environmental and socioeconomic damages in different parts of the world [4][5][6][7][8]. An integrated flood risk management approach with a focus on reducing the vulnerability of the societal system is required to reduce future flood risks [7].
Global climate change and unsuitable urbanization are considered as the significant factors of changing flood characteristics [9]. The frequency, scale, and extent of flood events are increasing worldwide [10]. For instance, flood events negatively affected more than 232,896 people and caused over 295 million USD damages in Morocco between 1995 and 2005 [11]. Likewise, the Department of Irrigation and Drainage (DID) in Malaysia reported that around 29,000 km 2 (8.79% of the area of Malaysia) land area and over 4.82 Remote Sens. 2021, 13, 4761 3 of 22 rainfall intensity, soil erodibility, flow accumulation, LULC, elevation, slope, and distance from drainage networks were combined.
Iran has been negatively affected by climate change over the past three decades. Prolonged droughts have changed many ecosystems in Iran, and significant parts of forests have disappeared [25]. Many regions in Iran are prone to floods due to their topographical and geographical characteristics. Moreover, a large number of rivers in Iran are susceptible to flash flooding during heavy rainfall [10]. For example, during spring 2019, almost all provinces in Iran experienced heavy precipitations, resulting in flood events in several parts of Iran [26]. To address flood-related issues at the basin scale, many studies have investigated the potential of remote sensing technology for flood mapping, monitoring, and risk assessment. For example, Shahabi et al. [27] developed a method to predict flood hazard zones using the Evidential Belief Function (EBF) model. They employed multiple remote sensing and GIS datasets, such as DEM, slope, curvature, Topographic Wetness Index (TWI), Stream Power Index (SPI), distance from rivers, rainfall, soil type, LULC map, and the Normalized Difference Vegetation Index (NDVI) to produce flood susceptibility maps in the Haraz Catchment in Mazandaran Province, Iran. Darand et al. [28] also employed daily precipitation and the Kriging method for flood monitoring in Iran for the period of 1962-2013. They found significant changes in precipitation trends especially near the Zagros Mountains, and therefore, more than 40% of Iran was determined as flood-prone zones. Vaghefi et al. [29] used five high-resolution climate models to estimate temperature, rainfall, and future floods in Iran between 1980 and 2004. They reported that the climate of Iran prolonged dry periods with irregular heavy precipitations, increasing the likelihood of floods between 2025 and 2049.
In the last few decades, the Iranian government has spent approximately 29 billion US dollars to mitigate various natural hazards. Drought, earthquakes, and floods have resulted in more government expenditures than other natural hazards, with more than 14, 6, 9, and 6.1 billion dollars incurred, respectively [30]. As discussed, a large portion of Iran is prone to flooding, and thus, it is necessary to develop an advanced system for flood hazard mapping. These flood hazard maps could effectively be employed for resilience actions and are considered a valuable criterion in urban growth planning [31]. This work aims to produce a flood hazard map of the Pol-e Dokhtar watershed located in the western part of Iran, which faced a massive flood in 2019. In this regard, open-access remote sensing and GIS datasets were integrated to produce an accurate flood hazard map that reduces the cost and increases the applicability of the implemented approach for other regions. To this end, nine parameters, including rainfall (i.e., satellite and in situ data), Digital Elevation Model (DEM), slope, distance from the main rivers, Topographic Wetness Index (TWI), LULC, soil type, Normalized Difference Vegetation Index (NDVI), and erosion rate within a fuzzy logic system were employed. Finally, the proposed approach was evaluated using Sentinel-1 SAR images captured over the study area before and on the day of the flooding event.

Study Area and Case Study
The study area is the Pol-e Dokhtar watershed located in Lorestan province, the western part of Iran, with an area of approximately 3738 km 2 and a total population of over 73,744 based on the 2015 census [32]. Pol-e Dokhtar is located between the latitude and longitude ranges of 32 • 28 and 33 • 28 N and 47 • 27 and 48 • 25 E, respectively, with an average elevation of about 700 m above the sea level (see Figure 1). The study area has a semi-arid climate with a noticeable temperature change within a year, ranging from −4 • C (around January) to 48 • C (in August and July). The annual rainfall value is approximately 358 mm. There are two main rivers in this area (i.e., Kashkan and Seymare). Moreover, in the vicinity of these rivers there are many illegal constructions, some of which have been destroyed by floods during the last two decades.
River is 30 cubic meters [13,33]. The scale of the economic damage of the spring flood in Pol-e Dokhtar was massive. The value of damages to production facilities, buildings and agricultural assets was estimated to be approximately 26 million USD [10]. About 57 bridges and culverts and 250 km of major roads were damaged in this province. Besides the financial damage in the province of Lorestan, over 46,000 families have been negatively affected and many people had to leave their properties and all their assets behind [10].

Flooded Areas and Rainfall History
Heavy rainfalls are one of the main reasons for flooding. Flooding is most frequently caused by excessive rainfall when natural streams are unable to control additional water. Therefore, one of the most important parameters affecting the flood is the amount of precipitation. The comparison between the Global Precipitation Measurement (GPM) rainfall data and the previous flood events in Pol-e Dokhtar (see Figure 2), showed that high precipitation rates caused several major floods in this area. For example, the highest amount of rainfall was observed in March 2019, causing one of the largest floods in this region. It is worth noting that we used an averaging method to combine the GPM precipitation and local rainfall datasets to improve the accuracy of the final rainfall map. Based on the provided reports of the United Nations Country Team (UNCT) in Iran, successive heavy rainfalls over the Pol-e Dokhtar watershed between 25 March and 1 April 2019, caused a devastating flood along the Kashkan River and its tributaries with a water flow rate of over 5000 cubic meters. The long-term average discharge of the Kashkan River is 30 cubic meters [13,33]. The scale of the economic damage of the spring flood in Pol-e Dokhtar was massive. The value of damages to production facilities, buildings and agricultural assets was estimated to be approximately 26 million USD [10]. About 57 bridges and culverts and 250 km of major roads were damaged in this province. Besides the financial damage in the province of Lorestan, over 46,000 families have been negatively affected and many people had to leave their properties and all their assets behind [10].

Flooded Areas and Rainfall History
Heavy rainfalls are one of the main reasons for flooding. Flooding is most frequently caused by excessive rainfall when natural streams are unable to control additional water. Therefore, one of the most important parameters affecting the flood is the amount of precipitation. The comparison between the Global Precipitation Measurement (GPM) rainfall data and the previous flood events in Pol-e Dokhtar (see Figure 2), showed that high precipitation rates caused several major floods in this area. For example, the highest amount of rainfall was observed in March 2019, causing one of the largest floods in this region. It is worth noting that we used an averaging method to combine the GPM precipitation and local rainfall datasets to improve the accuracy of the final rainfall map.

Datasets
In this study, nine parameters were used to produce a flood hazard map for the study area. These parameters include DEM, slope, rainfall, distance from the main river, TWI, LULC map, soil type map, NDVI, and erosion rate. These parameters are explained and illustrated in Table 1 and Figure 3. More details of each parameter are also briefly described in the following subsections. Additionally, bitemporal Sentinel-1 SAR images, acquired before and on the day of the flood event, were employed to assess the applicability of the produced flood hazard map.

Datasets
In this study, nine parameters were used to produce a flood hazard map for the study area. These parameters include DEM, slope, rainfall, distance from the main river, TWI, LULC map, soil type map, NDVI, and erosion rate. These parameters are explained and illustrated in Table 1 and Figure 3. More details of each parameter are also briefly described in the following subsections. Additionally, bitemporal Sentinel-1 SAR images, acquired before and on the day of the flood event, were employed to assess the applicability of the produced flood hazard map. Sentinel-2 was not used directly; only the products of the previous article were used [34]. Sentinel-2 10 m × 10 m Copernicus, the European Commission's (EC) Earth Observation Program Sentinel-2 was not used directly; only the products of the previous article were used [34].

Digital Elevation Model (DEM)
The elevation is among the most important parameters in flood risk studies. Areas with lower elevations have higher possibilities of being affected by rainfalls and are more flood-prone than regions with higher elevation [35,36]. This study used the Shuttle Radar Topography Mission (SRTM) DEM data with a spatial resolution of ~30 m (see Figure 3a). These data were published in September 2014 and are currently the best open-access elevation data over our study area [37].

Digital Elevation Model (DEM)
The elevation is among the most important parameters in flood risk studies. Areas with lower elevations have higher possibilities of being affected by rainfalls and are more flood-prone than regions with higher elevation [35,36]. This study used the Shuttle Radar Topography Mission (SRTM) DEM data with a spatial resolution of~30 m (see Figure 3a). These data were published in September 2014 and are currently the best open-access elevation data over our study area [37].

Slope
Low slopes increase the quantity of water soaked into the soil and, thus, cause a higher risk of flooding than higher slopes [38]. This is because more water is drained by higher slopes, which are generally located at high elevations, and much water gathers in areas with lower slopes. Various studies have adopted different categories for the slope map. However, since Pol-e Dokhtar is a mountainous area (Zagros Mountain), we categorized Remote Sens. 2021, 13, 4761 7 of 22 the slope map based on its location and topographic features. The slope map was generated using the (SRTM) DEM data and GIS tools in this study (see Figure 3b).

Rainfall
Rainfall is another important parameter that has been extensively employed in flood hazards and susceptibility mapping [39,40]. In this study, the Global Precipitation Measurement (GPM) precipitation data, along with in situ rainfall data from 15 synoptic stations, were employed to produce a highly accurate rainfall map. The GPM precipitation provides rainfall data with a spatial resolution of 0.1 • × 0.1 • . GPM usually provides higher accuracy than other precipitation datasets, such as those provided by the Tropical Rainfall Measuring Mission (TRMM; i.e., 0.25 • × 0.25 • ) [41]. An averaging method was used in this study to combine local rainfall data generated from 15 weather stations with GPM precipitation. The GPM is an international mission that integrates precipitation measurement from different sensors to provide sub-daily rainfall maps. Satellite precipitation estimates are less accurate, but cover a larger geographical area. On the other hand, the in-situ measurements are very accurate [42] with a sparse distribution over specific locations. To resolve each of these limitations, we combined these two datasets to produce an accurate rainfall map over the entire study area [43]. To this end, the in situ data, obtained from the meteorological organization of the province of Lorestan (see Table 2), were first used to produce a rainfall raster dataset using the Kriging interpolation approach. Kriging is an interpolation technique utilizing an estimator that is non-biased and has a linear minimum variance. A kriging method minimizes variance, which leads to the smoothing of the field, especially near synoptic stations [44]. Using the Kriging method is often beneficial, not only for estimating uncertainty but also for improving forecasts [45]. Compared to other interpolation methods, the Ordinary Kriging (ORK) was considered to provide the lowest RMSE [46]. Therefore, in this study, we interpolated in situ rainfall data from 15 synoptic stations using the ORK method. Finally, rainfall maps (i.e., GMP rainfall map and the map produced through in situ data and ORK method) were spatially overlaid and averaged to generate the final rainfall map (see Figure 3c).

Distance from the Main Rivers
The distance from the rivers significantly affects flood distribution. In fact, flood extent and intensity are more severe in regions located closer to the main rivers [47]. In this study, stream networks were delineated from DEM using the Flow Accumulation tool in ArcMap [48,49]. Finally, a raster with a spatial resolution of 30 m was produced using the multi-ring buffer analysis by assigning the distance from the main rivers to each pixel (see Figure 3d).

Topographic Wetness Index (TWI)
The TWI index indicates how slope affects hydrological processes. TWI describes the water accumulation trend at a particular location, and the local slope shows the influence of gravitational forces on the water flow [50,51]. TWI is commonly used to forecast the amount of soil moisture and describes the tendency of an area to accumulate water [52]. In this study, the slope and DEM datasets were applied to calculate TWI (see Figure 3e) using Equation (1).
where AS is the specific catchment area, and β is the local slope angle in degrees. The specified catchment area is an integration of surface and subsurface drainage from an upslope region per unit contour width [53,54].

Land Use/Land Cover (LULC)
Different types of land cover have different infiltration rates and debris flow, and thus, the LULC map plays an essential role in flood hazard mapping [55]. For instance, forests and farms have a higher infiltration capacity, while urban areas support the overland flow of water. In this study, we used the Iran-wide LULC map produced by Ghorbanian et al. [34]. This map was generated by processing time-series Sentinel-1 and Sentinel-2 images within the Google Earth Engine (GEE) cloud computing platform. This LULC map includes 13 LC classes of Urban, Water, Wetland, Kalut, Marshland, Salty land, Clay, Forest, Outcrop, Uncovered plain, Sand, Farmland, and Rangeland. In this study, the LULC map was reclassified using the criteria demonstrated in Table 3. For example, Urban areas have a very high risk for flooding, however Water and Marshland have very low risk [17]. Finally, this reclassified map ( Figure 4) was used as one of the inputs for flood hazard mapping.

Soil Type
One of the most important aspects of soil is its texture, which has a considerable impact on flooding [38]. In this study, soil types were extracted from the Geological Survey and Mineral Explorations of Iran. The study area mainly includes Rock Outcrops/Inceptisols, Bad Lands, Rock Outcrops/Entisols, Inceptisols/Vertisols and Inceptisols (see Figure 3g). The soil type parameter allowed determining the amount of water penetration in soil. For example, soil made from Entisols increases water retention capacity even more than the roots of plants [56]. Therefore, there is a very low risk of flooding in soils containing Entisol. However, the low water infiltration rate is one of the physiochemical characteristics of Badlands soil. Therefore, most rainfall in areas with Badland soil drains into surface runoff. Additionally, it has been reported that a concentration of sodium and chlorine near the surface of badland soils occurs during the dry season, which creates a repulsive force between the soil particles, causing the soil to become extremely hard. Badland ecosystems are also characterized by a high sodium and chloride concentration, which increases soil electric conductivity and makes it unsuitable for plant growth, resulting in bare landscapes. In areas with this type of soil, a short period of rainfall may cause severe flooding due to Remote Sens. 2021, 13, 4761 9 of 22 the inability of the soil to penetrate water and the lack of plants that retain water [57]. The soil type map was reclassified using the criteria reported in Table 4. Finally, this reclassified map ( Figure 5) was used as one of the inputs for flood hazard mapping.

Soil Type
One of the most important aspects of soil is its texture, which has a considerable impact on flooding [38]. In this study, soil types were extracted from the Geological Survey and Mineral Explorations of Iran. The study area mainly includes Rock Outcrops/Inceptisols, Bad Lands, Rock Outcrops/Entisols, Inceptisols/Vertisols and Inceptisols (see Figure  3g). The soil type parameter allowed determining the amount of water penetration in soil. For example, soil made from Entisols increases water retention capacity even more than the roots of plants [56]. Therefore, there is a very low risk of flooding in soils containing Entisol. However, the low water infiltration rate is one of the physiochemical characteristics of Badlands soil. Therefore, most rainfall in areas with Badland soil drains into surface runoff. Additionally, it has been reported that a concentration of sodium and chlorine near the surface of badland soils occurs during the dry season, which creates a repulsive force between the soil particles, causing the soil to become extremely hard. Badland ecosystems are also characterized by a high sodium and chloride concentration, which increases soil electric conductivity and makes it unsuitable for plant growth, resulting in bare landscapes. In areas with this type of soil, a short period of rainfall may cause severe flooding due to the inability of the soil to penetrate water and the lack of plants that retain water [57]. The soil type map was reclassified using the criteria reported in Table 4. Finally, this reclassified map ( Figure 5) was used as one of the inputs for flood hazard mapping.

Normalized Difference Vegetation Index (NDVI)
Vegetation cover on any surface reduces the water movement velocity and increases the water infiltration of the corresponding surface [58,59]. The conversion of vegetation cover to urban or barren areas increases the flooding risk since vegetation cover can mitigate flood intensity and velocity [29]. Since vegetation can influence local hydrodynamics, it can affect flood risk. In this regard, remote sensing data can be employed to provide new insights to support the management of high waterways [60]. For example, NDVI is a spectral index indicating the presence and vigor of vegetation cover. It can be calculated using the red and near-infrared bands of satellite images. In this study, NDVI values were computed using Landsat-8 satellite images with a 30 m spatial resolution and Equation (2) (see Figure 3h). We should note that the Landsat 8 image used in this study was captured on 7 March, which was very close to the date of the flood.

Normalized Difference Vegetation Index (NDVI)
Vegetation cover on any surface reduces the water movement velocity and increases the water infiltration of the corresponding surface [58,59]. The conversion of vegetation cover to urban or barren areas increases the flooding risk since vegetation cover can mitigate flood intensity and velocity [29]. Since vegetation can influence local hydrodynamics, it can affect flood risk. In this regard, remote sensing data can be employed to provide new insights to support the management of high waterways [60]. For example, NDVI is a spectral index indicating the presence and vigor of vegetation cover. It can be calculated using the red and near-infrared bands of satellite images. In this study, NDVI values were computed using Landsat-8 satellite images with a 30 m spatial resolution and Equation (2) (see Figure 3h). We should note that the Landsat 8 image used in this study was captured on 7 March, which was very close to the date of the flood.
where and are the surface reflectance values of the near-infrared and red bands, respectively. The NDVI map was used to categorize flooding hazards based on the numerical values. The classes with values more than 0.40 and less than −0.11 were considered less vulnerable and more susceptible to flooding, respectively.

Erosion Rate
Erosion rate is a crucial parameter that affects the possibility of flooding during heavy rainfalls. This is because the dissolution of soil in water affects the water flowing on the Earth's surface and significantly increases the water volume runoff in areas with a high erosion rate [61]. In this study, the erosion rate (Figure 3i) was determined using the where ρ N IR and ρ Red are the surface reflectance values of the near-infrared and red bands, respectively. The NDVI map was used to categorize flooding hazards based on the numerical values. The classes with values more than 0.40 and less than −0.11 were considered less vulnerable and more susceptible to flooding, respectively.

Erosion Rate
Erosion rate is a crucial parameter that affects the possibility of flooding during heavy rainfalls. This is because the dissolution of soil in water affects the water flowing on the Earth's surface and significantly increases the water volume runoff in areas with a high erosion rate [61]. In this study, the erosion rate (Figure 3i) was determined using the Universal Soil Loss Equation (USLE) formula (Equation (3)) to measure annual soil loss associated with sheet and flow erosions [32].
where A is the maximum amount of land loss (ton/year), R is the monthly average rainfall erosivity, K is soil erodibility index factor (ton/pear), L is the length factor of topography, S is the slope factor of the topography, C is crop management index factor, and P is soil conservation index factor. Soil erosion has a key role in floods. Most articles have not directly used erosion rate as a parameter for generating flood hazard maps because soil type and NDVI can be used to consider the impacts of soil erosion in flooding. However, in this study, we have directly used the soil erosion parameter to obtain more accuracy in determining flood hazard locations soils with a high potential for erosion can increase flood potential because soluble soil contributes to flooding. Besides increasing the volume of water, erosion is assumed to cause the transport of silt together with river water. A considerable amount of silt is deposited at the riverbed and the banks, resulting in a rise in the bed level and congestion of the water flow [62]. It is worth noting that erosion rates were provided by the Soil and Water Research Institute in Tehran, Iran.

Sentinel-1 Images
In this study, two Sentinel-1 Ground Range Detected (GRD) images, captured on 2 March 2019 (before the flood event) and 1 April 2019 (the date of flooding) were used to evaluate the accuracy of the produced flood hazard map. Several preprocessing steps were applied to each image to prepare them for flood detection. For example, orbit file correction was first applied to both images, followed by a thermal noise removal step. Subsequently, both images were converted to gamma-nought (γ • ) values through a radiometric calibration analysis, and the corresponding SRTM data were employed for Range-Doppler terrain correction [38]. Finally, a refined Lee filter with a 7 × 7 kernel size was implemented to reduce the undesirable speckle noise and to improve the signal-tonoise ratio [39].

Methodology
After preparing different input datasets, the proposed approach comprises three main steps of (1) input dataset fuzzification, (2) weight value determination of each dataset through an AHP technique, and (3) combination of different datasets using a fuzzy overlay analysis. These steps are separately explained in the following subsections.

Fuzzification
All input datasets (Section 3) were converted to fuzzy values using a linear fuzzy membership function. As illustrated in Figure 6, a linear fuzzy membership function converts the input values of each dataset to the range of 0-1 based on Equation (4). In this transformation, zero and one correspond to pixels with the very low and very high potential risk of flooding, respectively. Accordingly, for two input parameters (i.e., Rainfall, TWI), a linear function with a positive slope was utilized, while for the remaining parameters (i.e., DEM, slope, distance from the main river, LULC, soil type, NDVI, and erosion rate), a negative slope linear function was considered.
where Max and Min are the maximum and minimum values of each parameter, and x is

Analytic Hierarchy Process (AHP)
AHP is a decision-making method that involves organizing multiple-choice criteria into a hierarchy, evaluating their relative values, comparing alternative solutions for each criterion, and deciding on an overall ranking of the alternatives based on cost, benefit, and risk [63,64]. The multi-criteria decision-making technique has extensively been used for solving complex problems through adjusting suitable weights for different input criteria [65]. In this study, AHP was used to calculate the weights of all input parameters for flood hazard mapping. Regarding the determination of weight values, the pairwise comparison matrix was initially generated based on the recommendations of local experts [66]. Then, the normalized matrix was calculated by dividing each element of the pairwise comparison matrix by the sum of each column. Next, the average value of each row was considered as the final weight value of the corresponding parameter [64]. To ensure the correctness and suitability of the weight value determination step, the Consistency Ratio (CR) was computed to investigate the degree of consistency between weight values of different parameters (see Equation (5)). The CR values less than or equal to 0.1 indicate a suitable pairwise comparison matrix, while values greater than 0.1 indicate that the matrix should be reconsidered. To calculate CR, λ max (the maximum eigenvector) and Consistency Index (CI) were computed using Equations (6) and (7), respectively.
where RI is the Random Inconsistency that was set to 1.45 [64], a ij is pairwise comparison matrix element, and w i is the weight value of each parameter. In this study, the CR value was equal to 0.0019, demonstrating the suitability and consistency of weight values. The final weight value of each parameter, determined using the AHP technique, is provided in Table 5.

Fuzzy Overlay
After the fuzzification step, the fuzzy values of each parameter (i.e., each pixel in the corresponding parameter) were multiplied by the associated weight value (i.e., determined in the second step using the AHP technique). Subsequently, these fuzzy parameters were integrated using a fuzzy gamma operator to produce the final fuzzy flood hazard map [67]. (8), the fuzzy gamma operator is the algebraic product of the fuzzy sum and fuzzy product, which are then raised to the power of gamma (γ = 0.90) [68].

As illustrated in Equation
where n is the number of input parameters, f (i) is the pixel value of each input parameter, and F(γ) is the fuzzy flood hazard map.

Jenks Natural Breaks Classification
The Jenks optimization method, also called the Jenks natural breaks classification method, was applied to generate the flood hazard map [69,70]. Natural class breaks are one of the data clustering methods designed to recognize the best possible method to assemble similar values. This method is accomplished through classification values into different classes following the natural breaks or differences in the data. This is performed by minimizing the amount of variance between details in the same class and maximizing the variance between other classes [70]. This technique is frequently used in Geographic Information Systems (GIS) applications for raster data classification [71]. In this study, this method was used for the final classification since the Jenks optimization method has great performance for data with high variance and this was done in the GIS environment [72].

Final Flood Hazard Map
After the production of all input parameters, the Fuzzy logic membership was employed to determine the likelihood of flooding. The AHP method was also used to assess the weights of each parameter (see Table 5). Subsequently, these weights were applied to each fuzzy raster. These fuzzy rasters were then combined using a fuzzy gamma operator to produce the final fuzzy flood hazard map. Finally, the Jenks optimization method was applied to classify the final flood hazard map into five classes based on susceptibility to flooding.

Validation
As discussed, two Sentinel-1 images acquired before and on the day of a flood event were used for validation of the produced map. To this end, after applying the preprocessing steps (see Section 3.10), both images were classified into two classes of water and nonwater by considering the gamma-nought (γ • ) value lower than 0.02 as water areas. Then, the algebraic difference of these classified images produced the flooded areas on 1 April 2019. This flooded map was finally used to validate the result of the flood hazard map. Additionally, the proposed method was applied to another flood event in April 2020 to investigate its robustness. Figure 7 shows the flood hazard map of the study area produced using the proposed methodology. The final fuzzy map was classified into five different classes by natural breaks classification (Jenks) to provide an explicit flood hazard map. Visually, the final classified map has a satisfactory representation of the flood hazard zonation, in which locations with very high flood risk are appropriately bounded by high flood risk regions.

Results
to investigate its robustness. Figure 7 shows the flood hazard map of the study area produced using the proposed methodology. The final fuzzy map was classified into five different classes by natural breaks classification (Jenks) to provide an explicit flood hazard map. Visually, the final classified map has a satisfactory representation of the flood hazard zonation, in which locations with very high flood risk are appropriately bounded by high flood risk regions.    Figure 8. Comparing flooded areas with the produced flood hazard map, it was realized that 74.98% of flood areas were in regions with very high flood risk, 20.18% was located in regions with high flood risk areas, 4.25% was in moderate risk areas, and less than 1% was in low and very low risk regions. Therefore, the obtained results demonstrate the high potential of the implemented approach for an accurate flood hazard map in the study area.

Results
As mentioned in Section 4.6, the proposed method was applied to another flood event in Pol-e Dokhtar which occurred in April 2020 to evaluate its accuracy and robustness (see Figures 9 and 10). This flood event caused many damages to agricultural and urban areas. For example, according to the Copernicus Emergency Management Service (CEMS) reports, 1404 buildings were negatively affected by this flood. During this flood, 20 properties were totally destroyed and 1384 buildings/properties were damaged (Table 6) [73,74]. The comparison of the affected areas in Mamulan and Pol-e Dokhtar during this flood event with the results of the proposed flood hazard map showed that 100% of the buildings which were totally destroyed were located in the very high flood hazard class, and 97.9% of the buildings which were damaged or possibly damaged were over the high and very high flood hazard regions (see Figures 9 and 10). tion between flooded areas and flood hazard map classes are calculated based on Figure  8. Comparing flooded areas with the produced flood hazard map, it was realized that 74.98% of flood areas were in regions with very high flood risk, 20.18% was located in regions with high flood risk areas, 4.25% was in moderate risk areas, and less than 1% was in low and very low risk regions. Therefore, the obtained results demonstrate the high potential of the implemented approach for an accurate flood hazard map in the study area.  ness (see Figures 9 and 10). This flood event caused many damages to agricultural and urban areas. For example, according to the Copernicus Emergency Management Service (CEMS) reports, 1404 buildings were negatively affected by this flood. During this flood, 20 properties were totally destroyed and 1384 buildings/properties were damaged (Table  6) [73,74]. The comparison of the affected areas in Mamulan and Pol-e Dokhtar during this flood event with the results of the proposed flood hazard map showed that 100% of the buildings which were totally destroyed were located in the very high flood hazard class, and 97.9% of the buildings which were damaged or possibly damaged were over the high and very high flood hazard regions (see Figures 9 and 10).  It was important to prove the accuracy and robustness of the proposed method over other study areas. Therefore, in this study, the proposed method was applied to an independent flood-prone area in the province of Lorestan, Iran, called the Nurabad watershed. Figure 11 illustrates the produced flood hazard map over this watershed using the proposed methodology. Figure 12   It was important to prove the accuracy and robustness of the proposed method over other study areas. Therefore, in this study, the proposed method was applied to an independent flood-prone area in the province of Lorestan, Iran, called the Nurabad watershed. Figure 11 illustrates the produced flood hazard map over this watershed using the proposed methodology. Figure 12

Discussion
There are a wide range of approaches for flood analysis and forecasting, such as hydraulic, mathematical, survey and interviews, comparison, statistical, GIS, and remote sensing methods. Most of these methods are very expensive and require special equipment for collecting data and experts to analyze. In this way, this paper offers an efficient approach for flood hazard management community members who do not have extensive knowledge in remote sensing.
It is suggested that the authorities should prevent constructing new buildings in areas that are prone to flooding. However, most urban areas of Pol-e Dokhtar and Mamulan are falling under a high and very high susceptibility flood hazard. Therefore, diversion canals, river defenses, self-closing flood barriers, contour trenching or even dams should be constructed to control floods, especially during heavy rain seasons.
In this study, we improved the accuracy of flood mapping by incorporating the important parameters for flood mapping and improving the accuracy of some of them. For example, an accurate LULC map (i.e., generated using the Sentinel 1 and 2 images) with a 10-m resolution was used in this study instead of outdated low-spatial resolution LULC maps. Likewise, our rainfall data were derived by integrating GPM and synoptic data, which was more accurate than using them individually. Using such high-resolution and enhanced parameters resulted in a highly accurate flood risk map (accuracy =~95%) For instance, although [75] utilized Sentinel-1 images to generate flood risk maps for the north of Iran, its accuracy negatively affected due to using low-resolution rainfall data. Moreover, although [67] used new flood parameters, such as check-dams, all important geospatial datasets for flood risk mapping were not considered.
The accuracy of flood hazard maps was mostly not reported (e.g., research works [6, [76][77][78][79]), and when an accuracy was reported, the accuracy values were based on the train-test split procedure of historically flooded/non-flooded areas. Although the train-test split method is a fast and efficient approach to estimate the performance of a flood model, it is not sufficient when there are a few historical flood points. Moreover, since the nature of the test data were almost the same as training data, the evaluation of models with such a strategy might lead to errors. In this study, we used two validation methods which were more realistic. First, the flooded regions generated by pre-and post-flooded Sentinal-1 images were compared. Second, the results were compared with the CEMS reports. Additionally, the information about damaged buildings was used in the validation process where it was observed that 100% of these buildings were located in areas of very high flood hazard in our produced map. Although the proposed method has several advantages, it also has several limitations as well. For example, many man-made structures such as roads, dams, levees and urban areas, can change flood patterns. Although urban areas are considered in LULC, it would be great if other man-made structures are considered in future studies to improve the accuracy of flood hazard mapping. Additionally, we mainly used open-access remote sensing datasets. For example, the SRTM DEM data with an approximately 30 m spatial resolution were used, which is currently the best free topographic product over the study area. However, we strongly recommend using a high-resolution DEM, particularly when a stream network should be generated by DEM. Furthermore, high-resolution data are optimal for post-flood detection. High resolution airborne or satellite images can be more effectively applied to evaluate the amount of damage to infrastructure and agricultural regions. Moreover, optical remote sensing systems cannot collect observation during nighttime and cannot see under clouds or vegetation canopy.

Conclusions
In this study, we proposed an efficient approach to generate a reliable flood hazard map using remote sensing and GIS datasets and algorithms. To this end, the potential of Sentinel-1 images, multiple-criteria decision analysis, and fuzzy methods were employed. Two case studies, including the Pol-e Dokhtar and Nurabad watersheds, were also considered to evaluate the efficiency and robustness of the proposed approach of the flood hazard mapping. The experimental results showed that the proposed method had a high potential to accurately identify flood risk areas. For example, the accuracy of the risk of the floodprone areas were~95% for these two case studies. It was concluded that using higher resolution and more accurate flood parameters improved the accuracy of the flood risk map.