Evaluating the Performance of Hydrological Models for Flood simulation in Wangchu River Basin, Bhutan

1. Introduction

Floods have become increasingly frequent and severe worldwide, largely driven by climate change, which disrupts precipitation patterns and amplifies extreme rainfall events. These changes have already had devastating impacts. For instance, floods claimed over 220,000 lives globally between 1980 and 2013, and in 2017 alone, weather-related disasters caused economic losses exceeding US$300 billion [1,2]. Projections further underscore the escalating risks, with climate change expected to increase population vulnerability to river floods by 20–80% by 2030 and 40–150% by 2080 [1]. The Intergovernmental Panel on Climate Change (IPCC) Sixth Assessment Report forecasts a rise in global mean surface temperature of 1.5℃ to 4℃ above pre-industrial levels by the end of the 21st century [3]. This temperature increase is projected to disrupt precipitation patterns, leading to more extreme rainfall events and heightened flood risks in many river basins [4,5,6] As a result, these changes will have profound impact on lives, livelihoods, economics, and properties worldwide.

Given these escalating risks, there is an urgent need for comprehensive climate change adaptation and disaster risk reduction strategies to mitigate the growing threat of river floods. Regions at particularly high risk according to the IPCC [3] and Nagamani et al. [6], include South Asia, Sub-Sahara Africa and parts of Latin America. These areas are expected to experience some of the most severe impacts, highlighting the need for immediate and targeted action.

In South Asia, Bhutan is particularly vulnerable to severe flood hazards due to a combination of fragile landscape factors. The country’s rugged terrain, with steep slopes and narrow river valleys, along with its location in the path of the seasonal monsoon system, makes it highly susceptible to flood events [7,8]. Climate-induced changes further exacerbate this vulnerability, leading to significant impacts on the hydrological system.

Over the past several decades, Bhutan has experienced recurrent devastating floods, with major events recorded in 1950, 1956, 1960, 1968, 1994, 2000, 2004, 2009, 2016 and most recently in 2023 [8,9]. These flood events have resulted in significant loss of life, livelihoods, and historical monuments. Projections indicate that river discharge in the region will change annually, with annual rainfall expected to increase by 12.1% by 2050 and 27.1% by 2099 [10,11,12]. Notably, Syldon et al. [12] have identified the Wangchu River Basin (WRB) as a highly vulnerable region in Bhutan, where river inundation could increase by up to 60%, posing a severe threat to low-lying agricultural land.

Given the growing frequency and severity of floods, timely implementation of flood warnings and advanced hydrological modeling techniques have become crucial to mitigate the substantial damages [13,14]. Flood mitigation can be most effective by combining structural measures, non-structural measures, and institutions capacity [15,16]. However, in Bhutan’s mountainous terrain, structural measures such as discharge stage and cable system become impractical during high flood events. In this context, flood simulation can be approached through software-based hydrological models (a non-structural measure), which allows for more accurate river discharge estimation at various segments, irrespective of the timing and magnitude of flood events [16]

Developed countries have made significant advancements in predicting floods by employing advanced hydrological models. For instance, the USA has developed models like Hydrologic Engineering Centre – Hydrologic Modeling System (HEC-HMS), Soil and Water Assessment Tool (SWAT) and National Water Model (NWM) [13], while the European Commission created the Global Flood Awareness System (GloFAS) for global flood forecasting [17,18]. In contrast, developing countries, including Bhutan, face significant challenges in flood prediction due to data deficiencies, financial constraints, technical limitations, and lack of tailored flood forecasting models [13], resulting in increased vulnerability and uncertainty in flood prediction and response [19,20].

In Bhutan, the lack of monitoring in small and medium-sized rivers exacerbates flood prediction uncertainties, making it critical to select appropriate models for accurate simulations. Various hydrological models have been employed in Asia, with the Integrated Flood Analysis System (IFAS) providing particularly advantageous for countries such as developing and data-scarce regions like Bhutan [21,22]. IFAS has been used in Japan, Vietnam, Malaysia, the Philippines, Myanmar and Bhutan [23,24,25,26,27,28], showing promising results. The model’s integration of satellite rainfall data, Digital Elevation Models (DEM), and other geographical datasets makes it ideal for flood simulation in remote areas.

Similarly, the HEC-HMS model, although globally recognized, has seen limited use in Bhutan’s mountainous river basins [14,29,30,31]. Dorji et al. [32] and Fakhruddin [19] employed the model for flood simulation in river basins of Bhutan, however, more research is needed to explore its application in these challenging terrains.

