Application of hydrological model to assess river flow in the transboundary cryosphere and data-scarce watershed, a case study: Chitral-Kabul River Basin (C-KRB) in Pakistan

Water scarcity can affect the ecosystem and its functions as well as the economy of countries like Pakistan and Afghanistan as intense precipitation leads to severe hydrologic consequences such as massive flooding and erosion (Solomon et al. 2007), while scarcity leads to drought phenomena and eventually challenges to food security and livelihood. There are 145 countries which are riparian to one or more of 263 international basins across the world (Wolf 2002). The western rivers of the Indus River Basin originate from Afghanistan i.e., Kabul River. The Kabul River Basin (KRB) is a wider second major supplementary basin to the Indus River (Thomas et al. 2016). The Chitral-Kabul River Basin (C-KRB) with an area of 14,782 km² is a major water supplier to the KRB and is located across the Durand-line (border between Pakistan and Afghanistan) with no bilateral agreements so far over the shared water resources for water allocation and utilization between Pakistan and Afghanistan (Yıldız 2015); for instance the Indus Waters Treaty in 1960 (Pakistan–India) (Ali & Bhargava 2021) and Helmand River Treaty in 1973 (Afghanistan–Iran) (Goes et al. 2016).

Water treaties, negotiations and dialogues require precise data and information. However, due to the widespread conquered insecurity and scarcity of data in the area, the assessment of the river flow of the Kabul River is extremely challenging. At the same time, the development of future legal declarations and joint management structures based on common interests is very important. The C-KRB flow is very important for Pakistan due to water shortage, which is expected to decrease to 800 m³ per capita in the year 2025 (WWF Pakistan 2007). The average temperature of Pakistan increased by 0.6 ^oC from 1901 to 2000 and caused 15% decrease in glacier area (Shakir et al. 2010). Global warming, climate change and drastic population growth may affect the contribution of the C-KRB to the downstream Kabul River (Afghanistan). Being already underlined as a physically water-scarce country, it is deemed that Afghanistan will get further exposed to extreme weather severities and droughts with eventual threats to food security and a renewable water resources availability which is less than the threshold (i.e. 1,500 m³/capita/year) (Yang et al. 2003; Ahmad & Wasiq 2004). Therefore, despite the challenges of the limited number of hydro-meteorological stations and scarcity of data, it is very important to predict the peak and cumulative flow of the C-KRB.

Hydrological modeling approaches are advanced techniques of a wide range of applications in water resources planning, development and management. Different hydrological models such as Inverse Laplace Transform Diffusive Flood Routing model (ILTDFR) (Cimorelli et al. 2015), Soil and Water Assessment Tool (SWAT) model (Rostamian et al. 2008), University of British Columbia (UBC) model (Loukas & Vasiliades 2014) and Snowmelt Runoff Model (SRM) (Tekeli et al. 2005) have been used effectively for nearly real-time flow monitoring, and flood prediction with respect to its applicability, limitations and aerial or ground data availability. Herein, we have applied the HEC-HMS model, which has a wide range of application in both cryosphere (snow-covered) watersheds and no snow-cover watersheds (Yusop et al. 2007; Stojković & Jaćimović 2016; Gao et al. 2017; Bouragba et al. 2019). The HEC-HMS model can be applied in both event-based (Jin et al. 2015; Shahid et al. 2017) and continuous-based modeling (Chu & Steinman 2009; Azmat et al. 2017), commutating runoff volume, direct runoff, base and channel flow which includes watershed precipitation and evaporation methods for surface runoff simulations.

