Identifying Flood Source Areas and Analyzing High Flow Extremes Under Changing Land Use, Land Cover, and Climate in the Gumara Watershed, Upper Blue Nile Basin, Ethiopia

1. Introduction

Recently, flooding has become the most common and leading environmental hazard, with diverse impacts on human life, property, socioeconomic development, and agricultural production [1,2]. It is aggravated mostly by factors such as land use land cover (LULC) and climate changes [3,4]. LULC changes such as deforestation, agricultural expansion, and urbanization result in an increase in surface runoff by changes in the hydrologic characteristics of basins and watersheds [5,6]. While climate change causes flooding by changing the frequency, magnitude, and duration of extreme rainfall events [7,8,9].

For projecting extreme rainfall events that would cause flooding, general circulation models (GCMs) available in the Coupled Model Intercomparison Project Phases 5 and 6 (CMIP5 and CMIP6) play a significant role [10,11]. However, these models are often associated with biases due to their coarse spatial resolution [12]. To address this issue, the use of high-resolution downscaled outputs, such as those from the NASA Earth Exchange-Global Daily Downscaled Projections (NEX-GDDP) dataset, is generally recommended [13,14]. The NASA NEX-GDDP-CMIP5 and CMIP6 datasets provide daily temporal resolution and a spatial resolution of 0.25° (approximately 25 km × 25 km), making them suitable for watershed and basin scale studies [15,16]. Furthermore, the biases associated with climate models can be mitigated by evaluating and selecting the best-performing models, applying bias correction techniques, and using multi-model ensemble means of bias-corrected models [17,18,19].

Historical LULC maps generated from remotely sensed images have been extensively used to assess past watershed and basin conditions [3,20,21,22]. Additionally, projected LULC maps under different scenarios are crucial for understanding future conditions. Consequently, several studies have focused on LULC predictions in different regions [23,24,25]. For example, in the Gumara watershed, Belay et al. [26] projected future LULC changes under a business-as-usual (BAU) and a governance (GOV) scenario. The BAU scenario assumes the continuation of current trends in land use change and socioeconomic development, whereas the GOV scenario incorporates the Green Legacy Initiative (GLI), a program currently being implemented by the Ethiopian government. Their study revealed significant expansion in cultivated land and settlements under the BAU scenario, whereas the GOV scenario predicted increases in forest, shrubland, and grassland. Similarly, a study by Gebresellase et al. [27] in the upper Awash basin of Ethiopia reported the expansion of agricultural lands and urban settlements under the BAU scenario and improvements in forest cover under the GOV scenario.

Flood occurrence can be studied by estimating peak floods for different return periods. This is a crucial aspect in designing safe and cost-effective flood protection structures to minimize the destructive impacts of flooding on the environment [28,29,30]. Return period-based peak floods can be estimated using statistical flood frequency analysis and hydrological rainfall-runoff models. The traditional statistical flood frequency analysis involves estimation of flood quantiles via a best-fitting probability distribution on the basis of the assumption that historical discharge time series are stationary and independent [31]. However, rainfall-runoff hydrological models simulate the transformation of design storms (rainfall) into design floods, taking into account the effects of changing LULC and climate [32,33,34].

Flooding mostly occurs at certain downstream reaches of the main river. Hence, before flood control structures are constructed, quantifying the contribution of upstream subwatersheds to the peak flood magnitude at the main outlet of the watershed is often recommended. This method is commonly known as flood source area identification (FSAI) [35]. FSAI can be performed via multicriteria evaluation (MCE) and the unit flood response (UFR) approach [36]. The MCE approach can be implemented via a geographic information system (GIS), which considers the linear combination of weights of several flood control factors [37,38]. However, this approach does not account for flow routing through rivers. The second approach, UFR, can be performed via hydrological models by accounting for flood routing through rivers for watersheds with diverse conditions [36]. This method was first applied to the Damavand watershed in Iran by Saghafian & Khosroshahi [35] and has since been used in various watersheds. For example, Abdulkareem et al. [39] employed the UFR approach to assess the impact of LULC changes on FSAI for the Kelantan basin in Malaysia. In their study, the authors employed the Hydrologic Engineering Center-Hydrologic Modeling System (HEC-HMS) model and stated its importance for ranking subwatershed units on the basis of their contribution to peak discharge at the main outlet. Similarly, Fiorillo & Tarchiani [40] applied this approach to a catchment in Niger and reported its suitability for identifying flood source areas.

Several studies conducted across various basins and watersheds have demonstrated a potential increase in high-flow extremes in future periods. These studies typically employed hydrological models to simulate future streamflow based on GCM outputs. For instance, Yang et al. [41] assessed future changes in high-flow extremes in the Jiulong River Basin of Southeast China by simulating streamflow under changing LULC and climate conditions. Their findings indicated an increase in high-flow extremes due to intensified extreme rainfall and LULC changes driven by human activities. Similarly, Mahato et al. [42] showed an increase in flood frequency in the Brahmani River Basin, India, by simulating streamflow using a semi-distributed hydrological model driven by outputs from CMIP6 models. In the Ethiopian context, Kassaye et al. [43] showed a projected increase in high-flow extremes in the Baro River Basin by simulating streamflow using six GCMs under different scenarios. Likewise, Yalcin [44] indicated a projected rise in flood frequency in Bitlis Creek, Turkey.

The downstream area of the Gumara watershed, particularly near Lake Tana, is one of the most flood-affected regions in the Upper Blue Nile (UBN) basin of Ethiopia [45]. Several studies have addressed flood risk in this area through different approaches, including flood hazard mapping using the MCE method [46], statistical flood frequency analysis [30], hydraulic flood inundation modeling [47], and flood mapping via Sentinel-1 images [48]. However, return period-based peak flood estimation, FSAI, and extreme flow estimations under varying LULC and climate scenarios remain underexplored in the upstream region of the watershed. Additionally, most flood simulation studies in the Gumara watershed have relied on ground-based rainfall data. However, these observations often contain uncertainties due to missing data, sparse rain gauge networks, and the influence of topography [49,50]. To overcome these limitations, recent studies [51,52,53] have suggested merging ground-based rainfall data with gridded rainfall products (GRFPs). This approach enhances the accuracy of flood simulations by reducing the uncertainties of ground-based rainfall measurements.

Hence, the objectives of this study are to (1) examine the impact of LULC change on return period-based peak floods and flood source areas (FSA) for the observation period (1981–2019) and (2) assess the combined impacts of LULC and climate change on high flow extremes, considering historical (1981–2005), near-future (2031–2055), and far-future (2056–2080) periods. To achieve these objectives, flood simulations were performed using a calibrated and validated HEC-HMS model.

For the observation period, the study utilized merged rainfall estimates from ground-based observed rainfall and three gridded rainfall products (GRFPs): Multi-Source Weighted-Ensemble Precipitation (MSWEP) [54], Climate Hazards Group InfraRed Precipitation with Stations (CHIRPS) [55], and the Enhanced Global Dataset for the Land Component of the Fifth Generation of European Reanalysis (ERA5-Land) [56]. The merging process was performed using a random forest merging algorithm (RF-MERGE). LULC maps for the years 1985, 2000, 2010, and 2019 were incorporated alongside these merged rainfall estimates for this period.

To assess changes in high flow extremes between historical and future periods under changing LULC and climate conditions, high flow simulations were conducted using bias-corrected ensemble means of historical and projected rainfall data from four CMIP5 and four CMIP6 models. For the projected rainfall, the study accounted for representative concentration pathways (RCP4.5 and RCP8.5) for CMIP5 models and shared socioeconomic pathways (SSP2-4.5 and SSP5-8.5) for CMIP6 models. Additionally, future LULC scenarios under the BAU and GOV scenarios for 2035 and 2065 were considered for the near-future and far-future periods, respectively, alongside the future climate scenarios.

In summary, this study offers valuable insights into prioritizing flood source areas in the Gumara watershed and proposes targeted runoff reduction mechanisms. The estimation of high-flow extremes provides critical input for developing effective flood risk protection strategies, while the study’s methodology serves as a framework for flood simulations in similar watersheds across the UBN basin and beyond.

2. Materials and Methods

2.1. Study Area