Additionally, the development of Group on Earth Observation Global Water Sustainability (GEOGloWS) model, offering global historical discharge data and daily forecasts for around one million sub-basins worldwide [33], has demonstrated high accuracy in flood simulation across various countries like Australia, Brazil, Colombia, the Dominican Republic, Nepal and Peru [18,34]. Despite its success elsewhere, it is yet to be applied in Bhutan, indicating the need for comprehensively evaluation of each model for enhanced flood forecasting and risk management in Bhutan.

While the use of hydrological models in large basins is widespread, systematic evaluation of their relative efficiencies and advancement remain limited [14,16,22,29]. These models have primarily been employed for long-term discharge simulation, fewer studies focusing on flood event simulation. Furthermore, the lack of comparative studies and data in small basin presents a considerable challenge, particularly in regions like Bhutan where small rivers and streams are poorly monitored. To address these gaps, this study evaluates three hydrological models—IFAS, HEC-HMS and GEOGloWS—focusing on their suitability for simulating discharge during flood events in Wangchu and Thimchu river basins of Bhutan. Additionally, the study examines whether model performance varies between the area of 3556km² and 658.86km², respectively, using key performance metrics such as Nash-Sutcliffe Efficiency (NSE), Percentage Bias (PBIAS) and Peak Flow Error (PFE) metrics. The findings of this research will provide crucial insights for flood forecasting and risk management, thereby enhancing disaster preparedness in the face of increasing flood risks.

2. Materials and Methods

2.1. Study Area Description

This study focused on Bhutan, a small mountainous kingdom nestled between India and China. Specifically, we examine the upper Wangchu River Basin (WRB) in Western Bhutan, a densly populated region with 240,012 inhabitants and two economically important hydropower plants (Tala - 1020MW and Chhukha - 336MW) [11,35]. The comparative analysis was conducted in WRB (3,556 km²) and the Thimphu River Basin (TRB, 658.86 km²), as illustrated in Figure 1. The WRB is also known as Large River Basin (LRB) while the TRB is referred to as Small River Basin (SRB). Geographically, the area spans three distinct agro-ecological zones: alpine (3,500-7,500 m), cool temperate (2,600-3,600 m), and warm temperate (1,800-2,600 m) [12]. The basin experiences an annual average temperature of 13.4°C and an average of rainfall of 575.1mm/month in July and August. Flow measurements at the Chimakoti station (WRB) from 2001 to 2010 revealed seasonal variations in river discharge, with the highest flow of 251 m³/s occurring in August and the lowest of 26 m³/s in January [35].

2.2. Research Flow

Figure 2 illustrated the research methodology, encompassing the models setup, simulation, calibration, and validation stages for three hydrological models: IFAS, HEC-HMS and GEOGloWS. While IFAS and HEC-HMS models require extensive preprocessing of datasets from satellite and observation stations, the GEOGloWS model offers a more streamlined approach through its pre-developed package. The GEOGloWS model is freely available and provides immediate access through a user-friendly web-service interface. All three models use global data sets to address local data limitations. Models performance were evaluated using three commonly employed performance metrics: NSE, PBIAS and PFE metrics [36,37,38].

2.3. Data Preparation

In three hydrological models, we selected observed river discharge from two distinct flood events for calibration. The calibration period spanned from May 15th to June 14th, 2014, while the validation period covered October 1st to October 31st, 2021. These periods were specifically chosen to coincide with notable flood events, providing a robust foundation for model calibration and validation. For both the IFAS and HEC-HMS models, we employed rainfall data corresponding to these same time frames, as rainfall data serves as a crucial input parameter that directly influences flow patterns within the catchment area. The observed discharge and rainfall data, archived by National Centre for Meteorology and Hydrology (NCHM), Bhutan were used in these models. Daily rainfall data from six meteorological stations and discharge data from two hydrological stations were utilized in this study (Table 1).

For the GEOGloWS model, a minimum of one-year daily observed data is required for bias correction of simulated data at the respective hydrological station (Lozano et al., 2025). For this purpose, observed discharge from 1st January to 31st December, 2021 was used in the model. The bias corrected discharge was extracted from 1st to 31st October, 2021 enabling direct comparison with the validation results from the IFAS and HEC-HMS models.

2.4. Description and Parameter Setting Selected Model

