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

Abstract

Flood has become a major hazard globally, and in Bhutan, with its steep terrain and erratic rainfall, it has caused significant economic damage in recent years. Given these challenges, there is a lack of accurate flood prediction and management strategies. In this study, therefore, we evaluated three hydrological models—Integrated Flood Analysis System (IFAS), Hydrologic Engineering Centre Hydrologic Modeling System (HEC-HMS), and Group on Earth Observation Global Water Sustainability (GEOGloWS)—and identified the most suitable model for simulating flood events in the Wangchu River Basin in Bhutan. Furthermore, we examined the models’ performance in a large and a small basin using the Nash–Sutcliffe Efficiency (NSE), Percent Bias (PBIAS), and Peak Flow Error (PFE) metrics. Overall, the GEOGloWS model demonstrated the highest accuracy in simulating flood in the large basin, achieving NSE, PBIAS, and PFE values of 0.93, 3.21%, and 4.48%, respectively. In the small basin, the IFAS model showed strong performance with an NSE value of 0.84. The GEOGloWS model provides simulated discharge but needs to be bias corrected before use. The calibrated parameters can be used in the IFAS and HEC-HMS models in future studies to simulate floods in the Wangchu River Basin and adjacent basins with similar geographical characteristics.

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 USD300 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 °C to 4 °C 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 a profound impact on lives, livelihoods, economies, 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-Saharan 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 river 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.

In recent 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 losses of lives, 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] identified the Wangchu River Basin as a highly vulnerable region in Bhutan, where river flow could increase by up to 60%, posing a severe threat to low-lying agricultural land.

Given the growing frequency and severity of high river floods, the 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 institution capacity [15,16]. However, in Bhutan’s mountainous terrain, structural measures such as the discharge stage and cable system become impractical during high flood events. In this context, flood discharge simulation can be approached through software-based hydrological models (a non-structural measure), which allow 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 the 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 the lack of tailored flood forecasting models [13], resulting in increased vulnerability and uncertainty in flood discharge 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) proving to be particularly advantageous for 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 (DEMs), and other geographical datasets makes it ideal for flood simulation in remote areas.

Similarly, the HEC-HMS model, although globally recognized, has been of 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 the river basins of Bhutan; however, more research is needed to explore its application in these challenging terrains.