The Gumara watershed is part of the Lake Tana subbasin, which is located in the UBN basin of Ethiopia (Figure 1). The watershed is primarily drained by the Gumara River, which originates from Guna Mountain. Geographically, it spans from 11° 35′ N to 11° 55′ N latitude and from 37° 35′ E to 38° 15′ E longitude. The drainage area of the watershed is approximately 1,256.7 km². The topography of the watershed is diverse, with elevations ranging from 1,792 to 3,708 m above mean sea level (amsl) (Figure 1c). It is characterized by steep gradients and mountainous terrain in its upstream areas and flat downstream floodplains. The watershed experiences a unimodal rainfall pattern, with the main rainy season occurring from June to September. The watershed receives an average annual rainfall of 1,445 mm, and its average annual temperature is about 20.5 °C [57]. LULC in the watershed is dominated by cultivated land, which accounts for more than 80% of the total area, while the remaining portion of the watershed comprises forest, shrubland, grassland, and settlement areas [26]. The predominant soil texture classes in the region include clay, clay loam, sandy clay, silty clay, and silty clay loam [48]. The majority of the watershed is characterized by scattered rural settlements, where the local livelihood is primarily based on a mixed farming system (crop cultivation and livestock raising) [58]. The downstream part of the watershed, which is adjacent to Lake Tana, is susceptible to frequent flood events due to high runoff generated from the upland watershed [46,59].

2.2. Data Used

For this study, various spatial and temporal data were collected from different sources. These include observed ground-based rainfall and streamflow, gridded rainfall products, climate model data, and spatial data such as elevation, soil, and LULC data.

2.2.1. Observed Ground-Based Rainfall and Streamflow Data

Ground-based daily rainfall data from eight meteorological stations located in and around the Gumara watershed were collected from the National Meteorological Agency (NMA) of Ethiopia, covering the long-term observation period from 1981 to 2019 (Table 1).

The observed daily streamflow time series data (1981–2019) for the Gumara River at the Wanzaye gauging station were obtained from the Ministry of Water, Irrigation, and Electricity (MoIR) of Ethiopia. Selected events from this dataset were utilized for hydrological model sensitivity analysis, calibration, and validation.

2.2.2. Gridded Rainfall Products (GRFPs)