2.4.1. IFAS Model

The IFAS model, developed by the International Centre for Water Hazard and Risk Management (ICHARM), Japan, is freely available and can be downloaded from https://www.pwri.go.jp/icharm/research/ifas/index.html. It serves the rainfall-runoff simulation model designed for the efficient computation of discharge and flood forecasting especially in ungauged river basins [4,39]. The model works based on the principle of tank model for flood simulation and kinematic wave model for routing [26,40]. The model need an input data of DEM, land-use, soil data, and ground and satellite rainfall (Table 2 and Figure S1). Additionally, it features Geographic Information System functionality for the automatic generation of catchment boundary, river channel network, sub-basins, default baseline parameters, and import of boundary river basin.

In this study, as presented in Figure 2, the DEM, land-use, soil, and rainfall data were imported in the model from project information manager tool. The observed rainfall data was prepared in .csv file format and forced in the IFAS model. Then, the IFAS model converted rainfall data to discharge. Two-layer tank model was selected in this study as recommended in the previous studies [21,23]. River discharge was calculated from river course tank model based on Manning’s equation and Kinematic wave method [40]. The input parameters for river course tank model include the breadth of the river course, coefficient of the resume law, Manning’s roughness coefficient, initial water level in the river course, infiltration from river tank to the aquifer tank, and coefficient for cross section shape.

The IFAS model has an automatic function to generate default parameters. Default parameters were used for the initial simulation, however, it failed to synchronize with observed discharge pattern. Therefore, the parameters were calibrated accordingly to align simulated discharge with observed discharge using trial and error method manually. Table S1, S2 and S3 shows all the initial parameters of surface, aquifer and river course tank.

We incorporated all 28 model parameters during the simulation process. However, 5 sensitive parameters were calibrated based on recommendation of the previous studies in the regional studies [16,41]. Notably, parameters were not calibrated separately for SRB, as it falls with the LRB. These identified parameters were also suggested to apply in river basins of Bhutan having similar topography as highlighted in reports and seminar of NCHM, Bhutan.

The parameters such as Final Infiltration Capacity (SKF) and Surface Runoff Coefficient (SNF) values were increased in this model to capture the high flow. The Runoff of confined aquifers (AGD) and Runoff of confined aquifers (AUD) parameter values were reduced to capture low flow (Table 3). The parameters were calibrated carefully to reproduce flow dynamics within the acceptable range. The coefficient for cross section parameter was not calibrated because the simulated high flow aligned with the observed flow. Overall, the simulation with calibrated parameters demonstrated good synchronization with observed river discharge.

2.4.2. Description and Processing of HEC-HMS Model

The HEC-HMS model, developed by the US Army Corps of Engineers, is a physically based semi-distributed rainfall-runoff model used for simulating discharge in dendritic watersheds [29]. The model is widely recognized for its accuracy and extensive use for both event-based and continuous discharge simulations. The primary inputs to the model are precipitation, DEM, Land-use and soil type. We used Alaska Satellite Facility radio-metrically corrected terrain 12.5 resolution Advanced Land Observing Satellite Phased Array type L-band Synthetic Aperture Radar DEM to create subbasins and its characteristics. The DEM was downloaded from https://search.asf.alaska.edu/#/. The soil data with the resolution of 250m was downloaded from https://wocatapps.users.earthengine.app/view/dss-bhutan. The land-use type data of 2016 prepared by National Land Commission of Bhutan was used.

In this study, the river flow was simulated using different meta-models and methods of HEC-HMS such as (1) loss model (Soil Conservation Service Curve Number (SCS-CN)), (2) transform model (SCS Unit hydrograph (SCS-UH)) and (3) routing model (Muskingum). Detailed explanations and processing approaches are explained subsequently.

• Loss model

The SCS-CN method was used to estimate runoff from the catchment [42,43,44]. The method was chosen for its versatility and wide applicability in estimating surface runoff from each subbasins, utilizing equation 1, 2, 3, and 4 [45]. This method required data such as CNs, initial abstraction, and percentage imperviousness of the basin. For this, CNs was prepared in ArcGIS 10.8 with HEC-GeoHMS extension tool using land-use and hydrologic soil group (HSG). The raster data were converted to polygon and created union raster file. The SCS lookup table was created (Table S4) to assign CNs to land-use and HSG following the Technical Report - 55 standard and the NRCS (Natural Resources Conservation Service) land-use table [46]. The equations 1 to 7 compute parameters such as average CNs, accumulated precipitation, maximum retention potential, Time of Concentration ($\mathbf{T}\mathbf{c}$)and Lag time (${\mathbf{L}\mathbf{a}\mathbf{g}}_{\left(\mathbf{t}\right)}$). The curve number was calculated using the following equation.