Additionally, the development of the Group on Earth Observations Global Water Sustainability Version 1 (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, Bangladesh, Nepal, and Peru [18,34,35]. Despite its success elsewhere, it is yet to be applied in Bhutan, indicating the need for a comprehensive evaluation of each model for enhanced flood forecasting and risk management in this country.

While the use of hydrological models in large basins is widespread, the systematic evaluation of their relative efficiencies and advancement remains limited [14,16,22,29]. These models have primarily been employed for long-term discharge simulation, with fewer studies focusing on flood event simulation. Furthermore, the lack of comparative studies and data in small basins presents a considerable challenge, particularly in regions like Bhutan where small rivers and streams are poorly monitored. To address these gaps, in this study, we evaluated three hydrological models, namely, IFAS, HEC-HMS, and GEOGloWS, focusing on their suitability for simulating discharge during flood events in the Wangchu and Thimchu River Basins of Bhutan. Additionally, we examined whether model performance varied in the large basin (3556 km²) and small basin (658 km²) using key performance metrics such as the Nash–Sutcliffe Efficiency (NSE), Percent 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 examined the upper Wangchu River Basin (WRB) in Western Bhutan, a densely populated region with 240,012 inhabitants and two economically important hydropower plants downstream (Tala (1020 MW) and Chhukha (336 MW)) [11,36]. A comparative analysis was conducted in WRB (3556 km²) and the Thimchu River Basin (TRB, 658 km²), as illustrated in Figure 1. WRB is also known as the large basin while TRB is referred to as the small basin. Geographically, the area spans three distinct agro-ecological zones: alpine (3500–7500 masl), cool temperate (2600–3600 masl), and warm temperate (1800–2600 masl) [12]. The basin experiences an annual average temperature of 13.4 °C and an average rainfall of 575.1 mm/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 [36].

2.2. Research Flow

Figure 2 illustrates the research methodology for the three hydrological models: IFAS, HEC-HMS, and GEOGloWS. The IFAS and HEC-HMS models require extensive preprocessing of datasets from satellite and observation stations. These models also require a traditional approach of calibration and validation, with parameter setting. In this study, we did not use data from the same sources, such as DEM, land use, and soil data. The IFAS model is restricted to global datasets while the HEC-HMS model is flexible, with locally available datasets such as land use and soil data. The GEOGloWS model offers a more streamlined approach through its pre-developed package and web-based service. This model is freely available and provides immediate access through a user-friendly web service interface. It does not require calibration and parameter setting; however, it has to be bias corrected with the observed discharge data of the respective hydrological stations. The same discharge data were used for the validation of the IFAS and HEC-HMS models and the bias correction of the GEOGloWS model. The models’ performances were evaluated using three commonly employed performance metrics: NSE, PBIAS, and PFE [23,37,38].

2.3. Data Preparation

In the three hydrological models, we selected the observed river discharges from two distinct flood events for calibration and validation. To ensure that significant flood dynamics were incorporated, initially, we analyzed time series discharge data from 2012 to 2022 and identified the highest flood events within this period. Consequently, we selected the most significant event for calibration and validation in both the large and small basins. The calibration period spanned from 15 May to 14 June 2014, while the validation period covered 1 October to 31 October 2021. For both the IFAS and HEC-HMS models, we employed rainfall data corresponding to these same time frames, as rainfall data serve as a crucial input parameter that directly influences flow patterns within the catchment area. The observed discharge and rainfall data, archived by the National Center for Meteorology and Hydrology (NCHM), Bhutan, were used in these models. Although there are four hydrological stations in the study area (Figure 1), we used only two hydrological stations’ discharge data and six meteorological stations’ rainfall data in this study (

Table 1. Meteorological and hydrological station data used in this study.

Station Name	Station ID	Latitude	Longitude	Elevation (m)
Meteorological stations
Semtokha	12700046	27.44	89.42	2504
Begana	-	27.57	89.64	2520
Paro	12510046	27.38	89.68	2402
Drugyel Dzong	12580046	27.50	89.33	2467
Chapchu	12390046	27.07	89.55	2620
Haa	12510046	27.39	89.28	2726
Hydrological stations
Lungtenphu	12800045	27.45	89.66	2260
Chimakoti	12350073	27.11	89.53	2678

).

For the GEOGloWS model, a minimum of one-year daily observed data are required for the bias correction of simulated data at the respective hydrological stations [18,39]. Therefore, the observed discharge from 1 January to 31 December 2021 was used for the bias correction of the simulated discharge. The bias-corrected discharge was extracted from 1 to 31 October 2021, enabling direct comparison with the validation results of the IFAS and HEC-HMS models.

2.4. Description of the Selected Models and Parameter Setting

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 (accessed on 20 February 2023). It serves the rainfall–runoff simulation model designed for the efficient computation of discharge and flood forecasting, especially in ungauged river basins [4,40]. The model works based on the principle of a tank model for flood simulation and a kinematic wave model for routing [26,41]. The model needs input data on DEM, land use, soil, and ground rainfall (

Table 2. Source of satellite data for integrated flood analysis system model.

Dataset	Elevation	Land Use	Soil
Product name	GTOPO30	Land use land cover	Soil water-holding capacity
Resolution	30 arc second (1 km)	30 arc second (1 km mesh)	1 degree
Format	Raster (Tiles)	Raster (bil)	bil
Coordinate	WGS84	WGS84	90° N and 180° W
Coverage	World	World	World
Data source	https://earthexplorer.usgs.gov/ (accessed on 13 May 2023)	https://earthexplorer.usgs.gov/ (accessed on 15 May 2023)	https://www.fao.org/soils-portal/data-hub/soil-maps (accessed on 15 May 2023)

and Figure S1). Additionally, it features Geographic Information System functionality for the automatic generation of catchment boundaries, river channel networks, sub-basins, default baseline regional parameters, and the import of boundary river basins.

In this study, the DEM data were initially used to delineate river boundaries and the drainage network. DEM data of 1 km × 1 km resolution were downloaded, and the cell size was changed to 500 m using the automatic function in the IFAS model, as suggested by Lee and Kawata [16]. The reduction in cell size enhances the data’s capability to capture detailed information about the river basin. The two-layer tank model was selected as surface tank and aquifer tank in the vertical direction, and the river course tank on the right side of the model, as recommended in previous studies [22,24]. The surface tank was parameterized based on land use data and the aquifer tank based on the soil water-holding capacity. Output from both the surface tank and the aquifer tank provided the flow in the river course tank.

The surface tank model was used to divide the rainfall into surface flow, intermediate flow, and infiltration flows. The overland flow is calculated based on Manning’s equation, while ground infiltration is based on Darcy’s law [41]. Discharge from the land surface is responsible for increasing or decreasing peak discharge in the river. River discharge at the outlet was calculated from the river course tank model that was routed using the kinematic wave method. The input parameters for the river course tank model include the river width, Manning’s roughness coefficient, the constant coefficient, the initial water level in the river course, infiltration from the river tank to the aquifer tank, and the coefficient for cross-section shape [41]. The default parameters of the surface, aquifer, and river course tank model are presented in Tables S1–S3.

For the model setup, the coordinates of WRB were extracted using ArcGIS 10.8 and later assigned to the IFAS model to create the basin boundary. The daily time interval was selected for flood simulation. The downloaded DEM, land use, and soil data were imported into the model from the project information manager tool (Figure S1). The basin boundary, cell type, river course model, and sub-basins were created using DEM. The IFAS model has a function to automatically generate the river channel network based on elevation data. The observed rainfall data are prepared in .csv file format and forced into the IFAS model. Then, the IFAS model converts rainfall data using the inverse distance-weighted method to calculate the discharge. We incorporated all 28 of the IFAS model parameters during the simulation process.

During calibration, we tuned only five sensitive parameters associated with the surface, aquifer, and river course model based on the recommendations of previous studies [41,42]. The parameters include final infiltration capacity (SKF), surface roughness coefficient (SNF), runoff coefficient of unconfined aquifer (AUD), runoff coefficient of confined aquifer (AGD), and coefficient for cross-section shape (RLCOF). These identified sensitive parameters have also been suggested to apply in river basins of Bhutan with similar topography, as highlighted in reports and seminars of the NCHM, Bhutan.

2.4.2. Description and Processing of the 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 Advanced Land Observing Satellite Phased Array type L-band Synthetic Aperture Radar (ALOS PALSAR) DEM data to create the sub-basins and their characteristics. The DEM data had a 12.5 m resolution and were processed by the Alaska Satellite Facility. These data were downloaded from https://search.asf.alaska.edu/#/ (accessed on 16 August 2023). Soil data with a resolution of 250 m were downloaded from https://wocatapps.users.earthengine.app/view/dss-bhutan (accessed on 13 January 2024). The land use data of 2016, prepared by the National Land Commission of Bhutan, were used.

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

Table 3. Meta-models, methods, and parameters used in the HEC-HMS model.

Model	Method	Parameters	Functions
Loss	SCS-CN	Initial abstraction (mm)	The initial amount of rainfall lost due to surface storage, interception, and infiltration.
		Curve number	It reflects the basin’s runoff potential based on land use, soil type, and moisture conditions.
		Imperviousness (%)	The percentage of the basin area that is impervious, which contributes to surface runoff without infiltration.
Transform	SCS-UN	Lag time (min)	The time delay between the peak rainfall and peak runoff. It affects the shape of the hydrograph.
Routing	Muskingum	Travel time	The time it takes for water to move through a stream. It affects how quickly runoff reaches the outlet.
		Attenuation flood wave	It accounts for storage effects in the channel, leading to a flattened hydrograph and reduced peak discharge.

). Detailed explanations and processing approaches are provided below.