Advancements in remote-sensing technologies have enabled the acquisition of spatio-temporal hydro-climatic data for better understanding the hydrological processes and minimizing the uncertainties in the predictions process (Ogden et al. 2001; Emerson et al. 2005). Presently, the satellite or remotely sensed data is rather easily available which can be used in hydrological studies and management practices especially in ungauged basins. Satellite-based precipitation estimates are the only information source for poorly gauged watersheds in the tropics where extreme rainfall events cause heavy floods every year (Hong et al. 2010). Tropical Rainfall Measurement Mission (TRMM)-based rainfall retrieval is supposed to be more reliable than any other remotely sensed product (Nicholson 2005; Collischonn et al. 2008). TRMM-based rainfall estimations are usually well correlated with the ground-based measurements (Collischonn et al. 2008; Nair et al. 2009; Duncan & Biggs 2012; Islam et al. 2012). Both sampling errors and temporal errors may result in erroneous applications if applied without local calibration (AghaKouchak et al. 2009; Gebremichael et al. 2010). Therefore, TRMM-data was calibrated using available ground data and then used as input to the HEC-HMS model. Unavailability of reliable spatio-temporal variations in snow cover in the C-KRB was due to its complex terrain features and limited measuring stations. Therefore, Moderate Resolution Imaging Spectrometer (MODIS) satellite data was organized and used, as it is a useful data source for mapping regional and basin-scale snow-cover area (SCA) (Liu et al. 2017; Liu et al. 2018).

In the present study, application of the HEC-HMS model for runoff assessment in the snow-covered, glacier-developed and highly mountainous Chitral River watershed was attempted, which has been used by many researchers for various watersheds across the world (Yimer et al. 2009; Gyawali & Watkins 2013; Munyaneza et al. 2014; Mahmood & Jia 2016; Bhuiyan et al. 2017; Koneti et al. 2018; Khasmakhi et al. 2020). The model was calibrated and validated for the period of 2010–2011 (66% of the data) and 2012 (33% of the data) respectively. The main objective of this research was to develop and evaluate the HEC-HMS hydrologic model for the assessment of stream flow using free and easily accessible remotely sensed data in a cryosphere (snow-covered) watershed. The stream flow data has a wide range of application for agriculture planning and management practices, natural disaster management and first-hand information for bilateral agreement between the shared countries. It will also provide input stream data for other hydrological models (such as HEC-RAS) to develop a flood hazard map for the downstream C-KRB in Afghanistan.

This research was conducted on a major tributary of Kabul River originating in Pakistan which is interchangeably called the Kunar River in Afghanistan and Chitral River in Pakistan (Figure 1). It originates in the far north glaciated Hindu Kush mountains, approximately 7,683m above the sea level, and comprises the isolated areas of the western end of the Himalayas in Chitral district, area-wise the first largest district of Khyber Pakhtunkhwa (KP) province of Pakistan with a population of 318,689. Geographically, it is located at 35°17'54.49″N to 36°54'54.72″N latitude and 71°11'43.50″E to 73°53′24.50″E longitude. The catchment of the river covers an area of about 14,782 km².

Chitral River Basin is very important for the agriculture sector where the source of irrigation water originates from the snow melting over the mountainous peaks of Chitral. The Chitral River catchment has a dry Mediterranean climate according to Köppen classifications with almost no rainfall during the very hot summers (May to June). The precipitation decreases from southwest to northeast e.g., the mean annual rainfall is more than 600mm at Drosh (downstream of Chitral city) attributed to more intensive forest distribution than other regions of the C-KRB to approximately 450mm at the Chitral. Influences of the marginal monsoon are a source of heavy precipitation in summer (May to October) but it varies significantly in the winter season (November to April) year by year (Dahri et al. 2018). The maximum temperature recorded in Chitral meteorological station is 30 °C for the month of July while the minimum temperature recorded so far is −3.8 °C and −0.9 °C for the months of January and February, respectively. In winter, the night temperature occasionally drops to −10 °C and snowfall in winter heavily accumulates up to half a metre and at higher elevation snow-depth sometimes reaches up to six metres. The mean annual temperature and precipitation are 15.6 °C and 418mm, respectively.