$${\mathrm{C}\mathrm{N}}_{\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{g}\mathrm{e}}={\displaystyle \frac{{{\left(\mathrm{C}\mathrm{N}\right)}_{\mathrm{i}}\left(\mathrm{A}\right)}_{\mathrm{i}}}{{\mathrm{A}}_{\mathrm{T}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}}}$$

(1)

where, $\mathbf{i}$ is the number of sub-basin, Ai is the area of the particular sub-basin and ${\mathbf{A}}_{\mathbf{T}\mathbf{o}\mathbf{t}\mathbf{a}\mathbf{l}}$is the total area of the basin.

The accumulated excess precipitation is calculated by equation 2

$${\mathrm{P}}_{\mathrm{e}}={\displaystyle \frac{{(\mathrm{P}-\mathrm{I}\mathrm{a})}^2}{\mathrm{P}-\mathrm{I}\mathrm{a}+\mathrm{S}}}$$

(2)

where ${\mathbf{P}}_{\mathbf{e}}$ is accumulated precipitation excess, $\mathbf{P}$ unaccumulated rainfall depth, Ia is the initial abstraction and S is maximum retention potential.

The $\mathbf{I}\mathbf{a}$ is calculated by equation 3.

$$\mathrm{I}\mathrm{a}=0.2\mathrm{S}$$

(3)

And $\mathbf{S}$ is calculated following equation 4

$$\mathrm{S}={\displaystyle \frac{25400-254\mathrm{C}\mathrm{N}}{\mathrm{C}\mathrm{N}}}$$

(4)

The HEC-HMS model used 10 sub-basins to simulated river discharge in both LRB and SRB (Figure 3a). Subbasin2 covers the largest area of 617.58km², while the subbasin9 has the smallest area of 148.01km². Subbasin7 has the longest flow path of 54.01km. Drainage density falls under low category where it ranges from 0.066km/km² to 0.118km/km². The slope of the subbasins ranges from 47.91 to 60.03 degree indicating steep slope. The steeper the slope, the faster the surface runoff, thus reaches the outlet more rapidly.

The basin have four HSG namely dystric cambisols, eutric cambisols, haplic acrisols, and haplic lixisols classified under HSG B, C, and D (Figure 3c). The predominant soil type is dystric cambisols, while haplic acrisols, and haplic lixisols are found in scattered patches throughout the basin. HSG B indicates moderately low runoff potential, whereas HSG D indicates high runoff potential. Eleven land-use types were identified dominated by forest (Figure 3d), indicating low runoff potential and significant forest canopy abstraction. The spatial distribution of average curve numbers ranges from 60 to 92 in the basin (Figure 3b). Subbasin3 has the highest CN of 86 while subbasin8 has the CN of 60.

The calculated T_lag and calibrated parameters used in this model are presented in Table 4 and Table 5. The longest T_lag of 4.18hours was recorded in subbasin7 while the shortest T_lag of 2.04hours was recorded in subbasin4. The initial abstraction for the study area was 0.2mm, meaning the 20% of the rainfall is prevented by forest canopy from reaching to the land surface. Imperviousness of the area was 1.9. The K and X values were calibrated at 7 and 0.1, respectively. These calibrated parameters were used in synchronizing simulated discharge with observed discharge at hydrological stations.

2.

Transform model

SCS-UH was used to transform excess rainfall into surface runoff. This method requires ${\mathbf{L}\mathbf{a}\mathbf{g}}_{\left(\mathbf{t}\right)}$ parameter in minute, which represents the time interval from the center of the excess rainfall to the peak of hydrograph. ${\mathbf{L}\mathbf{a}\mathbf{g}}_{\left(\mathbf{t}\right)}$ is determined for each subbasin based on $\mathbf{T}\mathbf{c}$_, which is duration required for rainwater to travel from the most distant point in the watershed to the outlet. The $\mathbf{T}\mathbf{c}$ was calculated using Kirpich formula in equation 5, which incoperates factors such as river length, slope, and maximum retention potential of each subbasin [29].