• Loss model

The SCS-CN method was used to estimate runoff from the catchment [43,44,45]. This method was chosen for its versatility and wide applicability in estimating surface runoff from each sub-basin, utilizing Equations (1)–(4) [45]. This method required data such as CNs, initial abstraction, and the percentage imperviousness of the basin. For this, CNs were prepared in ArcGIS 10.8 with the HEC-GeoHMS extension tool using land-use and hydrologic soil group (HSG). The raster data were converted to polygons and a union raster file was created. An 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]. Equations (1)–(7) compute parameters such as average CN, accumulated precipitation, maximum retention potential, time of concentration ($\mathrm{T}\mathrm{c}$), and lag time (${\mathrm{L}\mathrm{a}\mathrm{g}}_{\left(\mathrm{t}\right)}$). The curve number is calculated as:

$${\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 ${\mathrm{C}\mathrm{N}}_{\mathrm{a}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{a}\mathrm{g}\mathrm{e}}$ is the averaged CN, $\mathrm{i}$ is the number of the sub-basin, ${\mathrm{A}}_{\mathrm{i}}$ is the area of the particular sub-basin, and ${\mathrm{A}}_{\mathrm{T}\mathrm{o}\mathrm{t}\mathrm{a}\mathrm{l}}$ is the total area of the basin.

The accumulated excess precipitation is calculated as:

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

(2)

where ${\mathrm{P}}_{\mathrm{e}}$ is the accumulated precipitation excess, $\mathrm{P}$ is the unaccumulated rainfall depth, $\mathrm{I}\mathrm{a}$ is the initial abstraction, and $\mathrm{R}$ is the maximum retention potential.

$\mathrm{I}\mathrm{a}$ is calculated as:

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

(3)

And $\mathrm{R}$ is calculated as:

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

(4)

The HEC-HMS model used 10 sub-basins to simulate river discharge in both large and small basins (Figure 3). The basin has four HSGs, 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). Sub-basin 3 has the highest CN of 86 while sub-basin 8 has a CN of 60.