In watersheds with diverse topographies, such as the Gumara watershed, sparsely distributed meteorological stations may not adequately capture rainfall characteristics. To address this issue, three GRFPs were accessed and merged with ground-based observed rainfall. The three GRFPs selected for this study were MSWEP, CHIRPS, and ERA5-Land. MSWEP is a global gridded rainfall product that integrates multiple data sources, including satellite observations, gauge data, and reanalysis products, with a spatial resolution of 0.1° × 0.1° (approximately 10 km × 10 km). CHIRPS combines satellite-based infrared imagery with ground station data to provide a high-resolution rainfall product with a spatial resolution of 0.05° × 0.05° (approximately 5 km × 5 km). ERA5-Land provides high-resolution reanalysis rainfall data produced by the European Centre for Medium-Range Weather Forecasts (ECMWF), with a spatial resolution of 0.1° × 0.1° (approximately 10 km × 10 km). These GRFPs were selected due to their superior performance in reproducing the rainfall characteristics of the UBN basin [60,61]. For this study, CHIRPS and ERA5-Land at each ground-based station were accessed via the Climate Engine (https://www.climateengine.org/), while the MSWEP dataset was extracted using a Python environment. For the three GRFPs, the time series data were accessed at a daily temporal scale for a common period (1981–2019). Table 2 summarizes the spatial resolutions, access, and references for these datasets.

2.2.3. General Circulation Models (GCMs)

Historical (1981–2005), near-future (2031–2055), and far-future (2056–2080) rainfall data from four selected CMIP5 and four CMIP6 models were accessed from the NASA NEX-GDDP dataset. The time series rainfall data for these models were extracted for each ground-based meteorological station using a Google Earth Engine (GEE) sample script, available at this repository: https://code.earthengine.google.com/74643b20363d1f0d4ee5875c456f5c14. For the CMIP5 models, future scenarios considered were RCP4.5 and RCP8.5, while for the CMIP6 models, SSP2-4.5 and SSP5-8.5 were utilized. The models were chosen on the basis of their better performance in reproducing historical extreme rainfall in the UBN basin following our recent work [62]. All the CMIP5 and CMIP6 models are statistically downscaled outputs with a daily temporal resolution and 25 km x 25 km spatial resolution. Table 3 summarizes detailed descriptions, including the resolutions of the climate models, modeling centers, and countries of the selected climate models. Details on the comparison and evaluation of the performance of the CMIP5 and CMIP6 models can be found in Belay et al. [62].

2.2.4. Elevation, Soil, and LULC Data

The elevation data of the study watershed were extracted from a 30-meter spatial resolution digital elevation model (DEM) of the Shuttle Radar Topography Mission (SRTM) [63]. The SRTM DEM data was accessed from USGS Earth Explorer via this link: https://earthexplorer.usgs.gov/. The data was primarily used to generate subwatershed units and other watershed characteristics, including watershed slope and stream network.

Hydrologic Soil Groups (HSGs) for the study area were extracted from a 250-meter spatial-resolution soil grid map established by the USDA. The data was accessed via this link: https://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=1566. As shown in Figure 2c, HSG-D (clay loam, sandy clay, silty clay, silty clay loam, or clay) covers approximately 89% of the study watershed. These soils have high runoff potential because of very slow infiltration rates [64]. The remaining portion of the study watershed is covered by HSG-C (sandy clay loam). These soils have a moderately high runoff potential due to slow infiltration rates [64].

Raster files of historical and future LULC maps of the Gumara watershed, with a 30-meter spatial resolution, were clipped from our previous work [26], which accounted for the entire Gumara watershed. In that work [26], the LULC maps were generated using composite images from Landsat-5/Thematic Mapper (TM) for the years 1985, 2000, and 2010, and Landsat-8/Operational Land Imager (OLI) for the year 2019. For these years, cloud-free images from January to March were accessed from the U.S. Geological Survey via the GEE data catalog. LULC classification was performed using a Random Forest (RF) algorithm [65] in the GEE environment, resulting in the identification of five dominant classes: forest, shrubland, grassland, cultivated land, and settlement areas, which describe the land use conditions of the study watershed. Specific information on the Landsat mission used, acquisition dates, and path/row details can be in the supplementary materials (Table 1S).

In our previous work [26], future land use prediction was performed via the Cellular Automata-Markov (CA-Markov) model available in the land change modeler (LCM) module of IDRISI Selva 17.0 [66]. For the prediction, historical LULC maps (baseline) and various land-use change driver variables were considered, including elevation, slope, distance from roads, distance from streams, distance from towns, and socio-economic development (evidence of the likelihood of change). The prediction was conducted under two "what-if" scenarios: the Business-as-Usual (BAU) scenario and the Governance (GOV) scenario. The BAU scenario considered the current trend of land use conversion and change drivers (transition probability), whereas the GOV scenario accounted for the massive tree seedling plantation program of the Ethiopian government, known as the Green Legacy Initiative (GLI) program [67]. Figures 2d-k show historical and predicted LULC maps of the Gumara watershed. Details on methods used for classification, accuracy assessment, LULC change drivers, scenario-based LULC prediction, and change detection for the entire Gumara watershed can be found in Belay et al. [26].

2.3. Methods

This study analyzed peak floods for different return periods (5-, 10-, 25-, 50-, and 100-year) and flood source areas during the observation period (1981–2019) under varying LULC conditions. Additionally, high flow extremes were assessed for historical (1981–2005) and future (2031–2080) periods under different LULC and climate scenarios. The analysis involved collecting and processing spatiotemporal data, including daily ground-based rainfall, streamflow, GRFPs (MSWEP, CHIRPS, and ERA5-Land), CMIP5 and CMIP6 models, elevation, LULC, and soil data.

Ground-based rainfall data were merged with the three GRFPs using the RF-MERGE algorithm [65] to enhance data accuracy. The merging process included four main steps: input data preparation, setting up the RF model, training the RF model, and performing the merging (prediction). For merging, ground-based rainfall data were used as the dependent variable, while GRFPs and other auxiliary data were used as independent variables. The resulting merged rainfall product was subsequently utilized in rainfall-runoff modeling, for extracting annual maximum (AM) 1-day rainfall series, and for correcting biases in climate model outputs. In the study, rainfall data from CMIP5 (MIROC5, MRI-CGCM3, CanESM2, and MPI-ESM1-2-LR) and CMIP6 (MPI-ESM1-2-LR, MRI-ESM2-0, CNRM-CM6-1, and BCC-CSM2-MR) models were also bias-corrected, and their ensemble means were predicted for historical and future periods to reduce uncertainties associated with climate models. Bias correction was performed using the quantile mapping (QM) method [68]. From the corrected daily rainfall data, AM 1-day rainfall series were extracted and subsequently used to estimate high flow extremes in the near-future and far-future periods.

Flood simulation was conducted using the HEC-HMS model (version 4.11) [69]. The process included subwatershed delineation, entering rainfall data, model selection, parameter estimation and sensitivity analysis, calibration and validation, and flood simulation under various climate and LULC scenarios. Using the calibrated model, peak floods for different return periods were determined, flood source areas prioritized, and high flow extremes for historical and future periods assessed. The overall methodology is outlined in Figure 3, with detailed explanations in subsequent sections.

2.4. Data Quality Control

In hydroclimatic studies, the use of long-term, homogeneous, and consistent data is generally recommended. However, the ground-based stations considered in this study have missing rainfall data and short-term rainfall records at the Hamusit, Wanzaye, and Arb Gebeya stations. The missing data were filled via the inverse distance weighting (IDW) technique, whereas short-term rainfall records were extended via the maintenance of variance extension, type 1 (MOVE.1) technique [70]. After filling in missing data and extending short-term data, a homogeneity test was employed to detect change points in the mean annual rainfall data via the nonparametric Pettitt test [71]. According to the Pettitt test statistics presented in the supplementary materials (Table 2S), no significant changes were detected in the annual rainfall data for the stations considered in this study.

Double-mass curve (DMC) analysis was also employed to check the consistency of the mean annual rainfall data. This analysis involved inspecting significant changes in slope by plotting the cumulative mean annual rainfall of all stations on the x-axis and the cumulative mean annual rainfall of each station on the y-axis. From the analysis, all the stations were found to be consistent (Figure 4a), since no significant changes in slope were detected in the DMC plot. Figures 4b, c, and d also present mean annual rainfall at each ground-based station, mean monthly rainfall, and mean monthly discharge data, respectively.

2.5. Random Forest Merging (RF-MERGE)

To improve the quality of the rainfall data, the random forest (RF) merging approach (RF-MERGE) was employed to combine daily ground-based observed rainfall with three GRFPs (MSWEP, CHIRPS, and ERA5-Land). RF is a supervised machine learning algorithm that can be used for both classification and regression [65]. Among the various merging approaches, RF-MERGE was selected based on our previous work [62], conducted in the study region. In that study, RF-MERGE demonstrated the best performance in reproducing observed rainfall compared to individual GRFPs and other merging approaches, as evaluated by both statistical and categorical performance measures. Additionally, RF-MERGE requires a minimal number of model parameters and is capable of handling high-dimensional data [52].

For the merging process, ground-based rainfall data, three GRFPs, and auxiliary data (including date, month, year, and elevation) were first prepared. Next, the RF model was set up to map the relationship between the dependent and independent variables. In this study, ground-based rainfall data served as the dependent variable, while the three GRFPs, along with the auxiliary data, were used as independent variables. Subsequently, the RF algorithm was trained using rainfall data recorded at stations located in the UBN basin (Table 3S in the supplementary materials). During training, the RF model learned the relationship between the dependent and independent variables. Based on this, the optimal parameter values of the RF algorithm, such as the number of trees ($N_t$), the minimum number of bootstrap samples ($N_{min}$), and the number of features ($N_f$) at the splitting node, were estimated. Finally, using the trained RF model, rainfall was predicted at stations located in and around the Gumara watershed (Table 2). The merging procedure was performed via the "randomForest" R package [72]. The analysis yielded estimated parameter values of 200, 6, and 3 for ${ N}_t$, $N_{min}$, and Nf, respectively. These parameters were determined through trial-and-error procedures and by comparison with parameters estimated in previous studies [73,74]. Using these parameters, the trained model was subsequently used to predict daily rainfall at the eight testing stations. The merged rainfall data were then used to extract the AM 1-day rainfall, correct the bias of the climate models, and calibrate/validate the hydrological model. A comprehensive review, including assumptions, procedures, and limitations of various merging approaches, such as RF-MERGE, can be found in our previous work [50].

2.6. Design Rainfall and Depth‒Duration‒Frequency (DDF) Curve

In this study, the annual maximum 1-day rainfall data series for the observation period (1981-2019) was extracted and fitted to commonly used probability distributions: normal (NOR), lognormal (LN3), gamma (GAM), Gumbel (GUM), Pearson type III (PE3), log Pearson type III (LPE3), and generalized extreme value (GEV) distributions. The parameters of these distributions were calculated via the L-moments method provided in Hosking & Wallis [31]. The Anderson‒Darling (AD) test [75] was employed to assess goodness-of-fit by comparing the CDFs of the observed data with those of the candidate distributions. This test is sensitive to deviations in the tails of the probability distribution, making it valuable for applications involving extreme values [76]. The best-fit distribution was selected on the basis of the lowest AD statistic ($A²$) at the 5% significance level and used for estimating design rainfall for different return periods.

To simulate flood hydrographs for different return periods (5-, 10-, 25-, 50- and 100-year), a short-duration (example with an hourly distribution) depth‒duration‒frequency (DDF) curve is often needed. For this, a short-duration hourly distributed DDF curve was developed using the following equation:

$$R_t=R_{24}\left(\frac{t}{24}\right){\left(\frac{b+24}{b+t}\right)}^n$$

(1)

where $R_t$ is the computed rainfall depth for a duration of t hours, $R_{24}$ is the design rainfall depth for a 24-hour (1-day) duration, $t$ is the rainfall duration in hours, $b$ is a parameter that controls the relationship between different durations and can be determined on the basis of regional studies, and $n$ is an exponent that adjusts the amount at which rainfall decreases with shorter durations and can be determined empirically on the basis of regional studies. For this study, the DDF curve for a 1-hour duration was developed using 0.3 for the parameter $b$ and 0.9 for the exponent $n$ following regional studies conducted in Ethiopia [77,78].

2.7. Bias Correction

It is common for climate model simulations to exhibit bias (uncertainties) due to their coarse spatial resolution, model assumptions, boundary conditions, spatial averaging, and other factors [79,80]. To address this issue, it is generally recommended that biases in climate models be corrected before their outputs are applied in impact assessment studies [81]. Among the different types of bias correction methods, the QM approach was used in this study because of its better performance in the study region [62]. The QM method uses the quantile‒quantile relationship to match the distribution function of climate model simulations with observed data [68]. The method assumes that the cumulative distribution functions (CDFs) of the observed and simulated datasets can be matched, meaning that the CDFs of the observed and simulated data should be similar enough for meaningful bias correction [82]. Based on this assumption, for a given simulated value $x_{sim}$ , QM bias correction adjusts it using:

$$x_{cor}=F_{obs}^{-1}\left[F_{sim}\left(x_{sim}\right)\right]$$

(2)

where $x_{cor}$ is bias corrected value, $F_{sim}\left(x_{sim}\right)$ is CDF of the simulated data, $F_{obs}\left(x\right)$ is CDF of the observed data, and $F_{obs}^{-1}$ is inverse CDF (quantile function) of the observed data.

In this study, prior to implementing QM bias correction, the CDFs of daily observed rainfall from RF-MERGE and simulated rainfall from historical CMIP5 and CMIP6 were computed and plotted, as shown in Figure 5, to assess their similarity. The graph indicates that the CDFs of the observed and simulated data are sufficiently similar, allowing the QM method to be applied.

Among the different types of QM approaches, the parametric transformation function (PTF) was chosen for this study because of its effectiveness in studying extreme rainfall events [83]. The bias correction was applied by matching the historical climate model data with the RF-MERGE data and then applying the correction to the future climate model data on the basis of historical relationships. A detailed explanation of the assumptions and formulations of the various QM bias correction methods can be found in Enayati et al. [82].

2.8. Rainfall-Runoff Modeling

Rainfall-runoff modeling can be conducted using various distributed and semi-distributed hydrological models. According to Ghonchepour et al. [84], factors such as modeling objectives, hydrologic conditions, limitations, assumptions, and data availability should guide model selection. In this study, the HEC-HMS model was chosen for its suitability in simulating rainfall-runoff processes under varying LULC and climatic conditions, aligning well with the study's objectives. HEC-HMS is a physically based, semi-distributed model widely recognized for its effectiveness in flood modeling [85,86,87]. It is compatible with the input data used in this study, including rainfall, streamflow, LULC, soil, and DEM data. Its semi-distributed structure facilitates accurate flow simulation at the subwatershed scale, which is critical for identifying flood source areas. Furthermore, HEC-HMS has been thoroughly tested and shown to perform well in modeling floods in similar watersheds within the UBN basin of Ethiopia [88,89].

The HEC-HMS model consists of several key components (managers), including basin models, meteorological models, time-series data, and control specifications [69]. The basin model defines the physical characteristics of the watershed, such as subwatersheds, river reaches (networks), junctions, outlets, and others, along with relevant hydrological parameters. The meteorological model specifies the type of rainfall data used in the simulation. The time-series data component is used to input rainfall and streamflow data, which are essential for driving the model. Lastly, control specification manages or sets the simulation time step, duration, and output intervals. These components work together to simulate hydrologic processes and provide flood hydrographs.

In this study, rainfall-runoff simulations in the HEC-HMS model involved subwatershed delineation (basin model), entering rainfall and streamflow data, model selection (for runoff volume, direct runoff, and routing), parameter estimation and sensitivity analysis, model calibration and validation, and flood hydrograph simulations under changing LULC and climate conditions as detailed in the subsequent sections.

2.8.1. Subwatershed Delineation

Subwatershed delineation was performed using a 30-meter resolution SRTM-DEM via the GIS tools in HEC-HMS 4.11, as this version is integrated with GIS components. The process involved setting up the coordinate system (projection), terrain reconditioning, filling sinks, drainage line processing, identifying streams, defining outlets, and delineating subwatersheds. Since the middle and upper portion of the Gumara watershed lies in steep gradient areas (slope greater than 15%), the 30-meter DEM resolution could introduce inaccuracies, particularly in calculating flow direction and accumulation. These limitations may result in misrepresentation of the river network of the watershed. To address this issue, DEM reconditioning techniques, such as stream burning and sink filling, were applied to improve the representation of hydrological features. Additionally, during subwatershed delineation, the homogeneity of various watershed characteristics, such as soil, LULC, elevation, and stream networks, was considered for merging similar smaller subwatersheds into one.

2.8.2. Rainfall and Streamflow Data

This study utilized daily rainfall data from RF-MERGE (1981–2019), rainfall DDF curves for return periods of 5, 10, 25, 50, and 100 years, and rainfall for both historical and future periods to simulate streamflow. Observed streamflow data for the Gumara River (1981–2019) were also employed to compare with simulated flows during model calibration and validation, focusing on selected flood events.

2.8.3. Model Selection

In the HEC-HMS, there are various models for computing runoff volume, direct runoff, and channel flow routing. For this study, the Soil Conservation Service Curve Number (SCS-CN) [90], SCS unit hydrograph (SCS-UH) [90], and Muskingum channel routing [91] models were selected for modeling runoff volume, direct runoff, and channel flow routing, respectively.

I.

Modeling runoff volumes (loss model)

Models available in the HEC-HMS generally compute runoff volume by subtracting the sum of various losses from the accumulated rainfall. Among these models, the Soil Conservation Service Curve Number (SCS-CN) method was selected for this study because it can compute runoff volume for a given rainfall event as a function of HSGs, LULC, and antecedent soil moisture conditions.

In the SCS-CN model, the excess rainfall ($R_e$) at a given time $t$ is the rainfall in excess of the infiltration capacity and other losses and can be computed as:

$$R_e=\frac{{(R-I_a)}^2}{R-I_a+S}$$

(3)

where $R$ is the accumulated rainfall depth at time $t$; $I_a$ is the initial abstraction (loss); and $S$ is the maximum retention potential of the watershed. Considering the analysis findings from many experimental watersheds, the SCS-CN provides an empirical relationship of $S$ and $I_a$, which is given as $I_a=0.2S$.Then, the expression of $R_e$ at a given time $t$becomes:

$$R_e={\displaystyle \frac{{\left(R-0.2S\right)}^2}{R+0.8S}}where(R\ge 0.2S)$$

(4)

The maximum retention ($S$) can be computed as a function of the CN via the following equation.

$$S={\displaystyle \frac{25400-254CN}{CN}}$$

(5)

In this study, 30-m spatial resolution gridded runoff CN maps were generated by considering soil and LULC characteristics (raster files), as well as the average antecedent moisture conditions (AMC-II). For the soil moisture condition, the AMC-II was accounted for because it represents a standard condition commonly used in hydrological modeling. The gridded CN maps were prepared via the ARC-GIS 10.4 environment [92]. The process involved preparing input data by aligning raster files (LULC and soil) to the same resolution, coordinate system, and extent; reclassifying LULC and HSG maps; combining them; and assigning curve numbers using a lookup table with a raster calculator. Figure 6 shows the gridded runoff curve numbers for both the historical and future periods generated based on these procedures. The results indicate that the CN values under the BAU scenario are relatively greater than those under the governance (GOV) scenario. This is attributed to the expansion of cultivated land and settlement areas under the BAU scenario. As a result, the maximum retention potential of the watershed will reduce, which may lead to an increase in peak flow occurrence due to higher surface runoff volumes. The composite CN for each subwatershed was estimated on the basis of the mean value of raster grids (pixels) within the boundaries of the subwatershed area. Similarly, the composite CN for the entire watershed was calculated as the mean value of raster grids within the boundary of the whole watershed.

II.

Direct runoff (transform model)

In the HEC-HMS model, several empirical and conceptual models are available to transform excess rainfall into direct runoff. Among these, the SCS unit hydrograph (SCS-UH) model [90] was selected for this study. The key parameter of this model is the lag time ($T_{lag}$), which represents the time difference between the center of mass of excess rainfall and the peak of the SCS-UH. Although the assumption of a single-peak hydrograph is a limitation of the SCS-UH, it was chosen because its parameters can be estimated from the physical characteristics of the watershed, specifically in relation to the time of concentration$T_c$. For this model, $T_c$ can be computed via the Kirpich equation [93], and $T_{lag}$ is then calculated as 0.6 times $T_c$:

$$T_C=0.0195\left({\displaystyle \raisebox{1ex}{$L^{0.77}$}\!\left/ \!\raisebox{-1ex}{$S^{0.385}$}\right.}\right)$$

(6)

where $T_C$ is the time of concentration (minutes), $L$ is the length of the longest flow path (meter), and $S$ is the average channel slope (m/m).

III.

Channel flow (routing model)

From the available flow routing models in the HEC-HMS model, the Muskingum routing model was selected for this study. This model uses a finite difference approximation of the continuity equation and is widely used for flood routing [94]. It accounts for flow attenuation resulting from the storage effect of the channel as flood flows travel through it.

The Muskingum routing model computes channel storage as follows:

$$S_t=K\left[X{I}_t+\left(1-X\right)O_t\right]$$

(7)

where ${ S}_t$, $I_t$, and $O_t$ are storage, inflow, and outflow at time $t$, respectively; $K$ is the travel time of the flood wave though the routing channel; and $X$ is a dimensionless weight factor ($0\le X\le 0.5).$ For this study, the Muskingum routing model parameters ($K$ and $X$) were determined through model calibration.

2.8.4. Model Sensitivity Analysis

Sensitivity analysis is a critical step in hydrological modeling aimed at reducing uncertainty and identifying model parameters that significantly influence the model's output [86,95]. In this study, sensitivity analysis was conducted for the parameters of the runoff volume (loss) and transform models by varying each parameter individually while keeping the others constant (local sensitivity analysis). Each parameter value was adjusted by ±25%, and the relative percentage change in peak discharge was assessed for each simulation. The sensitivity of each parameter was evaluated on the basis of the magnitude of the change in peak flow. The relative change in peak flow as a percentage was computed via the following equation:

$$Relative change\left(\mathrm{\%}\right)={\displaystyle \frac{Q_{new}-Q_{base}}{Q_{base}}}\times 100$$

(8)

where $Q_{new}$ is the peak flow after changing the parameter and where $Q_{base}$ is the peak flow from the baseline simulation. In this study, model sensitivity analysis was performed for selected events that occurred from 13 to 31 August 2010.

2.8.5. Model Calibration and Validation

Model calibration and validation were performed to match the observed and simulated data for selected events. Both automated and manual calibrations were employed to obtain better agreement between the simulated and observed flow events through iterative adjustments of the model parameters. The performance of the model during calibration and validation was evaluated through graphical comparison (visual comparison) and statistical performance measures such as the coefficient of determination ($R^2$), Nash–Sutcliffe efficiency ($NSE$), root mean square error ($RMSE$), and percent bias ($PBIAS$). These performance measures were selected from different categories: $R^2$ from standard regression-based methods, $NSE$ from dimensionless methods, and $RMSE$ and $PBIAS$ from error-based statistical methods [96].

$R^2$ expressed in equation 9 describes the proportion of the variance in the observed data explained by the model-simulated data.

$$R^2=\frac{{\sum }_{i=1}^N(O_i-\overline{O})(S_i-\overline{\mathrm{S}})}{\sqrt{{\sum }_{i=1}^N(\mathrm{O}-\overline{O}{)}^2}\sqrt{{\sum }_{i=1}^N{(S_i-\overline{S})}^2}}$$

(9)

where $O_i$ is the $ith$ observed discharge value for the event being evaluated, $S_i$ is the $ith$ simulated discharge value for the event being evaluated, $\overline{O}$ is the mean of the observed discharge values, $\overline{S}$ is the mean of the simulated discharge values, and $N$ is the total number of observations. $R^2$ values range from 0 to 1 (inclusive), with higher values indicating less error variance, and values higher than 0.5 are generally considered acceptable values [97].

The NSE is a dimensionless statistical performance measure that describes the relative magnitude of the residual variance (noise) compared with the observed data variance (information) [98]. $NSE$ is expressed as shown in equation 10:

$$NSE=1-\frac{{\sum }_{i=1}^N{(0_i-S_i)}^2}{{\sum }_{i=1}^N{(O_i-\overline{O})}^2}$$

(10)

$NSE$ranges from $-\infty$to $1$ $(-\infty <NSE\le 1)$, where $NSE=1$ indicates good agreement between the observed and simulated values, $0<NSE\le 1$ is generally considered acceptable, and $NSE\le O$ indicates poor or unacceptable model performance.

$RMSE$ expressed in equation 11 describes the standard deviation of the model prediction error.

$$RMSE=\sqrt{\frac{1}{N}{\sum }_{i=1}^N{(O_i-S_i)}^2}$$

(11)

The optimal value of $RMSE$ is 0, which indicates perfect agreement between the observed and simulated values.

where $PBIAS$ indicates the average tendency of the model-simulated values to be higher or lower than the observed values [99]. It is given as shown in equation 12:

$$\mathrm{PBIAS}(\%)=\frac{\sum _{i=1}^N(o_i-S_i)}{\sum _{i=1}^N{o}_i}\times 100$$

(12)

The optimal value of $PBIAS$ is 0, with low values indicating accurate model simulation. Positive $PBIAS$ values indicate overestimation, whereas negative values indicate model underestimation. According to Moriasi et al. [100], $-25\%<PBAS<25$ is an acceptable range for discharge simulation.

In this study, model calibration was performed for two selected events: Event-1, which occurred from July 1 to August 31, 1996, and Event-2, which occurred from July 5 to July 31, 2008. Model validation was conducted via Event-3, which occurred from August 2 to August 27, 2014. These events were specifically chosen to account for periods with high flood peaks.

2.9. Flood Source Area Identification (FSAI)

Flood source area identification (FSAI) is an important approach for identifying the main sources of downstream flooding. It helps rank subwatersheds on the basis of their contribution to peak flood magnitude at the main watershed outlet. In this study, FSAI was performed using the UFR, which was introduced by Saghafian & Khosroshahi [35]. In the UFR approach, the contribution of each subwatershed can be computed as:

$$F{I}_n=\frac{Q_{\mathit{all}}-Q_{\mathit{all}-n}}{Q_{\mathit{all}}}\times 100$$

(13)

$$f{i}_n=\frac{Q_{all}-Q_{all-n}}{A_n}$$

(14)

where ${FI}_n$ is the gross flood index of the $n^{th}$ subwatershed in percent, $Q_{all}$ is the outlet peak discharge in m³/s when all the subawtersheds are included in the base simulation, $Q_{all-n}$ is the outlet peak discharge in m³/s when $n^{th}$ subawatersheds are excluded from the simulation, ${fi}_n$ is the flood index of the $n^{th}$ subawtershed (m³/s/km²) based on the area of the subwatershed, and $A_n$ is the area of the $n^{th}$ subwatershed in km².

In the UFR approach, the subwatershed with the highest flood index value is ranked first (indicating the greatest contribution to peak discharge), whereas the subwatershed with the lowest flood index value is ranked last, indicating the least contribution. In this study, a calibrated HEC-HMS model incorporating all subwatershed components was first developed. The flood indices FI and fi for each subwatershed were subsequently calculated by disconnecting or eliminating a chosen subwatershed unit while maintaining the connections of the other subwatersheds. This process was repeated for all subwatersheds, allowing for a comprehensive analysis. In this regard, the analysis was also conducted under different conditions to investigate the impacts of LULC changes on the flood index values.

2.10. Combined Impacts of LULC and Climate Change on High Flows

The combined impacts of LULC and climate change on high flows, specifically on the annual maximum (AM) 1-day flows, were analyzed by simulating historical and future conditions using the calibrated HEC-HMS model, driven by AM 1-day rainfall. To assess these impacts on high flow extremes, flow duration curves (FDCs) were constructed from the AM flow series. The FDC construction process involved (1) sorting AM flows in descending order, (2) ranking the AM flows, (3) calculating the probability of exceedance for each flow value, and (4) plotting the AM flow series on the y-axis against the probability of exceedance on the x-axis. In this study, the probability of exceedance was calculated using the Weibull formula:

$$P=\frac{m}{N+1}\times 100$$

(15)

where $P$ is the probability of exceedance (%), m is the rank of the flow value, and N is the total number of flow data points.

3. Results

3.1. Rainfall Climatology

3.1.1. Mean Monthly and Annual Rainfall

Figure 7 presents the mean monthly rainfall of the study watershed, estimated using both observed data and GCM-simulated rainfall data. Observed rainfall estimates were derived by merging ground-based observations with three GRFPs (MSWEP, CHIRPS, and ERA5-Land) using the RF-MERGE approach, which enhances accuracy by integrating diverse sources. Historical and future GCM-simulated rainfall estimates were derived from statistically downscaled and bias-corrected ensemble means from CMIP5 (MIROC5, MRI-CGCM3, CanESM2, and MPI-ESM1-2-LR) and CMIP6 (MPI-ESM1-2-LR, MRI-ESM2-0, CNRM-CM6-1, and BCC-CSM2-MR) models.

Results indicate that the GCMs effectively captured observed rainfall patterns, despite projected increases in rainfall under RCP and SSP scenarios. Approximately 80% of the watershed’s rainfall occurs during the main rainy season (June to September), consistent with the regional monsoon. For the historical period (1981–2005), the mean annual rainfall was estimated to be 1336 mm based on RF-MERGE data, while CMIP6 models estimated a slightly higher mean annual rainfall of 1379 mm. For future simulations, the highest projected mean annual rainfall was 1495 mm under the SSP5-8.5 scenario (2056–2080), indicating an expected increase in mean annual rainfall in the future period. These findings align with previous studies [101,102,103], which similarly reported projected increases in rainfall in the region.

Figure 8 illustrates the spatial distribution of mean annual rainfall (MARF) in the Gumara watershed for both the historical and future periods. The spatial distribution patterns are similar for both periods, with the primary difference being an increase in rainfall during the future period. Overall, the mean annual rainfall in the watershed ranges from approximately 1320 mm to 1560 mm, with the highest rainfall values concentrated in the northeastern part of the watershed and the lowest rainfall occurring in the northern part (Figure 8).

3.1.2. Design Rainfall and Depth‒Duration‒Frequency (DDF) Curve

Design rainfall for the study watershed for various return periods was estimated by fitting the AM 1-day rainfall extracted from RF-MERGE estimates (Figure 9a) to candidate distributions. The goodness-of-fit test identified the GEV distribution as the best fit. Consistent with this finding, similar studies [104,105] have reported the reliability of the GEV distribution for predicting extreme rainfall indices, such as the AM 1-day rainfall. The estimated design rainfall values for a 1-day duration using the GEV distribution were 70 mm, 75 mm, 81 mm, 85 mm, and 88 mm for the 5-, 10-, 25-, 50-, and 100-year return periods, respectively. However, to analyze flood occurrence, shorter-duration rainfall depths are often required for hydrological models. For this, a DDF curve for a 1-hour duration was developed using a regional equation. Figure 9b presents this DDF curve, which predicts extreme rainfall intensity over various durations and return periods.

3.1.3. Annual Maximum (AM) 1-Day Rainfall

The historical and projected annual maximum 1-day rainfall series were extracted from CMIP5 and CMIP6 multimodal ensemble means. The extracted 1-day rainfall was converted to hourly distributed rainfall using a regional equation and subsequently used to simulate the annual maximum flow series. The box plot in Figure 10 illustrates the distribution of the extracted data. As shown, the data exhibits high variability across all periods, with an increase in future periods relative to historical values. Across all periods, the extracted maximum 1-day rainfall ranges from 23 mm to 83 mm, with the highest values observed under the SSP climate scenarios.

3.2. Land Use Land Cover Change Analysis

LULC maps of the study watershed for the historical (1985, 2000, 2010, and 2019) and future (2035 and 2065) years were extracted from our previous work [26], which analyzed the entire Gumara watershed. Figures 11a and b present the percentage coverage of different LULC classes for both the historical and future years, respectively. During the historical period, a significant increase in cultivated land and settlements was observed at the expense of a decrease in other land use classes. From 1985–2019, the percentage area of cultivated land and settlements increased from 63.70% (800.07 km²) to 88.39% (1110.18 km²) (Figure 11a). In contrast, the percentage areas of forest, shrubland, and grassland decreased from 4.98%, 28.21%, and 2.35% to 1.97%, 7.88%, and 0.59%, respectively (Figure 11a). These findings align with previous studies by Anteneh & Mohammed [106] and Chakilu & Moges [107], which also reported an increase in agricultural cultivation in the region.

Figure 11b shows the predicted percentage area of each LULC class in the study watershed for the years 2035 and 2065 under the BAU and GOV scenarios. Under the BAU scenario (which is based on the current trend of LULC change), cultivated land is expected to dominate, covering 91.34% and 91.78% of the area in 2035 and 2065, respectively, at the expense of decreases in other LULC classes (forest, shrubland, and grassland). Conversely, under the GOV scenario (which is based on the green legacy initiative program), forest cover is predicted to increase, reaching 50.98 km² (4.06%) in 2035 and 66.06 km² (5.26%) in 2065. In this scenario, the expansion of cultivated land will be more controlled than in the BAU scenario. These findings are consistent with a study conducted in the upper Awash basin of Ethiopia by Gebresellase et al. [27], which indicated an increase in agricultural lands and settlements under the BAU scenario and improvements in forests under the GOV scenario.

3.3. Hydrological Modeling

3.3.1. Subwatershed Delineation

Accounting for the homogeneity of various watershed characteristics, such as soil, LULC, elevation, and stream networks the Gumara watershed was delineated into eight subwatersheds, labeled W1 to W8 (Figure 12). The area of the subwatersheds ranges from 99.9 km² (W6) to 215 km² (W4) (Figure 12). Table 4 also presents the computed subwatershed (area and slope) and stream characteristics (length and slope).

3.3.2. Model Sensitivity Analysis

Sensitivity analysis was performed to prioritize parameters for calibration and improve model efficiency. In this study, sensitivity analysis of the HEC-HMS model was conducted on key parameters, including $CN$, which controls runoff volume, $I_a$ that influences the proportion of rainfall lost before runoff begins, and $T_{lag}$ that affects the timing of peak discharge. For this analysis, flow events from August 13 to August 31, 2010, were considered (Figure 13). The results indicate that the $CN$ was the most sensitive parameter. A direct relationship between flow and the $CN$ was observed: increasing the $CN$ resulted in higher discharge, while decreasing the $CN$ led to lower discharge. Specifically, a 25% increase in the $CN$ led to a 9.79% rise in peak discharge, while a 25% decrease in the $CN$ caused a 12.71% reduction in peak discharge (Table 5). These findings are consistent with studies by Tassew et al. [89] in the Gilgel Abay catchment and Zelelew & Melesse [108] in a small watershed in northwest Ethiopia, both of which also identified the $CN$ as the most sensitive parameter. In contrast, $I_a$ and $T_{lag}$ were less sensitive, with minimal impact on peak discharge. Furthermore, since the SCS-CN method assumes uniform soil, land use, and rainfall conditions within a subwatershed, the model output (runoff volume) may be subject to uncertainties in watersheds with diverse topography, such as varying slopes and elevations. To address this limitation, it is crucial to adjust CN to local conditions based on field surveys and topographic data.

3.3.3. Model Calibration and Validation

Model calibration and validation were performed for selected events. These events were chosen on the basis of the availability of gap-free discharge data and discharge series with high flood peaks. Figure 14 presents a comparison of the observed and simulated discharge hydrographs for these selected events.

In flood modeling, the most critical aspect of the hydrograph is the peak discharge, as it is the peak that has the potential to cause flooding downstream. In this study, the model adequately captured the timing of the hydrograph peaks, although it tended to slightly overestimate peak discharge during calibration and underestimate during validation (Figure 14 and Table 6).

Table 6 summarizes the observed and simulated peak discharge and model performance evaluation results during the calibration and validation periods in terms of ${\mathrm{R}}^2$, $\mathrm{N}\mathrm{S}\mathrm{E}$, $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}$, and $\mathrm{P}\mathrm{B}\mathrm{I}\mathrm{A}\mathrm{S}$. The observed and simulated discharges were in close agreement. For example, ${\mathrm{R}}^2$ value of 0.89, $\mathrm{N}\mathrm{S}\mathrm{E}$ value of 0.82, $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}$ value of 27.1 m³/s, and $\mathrm{P}\mathrm{B}\mathrm{I}\mathrm{A}\mathrm{S}$ value of -5.19% were found for Event 1 (calibration). Like during the calibration period, the observed and simulated discharges closely agreed during the validation period (Event 3), with ${\mathrm{R}}^2$ value of 0.87, $\mathrm{N}\mathrm{S}\mathrm{E}$ value of 0.85, $\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}$ value of 26.7 m³/s, and $\mathrm{P}\mathrm{B}\mathrm{I}\mathrm{A}\mathrm{S}$ value of -3.99%. According to Moriasi et al. [100], the statistical performance measures computed in this study are within an acceptable range, as the $R²$ and $NSE$ values exceed 0.5, and the $PBIAS$ values fall within ±25% for all events, indicating minimal error. Similarly, previous studies have demonstrated the effectiveness of the HEC-HMS model in simulating discharge in the UBN basin of Ethiopia. For instance, Tassew et al. [89] applied HEC-HMS in the Gilgel Abay catchment, UBN basin of Ethiopia, and achieved R² values between 0.68 and 0.95 and NSE values ranging from 0.63 to 0.92 during calibration and validation for selected flow events, indicating acceptable model performance. Similarly, Yilma & Kebede [109] used HEC-HMS in the Dabus subbasin, achieving calibration results with R² = 0.81 and NSE = 0.78, while validation yielded R² = 0.87 and NSE = 0.79, indicating the reliability of the HEC-HMS model for rainfall-runoff simulations. These findings suggest that the HEC-HMS model can reliably simulate discharge in watersheds with topographic, soil, and LULC characteristics similar to the Gumara watershed. Table 7 also presents the calibrated (optimal) parameter values of the HEC-HMS model for the runoff volume (loss), direct runoff (transform), and routing components. These calibrated parameter values fall within the acceptable ranges specified in the HEC-HMS Technical Reference Manual [110].