$$\mathrm{T}\mathrm{c}=0.0195\ast {\mathrm{L}}^{0.77}\ast {\mathrm{S}}^{-0.385}$$

(5)

$\mathbf{T}\mathbf{c}$ is the time of concentration (in min), S is the average watershed slope, and L is the longest flow length (in meters) of a basin. The ${\mathbf{L}\mathbf{a}\mathbf{g}}_{\left(\mathbf{t}\right)}$ is calculated using equation 6.

$${\mathrm{L}\mathrm{a}\mathrm{g}}_{\left(\mathrm{t}\right)}=0.6\mathrm{T}\mathrm{c}$$

(6)

3.

Routing model

Muskingum method was used to calculate the river discharge at the outlet. This method developed in the 1930s is widely used in natural river channels and river engineering practices due to its simplicity and effectiveness [29,45]. It accounts for the gradual reduction of runoff as it moves along the river channel, primarily due to the channel storage effect. The method requires two key parameters, flood travel time ($\mathbf{K}$), measured in hours, and represents the time it takes for the flood wave to travel through the river reach. The value of $\mathbf{K}$ typically ranges from 0.1 to 150 hours, and is estimated as the ratio of the river length to the flow velocity. Secondly, the attenuation flood wave ($\mathbf{X}$) is a dimensionless parameter that reflects the influence of channel storage on the flood wave. The value of$\mathbf{X}$ ranges from 0.1 to 0.5, where a higher value indicates a steeper channel slope and a larger outflow volume. The basin storage was computed using equation 7.

$$\mathrm{S}=\mathrm{K}\left[\mathrm{X}\mathrm{I}+\left(1-\mathrm{X}\right)\mathrm{Q}\right]$$

(7)

In which $\mathbf{K}$ is flood wave traveling time (0 ≤ K ≤ 150), $\mathbf{X}$ is a weighting factor, $\mathbf{I}$ is inflow, $\mathbf{Q}$ is outflow, and $\mathbf{S}$ is storage.

2.4.3. Description of GEOGloWS Model

The third model, the GEOGloWS model is an open source, web-based software that access ECMWF forecast and historical discharge services (https://apps.geoglows.org/apps/geoglows-hydroviewer/). The model was created in 2017 by Brigham Young University and European Centre for Medium Range Weather Forecast (ECMWF). The model generates 15-day ensemble forecast and historical simulation since 1979 for watersheds having area greater than 150km² worldwide [18,33]. The model simulates discharge using the GEOGloWS ECMWF Streamflow hydroviewer, forced with the ECMWF Reanalysis version 5 (ERA-5) datasets [47].

In this study, the historical discharge for respective river basins was simulated using geographical coordinates and the GEOGloWS Reach ID. The Reach ID of the river basins were identified using the GEOGloWS hydroviewer. Any discrepancies between the coordinates and the Reach ID were corrected to ensure both represent the same river basin. Analysis was done in the google colab notebook environment. The GEOGloWS and Hydrostats python packages serve a client for the GEOGloWS model, facilitating programmatic access to the data service. We simulated historical flow and bias-corrected with observed discharge data at Chimakoti and Lungtenphu stations.

For bias correction, the model used the method developed by Farmer et al. [48] and Lozano et al. [18]. The model was developed into python package called GEOGloWS, which used geoglows.bias tool for effectively reducing bias, enhancing correlation and aligning flow variability. It uses a simulated flow duration curve to simulated streamflow data to a non-exceedance probability on the hydrograph. Then, using the observed flow duration curve, the non-exceedance probability estimated in the previous step is used to determine the corresponding observed streamflow. This observed streamflow was substituted into the sequence of simulated streamflows, resulting in a bias-corrected hydrograph.

2.5. Model Performance Evaluation Metrics

To evaluate performance of three hydrological models, numerical metrics namely NSE, PBIAS and PFE as recommended by Moriasi et al., [37] was employed. Briefly, firstly, NSE compares the variance of simulated data to the variance of observed data, indicating how well the plot of observed versus simulated data fits the 1:1 line [36]. The value ranges from $-\mathbf{\infty }$ to 1.0, with values above 0.6 indicating acceptable during calibration, while values exceeding 0.8 are considered excellent. For validation, NSE value exceeding 0.5 is considered acceptable and NSE value greater than 0.7 is considered very good [36,49]. NSE is computed following equation 8.

$$\mathrm{N}\mathrm{S}\mathrm{E}=1-{\displaystyle \frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{(\mathrm{Y}}_{\mathrm{i}}^{\mathrm{o}\mathrm{b}\mathrm{s}}-{\mathrm{Y}}_{\mathrm{i}}^{\mathrm{s}\mathrm{i}\mathrm{m}}{)}^2}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{(\mathrm{Y}}_{\mathrm{i}}^{\mathrm{o}\mathrm{b}\mathrm{s}}-{\mathrm{Y}}^{\mathrm{o}\mathrm{b}\mathrm{s},\mathrm{m}\mathrm{e}\mathrm{a}\mathrm{n}}{)}^2}}$$