Table 4. Sub-basin characteristics of the HEC-HMS model.

Sub-Basins	Area (km2)	Length (km)	Centroidal Flow Path (km)	Slope (Degree)	Drainage Density (km/km2)	Curve Numbers	Lag Time (hour)
Sub-basin1	418	34	12	60.03	0.066	83	2.55
Sub-basin2	617	48	19	57.65	0.094	78	3.54
Sub-basin3	432	42	21	57.65	0.089	86	3.67
Sub-basin4	216	26	8	50.48	0.093	73	2.04
Sub-basin5	533	39	15	50.96	0.088	65	3.84
Sub-basin6	221	30	12	47.91	0.099	69	2.36
Sub-basin7	544	54	23	54.39	0.084	78	4.18
Sub-basin8	227	34	14	52.73	0.118	60	3.16
Sub-basin9	148	20	6	52.21	0.109	68	2.16
Sub-basin10	160	34	15	56.99	0.088	72	2.38

shows the sub-basin characteristics of the model. Sub-basin 2 covers the largest area of 617 km², while sub-basin 9 has the smallest area of 148 km². Sub-basin 7 has the longest flow path of 54 km. Drainage density falls under the low category, where it ranges from 0.066 km/km² to 0.118 km/km². The slopes of the sub-basins range from 47.91° to 60.03°, indicating a steep slope. The steeper the slope, the faster the surface runoff, and thus, it reaches the outlet more rapidly.

2.

Transform model

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

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

(5)

$\mathrm{T}\mathrm{c}$ is the time of concentration (in minutes), $\mathrm{S}$ is the average watershed slope, and $\mathrm{L}$ is the longest flow length (in meters) of a basin. ${\mathrm{L}\mathrm{a}\mathrm{g}}_{\left(\mathrm{t}\right)}$ is calculated as:

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

(6)

3.

Routing model

The 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,47]. It accounts for the gradual reduction in runoff as it moves along the river channel, primarily due to the channel storage effect. The method requires two key parameters. The flood travel time ($\mathrm{K}$), measured in hours, represents the time it takes for the flood wave to travel through the river reach. The $\mathrm{K}$ value ranges from 0.1 to 150 h and is estimated as the ratio of the river length to the flow velocity.