3.4. Return Period-Based Peak Discharge and Runoff Volume Simulation

Flood and runoff volume hydrographs were simulated in the calibrated HEC-HMS model. Figures 15a and b compare the peak discharge and total runoff volume of simulated hydrographs under different LULC conditions based on DDF curves of different return periods. From the results, increases in peak discharge and runoff volume were observed in the study watershed from 1985 to 2019. Among the simulation results, maximum peak discharge values of 211.7, 278.5, 359.5, 416.7, and 452.7 m³/s were found for the 5-, 10-, 25-, 50-, and 100-year return periods, respectively, under LULC-2019 (Figure 15a). For the same period, runoff volume values of 5.9, 7.6, 9.6, 11.05, and 11.7 mega cubic meters (MCM) were found for the respective return periods (Figure 15b).

Historical observations from 1981 to 2019 indicate that the average peak flood magnitude of the Gumara River was approximately 260.0 m³/s, closely matching the simulated peak flood magnitudes for the 5- and 10-year return periods. The simulated results also align with historical flood events recorded during the main rainy months (July and August) in 1996, 2006, 2008, 2014, and 2019, which caused significant inundation of cultivated lands and rural settlements downstream of the Gumara watershed [45,111]. Furthermore, previous studies [30,47] that employed at-site flood-frequency analysis within the Gumara watershed reported peak flood magnitudes comparable to the simulated results of this study.