Meteorological data (maximum temperature, minimum temperature and precipitation) at three different stations in Chitral district were collected over 2010–2012 from Pakistan Meteorological Department (PMD). These stations are located at Chitral city (1,497.8m AMSL, longitude 71° 50′ E, latitude 35° 51′ N), Drosh (1,463.9m AMSL, longitude 71° 47′ E, latitude 35° 34′ N) and Mir Khani (1,250m AMSL, longitude 74° 42′ E, latitude 35° 30′ N) as shown in Figure 2. These data were used for calculating different initial input parameters for the HEC-HMS model, such as lapse rate for temperature data interpolation at the above-mentioned stations in the region. Daily stream flow gauge data at two-gauge stations of the C-KRB, i.e., Chitral city and Arandu gauge station (outlet of the basin) was provided by the Water and Power Development Authority (WAPDA) Lahore, Pakistan.

In this study, the TRMM satellite rainfall product 3B43 version 7 with monthly temporal and spatial (0.25°×0.25°) resolution was used and calibrated against the above-mentioned meteorological stations to check the accuracy of TRMM rainfall data in the C-KRB from 1 January 2000 to 31 December 2016. The TRMM satellite daily product 3B42 was also used for snow-water equivalent calculation. The hypothetical area averaged gauge-station was considered for each of the 11 sub-basins (constant precipitation over the entire sub-basin) of the entire C-KRB during downloading of the precipitation data from the Giovanni website [https://giovanni.gsfc.nasa.gov/giovanni/].

Therefore, TRMM three-hourly product 3B43 for the period of 1 January 2010 to 31 December 2012 were downloaded, projected to WGS 1984 UTM ZONE 43 and used in the HEC-HMS model for each sub-basin in which constant precipitation was considered over the entire sub-basin.

Land surface elevation is digitally represented by DEM projected to WGS 1984, UTM ZONE 43. The DEM (30m spatial resolution), generated from data provided by the US Geological Survey, was used in the HEC-HMS model for extraction of terrain features such as watershed delineation or drainage basin, sub-basins and river channel networks. The elevation of the C-KRB varies from 1,051 to 7,683m above sea level (AMSL) (Figure 3).

Soil data represents one of the main parameters in hydrological modeling and specifically it is used for generating the runoff curve number (CN) of each sub-basin of the study area. Moreover, it can also be used for finding the soil water holding capacity, hydraulic conductivity and soil textural information, etc. The HEC-HMS model requires soil information such as soil texture and physico-chemical properties. The World Digital Soil Map developed by the Food and Agriculture Organization (FAO) was used to show the soil texture in the study regions (source: http://www.fao.org/soils-portal/soil-survey/soil-maps-and-databases/faounesco-soil-map-of-the-world/en/). The soil map showed that sandy clay loam is the main dominant soil type (Figure 4).

Determination of the curve number (CN) values of each watershed was estimated using the land-cover data (Globe Cover-2009) (Defourny et al. 2009) and other parameters such as surface runoff, erosion and evapotranspiration. The union of land-cover and soil map was used for computing the CN values and initial values for HEC-HMS parameters such as for the SCS (Soil Conservation Service) Curve Number loss method and time of concentration and storage coefficient for transform method. Different land-cover classes of the study area are shown in Figure 5.

In cryosphere watershed (snow/ice covered) modeling applications, calibration of the temperature index snowmelt method is recommended (Daly et al. 1999). The temperature index snowmelt method in HEC-HMS and LBRM has recently been applied to the Upper Euphrates River Basin in Turkey (Yilmaz et al. 2012). Hence, snow cover data was utilized to obtain different temperature index snowmelt method input parameters such as melt rate and snow water equivalent. There were limitations of the snow water equivalent (SWE) data. Recently, snow cover information has also been proposed by Pistocchi et al. (2017) for obtaining snow water equivalent and snowmelt rate. Therefore, MODIS images from the Terra satellite were used for extracting the snow cover information. MOD10A2 product (https://nsidc.org/data/modis/data_summaries) (Hall et al. 2006) of Snow Cover 8-Day L3 Global 500m SIN Grid V005 images were downloaded covering a temporal range (1 January 2010, 31 December 2012) of eight-day snow cover and projected to WGS 1984, UTM Zone 43.

The snowmelt that occurs from beneath the snowpack is defined as the ground melt. For relatively shallow seasonal snow cover (SWE<12in.), it is set to zero (US Army corps of Engineers, Hanover, NH, USA, unpublished internal report, 2008).

The HEC-HMS is a numerical model and a computer program that includes a large set of methods to simulate watershed, water channel and water control structure behavior, and predict flow, timing and storage. The computation of runoff volume, direct runoff, base and channel flow, watershed precipitation and evaporation are included in the simulation methods. The HEC-HMS model can be applied for both event-based simulation (Jin et al. 2015; Shahid et al. 2017) and continuous-based simulations, as used by many researchers (Chu & Steinman 2009; Azmat et al. 2017). For performing data analysis in modeling, the following software is implemented.

The Geospatial Hydrologic Modeling System (HEC-Geo-HMS) is public domain extension software programmed as executable codes developed with a Cooperative Research and Development Agreement by the Hydrologic Engineering Center (HEC) and Environmental System Research Institute, Inc. (ESRI). It is supported by United States Army Corps of Engineers (USACE) Research Development funding. The HEC-Geo-HMS visual interface allows the user to visualize spatial information, document watershed characteristics, perform spatial analysis, delineate sub-basins and river networks, assist with report preparation and create input parameter files of the HEC-HMS model (Fleming & Doan 2013). It was developed as a geospatial hydrology toolkit used with ArcGIS for hydrologists and engineers with limited GIS experience. The preprocessing of DEM was done through GIS techniques. It provides the hydrological boundary of the catchment. Different characteristics of the catchment, such as length of river, slope of river, longest flow path and the basin's centroid, were estimated from the characteristics toolbar. The HEC-HMS project was created through the HMS pull bar menu of the HEC-Geo-HMS toolbar. The generated project then was imported in the HEC-HMS (3.5 version) for further processing and utilized for simulations.

The conceptual framework (Figure 6) represents the work procedures to develop the HEC-HMS model for three-hour-based rainfall–runoff simulation of the C-KRB at Arandu hydrometric station (outlet of the C-KRB). ArcGIS and the HEC-HMS model were used as the main processing engines for simulating flow. The whole process was carried out in five steps: data collection, initial input parameter estimations, HEC-HMS model processing, calibration and validation.

Delineation of rivers and their sub-basins are the substantial initial inputs required in model set-up. The DEM was used to divide the entire watershed into several sub-basins, the topography, contour and slope of which are derived using HEC-GeoHMS tools, such as, fill sinks, flow direction, flow accumulation, stream definition, stream segmentation, catchment grid delineation and catchment polygon processing etc. Through the automatic delineation of the entire watershed (Figure 7), 11 sub-basins were created in the C-KRB namely as W1, W2, W3, W4, W5, W6, W7, W8, W9, W10 and W11, which represent 21.21%, 12.46%, 8.83%, 2.96%, 8.64%, 2.16%, 10.81%, 6.89%, 7.27%, 2.52% and 16.25% of the area of the entire C-KRB respectively.

CN is a non-dimensional hydrological parameter which varies from 0 to 100 to predict the response runoff or infiltration from excess rainfall. Consequently, the Soil Conservation Services–Curve Number (SCS-CN) method has been used on small basins for long-term hydrological simulation and several models have been developed in the past two decades. It depends upon the soil type, land use, hydrological condition and antecedent moisture condition. In this study, an SCS-CN grid was developed with soil and land-cover data of the study area, following the proposed methodology of Merwade (2012) using HEC-GeoHMS and GIS. CNs were generated based on preparation of land-cover data through downloading the Global Land Cover Map (Globe Cover 2009 version 2.3) developed in 2010 by the European Space Agency (ESA) (Defourny et al. 2009). Each grid of the data has a unique value or number representing a specific land-cover class defined by the United States Geological Survey (USGS), Land Cover Institute (LCI). The 15 different land-cover classes were reclassified into four broad classes illustrated in Figure 5, namely, water body and wetlands, urban/residential areas, and forest and agricultural lands as described by Shahid (2015); but the C-KRB has negligible urban areas in the map extracted (Figure 5) from the Globe Cover map and hence only three classes were obtained, excluding the urban areas.

Preparation of the soil information is a basic input for generating CN. For this purpose, the study area was extracted from the digital soil map of the world developed by the FAO with specific soil codes which represent different soil hydrologic groups, using guidelines provided by the National Engineering Handbook of Hydrology (USDA-NRCS 2007) and described in Table 1. In the final step, processed soil and land-cover information were merged. The CN look-up table was created with six fields, namely, Land Use (LU) value, description of LU, and Soil Hydrologic Groups i.e., A, B, C and D, as shown in Table 2.

The table values for each land-cover and soil group were obtained from SCS TR55 of the United States Department of Agriculture Natural Resources Conservation Service (USDA-NRCS) (USDA-NRCS 1986). The CN look-up table along with the merged soil and land-cover feature classes were employed using HEC-Geo-HMS extension tools in ArcGIS for generating the CN grid (Figure 8).

Finally, the required mean CN values for each sub-basin of the study area were obtained (Figure 9), which were used for estimation of initial input values of various parameters such as lag time, time of concentration and storage coefficient for all the sub-basins.

Snow Water Equivalent (SWE) is a common snowpack measurement which determines the amount of water contained within the snowpack. It can be expressed as the depth of water that would theoretically result if the entire snowpack were melted. SWE can also be described as the equivalent amount of liquid water stored in the snowpack, mathematically defined as the product of the snow layer's depth and its density. Snow-cover area is an important input variable for computing the different HEC-HMS parameters for a snowmelt model such as SWE and melt-rate etc. MODIS Terra Snow Cover 8-Day L3 Global (MOD10A2) images from the website of Earth data during 2010–2012 were downloaded and projected to WGS 1984 UTM ZONE 43. Using GIS, the study area was extracted by using the mask technique.

In Equation (3), C_o is a constant (mm/^oC day), R is rainfall (mm/day), di is the day of the year (1=1 January) and default values are of α=0.25 and β=0.01 which is stated by Burek et al. (2013). Similarly, in Equation (4), Y, X1, X2 and X3 are constants for computing C and C₀, which are computed as a single value for all the days of the dry/winter season (November to May) using Equations (5)–(8) respectively.

For time series data, a hypothetical area-averaged gauge-station (constant precipitation over the entire sub-basin) was considered for each of the 11 sub-basins of the entire C-KRB using the specified hyetograph model. Observed daily discharge data from the Arandu hydrometric station (Chitral River outlet to downstream Afghanistan) was utilized which was obtained from the Surface Water Hydrology Project (SWHP), Water & Power Development Authority (WAPDA) Lahore, Pakistan. The temperature data for each sub-basin was obtained by lapse rate using meteorological data. SWE is also an important parameter in simulation of flow magnitude in snow-cover areas, which was computed for each sub-basin as described by Pistocchi et al. (2017).

The calibration or optimization of the HEC-HMS model for initial input parameters (Table 3) was carried out to produce the best suitable input variables according to the simulated and observed stream flow for the C-KRB over two subsequent hydrological years (2010–2011) independently. The model was first optimized for 2010 and then again applied for 2011 using the optimized results of 2010. The Univariate Gradient method with Sum of Squared Residuals objective function was used in the optimization process, with a maximum number of iterations of 300 at the 60 input parameters and tolerance value of 0.01. The optimized parameters along with known parameters of discharge and rainfall data were used for validation by running the model for the year 2012.

The optimization of the model was carried out for two independent hydrological years, i.e., 2010 and 2011, respectively. The results of the optimization trail for the year 2010 show a close relation between the simulated and observed stream flow at the Arandu hydrometric station (outlet of the C-KRB) as shown in Figure 10. The simulated stream flow was fluctuating and overestimated from February until May, and underestimated from June until December, while overall underestimation of stream flow was found to be owing to uncertainties in TRMM rainfall data because the ground rain-gauges and soil information were limited, which affected the model's performance. All the initial input parameters optimized in calibration are shown in Table 3, except initial and maximum deficit.

The model performance was checked by the NSE resulting in a value of 0.52, the deviation of runoff volume D_v was −28.84%, difference in peak flow (D_p) was 42.77%, and R² was 0.60. From the scatter graph (Figure 11) of the simulated vs observed stream flow, it is clear that the model simulations were closer to the observed discharge during the low-rate (below 1,000 cumec) period but underestimated during the increasing period (above 1,000 cumec).

Calibration was further conducted using the obtained parameters from the calibration year 2010 (Table 4) as input parameters for 2011, and the renewed calibration parameters were obtained as shown in Table 5. The simulated versus the observed stream flow for calibration over 2011 is illustrated in Figure 12. The simulated stream flow was found to be fluctuating in comparison with the observed data. Modeling-uncertainty and less accurate TRMM satellite rainfall data as well as limited soil information for each sub-basin had an effect on the model outputs.

Nevertheless, the developed model produced a satisfactory result with an NSE of 0.56 according to Moriasi et al. (2007), deviation of runoff volume D_v of −28.95%, difference in peak flow (D_p) of −33.46, and R² of 0.78 as presented in Figure 13. The scatter graph showing the simulated versus the observed stream flow with a 1:1 line for the calibration year (2011) in Figure 13 also supports the feasibility of the model for stream flow prediction in the study area.

The model was validated for the hydrological year 2012 utilizing the optimized parameters derived from the calibration year (Table 5). The performance of the validated model was evaluated by NSE and R² values, which were about 0.70 and 0.89, respectively, which is acceptable according to Moriasi et al. (2007). The results stated that significant fluctuations in magnitude of the simulated stream flow were observed (Figure 14), which was probably due to failure of the remotely sensed rainfall in capturing certain types of precipitation, along with the low accuracy of spatial and temporal resolution of the remotely sensed rainfall in not representing the true ground conditions in exact time (Ward et al. 2011). If the catchment area has highly accurate and a greater number of rain gauges, the model will more accurately simulate river flow.

Figure 14 illustrates that the simulated surface peak outflow for the validation year 2012 was 850 m³/s on 1 August, which was quite close to the observed peak flow of 861 m³/s on the same day. The generated peak stream-flows during July and August are probably attributable to the enhanced melting of snow and glaciers, due to the increased temperature of up to 38 °C that accelerates the snow melting rate. Although the model underestimated stream flow over the validation year with the simulated total outflow of about 430mm compared with the observed outflow of about 588mm, the model efficiency generated an acceptable output, i.e., NSE of 0.70 and R² of 0.89. In Figure 15, the scatter graph reveals that the model-simulated stream flow approximates to the observed flows at low-water periods but differs more at high-water periods of the study area. It is believed that the deviation of the simulated stream flow from the observed flows was mainly because of the uncertainty in TRMM rainfall data, since the model was more sensitive to rainfall data. Overall, the model accuracy was in the acceptable range (R²=0.89). Therefore, the model can be applied in future in the C-KRB and other similar watersheds.

This study demonstrated the feasibility of the HEC-HMS model in efficiently and effectively simulating rainfall–runoff with either daily or three-hour series of precipitation data for a highly mountainous, elevated and scarcely gauged basin like the C-KRB. The HEC-HMS model as a useful tool can provide the detailed hydrological processes at each sub-basin. Moreover, TRMM precipitation data is a useful alternative to local rain-gauge data and therefore can be utilized for HEC-HMS modeling in the absence of on-site meteorological observations.

The model was first calibrated for the year 2010 at the Chitral river-gauge station (hydrometric station at Arandu), which yielded an NSE of 0.52, deviation of runoff volume D_v of −28.84%, difference in peak (D_p) of 42.78%, and R² of 0.60, whereas for the calibration year 2011, the NSE-value increased to 0.57, R² value to 0.78 while D_p and D_v decreased to −33.46% and −28.95%, respectively. The calibration was performed using the parameters obtained from the model to obtain more refined and optimized values of initial input parameters for simulating more representative runoff in the validation year. Furthermore, the model was validated for the year 2012, which resulted in an NSE-value of 0.70, R² value of 0.89, D_p of −4.45% and D_v of −26.95%, as shown in Table 6. The HEC-HMS model does not use precipitation data as input in both forms (i.e., liquid rainfall and snowfall). Therefore, the simulated runoff peaks were poorly captured. Moreover, snow-cover information was not used directly in HEC-HMS; therefore, snowmelt runoff is potentially translated into simulated runoff errors. The HEC-HMS simulated runoff responded significantly to the snow cover during the months of July–August while the fluctuation during September–October was due to the snowmelt runoff contribution from the summer months (i.e., June–August) of the year which caused undulation in the simulated runoff magnitude, although the results of model validation (i.e. NSE=0.70) revealed better performance than other relevant studies of HEC-HMS application in the snow-cover watersheds (Gyawali & Watkins 2013; Khasmakhi et al. 2020). It is noteworthy that the HEC-HMS model is applicable for future flow prediction and performance is significantly dependent upon precisely estimated initial parameters and accuracy of rainfall data, which play a substantial role in improving the model runoff simulation and its efficiency.

In this study, the HEC-HMS model developed can be utilized for future stream flow prediction using the resultant calibrated parameters for the year 2011. Moreover, model validation results for the year 2012 showed that the NSE (i.e. 0.70) and R² (i.e. 0.89) values were very satisfactory and found an acceptable performance of the model according to the criteria of Ritter & Muñoz-Carpena (2013). In addition, the performance efficiency was in the range reported by other similar studies which declared that the HEC-HMS performance is closer to the approximately similar characteristics of nearby watershed study areas such as the Indus River Basin, Jhelum River (NSE value 0.70) studied by Azmat (2015) and Azmat et al. (2017) and as other researchers also used it for snow-covered catchments (Daly et al. 1999; Yilmaz et al. 2012). The slight uncertainty of the model and relatively poor correlations between the simulated and the observed stream flow was mainly attributed to the limited number of hydro-meteorological stations and consideration of average precipitation at each entire sub-basin. Also, the unevenly distributed locations of these hydro-meteorological stations could not capture the actual ground meteorological spatial variability caused by different physiographic attributes. The coarser spatial resolution of the TRMM rainfall data in such a complex river basin (i.e., C-KRB) also significantly affected the model performance.

The modeling experiments highlighted that in the case of more hydro-meteorological stations being available in the study watershed, better performance of the model could be expected. Moreover, in the present study, stream flow of the C-KRB was assessed using the present rainfall data (i.e., TRMM rainfall product 3B42), but the calibrated model could also be used to project future stream flow by using the representative concentration pathways (i.e., global climate models (GCMs) and regional climate models (RCMs)). This will help in understanding the future expected stream flow against the increasing future water demand for agricultural, municipal and industrial sectors, which may trigger the move of such simulation towards sustainable surface and groundwater resources management. The model performance can be further improved by using datasets from several rain gauges, stream flow gauges and snow-water-equivalent gauges. The present study will add value to cumulative stream flow simulation using HEC-HMS in a snow-cover area which may assist in shaping future river basin management plans.

The authors highly appreciate the Water & Power Development Authority (WAPDA) and Pakistan Meteorological Department (PMD), Pakistan, for providing the required data and making this research possible. The authors also thank Wanchang Zhang's research group students for their review and comments. This work was jointly supported by the Key R&D and Transformation Program of Qinghai Province [Grant No. 2020-SF-C37] and the National Key R&D Program of China [Grant No. 2016YFA0602302].

The authors declare that they have no conflict of interest or personal relationships that could have appeared to influence the work reported in this paper.

All relevant data are included in the paper or its Supplementary Information.