Secondly, the attenuation flood wave ($\mathrm{X}$) is a dimensionless parameter that reflects the influence of channel storage on the flood wave. The parameter $\mathrm{X}$ is a weighting factor that defines the non-linearity of the routing process, influencing how the flow is distributed over time. Typically, $\mathrm{X}$ values range between 0 and 0.5, where 0 represents pure translation and 0.5 represents maximum diffusion. The basin storage was computed as:

$$\mathrm{S}=\mathrm{K}\left[\mathrm{X}{\mathrm{Q}}_{\mathrm{i}\mathrm{n}}+\left(1-\mathrm{X}\right){\mathrm{Q}}_{\mathrm{o}\mathrm{u}\mathrm{t}}\right]$$

(7)

where $\mathrm{K}$ is the flood wave traveling time (0 ≤ $\mathrm{K}$ ≤ 150), $\mathrm{X}$ is a weighting factor, ${\mathrm{Q}}_{\mathrm{i}\mathrm{n}}$ is the inflow, ${\mathrm{Q}}_{\mathrm{o}\mathrm{u}\mathrm{t}}$ is the outflow, and $\mathrm{S}$ is the storage.

2.4.3. Description of GEOGloWS MODEL

The third model, GEOGloWS, is an open-source, web-based software that accesses the European Centre for Medium Range Weather Forecast (ECMWF)’s forecast and historical discharge services (https://apps.geoglows.org/apps/geoglows-hydroviewer/, accessed on 13 June 2023). This model was developed in 2017 by ECMWF [33]. It has been generating 15-day ensemble forecasts and historical simulations since 1979 for watersheds with an area greater than 150 km² worldwide [18,34]. The model simulates discharge using the GEOGloWS ECMWF Streamflow hydroviewer, forced with the ECMWF Reanalysis version 5 (ERA-5) datasets [48].

The GEOGloWS model needs two input datasets such as runoff estimates from the Hydrology-Tiled ECMWF Scheme for Surface Exchanges over Land (HTESSEL) and the delineated network of watershed boundaries and streamlines. The runoff in each sub-basin is computed as the sum of each cell’s contributing area (km²) multiplied by its forecasted runoff with the help of Equation (8) [35]. This computation is repeated for all the sub-basins at each time step of the runoff data to produce runoff. The model calculates a cumulative runoff volume at each time step as an incremental contribution from each basin using the Muskingum routing method (35).

$$\mathrm{V}\mathrm{o}\mathrm{l}\mathrm{u}\mathrm{m}\mathrm{e}(\mathrm{r}\mathrm{u}\mathrm{o}\mathrm{f}\mathrm{f})=\sum _{\mathrm{c}\mathrm{e}\mathrm{l}\mathrm{l}\mathrm{s}}\mathrm{A}\mathrm{r}\mathrm{e}\mathrm{a}\left(\mathrm{i}\right)\ast \mathrm{R}\mathrm{u}\mathrm{n}\mathrm{o}\mathrm{f}\mathrm{f}\left(\mathrm{i}\right)$$

(8)

where $\mathrm{A}\mathrm{r}\mathrm{e}\mathrm{a}\left(\mathrm{i}\right)$ is the area of the sub-basin and $\mathrm{R}\mathrm{u}\mathrm{n}\mathrm{o}\mathrm{f}\mathrm{f}\left(\mathrm{i}\right)$ is the runoff from the sub-basin.

In this study, the historical discharge for the respective river basins was simulated using geographical coordinates and the GEOGloWS Reach ID. The Reach IDs of the river basins were identified using the GEOGloWS hydroviewer. Any discrepancies between the coordinates and the Reach IDs were corrected to ensure that both represent the same river basin. Analysis was performed in the Google Colab notebook environment. The GEOGloWS version 0.27.1 and Hydrostats Python packages serve as clients for the GEOGloWS model, facilitating programmatic access to the data service. We simulated daily historical flow and bias corrected it with the observed discharge data at the Chimakoti and Lungtenphu stations.

For bias correction, the model used the method developed by Farmer et al. [49] and Lozano et al. [18]. The model was developed into a Python package called GEOGloWS, which used the geoglows.bias tool to effectively reduce bias, enhance correlation, and align flow variability. It used a simulated flow duration curve to simulate 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 was used to determine the corresponding observed streamflow. This observed streamflow was substituted into the sequence of simulated streamflows, resulting in a bias-corrected hydrograph. The model does not have to be calibrated separately, and it works based on the bias correction technique implemented in the GEOGloWS package.

All analyses of the GEOGloWS model were performed in the Google Colab environment. The discharge from the model was downloaded using the GEOGloWS Representational State Transfer Programming Interface (REST API). The model simulated and bias corrected all the discharge data since 1979. However, for the comparison with the IFAS and HEC-HMS models, we extracted the required timeframe datasets for the calculation of metrics and analysis.

2.5. Model Performance Evaluation Metrics

To evaluate the performance of the three hydrological models, numerical metrics, namely, NSE, PBIAS, and PFE, as recommended by Moriasi et al. [38] and Chen et al. [23], were employed. Briefly, first, 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 [37]. The value ranges from $-\mathbf{\infty }$ to 1.0, with values above 0.6 being acceptable during calibration, while values exceeding 0.8 are considered excellent. For validation, NSE values exceeding 0.5 are considered acceptable and NSE values greater than 0.7 are considered very good [37,50]. NSE is computed as:

$$\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}}$$