Table 8 also shows the relative changes in peak discharge and runoff volume between two different LULC conditions (LULC-1985 to LULC-2019). The results indicate that the relative effect of LULC change decreases in flood peaks and runoff volumes of higher return periods, which is consistent with previous studies [34,112,113]. For example, the LULC change from 1985–2019 increased the peak discharge by 38.7% for a 5-year return period but only by 22.4% for a 100-year return period. Similarly, the LULC change from 1985–2019 increased the runoff volume by 24.9% for a 5-year return period but only by 17.6% for a 100-year return period. This could be due to the diminishing impact of land use changes for higher return periods, such that further increases in rainfall do not significantly affect peak discharge [39]. The analysis results also show that the impact of LULC change is more noticeable on peak discharge than on runoff volume.

Historical observations from 1981 to 2019 indicate that the average peak flood magnitude of the Gumara River was approximately 260 m³/s, which aligns closely with the simulated peaks for the 5- and 10-year return periods of LULC-2019. The simulated the return period-based peak flood magnitudes simulated using are also consistent with historical flood events recorded during the main rainy months (July and August) in the years 1996, 2006, 2008, 2014, and 2019. These events resulted in substantial inundation of cultivated lands and rural settlements downstream of the Gumara watershed.

3.5. Flood Source Area Identification