(8)

where, ${\mathbf{Y}}_{\mathbf{i}}^{\mathbf{o}\mathbf{b}\mathbf{s}}$and ${\mathbf{Y}}_{\mathbf{i}}^{\mathbf{s}\mathbf{i}\mathbf{m}}$ represent the observed and simulated discharge data respectively, ${\mathbf{Y}}^{\mathbf{o}\mathbf{b}\mathbf{s},\mathbf{m}\mathbf{e}\mathbf{a}\mathbf{n}}$ is the mean of observed discharge data, $\mathbf{n}$ is the total number of observations, and $\mathit{i}$ represent the time series of the observed data.

Further, the models were evaluated using PBIAS to quantify river discharge errors [37]. It is the deviation of data being evaluated, expressed as a percentage, where its optimal value of 0.0 indicates accurate model representation. Positive PBIAS value indicates the underestimation of the simulated data whereas the negative value indicates the overestimation of simulated discharge. PBIAS values within the range of -15% and +15% is considered acceptable according to Rizwan et al. [49]. PBIAS is computed as shown in equation 9:

$$\mathrm{P}\mathrm{B}\mathrm{I}\mathrm{A}\mathrm{S}={\displaystyle \frac{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}\left({\mathrm{Y}}_{\mathrm{i}}^{\mathrm{o}\mathrm{b}\mathrm{s}}-{\mathrm{Y}}_{\mathrm{i}}^{\mathrm{s}\mathrm{i}\mathrm{m}}\right)\ast 100}{{\sum }_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{Y}}_{\mathrm{i}}^{\mathrm{o}\mathrm{b}\mathrm{s}}}}$$

(9)

The PFE, which is flood specific metric was used to measure the difference between simulated and observed peak flow. This metric measures the accuracy of the models in capturing peak flow during the flood events. The error closer to 0, indicates the model’s effectiveness in capturing the peak flow [38]. The PFE was calculated using equation 10.

$$\mathrm{H}\mathrm{F}\mathrm{E}={\displaystyle \frac{{(\mathrm{Y}}^{\mathrm{o}\mathrm{b}\mathrm{s},\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}-{\mathrm{Y}}^{\mathrm{s}\mathrm{i}\mathrm{m},\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}})}{{\mathrm{Y}}^{\mathrm{o}\mathrm{b}\mathrm{s},\mathrm{p}\mathrm{e}\mathrm{a}\mathrm{k}}}}\ast 100$$

(10)

In this equation, ${\mathbf{Y}}^{\mathbf{s}\mathbf{i}\mathbf{m},\mathbf{p}\mathbf{e}\mathbf{a}\mathbf{k}}$ and ${\mathit{Y}}^{\mathit{o}\mathit{b}\mathit{s},\mathit{p}\mathit{e}\mathit{a}\mathit{k}}$ represent simulated high flow and observed high flow, respectively.

3. Results

3.1. Models Performance Result

The simulated discharge was calibrated and validated at the outlets of Chimakoti (LRB) and Lungtenphu (SRB) stations. The default regional parameters failed to align the simulated discharge with observed values, necessitating calibration to improve accuracy. During validation, the IFAS model achieved NSE values of 0.87 and 0.84 at LRB and SRB, respectively (Figure 4), which are considered very good [49]. However, the NSE value during validation at LRB was reduced compared to calibration, whereas at SRB, the validation NSE exceeded the calibrated value.

The PBIAS remained within the acceptable range, with values of -8.51 for LRB and 13.04 for SRB. Despite this, the validated discharge was overestimated and underestimated in LRB and SRB, respectively compared to observed discharge. Additionally, the validated peak flow fell below observed peak at both the stations, indicating underperformance in capturing extreme events. The IFAS model demonstrated high efficiency in simulating river discharge, leveraging freely available satellite datasets to enhance performance.