(9)

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

Further, the models were evaluated using PBIAS to quantify river discharge errors [38]. PBIAS is the deviation of data being evaluated, expressed as a percentage, where its optimal value of 0.0 indicates accurate model representation. A positive PBIAS value indicates an underestimation of the simulated data whereas a negative value indicates an overestimation of simulated discharge. PBIAS values within the range of −15% and +15% are considered acceptable, according to Rizwan et al. [50]. PBIAS is computed as:

$$\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}}}}$$

(10)

PFE, which is a flood-specific metric, was used to measure the difference between the simulated and observed peak flow. This metric measures the accuracy of the models in capturing peak flow during flood events. An error closer to 0 indicates the model’s effectiveness in capturing the peak flow [23]. The PFE was calculated using Equation (11).

$$\mathrm{P}\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$$

(11)

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

3. Results

3.1. Results of Model Performance

The calibrated parameters of the IFAS and HEC-HMS models are presented in Table 4 and

Table 5. Calibrated parameters of IFAS model in WRB.

Tank Model	Parameters	Notations	Units	Default Parameter	Calibrated Parameters
Surface tank	Final infiltration capacity	SKF	cm/s	0.0005	0.005
				0.00002	0.002
				0.00001	0.001
				0.000001	0.0001
				0.00001	0.0001
	Surface roughness coefficient	SNF	m−1/3/s	0.70	1.50
				2.00	2.00
				2.00	2.00
				0.10	1.00
				2.00	2.00
Aquifer tank	Runoff coefficient of unconfined aquifer	AUD	(1/mm/day)1/2	0.10	0.02
	Runoff coefficient of confined aquifer	AGD	1/day	0.003	0.001
River course tank	Coefficient for cross-section shape	RLCOF	Non-dimensional	1.4	1.4

. These parameters were used for the validation of the models. In the case of the GEOGloWS model, the bias-corrected discharge was used for the comparison with the IFAS and HEC-HMS models. The simulated discharge was calibrated and validated at the outlets of the Chimakoti (large basin) and Lungtenphu (small basin) stations. The details of the performance metrics and descriptive statistics of the models are presented in Figure 4 and Table S5.

3.1.1. IFAS Model

The sensitive parameters of such final infiltration capacity and surface roughness coefficient values were increased in this model to capture the high flow. The runoff coefficient of the confined aquifer and the runoff coefficient of the unconfined aquifer parameter values were reduced to capture low flow (Table 5). The parameters were calibrated carefully to reproduce flow dynamics within the acceptable range. The coefficient for the cross-section shape parameter was not calibrated because the simulated high flow aligned with the observed flow. Overall, the simulation with calibrated parameters demonstrated good synchronization with the observed river discharge.