The ranking of subwatersheds on the basis of their contribution to peak discharge at the main outlet of the Gumara watershed was determined via the UFR approach [35]. In this analysis, peak discharge was simulated in the calibrated HEC-HMS model using the DDF curve of the 50-year return period. Table 9 presents the peak discharge contributions and rankings of subwatersheds under different LULC conditions. According to column 3 of Table 9, subwatershed W1 generated the highest peak discharge values under LULC-1985 and 2000 (74.9 and 81.8 m³/s, respectively). In contrast, subwatersheds W2 and W6 generated the smallest peak discharge values under the respective LULC conditions (37.5 and 39.9 m³/s, respectively). Under LULC-2010 and 2019, subwatershed W4 generated the largest peak discharge values (84.6 and 90.7 m³/s, respectively), whereas W6 (42.9 m³/s) and W2 (51.4 m³/s) generated the smallest peak discharge values for the respective LULC conditions.

Additionally, subwatershed units or (flood source areas) were ranked on the basis of flood index values such as $FI$ and $fi$ (Table 9). According to the results, different ranking index values were found under different conditions. For example, under LULC-1985, the W2 subwatershed, which is seventh in terms of area and eighth in terms of peak discharge generation, acquired seventh and sixth ranks in terms of the $FI$ and $fi$ indices, respectively (Table 9, columns 2, 3, 6, 8). While it maintains the same $FI$ rank under LULC-2000, it jumps to seventh in terms of the $fi$ index. Similarly, under LULC-2010, the W8 subwatershed, which is sixth in terms of area and peak discharge generation, acquired seventh rank in the $FI$ index and sixth rank in the $fi$ index (Table 9, columns 2, 3, 6, 8). However, it jumps to eighth and seventh ranks in terms of the $FI$and $fi$ indices, respectively, under LULC-2019.