The performance of the HEC-HMS model was lower than that of the IFAS model in both LRB and SRB. The NSE values obtained were 0.74 and 0.61 in LRB and SRB, respectively, which fall within the acceptable range but indicate weaker performance. The discharge was significantly underestimated with PBIAS values of 28.67 in LRB and 2.54 in SRB. Additionally, peak flows varied between stations; peak discharge was underestimated in LRB (PFE = 11.80) and overestimated in SRB (PFE = -15.58). The overall lower performance metrics can be attributed to the significant underestimation of low flows at both stations, affecting the model’s accuracy in simulating discharge dynamics (Figure 5).

The GEOGloWS model was bias corrected using observed discharge at Chimakoti and Lungtenphu stations. After bias correction, the model achieved an NSE value of 0.93 in LRB, indicating very good performance, whereas in SRB, the NSE value showed low performance with the value of 0.39 (Figure 4). This substantial variation in performance between two stations highlights the model’s inconsistent accuracy across different basin scales.

Despite this discrepancy, the PBIAS values were well within the acceptable range of 3.21 and 3.06 in LRB and SRB, respectively, although the discharge remained underestimated at both stations. In terms of PFE, the model underestimated peak discharge in LRB (PFE = 4.48) while overestimated SRB (PFE = -9.85). While the validated model effectively captured low flow conditions at both stations, it struggled with peak flow simulations, particularly at SRB. Overall, the GEOGloWS model demonstrated strong performance in SRB but delivered unsatisfactory results SRB.

3.2. Performance Based on Basin Size

Among the three models, the IFAS model estimated the highest discharge of 764.86m³/s, while the HEC-HMS model produced the lowest estimate of 40m³/s LRB. In SRB, the HEC-HMS model exhibited the wide range of discharge estimates, with a maximum of 115.5m³/s and minimum of 12m³/s during validation. The average NSE values of all models were 0.84 (very good) in LRB and 0.61 (acceptable) performance in SRB. However, the average discharge was consistently underestimated at both stations. Table S5 provides a detailed summary of descriptive statistics and model performance metrics.

At the individual basin level, the GEOGloWS model performed best in LRB achieving NSE value of 0.93, which reflects very good performance. Although the validated discharge was underestimated, both PBIAS and PFE remained within ±5%, indicating minimal bias. Meanwhile, the IFAS model demonstrated strong performance in SRB, with an NSE value of 0.84 significantly higher than the other models. The HEC-HMS model achieved acceptable NSE values, whereas the GEOGloWS model failed to meet the acceptable range in SRB. Despite this, the PBIAS and PFE values for all models were within the acceptable range in LRB.

Overall, the GEOGloWS model was best suited for flood simulation in LRB while the IFAS model performed more effectively in SRB. These findings highlight the importance of model selection based on basin characteristics to ensure reliable flood events simulation.

4. Discussion

The simulation of hydrological models remotely in a river basin saves time and resources especially during flood events. However, the selection of most appropriate hydrological model for river flow in varying eco-climatic conditions of the basin is a daunting and time consuming task [13]. Therefore, this study applied IFAS, HEC-HMS and GEOGloWS model for flood simulation in WRB. Our study reveals that the GEOGloWS model simulated river flow effectively in large river basins, while the IFAS model is effective in small river basins.

The calibrated parameters of IFAS model demonstrated a notable synchronization of simulated discharge patterns with observed discharge in both LRB and SRB. However, default parameters failed to accurately synchronize the hydrograph shape compared to observed discharge. This discrepancy might have arisen because the selected parameters were based on a regional study, necessitating further calibration at individual river basin. Similar findings were highlighted where default parameters failed to capture the high and low flow in the river [50]. The sensitive parameters that significantly impacted discharge simulation in this study are SNF, SKF, AUD and AGD. These parameters were also identified through sensitivity analysis [16] and calibration process [41,51,52], which played a crucial role in synchronizing simulated discharge with observed discharge. The IFAS model demonstrated higher NSE accuracy in LRB compared to SRB, which is consistent with previous study [24]. Smaller area will have courser resolution of satellite data, making it challenging to detect minute soil and land-use characteristics of the basin.

