Flood vulnerability factors
Slope
Slope plays a vital role in flood vulnerability mapping29,30 as it directly affects surface runoff, precipitation dynamics, and the speed at which water flows over the land surface31,32. Areas with low slopes, typically found in gentle floodplains, are more susceptible to flooding and are therefore assigned higher ratings, while steeper slopes are less prone to flooding and receive lower ratings33,34.
To create the slope map, a 30 m resolution DEM was processed in ArcGIS using the Spatial Analysis Tool and the Surface Analysis and Slope functions. The resulting slope raster layer was classified into five categories using natural breaks (Junk), as shown in Fig. 3 and Table 2 (ArcGIS Desktop version 10.8 developed by Esri (https://www.esri.com)). The predominant slope categories are 0.0°–1.7°, 1.8°–4.9°, and 5.0°–8.1°, which collectively cover 31.95%, 26.47%, and 25.28% of the study area, respectively. Steeper slope classes, categorized as steep (8.2°–11.7°) and very steep (11.8° and above), constitute 12.82% and 3.48% of the total area, respectively.
(a) Reclassified raster map of slope and (b) slope severity classification of the study area.
Soil type
Soil types vary in their capacity to absorb water35,36, making infiltration a critical factor in managing stormwater runoff. Infiltration significantly influences the amount of water contributing to surface runoff during rainfall events. Key factors affecting infiltration rates include soil porosity (total pore volume), grain size distribution, pore structure, soil aggregation, and organic matter content37,38.
According to the Ethiopian Ministry of Agriculture’s soil classification, the Wuseta Watershed contains diverse soil types with differing levels of flood vulnerability (Fig. 4 and Table 3) (ArcGIS Desktop version 10.8 developed by Esri (https://www.esri.com)). Pellic vertisols, which are highly prone to flooding, occupy the largest share of the watershed, covering 7.40 km2 (41.78% of the area). Haplic Alisols, located in the high flood vulnerability zones, account for 8.84 km2 (49.91% of the area). Meanwhile, Eutric Vertisols, which are associated with medium flood vulnerability levels, span 1.47 km2 (8.31% of the watershed).

(a) Reclassified raster map of soil type and (b) soil type severity classification of the study area.
Land use/land cover
The classification and analysis of land use/land cover (LU/LC) type were performed using Landsat 8 OLI satellite imagery. The images were sourced from the United States Geological Survey (USGS) Earth Explorer platform (https://earthexplorer.usgs.gov/) and selected based on image quality, prioritizing those with minimal or no cloud cover26,39. The images were geo-referenced to the WGS 84 datum and the Universal Transverse Mercator Zone 37 North coordinate system. Vegetation cover type and density significantly influence soil infiltration capacity, as established in several studies 40,41. These studies indicate that increased vegetation and higher organic matter content enhance soil infiltration while simultaneously reducing bulk density. Conversely, desert environments or areas recently impacted by deforestation due to fire or human activity experience reduced infiltration capacity, resulting in higher surface runoff and shorter lag times12,42.
In the Wuseta Watershed, LULC types (Fig. 5 and Table 4) (ArcGIS Desktop version 10.8 developed by Esri (https://www.esri.com)) were reclassified based on their ability to absorb stormwater. These classifications were standardized to a uniform scale to facilitate flood vulnerability analysis and integrated into the overall flood vulnerability scoring for the land cover factor. Eighty-five (85) sample points were utilized to evaluate classification accuracy, resulting in an overall accuracy of 96.67% and a kappa coefficient of 95.87, indicating a high level of reliability in the classification results.

(a) Reclassified raster map of LU/LC and (b) LU/LC severity classification of the study area.
Elevation
The DEM was converted into an elevation raster layer using the ArcGIS conversion tool. This raster layer was derived from the generated TIN and subsequently classified into five subgroups using the standard equal-interval classification method43,44. This method divides the range of elevation values into intervals of equal size, facilitating the identification of distinct categories. The ArcMap tool automatically determined the interval breaks, and the values were reassigned based on flood risk classification. In this system, the lowest elevations, which are most susceptible to flooding, were assigned a rank of five, while higher elevations, less prone to flooding, received a rank of one43,44.
The procedure for developing the elevation factor mirrored the steps used for slope classification. The study area raster layer was reclassified to reflect its vulnerability to flooding (Fig. 6 and Table 5) (ArcGIS Desktop version 10.8 developed by Esri (https://www.esri.com)). Lower elevations correspond to flatter terrain, which increases flood susceptibility7,33. Consequently, elevation was categorized into five distinct classes based on flood risk.

(a) Reclassified raster map of elevation and (b) elevation severity classification of the study area.
Drainage density
The drainage characteristics of the Wuseta Watershed were extracted from the digitized drainage network (Fig. 7 and Table 6) (ArcGIS Desktop version 10.8 developed by Esri (https://www.esri.com)), with drainage density calculated using the Line Density tool in the Spatial Analyst extension. This tool computes the total length of drainage lines per unit area, offering insight into the efficiency of surface water conveyance33. The resulting drainage density layer was classified into five categories using the natural breaks (Jenks) method, which objectively identifies groupings and patterns in the data. Each category was then assigned a rank from one to five, where a rank of five indicates very high drainage density (and hence greater flood risk), and a rank of one denotes very low drainage density45.

(a) Reclassified raster map of drainage density and (b) drainage density severity classification.
Areas with high drainage density tend to have a greater concentration of stream channels, which accelerates the movement of runoff and reduces infiltration time. This increases the potential for surface water accumulation in downstream low-lying areas, particularly during high-intensity rainfall events. In such zones, the landscape is typically dissected by numerous small channels that contribute to rapid flow convergence and localized flooding. According to46, high drainage density is often associated with impervious or compacted soils and steep terrain, both of which enhance runoff and reduce groundwater recharge. Conversely, areas with low drainage density reflect sparse stream networks and often indicate higher infiltration capacity or flatter terrain, where water tends to spread out rather than rapidly concentrating into flow channels. These regions exhibit reduced runoff generation and typically experience lower flood vulnerability.
In the Wuseta Watershed, the spatial analysis showed that southern and southeastern zones, where drainage density was highest, aligned with the zones classified as highly flood-vulnerable. Field observations also confirmed that these areas experienced frequent surface waterlogging and flash flooding during the Kiremt season, supporting the modeled results. Therefore, the inclusion of drainage density as a key factor in the flood vulnerability model is not merely a generic assumption but is substantiated by both hydrological theory and local empirical evidence. The ranking system effectively reflects the differential impact of drainage network concentration on flood dynamics across the watershed.
Rainfall factor
Rainfall is a critical determinant in evaluating flood severity, as regions receiving higher precipitation are generally more prone to flooding, while areas with lower rainfall face reduced risk. However, in the study area, the Ethiopian Meteorological Agency has established only a single weather station located within the Wuseta watershed. Due to this limited data coverage, rainfall is assumed to be uniformly distributed across the entire region. This assumption significantly undermines the accuracy and reliability of flood severity assessments, as it fails to capture spatial variability in precipitation patterns.
Pairwise comparisons of the factors
Flood vulnerability areas were determined using weighted overlapping factors derived from slope, elevation, drainage density, land use/land cover, and soil type33,47. The weights assigned to each factor were established through consultations with stakeholders and supported by relevant literature. The analysis was conducted within GIS software using the analytical hierarchy process (AHP), a multi-criteria decision-making approach introduced by28. This method relies on pairwise comparisons to derive factor weights, a widely recognized technique in decision analysis48,49. The AHP approach involves creating a square reciprocal matrix, where the principal eigenvectors are used to calculate the weights as shown inTables 7, 8, 9 and 10. Eigenvectors, also known as characteristic or latent vectors, are solutions to matrix equations that represent the relative importance of each factor50. The factor weight coefficients were computed using the decision support/weight module in GIS software, which applied a 9-point continuous scale for pairwise comparisons28. The resulting principal eigenvector represented the combined influence of the flood vulnerability factors. Finally, these weights were normalized to sum up to one, ensuring a consistent weighted linear combination for flood vulnerability mapping44. The thematic values used in the analysis are presented in Tables 7 and 8, with weights determined through a combination of stakeholder engagement and expert judgment. As detailed in Sect. “Field”, the first set of weights was derived from structured questionnaires completed by 42 individuals across seven kebeles (local administrative units) in Debre Markos town, along with insights gathered from seven focus group discussions involving both community members and administrative staff. This participatory approach ensured that local knowledge, lived experiences, and perceptions were directly incorporated into the assessment process. The second set of weights was informed by the researchers’ professional expertise and guided by findings from previous studies on flood vulnerability51,52. These two combinations provide the final values for each thematic factor.
The Consistency Index (CI), which is a measure of departure from consistency, was calculated using the formula:
$$CI = frac{lambda max – n}{{n – 1}}$$
(1)
where, n = number of factors (i.e. 5) and λ = average value of the consistency vector determined in above table, λmax = 5.374.
Based on the above equation, (text{CI}) =5.374–5/5–1 = 0.093.
In order to assess the robustness of the expert view the Consistency Ratio (CR) was calculated using Eq. (2).
$${text{Consistency ratio }}left( {{text{CR}}} right),{ }CR = frac{CI}{{RI}}$$
(2)
where, RI is the random inconsistency index whose value depends on the number (n) of factors being compared; for n = 5, RI = 1.12 as illustrated in Table 9.
Using Eq. (3):
$${text{CR}} = frac{0.09}{{1.12}} = 0.083$$
(3)
Since the calculated consistency ratio of 0.083 is below the threshold of 0.1, it indicates an acceptable level of consistency in the pairwise comparison. As a result, the following weights were assigned to the factors: slope (0.503), elevation (0.26), drainage density (0.134), soil type (0.0678), and land use/cover (0.0348). The derived eigenvector values were used as coefficients for the various flood vulnerability factors: land use/cover, slope, elevation, drainage density, rainfall, and soil type. These weighted factors were then integrated into a weighted overlay analysis in ArcGIS 10.8 to produce the final flood risk map for the Wuseta watershed. The analysis was conducted using the equations outlined below:
$$begin{aligned} {text{Flood}};{text{vulnerability }} & = , 0.{5}0{28 } times , left[ {{text{Slope}}} right] , + , 0.{26 } times , left[ {{text{Elevation}}} right] , \ & quad + , 0.{1344 } times , left[ {{text{Drainage}};{text{density}}} right] , + , 0.0{678 } times , left[ {{text{Soil}};{text{type}}} right] , \ & quad + , 0.0{348 } times , left[ {{text{Land}};{text{use}}/{text{cover}}} right] \ end{aligned}$$
Flood vulnerability map
The flood vulnerability assessment map for the Wuseta watershed (Fig. 8 and Table 11) (ArcGIS Desktop version 10.8 developed by Esri (https://www.esri.com)) was created by integrating flood-generating factors such as slope, rainfall, drainage density, elevation, land use/cover, and soil type using AHP techniques and weighted overlays. The analysis revealed that the watershed contains areas classified as very low (0.001 km2), low (4.41 km2), moderate (7.86 km2), high (4.65 km2), and very high (0.07 km2) flood vulnerability. These classifications were determined based on expert evaluations of the contributing factors. The results indicate that high and very high flood vulnerability zones are predominantly located in the middle and downstream areas of the watershed, which are low-lying and flat. These regions should be prioritized for flood management and mitigation efforts before addressing other parts of the watershed. Notably, one drainage outlet crosses the main road near Debre Markos University, while another outlet with significant flow accumulation is situated close to the Enterprises Cooperative Society near the university. The largest portion of the watershed, comprising 46.28%, falls under the moderate flood vulnerability category, while the upstream areas exhibit slight to very slight flood vulnerability.

(a) Flood vulnerability ranking and (b) flood vulnerability severity in the study area.
Verification of identified flood vulnerability areas
To validate the developed flood vulnerability area, a combination of field surveys, GPS-based ground truth points, and community-sourced flood history data were employed. Site visits provided valuable insights into areas that have historically experienced recurrent flooding (Fig. 1). The comparison revealed a strong spatial correlation between these community-identified flood-prone zones and GPS-based ground truth points with the high-vulnerability areas delineated by the model, thereby supporting the model’s credibility.
The successful verification through field visits and GPS data significantly enhanced the confidence in the accuracy of the flood vulnerability maps, reinforcing their effectiveness in identifying flood-prone regions and informing decision-making processes for mitigation strategies45. In this study, flood vulnerability maps generated using ArcGIS were validated through overlapping of maps developed using GPS-based ground truth points with maps developed by models. This was analyzed as a total of 75 GPS readings were collected from various flood-prone areas, with 62 of these readings (83%) aligning with the identified flood vulnerability zones as statistically analyzed as shown in Fig. 9 (ArcGIS Desktop version 10.8 developed by Esri (https://www.esri.com)). The field validation demonstrated that the flood vulnerability maps generally corresponded well with the actual flood-affected regions33,47. The comparison between the predicted and actual flood-prone areas confirmed the reliability of the maps. Additionally, the validation process revealed areas with higher flood vulnerability that were not initially evident in the maps, highlighting potential gaps in flood risk assessment18,19. This insight emphasizes the importance of iterative validation and community engagement in refining flood risk models. These findings provided valuable insights into regions that require further attention in flood risk management. Overall, the results underscore the reliability and utility of the generated flood vulnerability and risk maps, demonstrating their potential for supporting future flood risk assessments and informing effective mitigation strategies.

Validation of identified flood vulnerability zones using GPS-based field data.
Likewise, the ROC analysis was performed to evaluate the performance of the flood vulnerability map based on classified vulnerability zones: very slight, slight, moderate, high, and very high. A total of 75 reference points were analyzed, among which 62 were actual flooded points and 13 were non-flooded. By progressively aggregating vulnerability zones using thresholds (from very high to very slight), true positives (TP), false positives (FP), false negatives (FN), and true negatives (TN) were calculated for each threshold.
At the highest threshold (≥ 5), which includes only the very high, the map correctly identified 13 flooded points (TP) with zero false positives (FP), yielding a true positive rate (TPR) of 0.10 and a false positive rate (FPR) of 0.00. As the threshold was lowered to include more zones, sensitivity (TPR) increased, but at the cost of reduced specificity (higher FPR). At threshold ≥ 4 (including high and very high), TPR rose to 0.38 with an FPR of 0.14. When zones moderate, high, and very high (≥ 3) were included, the TPR reached 0.721, while the FPR increased to 0.28. At threshold ≥ 2, covering all except very slight, the map achieved perfect sensitivity (TPR = 1.00) but with an FPR of 0.43. Including all zones (≥ 1) resulted in both TPR and FPR reaching 1.00, indicating no discriminatory power. Using the trapezoidal rule, the Area Under the Curve (AUC) was calculated to be 0.81, as shown in Fig. 10, which suggests good discriminatory ability of the flood vulnerability map in distinguishing between flooded and non-flooded areas.

ROC curves illustrating the performance evaluation of flood vulnerability maps.