Furthermore, previous studies [39,113] recommend the $fi$ index for ranking because it considers changes in peak discharge per unit area, unlike the $FI$ index, which considers only changes in peak discharge. Therefore, more attention needs to be given to subwatersheds with high $fi$ index values during prioritization to reduce runoff generation in these subwatersheds. Figure 16 shows the $fi$ index of each subwatershed of the Gumara watershed for the baseline period. According to the results, W6 was the most runoff-generating subwatershed under LULC-1985, with a $fi$ index value of 0.303, whereas W4 was the most runoff-generating subwatershed under LULC conditions of 2000, 2010, and 2019, with $fi$index values of 0.328, 0.34, and 0.368, respectively. These show the combined effects of LULC and other related factors such as the location of subwatersheds within the entire watershed, river routing, elevation, slope, LULC, stream network, reach length and slope, rainfall distribution, and other physical features of the watershed under study in prioritizing subwatersheds [36,113,114]. The combined effects of these factors can be well simulated via hydrologic models such as the HEC-HMS model used in this study.

3.6. Combined Impacts of LULC and Climate Change on High Flows

Historical and projected high flows of the study watershed were simulated using the calibrated HEC-HMS model, based on the AM 1-day rainfall. Table 10 below presents the mean AM flows, along with relative percentage change values. The results indicate that greater high-flow magnitudes were observed when future climate scenarios were combined with the BAU land use scenario compared to the GOV scenario. Specifically, the highest flow value of 254.8 m³/s was found under the combined BAU-2065 land use and SSP5-8.5 (2056–2080) climate scenario, while the highest value of 230.6 m³/s was found under the combined GOV-2065 land use and SSP5-8.5 (2056–2080) climate scenario (Table 10).

Considering the relative change values, higher positive change values were observed when future climate scenarios combined with the BAU scenario ranged from 11.3% to 20.5% (Table 10). Among the computed changes under the BAU scenario, the greatest increase is projected under SSP5-8.5 (2056–2080) with a 20.5% increment. These significant increases in high flows in the future suggest a higher likelihood of more frequent flooding in the region.

In contrast, the change values are smaller when future climate scenarios are combined with the GOV land use scenario (Table 10), ranging from 5.9% to 9.0%. Among the computed changes under the GOV scenario, the greatest increase is projected under SSP5-8.5 (2056–2080) with a 9.0% increment. These lower change values can likely be attributed to afforestation efforts and controlled cultivation under the GOV scenario, which could enhance infiltration, reduce surface runoff, and contribute to more stable streamflow, especially during extreme rainfall events.

Furthermore, flow duration curves (FDCs) were constructed to assess the combined impact of LULC and climate change on high flows (Figures 17a and b). The results showed significant variations of high flows, particularly when different climate scenarios combined with the BAU land use scenario (Figure 17a) than combined with the GOV scenario (Figure 17b). The change in the high flow of the study watershed was also assessed by classifying high flow magnitudes into four categories depending on different ranges of probability of exceedance (Q0–Q25, Q26–Q50, Q51–Q75, and Q76–Q100) (Table 11). Based on these categories, the change in high flow is more visible in the Q0–Q25. In this category (Q0–Q25), the SSP5-8.5 (2056–2080) scenario showed greater change in high flow when combined with the BAU-2065 and GOV-2065 land use scenarios as shown in Figures 17a and b, respectively.

In summary, the use of FDCs provides a robust approach to capture and compare the impacts of different LULC and climate scenarios on high flow extremes. The increases in high flow extremes under future scenarios, particularly in the BAU case, emphasize the urgent need for sustainable land use practices and climate adaptation strategies to mitigate future flood risks in the study watershed.

4. Discussion

4.1. Flood Source Area Identification and Analysis of High Flow Extremes