The combination of SCS CN, SCS UH, and Muskingum method showed results within the acceptable range. However, graphical observations and metrics illustrated that the HEC-HMS model does not performed well in both LRB and SRB. Further, slight variation was observed in the simulated discharge in recession limb during both calibration and validation phase. This discrepancy is attributed to runoff reaching the river channel quickly and sharp rise in peak discharge due to reduced infiltration and fast runoff. Similar results were observed in hilly river basin, where variations in recession and rising limb were observed [53]. Furthermore, the peak discharge was overestimated in small basin and underestimated in large basin during validation phase. The HEC-HMS model failed to capture the high flow which is an important performance indicator for event-based modeling. Therefore, this model can be applied in WRB, however, it needs further research on capturing the peak discharge of river flow.

In case of GEOGloWS model, there was large variation between LRB and SRB. The NSE value was below the acceptance value in small basin, while the NSE value was high in large basin. This is likely due to the use of ERA-5 rainfall with course resolution data which cannot capture the detailed rainfall patterns. Small basins are more sensitive to local climate and topographical features. In contrast, large basins benefit from the averaging effect of spatial variability and align better with global parameterizations. Additionally, the global bias-correction of ERA-5 rainfall may not address localized biases, further impacting small basin’s accuracy. Small basins model performance could be improved with higher resolution data, regional calibration, and better representation of localized processes.

The GEOGloWS model demonstrated more accurate simulation of low flow than high flow after bias correction in both LRB and SRB. This observation is consistent with the findings of Lozano et al. [18] in Dominican Republic, where significant improvements were observed in low flow simulation following bias correction. Hales et al., [33] also identified several problems of GEOGloWS model such as seasonally and spatially varying bias in flow magnitude and failure to capture short duration flood events. In this study, the flood peak has arrived before time and failed to capture the high flow in SRB.

Based on the findings of this study, it is recommended that future research should incorporate hourly rainfall and discharge data to refine the calibration of IFAS and HEC-HMS models, particularly to improve the accuracy in predicting the peak flows. The GEOGloWS model can be applied only in gauged basins for bias correction and validation of simulated discharge. This limitation highlights the critical importance of acquiring at least minimal observed data for one year. Alternatively, discharge data from other hydrological models can be substituted for bias correction. Although bias correction cannot be directly implemented in fully ungauged basins, collecting minimal observed discharge data is a more cost-effective solution than establishing long-term hydrological stations. The current study utilized PBIAS, NSE and PFE metrics to evaluate model performance, these measures do not capture the uncertainties associated with long-term simulations. Therefore, additional research is necessary to explore these uncertainties and develop more robust methods for long-term forecasting.

Furthermore, given the unique challenges posed by fast flowing mountainous rivers, this study underscores the need for expanded research across multiple small river basins within such environments. Investigating a broader range of basins would provide a more comprehensive understanding of model performance variability and enhance flood prediction and management strategies in these critical areas.

5. Conclusions

As the impact of climate change continues to manifest, particularly in regions like South Asia, the frequency and severity of river flood events are expected to increase. The use of software-based hydrological models for estimating discharge at different segments of a river becomes crucial to mitigate the growing threat. Despite this pressing need, Bhutan, like many developing and mountainous region, has seen limited progress in flood prediction through hydrological modeling. Recognizing this gap, this study addressed the need for improved flood management in Bhutan by evaluating models, namely IFAS, HEC-HMS and GEOGloWS. This study identified the most suitable hydrological model for simulating river discharge during flood events in Wangchu Basin using NSE, PBIAS, and PFE metrics.

Among three models, GEOGloWS exhibited commendable performance in simulating discharge in large basin, with NSE value of 0.93. Although the discharge was underestimated, the PBIAS value was 3.21, indicating minimal bias. Conversely, the IFAS model demonstrated higher efficiency in simulating discharge in the small basin, with NSE and PBIAS values of 0.84 and 13.04, respectively. However, peak discharge was underestimated, with PFE of 15.30. The HEC-HMS model underestimated discharge, showing inaccuracies in estimating low flows in large river basin and peak flows in small river basin.

The study concluded that the calibrated parameters performed well for flood simulation in Wangchu river basin, although further improvement is needed to increase model accuracy. These calibrated parameters hold potential for broader application in Wangchu river basin and similar geographical regions. The bias corrected GEOGloWS model demonstrated its usefulness for decision-makers, requiring fewer technical resources and expertise. This capacity lessens the financial and technical burdens on government. The majority of the work is manageable remotely using observed discharge and satellite datasets, offering a promising approach for flood simulation in the region.