For the IFAS model, the default regional parameters failed to align the simulated discharge with the observed discharge, necessitating calibration to improve accuracy. During calibration, the IFAS model achieved an NSE value of 0.88 in the large basin and 0.76 in the small basin (Table S5). However, during validation, the model obtained NSE values of 0.87 and 0.84 in the large and small basins, respectively (Figure 4), which are considered very good [50]. The NSE value during validation in the large basin was decreased compared to calibration, whereas the NSE value exceeded the calibrated value in the small basin.

PBIAS, which calculates the overestimation and underestimation of river discharge, remained within the acceptable range during both calibration and validation in the large basin. The model obtained PBIAS values of −8.51% in the large basin and 13.04% in the small basin during validation. Despite this, the validated discharge was overestimated and underestimated in the large and small basins, respectively, compared to the observed discharge. Additionally, the validated peak flow was below the observed peak in both 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.

3.1.2. HEC-HMS Model

Table 6. Calibrated parameters used in the HEC-HMS model.

Parameters	Calibrated Values
Initial abstraction (mm)	0.2
2. Imperviousness (%)	1.9
3. Travel time (K)	7
4. Flood wave (X)	0.1

presents the calibrated parameters used in this model. The initial abstraction for the study area was 0.2 mm, meaning that 20% of the rainfall is retained by the forest canopy and does not contribute to runoff. The imperviousness of the area was 1.9. The K and X values were calibrated at 7 and 0.1, respectively. The K value indicates that the basin has a steep topography, which results in a shorter time for floodwaters to travel. These calibrated parameters were used in the synchronization of the simulated discharge with the observed discharge at hydrological stations.

The performance of the HEC-HMS model was lower than that of the IFAS model in both the large and small basins. Compared to the validated NSE values, the calibrated NSE values were high in the large basin and low in the small basin (Table S5). During validation, the NSE values obtained were 0.74 and 0.61 in the large and small basins, respectively, which falls within the acceptable range but indicates weaker performance compared to other models. The discharge was significantly underestimated, with PBIAS values of 28.67% in the large basin and 12.54% in the small basin. Additionally, peak flows varied between stations; peak discharge was underestimated in the large basin (PFE = 11.80%) and overestimated in the small basin (PFE = −15.58%).

The overall lower performance metrics can be attributed to the significant underestimation of low flows in both stations, affecting the model’s accuracy in simulating discharge dynamics (Figure 5). This could be attributed to the general use of an initial abstraction value of 0.2 mm rainfall and the poor estimation of curve numbers. A lower curve number will decrease the surface runoff and increase infiltration, leading to the underestimation of low flows. As illustrated in Figure 3d, most of the area is covered by forests and shrubs, alongside a few patches of built-up areas.

3.1.3. GEOGloWS Model

The simulated discharge from the GEOGloWS model was bias corrected using the observed discharge in both the large and small basins. The bias-corrected result was used for the comparison with other models. All the performance metrics were out of range during the simulation of the model in both basins. However, after bias correction, the model had an NSE value of 0.93 in the large basin, indicating very good performance, whereas in the small basin, the NSE value showed low performance with a value of 0.39 (Figure 4). During simulation, the discharge was overestimated, with PBIAS values of −59.06% and −27.65% in both the large and small river basins, respectively. Despite this discrepancy, the PBIAS values were well within the acceptable range, with PBIAS values of 3.21% in the large basin and 3.06% in the small basin after bias correction. The overall discharge was overestimated during simulation and underestimated after bias correction in both stations.

In terms of PFE, the model underestimated peak discharge in the large basin (PFE = 4.48%) while overestimating it in the small basin (PFE = −9.85%). The bias-corrected discharge effectively captured low flow conditions in both stations; however, it struggled with peak flow correction, particularly in the small basin. Overall, the GEOGloWS model demonstrated excellent performance in the large basin but delivered unsatisfactory results in the small basin. The significant variation in performance between the large and small basins underscores the model’s inconsistent accuracy across varying basin scales and the challenges of applying a global model at a local scale.