The downstream part of the Gumara watershed is the most frequently flood-impacted area due to the high runoff contributed by the upstream watershed. This area, which includes cultivated land, rural and urban settlements, and infrastructure, is adversely affected by frequent floods, particularly during the main rainy months of July and August [46,47]. Land use and land cover (LULC) changes, such as deforestation and the expansion of cultivated land, combined with increases in extreme rainfall due to climate change in the upstream watershed, are key factors driving the increased frequency of floods in the region. For these reasons, this study focused on FSAI and the analysis of high-flow extremes under changing LULC and climate conditions in the upstream part of the watershed.

From the FSAI analysis, subwatersheds labeled from W1 to W8 showed varying flood index values under different LULC conditions. Consistent with these findings, a study by Saghafian et al. [113] in the Golestan watershed, located northeast of Iran, and Abdulkareem et al. [39] in the Kelantan River Basin, Malaysia, showed the significant contribution of LULC in identifying flood source areas. Overall, for the Gumara watershed, subwatersheds labeled W4 and W6 were identified as the most significant runoff-contributing subwatersheds. This suggests that these two subwatersheds play the largest role in contributing to flooding in the downstream area, making them priority areas for intervention. To mitigate runoff and reduce downstream flooding, it is essential to implement measures in these subwatersheds that enhance infiltration. These measures could include effective land use management strategies (such as afforestation, reforestation, terracing, and contour farming), soil and water conservation practices (e.g., check dams and gabions), and green infrastructure solutions.

The projected high flow extremes for the future period in this study consistently showed increases relative to the historical period. Notably, the percentage increases were higher when future climate scenarios were combined with the BAU land use scenario. This can be attributed to the fact that, under the BAU scenario, the clearing of forests and other vegetative cover for agricultural expansion leads to decreased infiltration rates and increased surface runoff, resulting in higher flows during extreme rainfall events. These findings are in line with previous studies [115,116], which also projected an increase in streamflow under future climate scenarios combined with the BAU land use scenario. In contrast, when different climate scenarios were combined with the GOV land use scenario, the percentage increases were lower. This is because land use management strategies under the GOV scenario, such as afforestation, reforestation programs, and controlled expansion of cultivated land and settlements, are expected to enhance infiltration and reduce runoff during extreme rainfall events.

Among all simulations of high flows for future periods, the most alarming increase in high flows was observed under the SSP5-8.5 (2056–2080) climate scenario combined with the BAU-2065 land use scenario consistent with recent [42,43,44]. These results indicate that the Gumara watershed will face more frequent floods during this period, which could significantly disrupt the socioeconomic conditions of the community, as their livelihoods depend on crop cultivation. To reduce the adverse impact of future flood events, this suggests a pressing need to implement various flood protection and adaptation strategies in the downstream part of the watershed. First, it is recommended that local governments implement strict land use management policies that promote afforestation and reforestation to counterbalance the impacts of climate change on flooding. Second, local communities, basin initiatives, and non-governmental organizations should actively participate in soil and water conservation activities. Third, the construction of flood protection infrastructure in the downstream areas, particularly near cultivated land and settlements, is essential to mitigate the adverse effects of flooding on socioeconomic activities and infrastructure. The feasibility of these flood protection and adaptation measures requires collaboration among local governments, non-governmental organizations, basin initiatives, and the local community.

Overall, this study investigated the impacts of LULC and climate changes on FSAI and high-flow extremes in the Gumara watershed. The identification of flood source areas provides valuable insights for prioritizing subwatersheds that contribute most significantly to peak discharges. This spatially explicit information can guide targeted interventions to mitigate flood risks. By integrating historical and future LULC scenarios with projected climate data, the findings indicate a substantial increase in high-flow extremes in future periods. These increases suggest the likelihood of more frequent flood events, which could adversely impact downstream areas. In this regard, the study provides critical insights into future flood risks and offers a foundation for developing effective flood risk management strategies to mitigate the adverse impacts of frequent flooding in downstream areas. Furthermore, the methodology developed in this study, particularly the application of merged rainfall data from multiple data sources in flood studies, offers a significant contribution to improving hydrological simulations.

4.2. Generalizability of the Study Findings

The Gumara watershed, located in the UBN basin, has unique geographic and socio-economic characteristics. Geographically, it is characterized by steep upstream slopes and flat downstream floodplains. Socioeconomically, the watershed is dominated by rural areas with scattered settlements, where the livelihood is primarily based on a mixed farming system (crop cultivation and livestock raising). While these characteristics make it a valuable case study for assessing the impacts of LULC and climate change, they may limit the generalizability of findings. The results of this study may be most applicable to watersheds with similar socio-economic conditions and geographic settings. Regions with different geographic and socio-economic characteristics, such as urbanized areas or those with advanced flood management infrastructures, may experience varying hydrological responses. However, the study’s methodological framework offers a transferable approach that can be adapted to other regions with comparable challenges.

4.3. Limitations of the Study and Future Research Directions

As mentioned in the methods section, historical and future climate scenarios (RCP and SSP) from GCMs are often associated with uncertainties [80,82]. To reduce these uncertainties, it is recommended to correct biases in climate models and use ensembles of multiple models instead of relying on a single model. In this study, bias correction of climate models and the ensemble mean of GCMs were employed to mitigate these uncertainties. Future research could explore various machine learning ensemble prediction (MLEP) approaches to further reduce uncertainties in climate models and use these outputs as inputs to hydrological models.

Both the BAU and GOV LULC scenarios considered in this study rely on assumptions about future socioeconomic and policy developments. However, these assumptions may not fully capture the uncertainties of future land use changes. For instance, the BAU scenario assumes continuity in current land use trends, while the GOV scenario assumes improvement of forest cover. However, sudden policy shifts and regional economic changes could significantly alter the trajectory assumed in these scenarios. In this context, future research could explore the uncertainties associated with these scenarios by considering more intermediate scenarios to further reduce the uncertainties associated with land use projections.

In this study, the antecedent soil moisture was assumed to be at an average condition, as it represents a standard condition for hydrological modeling. Future research could integrate rainfall data with satellite-based soil moisture data to further enhance flood simulations using hydrological models.

Additionally, future research could extend this study by incorporating hydraulic flood modeling for downstream areas of the Gumara watershed. Using the outputs of this study as inputs for detailed hydraulic modeling would enable a more comprehensive analysis of flood inundation and hazard mapping in the region. This approach could provide valuable insights for flood risk management and inform targeted mitigation strategies.

5. Conclusions

This study assessed the impacts of LULC and climate changes on high flow extremes in the Gumara watershed using the HEC-HMS model. Merged rainfall data (1981–2019) and multi-model ensemble means from CMIP5 and CMIP6 models were utilized for historical (1981–2005), near-future (2031–2055), and far-future (2056–2080) periods, with future climate scenarios projected under RCP4.5, RCP8.5, SSP2-4.5, and SSP5-8.5 pathways. Historical LULC maps (1985, 2000, 2010, and 2019) and projected LULC scenarios (BAU-2035, BAU-2065, GOV-2035, and GOV-2065) were incorporated to analyze the effects of LULC change on flood risk.

The long-term observation period (1981–2019) showed a significant increase in peak discharge under recent LULC conditions, with peak values rising by 24.9% for the 5-year return period and by 17.6% for the 100-year return period, indicating a decreasing influence of LULC changes with higher return periods. The flood source areas identified using the UFR approach revealed subwatersheds W4 and W6 as consistent contributors to high runoff across various LULC conditions, making them priority areas for flood mitigation efforts.

Future projections indicate a further increase in high flows, especially under the SSP5-8.5 (2056–2080) climate scenario combined with BAU-2065 land use, where flows are projected to increase by up to 20.5%, suggesting an increased flood risk in the region. In contrast, the GOV scenario, which promotes afforestation, shows a smaller flow increase of up to 9.0%, demonstrating that sustainable land management practices can help mitigate climate change impacts on extreme flows.

In summary, this study integrated RF-MERGE, GIS, and HEC-HMS models to investigate the effects of LULC and climate changes on high-flow extremes in the Gumara watershed. The findings provide valuable insights for watershed managers and land use planners to prioritize subwatersheds and develop effective flood risk management practices. The study suggests the coordinated participation of local governments, non-governmental organizations, basin initiatives, and communities to collaboratively mitigate the adverse effects of future flood events. Future research could expand on this work by using climate model ensembles based on MLEP approaches, integrating satellite-based soil moisture data, and conducting hydraulic flood modeling downstream of the Gumara watershed. These efforts would further inform flood mitigation strategies across the region.