3.2. Performance Based on Basin Size

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

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

Overall, the GEOGloWS model was best suited for flood simulation in the large basin, while the IFAS model performed more effectively in the small basin. These findings highlight the importance of model selection based on basin characteristics to ensure reliable flood event 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 the most appropriate hydrological model for river flow in the varying eco-climatic conditions of the basin is a daunting and time-consuming task [13]. Therefore, this study applied the IFAS, HEC-HMS, and GEOGloWS models for flood simulation in WRB. Our study revealed that the GEOGloWS model simulated river flow effectively in the large river basin, while the IFAS model was effective in the small river basin.

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

The combination of SCS CN, SCS UH, and the Muskingum method showed results within the acceptable range. However, the graphical observations and metrics illustrated that the HEC-HMS model did not perform well in both the large and the small basins. Furthermore, a slight variation was observed in the recession limb of the simulated discharge during both the calibration and validation phases. This discrepancy is attributed to runoff reaching the river channel quickly and the sharp rise in peak discharge due to reduced infiltration and fast runoff. Similar results were observed in hilly river basins, where variations in recession and rising limbs were observed [54]. Furthermore, the peak discharge was overestimated in the small basin and underestimated in the large basin during the 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, further research is needed on capturing the peak discharge of the river flow.

In the case of the GEOGloWS model, there was a large variation between the large and small basins. The NSE value was below the acceptance value in the small basin, whereas it was high in the large basin. This is likely due to the use of ERA-5 rainfall with coarse-resolution data, which cannot capture 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 the accuracy in small basins. The small basin model performance could be improved with higher resolution data, regional calibration, and a better representation of localized processes.

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

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 the IFAS and HEC-HMS models, particularly to improve their accuracy in predicting peak flows. Additionally, the GEOGloWS model could be more effective if it were tailored for simulating hourly discharge data. The GEOGloWS model can be applied only in gauged basins for the bias correction of simulated discharge. This limitation highlights the critical importance of acquiring at least one year of observed data for bias correction. Alternatively, discharge data from other hydrological models can be substituted for bias correction. The current study utilized PBIAS, NSE, and PFE metrics to evaluate model performance, but 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.

The calibrated parameters from the IFAS and HEC-HMS models can be effectively applied to ungauged basins in Bhutan due to the country’s consistent mountainous topography. The hydrological characteristics across Bhutan’s river basins exhibit similar responses to precipitation events, making the transferability of calibrated parameters feasible. This allows for reliable model performance in ungauged basins when using parameters derived from gauged basins with comparable catchment features. Consequently, IFAS and HEC-HMS can serve as valuable tools for hydrological analysis and flood forecasting in data-scarce regions within Bhutan. In the case of the GEOGloWS model, bias correction cannot be directly implemented in fully ungauged basins. Therefore, collecting a minimum of one year of observed discharge data is a more cost-effective solution than establishing long-term hydrological stations for flood discharge simulation.

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 regions, 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 three such models, namely, IFAS, HEC-HMS, and GEOGloWS. We then identified the most suitable hydrological model for simulating river discharge during flood events in the Wangchu River Basin using NSE, PBIAS, and PFE metrics.

Among the three models, GEOGloWS exhibited commendable performance in simulating discharge in the large basin, with an NSE value of 0.93. Although the discharge was underestimated, the PBIAS value was 3.21%, indicating minimum 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 a PFE of 15.30%. The HEC-HMS model underestimated discharge, showing inaccuracies in estimating low flows in the large basin and peak flows in the small basin.

We concluded that the calibrated parameters performed well for flood simulation in the Wangchu River Basin, although further improvement is needed to increase model accuracy. These calibrated parameters hold potential for broader application in the 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 governments. Most of the work is manageable remotely using observed discharge and satellite datasets, offering a promising approach for river flood discharge simulation in the region.