Abstract

Discharge over the Narayani river catchment of Nepal was simulated using Statkraft's Hydrologic Forecasting Toolbox (Shyft) forced with observations and three global forcing datasets: (i) ERA-Interim (ERA-I), (ii) Water and Global Change (WATCH) Forcing Data ERA-I (WFDEI), and (iii) COordinated Regional Climate Downscaling EXperiment with the contributing institute Rossy Centre, Swedish Meteorological and Hydrological Institute (CORDEX-SMHI). Not only does this provide an opportunity to evaluate discharge variability and uncertainty resulting from different forcing data but also it demonstrates the capability and potential of using these global datasets in data-sparse regions. The fidelity of discharge simulation is the greatest when using observations combined with the WFDEI forcing dataset (hybrid datasets). These results demonstrate the successful application of global forcing datasets for regional catchment-scale modeling in remote regions. The results were also promising to provide insight of the interannual variability in discharge. This study showed that while large biases in precipitation can be reduced by applying a precipitation correction factor (p_corr_factor), the best result is obtained using bias-corrected forcing data as input, i.e. the WFDEI outperformed other forcing datasets. Accordingly, the WFDEI forcing dataset holds great potential for improving our understanding of the hydrology of data-sparse Himalayan regions and providing the potential for prediction. The use of CORDEX-SMHI- and ERA-I-derived data requires further validation and bias correction, particularly over the high mountain regions.

INTRODUCTION

The Himalayan and adjacent Tibetan Plateau water supply is intricately linked to the livelihoods, economic and social, to millions of people (Bookhagen & Burbank 2010; Immerzeel et al. 2010; Remesan et al. 2019; Zhang et al. 2019). This region provides for drinking water, hydropower generation, agricultural demands, as well as water-powered grain mills and other agro-economic activities (Ménégoz et al. 2013). Large Asian rivers, such as the Indus, Sutlej, Ganges, and Brahmaputra, are widely acclaimed for their great cultural, spiritual, economic, and ecological significance (Kumar 2017). The Ganges River and its tributaries alone fulfill significant water demands for more than 250 million people. Originating in the Himalayas, the river travels over 2,500 km, aggregating water from its tributaries to become the third largest freshwater delta (Chowdhury & Ward 2004).

Due to the high precipitation rates (1,500–2,500 mm/year; Dahal & Hasegawa 2008) and steep elevation gradients, hydropower has significant potential in the region and, if well managed, could provide a resource for economic growth. Hydropower potential depends on climatic conditions, particularly on precipitation, evaporation, temperature, and snow/ice in the catchment (Edenhofer et al. 2011). However, the barriers for sustainable water resource management in this region are manifold. First, there are numerous socio-political issues related to trans-border management issues (Biswas 2011). Despite efforts, a comprehensive framework to guide the development and international cooperation has yet to be ratified (Biswas 2011). Second, portions of the region are heavily glaciated and climate change is exacting an immediate and tangible impact (Immerzeel et al. 2010; Dehecq et al. 2019).

In this paper, we focus on a third significant barrier: providing robust and reliable analysis of the water resource potential and impacts from climate change. Data scarcity and the complexity of the terrain creating large spatial gradients and variability in weather and climate make such analyses a challenging task. The extremely heterogeneous topography in the region presents a great challenge for the accurate measurement of meteorological variables, giving rise to data scarcity (Pellicciotti et al. 2012). Although there are some meteorological stations, they are not well distributed in space. To overcome these challenges, regional climate models and reanalysis data offer gridded datasets of many meteorological variables, although with a rather coarse spatial resolution (Guo & Su 2019). The Coordinated Regional Climatic Downscaling Experiment (CORDEX) (Giorgi et al. 2009) provides extensively used regional datasets for past and future climate. For a historical perspective, the reanalysis dataset ERA-Interim (ERA-I) from the European Center for Medium-Range Weather Forecasts (ECMWF) (Dee et al. 2011) is widely used (Li et al. 2013, 2018; Xu et al. 2016). Bharti & Singh (2015) reported that the ERA-I precipitation is largely overestimated over the Indian Himalayan region. To address the bias inherited in ERA-I, bias- and elevation correction-based monthly observations were carried out when creating the Water and Global Change Forcing Data ERA-I (WFDEI) datasets (Weedon et al. 2014; Kim et al. 2019).

Various factors make hydrological modeling in the region challenging, including a large spatial variability in hydrometeorological variables, steep gradients, marked seasonality driven by the Indian Monsoon, and contrasting moisture regimes between the high elevation Tibetan Plateau and regions in the vicinity of the Indian Ocean. Changes in the climatic condition in the region may lead to changes in regional water balance components impacting hydrologic regimes (Koch et al. 2011). Moreover, snow and glacier storage and melt play an important role for the river discharge generation (Radić & Hock 2014; Li et al. 2015, 2016). The effective assessment and implementation of hydropower projects and other water resource management projects depend on a thorough analysis of the discharge and hydrological storages (Bakken et al. 2013).

Hydrologic simulation of discharge and other water balance components (evapotranspiration, snow, and groundwater storage) (Bhattarai & Regmi 2015; Matt & Burkhart 2018; Li et al. 2019) is used for the analysis of available water resources both in the past and future analyses, but the quality of discharge simulation remains a challenge (Rochester 2010; Engeland et al. 2016; Kauffeldt et al. 2016). The choice of a suitable hydrological model and appropriate forcing data are critical for any analysis and will greatly affect the outcome (Kauffeldt et al. 2016). Currently, no standardized or community modeling framework within hydrology and prior studies of the water resources of the Himalayan region have been conducted using different types of hydrologic models ranging from simple conceptual models (e.g. Pradhananga et al. 2014; Bhattarai & Regmi 2015; Skaugen & Weltzien 2016; Bhattarai et al. 2018) to more advanced, distributed models (e.g. Pellicciotti et al. 2012; Jain et al. 2017). Hydrological modeling and water balance studies in Himalayan regions have taken a range of approaches and used different models to focus on topics such as glacier melt and retreat, water balance, flooding, and the impact of climate change (Bookhagen & Burbank 2010; Pellicciotti et al. 2012; Wortmann et al. 2014; Khanal et al. 2015). Recent studies have addressed the performance of large-scale forcing datasets, including the WFDEI and ERA-I, for discharge simulation in various regions across the world (e.g. Li et al. 2013, 2016; Weedon et al. 2014; Essou et al. 2016; Nkiaka et al. 2017), but, to date, the use of these global forcing datasets to simulate discharge in a Himalayan catchment is limited.

In this study, we use a distributed, conceptual hydrologic framework with an energy balance-based snow routine that has been demonstrated to perform well in the region (Xu et al. 2015; Hegdahl et al. 2016; Matt & Burkhart 2018). The main goal of this work is to evaluate the impact of four different forcing datasets on the simulated discharge and associated water balance assessment for the Narayani river catchment in Nepal. The forcing datasets include WFDEI, ERA-I, CORDEX, and a hybrid dataset referred to as Observed + WFDEI (all datasets are described in the Data and Methods section). Our specific objectives are to evaluate: (i) the discharge sensitivity to the forcing dataset for hydrologic modeling inflow, (ii) model performance using bias-corrected versus non-bias-corrected forcing as input to a (calibrated) hydrological model, and (iii) the sensitivity to forcing data for water balance assessment. Furthermore, we assess the overall ability of the Statkraft Hydrological Forecasting Toolbox (Shyft) (Burkhart et al. 2016) to simulate discharge in this complex and data-sparse regions.

STUDY AREA

The Narayani river catchment lies in the central part of Nepal. About 13% of the total area lies in China (Figure 1), and thus, it is a transboundary catchment. The main Narayani river gauging station is located in Narayanghat (27°42′30′′N, 84°25′50′′E) and is operated by the Department of Hydrology and Meteorology, Government of Nepal (DHM, GoN). The catchment area is 31,692 km2 and is partly glacier covered (∼8%; Omani et al. 2017) with elevation ranging from 175 m a.s.l. in the south to 8,148 m a.s.l. in the north. The hydrological regime is heavily influenced by the season. The seasons are defined as monsoon (June to September), post-monsoon (October to November), winter (December to February), and pre-monsoon (March to May) (Shrestha et al. 1999; Bhattarai et al. 2019). Almost 80% of total annual precipitation occurs during the monsoon period (Nayava 1974; Kripalani & Sontakke 1996). Tributaries to the Narayani river are either monsoon fed (those originating in middle and high mountain regions) or glacial and snow melt fed (those originating in a higher Himalayan region). The Narayani river catchment lies also in the primary hydropower development region of Nepal, providing 44% of the total electric generation (Adhikari 2006).

Figure 1

Narayani river catchment with the location of the discharge station and meteorological stations (72 precipitation and 30 temperature) in the region (both within and outside the catchment boundary).

Figure 1

Narayani river catchment with the location of the discharge station and meteorological stations (72 precipitation and 30 temperature) in the region (both within and outside the catchment boundary).

DATA AND METHODS

Meteorological forcing data

Forcing variables used in the study include air temperature (T), precipitation (P), relative humidity (RH), wind speed (WS), and shortwave incoming solar radiation (S) at a daily resolution. Details of each dataset are described in the following sections.

Observed data

The DHM, GoN is responsible for collecting and disseminating hydrological and meteorological information for water resources, agriculture, energy, and other development activities in the country. Most of the measured data are based on the conventional manual observation. Measurements from the automatic weather station (AWS) are available only after the year 2000 (Karki 2010). The geographic location and installation date for each station are available from http://dhm.gov.np/meteorological-station/.

Observed temperature and precipitation station data which are located inside or in the vicinity of the catchment (Figure 1) were collected from the DHM. Data from the DHM were manually plotted for each station and screened for data quality. Clearly, erroneous departures from the historical pattern were removed manually from the datasets. The highest temperature station elevation is 3,870 m a.s.l., located in Chhoser (DHM st. no 633), and the highest precipitation station elevation is 3,705 m a.s.l., located in Mustang (DHM st. no 612).

To accommodate missing data, stations with less than 10 years of record or missing more than 15% of the observations were removed from the datasets. Ten precipitation stations (out of 82) did not meet this criterion and were removed from the datasets. Similarly, seven (out of 30) temperature stations were removed from datasets. This resulted in 72 precipitation stations and 23 temperature stations that were used for further analysis. Maximum numbers of stations with temperature and precipitation data were observed in the 2000–2009 period. Therefore, in this study, data for the period 2000–2009 were used.

All the stations are located below 4,000 m a.s.l., covering only 60% of the total catchment area (Figure 2(a)). Generally, temperature and precipitation data from higher elevations are missing, which add uncertainty to the current study. Daily mean temperature and mean annual precipitation for each station are plotted versus station elevation in Figure 2(b). Temperature shows as expected a consistent decrease with elevation, whereas precipitation shows a more mixed picture, with no clear trend below approximately 2,000 m a.s.l., but a notable decrease in precipitation above, i.e. from about 2,000 to 3,500 m a.s.l. Normally precipitation increases with elevation (Daly et al. 1994) in the mountainous region due to the orographic effect, but in the Himalayan region, an opposite pattern has been reported after a certain elevation level (e.g. Nayava 1980; Kansakar et al. 2004), in agreement with what is observed for the Narayani river catchment. From the analysis, station average daily temperature and mean annual precipitation for the period 2000–2009 were 10.7 °C and 1,292 mm/year, respectively.

Figure 2

(a) Station altitude (temperature and precipitation) plotted against the hypsographic curve for the Narayani river catchment and (b) station average daily mean temperature and annual mean precipitation (for the period 2000–2009) plotted against altitude. Station name, number, and data availability are provided in the Supplementary Data.

Figure 2

(a) Station altitude (temperature and precipitation) plotted against the hypsographic curve for the Narayani river catchment and (b) station average daily mean temperature and annual mean precipitation (for the period 2000–2009) plotted against altitude. Station name, number, and data availability are provided in the Supplementary Data.

Discharge observed at the Narayanghat station is available from the year 2000–2009 with no missing values, and this is also the period used for the hydrological model calibration and validation. Daily average discharge for the period 2000–2009 was 1,482.5 m3/s.

Reanalysis and regional climate model data

A summary of the available time periods and the resolution of the forcing datasets are provided in Table 1. All datasets are well documented and are freely available.

Table 1

Summary of selected forcing datasets (n refers to the number of grid cells in the Narayani river catchment)

Forcing datasetData periodSpatial resolution (degrees)Temporal resolutionnReference
ERA-I 1979–2015 0.75 × 0.75 3 h, 6 h 12 Dee et al. (2011)  
WFDEI 1979–2016 0.5 × 0.5 Daily 20 Weedon et al. (2014)  
CORDEX-SMHI 1980–2010 0.44 × 0.44 Daily 25 Sanjay et al. (2017)  
Observed + WFDEI 1999–2010 – Daily –  
Forcing datasetData periodSpatial resolution (degrees)Temporal resolutionnReference
ERA-I 1979–2015 0.75 × 0.75 3 h, 6 h 12 Dee et al. (2011)  
WFDEI 1979–2016 0.5 × 0.5 Daily 20 Weedon et al. (2014)  
CORDEX-SMHI 1980–2010 0.44 × 0.44 Daily 25 Sanjay et al. (2017)  
Observed + WFDEI 1999–2010 – Daily –  

ERA-I is a reanalysis global forcing dataset available from 1979, produced by the ECMWF. ERA-I temperature results from the assimilated surface temperature (Essou et al. 2016), while precipitation data are based on a reanalysis of precipitation fields generated with a meteorological model (Berrisford et al. 2011; Dee et al. 2011). The obtained precipitation data are not scaled using observation data. ERA-I is continuously updated once per month, with a delay of 2 months, and is freely available from http://apps.ecmwf.int/datasets/. Monthly mean meteorological variables averaged over the Narayani river catchment using the Shyft interpolation routines (see the section ‘Spatial-temporal interpolation of forcing data’) are shown in Figure 3. The catchment average mean temperature for the period (2000–2009) was 7.6 °C, while catchment average annual precipitation was 4,660 mm/year.

Figure 3

Catchment average monthly mean Shyft-interpolated meteorological variables (for the period 2000–2009) for each forcing dataset.

Figure 3

Catchment average monthly mean Shyft-interpolated meteorological variables (for the period 2000–2009) for each forcing dataset.

The WFDEI is a global forcing dataset obtained by downscaling and bias-correcting ERA-I data (Weedon et al. 2014; Raimonet et al. 2017). The temporal and horizontal resolution of the dataset is shown in Table 1. WFDEI has two sets of rainfall generated by using either Climate Research Unit (CRU) or Global Precipitation Climatology Centre (GPCC) precipitation correction methods (Weedon et al. 2011). In this study, GPCC-corrected data were preferred to CRU because of their higher resolution and data quality (Weedon et al. 2014; Raimonet et al. 2017). The WFDEI dataset is freely available online from ftp.iiasa.ac.at. Catchment average mean daily temperature and annual precipitation (for the period 2000–2009) were found to be 7.9 °C and 1,764 mm/year, respectively.

CORDEX is a program sponsored by the World Climate Research Programme (WCRP), to produce an improved generation of regional climate change projections (Giorgi et al. 2009). CORDEX has two datasets, referred to Evaluation and Historical. Evaluation is run within reanalysis and is used to ‘mimic’ observations (i.e. represent real weather), whereas Historical is run within a climate model, and the results can only be used in a climatological sense. Daily datasets of the Evaluation product over the study area are available for the South Asia CORDEX domain (Region 6) from two institutional runs: Indian Institute of Tropical Meteorology (CORDEX-IITM) and Swedish Meteorological and Hydrological Institute (CORDEX-SMHI) (Sanjay et al. 2017). Based on a comparison of seasonal precipitation patterns over the study area, the CORDEX-SMHI was selected for this study. The horizontal resolution and available periods for the CORDEX-SMHI data are shown in Table 1. CORDEX-SMHI data can be freely downloaded on request from http://cccr.tropmet.res.in/home/index.jsp. Catchment average mean daily temperature and annual precipitation (for the period 2000–2009) are found to be 6.1 °C and 5,431 mm/year, respectively.

Topographical and land cover datasets

In this study, a digital elevation model (DEM) of 90 m spatial resolution from NASA's Shuttle Radar Topography Mission (NASA-SRTM) was used. The NASA-SRTM DEM is freely available for download from https://eros.usgs.gov/. The hydrology tool available in the software package, System for Automated Geoscientific Analysis (SAGA) (Conrad et al. 2015) under Qgis-2.18 (QGIS Development Team 2016), was used for automatic catchment delineation. Catchment delineation was based on the gauge point at Narayanghat and NASA-SRTM DEM. Catchment slope and aspect data for Shyft (see the Hydrological Model section) were generated during the catchment delineation sub-processes in Qgis. The calculated catchment domain was gridded into 2 km × 2 km cells for input to Shyft. The centroid of each vector grid cell was calculated using Qgis as further input to the model.

A land cover map (0.5° × 0.5°, resolution) providing forest, lake, glacier, and reservoir cover of the study area was extracted from the Moderate Resolution Imaging Spectroradiometer (MODIS) land cover data (Channan et al. 2014). MODIS land cover data are available free online at http://glcf.umd.edu/data/lc/. Land cover fractions of forest, lake, glacier, and reservoir cover area at the centroid point per each grid cell within the catchment were also calculated. The python code and algorithm for calculating land cover fraction per each vector grid cell at centroid are available online at https://github.com/felixmatt/shyft-gis.

Spatial interpolation of observed and gridded forcing data

Shyft is able to ingest daily spatially distributed meteorological data as input variables. Here, we interpolate both observations (station data) and gridded forcing data (WFDEI, ERA-I, and CORDEX-SMHI).

The gridded forcing data have a coarser resolution than the simulation domain in the hydrological modeling framework of Shyft (Table 1). Accordingly, Shyft provides interpolation routines used for downscaling the forcing datasets to the model grid cell scale. Data interpolation methods were selected based on prior studies by Sapkota (2016), Lombraña (2017), Matt et al. (2018), Matt & Burkhart (2018), and Teweldebrhan et al. (2018). Observed temperature and gridded data were extrapolated/interpolated by three-dimensional Kriging where elevation is the third dimension (Oliver & Webster 1990). Observed precipitation, gridded precipitation, global radiation, WS, and RH data were extrapolated/interpolated by inverse distance weighting (IDW) (Shepard 1968). A precipitation gradient of −0.07 mm/100 m (Table 2) and a temperature gradient of −0.6 °C/100 m are used during interpolation. The precipitation gradient was calculated from the observed precipitation data. The different interpolation methods along with their parameter values are presented in Table 2.

Table 2

Parameter values used for IDW and Kriging interpolation

 IDW
Kriging
ParameterValueParameterValue
Maximum distance 400,000 m Temperature gradients −0.6 °C/100 m 
Maximum number 10 Range 20,000 m 
Distance measure factor z-scale 20 
Resolution 2 km × 2 km Resolution 2 km × 2 km 
– – Sill 25 °C 
Precipitation gradient −0.07 mm/100 m   
 IDW
Kriging
ParameterValueParameterValue
Maximum distance 400,000 m Temperature gradients −0.6 °C/100 m 
Maximum number 10 Range 20,000 m 
Distance measure factor z-scale 20 
Resolution 2 km × 2 km Resolution 2 km × 2 km 
– – Sill 25 °C 
Precipitation gradient −0.07 mm/100 m   

Hydrological model

The Shyft modeling framework was used to simulate daily discharge from the catchment. The modeling framework has three main hydrologic modeling routines (https://gitlab.com/shyft-os). These three models are different in the way to calculate evapotranspiration and snow estimation and melt.

In this study, the PT_GS_K model was used for discharge simulation. The PT_GS_K model uses the Priestley–Taylor (PT) method (Priestley & Taylor 1972) for estimating potential evapotranspiration, the Gamma snow routine (GS; described in the Gamma Snow section) for snowmelt, sub-grid snow distribution and mass balance calculation, and a simple storage–discharge function (K) for catchment response calculation (Lambert 1972; Kirchner 2009). Actual evapotranspiration (AE) was assumed to take place from the snow-free area and was estimated as a function of potential evapotranspiration and a scaling factor.

The catchment response function ‘K’ is based on the storage–discharge relationship concept described in Kirchner (2009) and represents the sensitivity of discharge to storage changes as given by the following equation. Generally, discharge is nonlinear and varies by many orders of magnitude, and it is recommended to use log-transformed discharge values to avoid numerical instability. 
formula
(1)
where P, Q, and E are the rate of precipitation (input from either snow melt or rain), discharge, and AE (from the snow-free area), respectively, in units of depth per time.
The idea behind this method is that the discharge sensitivity to changes in storage, i.e. g(Q), can be estimated from the time series of the discharge alone through fitting empirical functions to the data, such as the quadratic equation (Kirchner 2009), which is given by: 
formula
(2)
where c1, c2, and c3 are the catchment-specific outlet parameters (hereafter called Kirchner parameters) obtained during the model calibration. There is no routing function in the model, but the Kirchner response function represents a delayed outflow from storage within the catchment.

Gamma snow

In Shyft, the Gamma snow routine, an energy balance approach (Equation (3)) for snow ablation, and the snow depletion curve following a gamma distribution are combined into a single routine. The energy balance approach in the Gamma snow routine is based on DeWalle & Rango (2008) and also briefly explained in Hegdahl et al. (2016). The net energy flux () at the surface available for snow ablation is expressed as follows: 
formula
(3)
where S is the net shortwave radiation, Lin and Lout are the incoming and outgoing longwave radiations, HSE and HL represent sensible and latent heat fluxes, and EG is the net ground heat flux calculated using a bulk-transfer approach. Two parameters defining the wind profile, intercept (wind constant), and slope (wind scale) are determined either by model calibration or as provided (Table 3). For a given time step (t), the snow albedo () at each cell depends on the minimum and maximum albedo as well as the albedo decay rate, temperature, and snowfall as described in Hegdahl et al. (2016): 
formula
(4)
Table 3

Model calibration parameters with upper and lower limits

ParameterDescription and unitParameter used in the submodelLower limitUpper limitSources
cOutlet empirical coefficient 1 (–) Equation (4) −8.0 0.0 Sapkota (2016); Lombraña (2017)  
cOutlet empirical coefficient 2 (–) Equation (4) −1.0 1.2 Sapkota (2016); Lombraña (2017)  
cOutlet empirical coefficient 3 (–) Equation (4) −0.15 −0.05 Sapkota (2016); Lombraña (2017)  
wind scale Slope in turbulent wind function (m/s) GS 1.0 6.0 Sapkota (2016); Lombraña (2017)  
tx Temperature threshold rain/snow (°C) GS −3.0 2.0 Sapkota (2016); Lombraña (2017)  
FDR Fast albedo decay rate during melt (days) GS 5.0 15.0 Hegdahl et al. (2016); Sapkota (2016)  
SDR Slow albedo decay rate during cold condition (days) GS 20.0 40.0 Hegdahl et al. (2016)  
p_corr_factor Scaling factor for precipitation (–) GS 0.4 2.0 Sapkota (2016)  
Prescribed parameters
ae scale factor Scaling factor for AE (–) AE 1.0 1.0 Sapkota (2016); Lombraña (2017)  
Surface magnitude Snow heat constant (mm SWE) GS 30.0 30.0 Sapkota (2016)  
wind constant Intercept in turbulent wind function (–) GS 1.0 1.0 Lombraña (2017)  
min albedo Minimum snow albedo used in snow routine GS 0.6 0.6 Sapkota (2016)  
max albedo Maximum snow albedo used in snow routine GS 0.9 0.9 Sapkota (2016)  
max water Frictional max water constant of snow (–) GS 0.1 0.1 Hegdahl et al. (2016)  
ParameterDescription and unitParameter used in the submodelLower limitUpper limitSources
cOutlet empirical coefficient 1 (–) Equation (4) −8.0 0.0 Sapkota (2016); Lombraña (2017)  
cOutlet empirical coefficient 2 (–) Equation (4) −1.0 1.2 Sapkota (2016); Lombraña (2017)  
cOutlet empirical coefficient 3 (–) Equation (4) −0.15 −0.05 Sapkota (2016); Lombraña (2017)  
wind scale Slope in turbulent wind function (m/s) GS 1.0 6.0 Sapkota (2016); Lombraña (2017)  
tx Temperature threshold rain/snow (°C) GS −3.0 2.0 Sapkota (2016); Lombraña (2017)  
FDR Fast albedo decay rate during melt (days) GS 5.0 15.0 Hegdahl et al. (2016); Sapkota (2016)  
SDR Slow albedo decay rate during cold condition (days) GS 20.0 40.0 Hegdahl et al. (2016)  
p_corr_factor Scaling factor for precipitation (–) GS 0.4 2.0 Sapkota (2016)  
Prescribed parameters
ae scale factor Scaling factor for AE (–) AE 1.0 1.0 Sapkota (2016); Lombraña (2017)  
Surface magnitude Snow heat constant (mm SWE) GS 30.0 30.0 Sapkota (2016)  
wind constant Intercept in turbulent wind function (–) GS 1.0 1.0 Lombraña (2017)  
min albedo Minimum snow albedo used in snow routine GS 0.6 0.6 Sapkota (2016)  
max albedo Maximum snow albedo used in snow routine GS 0.9 0.9 Sapkota (2016)  
max water Frictional max water constant of snow (–) GS 0.1 0.1 Hegdahl et al. (2016)  

In the table, ‘K’ is the catchment response function; ‘GS’ is the Gamma snow; ‘PT’ is the Priestley and Taylor; and ‘AE’ is the actual evapotranspiration.

Table 4

Definitions of the evaluation criteria

CriteriaMathematical expressionDescriptionBest value
  Volume difference in percentage 
NSE  Nash–Sutcliffe efficiency 
  Squared-root transformed Nash–Sutcliffe efficiency 
 

 
Modified Kling–Gupta efficiency 
 
 
Goodness of fit with respect to the benchmark series 
CriteriaMathematical expressionDescriptionBest value
  Volume difference in percentage 
NSE  Nash–Sutcliffe efficiency 
  Squared-root transformed Nash–Sutcliffe efficiency 
 

 
Modified Kling–Gupta efficiency 
 
 
Goodness of fit with respect to the benchmark series 

n: total number of days in the evaluation period; and : simulated and observed discharges; and : standard deviation for simulated and observed discharges; and : arithmetic mean for simulated and observed discharges; and : total discharge volume for simulated and observed discharges; Qbench: monthly long-term average observed discharge.

In Equation (4), FDR and SDR denote fast and slow snow cover decay rates, respectively. In this study, and are prescribed (refer Table 3).

Within the Gamma snow routine, precipitation falling in each cell is classified as solid or liquid depending on a threshold temperature (tx) and the actual cell temperature. Snow distribution within each cell is estimated by using a three-parameter gamma probability distribution. The third parameter in the gamma probability distribution represents the bare ground fraction in the cell. Finally, snowmelt depth (mm/day) is calculated by multiplying (available energy) with the latent heat of fusion for water.

A temperature index model which does not require glacier ice albedo was used to calculate glacier melt (see Hock 2003). The glacier reservoir was assumed to be infinite, and the glacier area was assumed to be constant throughout the simulation periods. Within a glacier-covered cell, glacier melt only happens from the snow-free fraction.

Parameters and calibration

Hydrological simulation from distributed models such as Shyft generally requires the estimation of model parameters through calibration with measured data (Madsen 2003). In the PT_GS_K routine, there are 14 parameters, out of which eight parameters (Table 3: top eight parameters) have been found to be significantly more sensitive than the rest (Teweldebrhan et al. 2018) and were selected for calibration in this study. The remaining six parameters were prescribed (Table 3: lower six parameters). The precipitation correction factor (p_corr_factor), which is used to correct bias, was also set as a calibration parameter.

Manual calibration can be time-consuming and subjective; therefore, an automatic calibration was carried out. Upper and lower limits for each parameter are shown in Table 3. CREST v.2.1, the Shuffle Complex Evolution University of Arizona (SCE-UA) (Duan et al. 1992), is used as the kernel algorithm in the automatic calibration process. Typically, the procedure involves the selection of samples in the parameter space through the use of competitive evolution schemes, such as the simplex scheme, to reproduce better the observations. After several iterations, either due to convergence or when the maximum number of iteration is reached, the best set of parameter values based on the Nash–Sutcliffe coefficient of efficiency (NSE) is determined (Chu et al. 2010).

In this study, the model was calibrated and validated for each forcing dataset independently (hereafter named ‘independent calibration mode’). The model calibration was based on observed discharge data from 2000 to 2004, and validation was based on data from 2005 to 2009.

Water balance estimation

A water balance analysis is a useful tool to describe the principal components of water in and out of a catchment (Rochester 2010), where the volume of water inflow should be balanced with water outflow, assuming no changes in storage . The following water balance components: mean annual precipitation, discharge, and AE were calculated. Water balance components were calculated using calendar years, so that the change in storage is mainly a difference in storage between first and last days of simulation.

The annual change in storage was calculated according to the following equation: 
formula
(5)
where P is precipitation, Q is the discharge, ET is actual evapotranspiration, G is the glacier melt, and the unit of measurement is mm/year. is the change in storage per time unit.

Model performance evaluation

To determine the agreements between observed and simulated discharges, NSE (Nash & Sutcliffe 1970), Kling–Gupta efficiency (KGE) (Gupta et al. 2009), square-root-transformed NSE (NSEsqrt) (Seiller et al. 2012), benchmark series (Gbench) (Seibert 2001), and percentage volume difference (Dv) (Martinec et al. 1998) were used (Table 4). The NSE, NSEsqrt, and KGE efficiency measures are used to access the predictive power of the hydrological model and can range from −∞ to 1. An efficiency of 1 corresponds to a perfect match of simulated and observed discharges. To improve the estimation of the performance error, KGE considers three components: bias (), variability (), and linear correlation coefficient () to overcome the problems associated with NSE, i.e. NSE results in the underestimation of the streamflow variability and the runoff peaks (Gupta et al. 2009). NSEsqrt gives more emphasis to the overall agreement between observed and simulated discharges (Seiller et al. 2012; Peña-Arancibia et al. 2015) as compared to NSE, so it is also included as an evaluation criterion. Percentage volume difference (Dv) gives the percentage bias between the simulated and observed series and can range from −∞ to +∞. Dv equals to zero indicates a perfect agreement between simulated and observed discharges. Since the observed discharge shows strong seasonal patterns, we also used a benchmark series (Gbench) for model evaluation. The monthly average discharge was used as a benchmark series. Gbench is negative if the model performance is poorer than the benchmark, zero if the model performs as well as the benchmark, and positive if the model is superior, with a highest value of one for a perfect fit. In this study, we aimed to achieve NSE > 0.7, KGE > 0.7, NSEsqrt > 0.7, Gbench > 0.5, and Dv within ±15%, during both the model calibration and validation periods.

RESULTS

Meteorological forcing data analysis

As discussed in the section ‘Reanalysis and regional climate model data’, the resolutions of gridded datasets are different. For the sake of comparison, Shyftinterpolated meteorological variables (for the period 2000–2009) for each forcing dataset were used. Temperature shows an overall agreement where all three datasets demonstrate similar seasonality and distinct cool or negative bias as compared to the observations during the winter and pre-monsoon seasons (Figure 3). The bias was the strongest in the CORDEX-SMHI data during the September–December period. Long-term spatially averaged temperature for July (Figure 4) shows that the CORDEX-SMHI was not able to capture higher daily average temperatures in the river valley. Daily average temperatures during July over the lower elevation region from ERA-I and CORDEX-SMHI were lower than the temperatures from the Observed and WFDEI interpolated datasets. What was most interesting is that negative catchment average temperatures, which were not at all captured by the interpolated observations (Figure 3), will have a significant impact on any simulation. We attribute the lack of negative catchment average temperatures to the lack of stations above 4,000 m a.s.l. combined with a too low-temperature lapse rate. This is a typical challenge when observations lack representativity to topography.

Figure 4

Distributed mean daily Shyft-interpolated temperature (°C) for July over the Narayani river catchment for each forcing dataset.

Figure 4

Distributed mean daily Shyft-interpolated temperature (°C) for July over the Narayani river catchment for each forcing dataset.

Precipitation varies notably among the forcing datasets. WFDEI precipitation agrees well with the observation (Figures 3 and 4). The strongest positive bias was observed for CORDEX-SMHI and ERA-I datasets (Figure 5). Long-term daily averaged precipitation for July–September shows the strongest deviation among forcing. The catchment precipitation for July was more than 150 mm/day in CORDEX-SMHI and ERA-I datasets, and the highest precipitation was found in the lower elevations. A seasonal comparison shows that the WFDEI was slightly higher than the observed precipitation during the monsoon season (Figure 6). Similar findings are reported for the Indus catchment by Dahri et al. (2016).

Figure 5

Distributed daily mean Shyft-interpolated precipitation (mm) for July over the Narayani river catchment for each forcing dataset.

Figure 5

Distributed daily mean Shyft-interpolated precipitation (mm) for July over the Narayani river catchment for each forcing dataset.

Figure 6

Seasonal precipitations from different forcing datasets.

Figure 6

Seasonal precipitations from different forcing datasets.

Monthly averaged wind speeds and incoming shortwave radiation from WFDEI and ERA-I were fairly similar but different for CORDEX-SMHI. Station observations for RH, WS, and shortwave incoming global radiation were not available for the catchment, but previous global studies by Iizumi et al. (2014) and Weedon et al. (2014) showed that the WFDEI is more representative for the observations.

Model parameters

Shyft was calibrated for each forcing dataset, and the calibrated parameters are shown in Table 5. Normal ranges for all calibrated parameters are given in Table 3. The highest range in values (greater than 0.15 variance in the normalized scale of −1 to 1) was seen for the parameters p_corr_factor (0.4–1.34), wind scale (2.01–4.84), threshold temperature (tx) (−1.68 to 0.94), and slow albedo decay rate (29.04–39.98). A relatively small variation in parameter values (less than 0.05 variance in the normalized scale of −1 to 1) was seen for the Kirchner coefficients (c1 (−6.64 to −5.51), c2 (0.14 to 0.44), and c3 (−0.05 to −0.07)) and fast albedo decay rate (9.94–12.03). The highest variance (0.27) in a normalized scale of −1 to 1 was found for p_corr_factor. The lowest and highest p_corr_factor values were observed for CORDEX-SMHI (0.4) and Observed + WFDEI (1.34) forcing dataset, respectively. The higher p_corr_factor for Observed + WFDEI is interesting and suggests an underprediction of precipitation in the catchment when gridded precipitation was based on Shyft-interpolated station observations. Unique to the CORDEX-SMHI dataset, the Gamma snow threshold temperature (tx) was found to be positive (0.94). Less sensitive parameters as reported by Teweldebrhan et al. (2018) were not calibrated but are listed with given values in Table 3.

Table 5

Estimated parameter values from the independent calibration mode for the four forcing datasets

ParametersUnitObserved + WFDEIWFDEIERA-ICORDEX-SMHI
c– −5.51 −6.09 −6.64 −5.99 
c– 0.44 0.25 0.14 −0.20 
c– −0.05 −0.05 −0.05 −0.07 
tx °C −1.68 −1.38 −0.36 0.94 
wind scale m/s 2.01 4.70 4.84 3.40 
FDR days 11.59 12.03 9.94 10.75 
SDR days 39.98 32.66 29.04 31.55 
p_corr_factor – 1.34 1.09 0.41 0.4 
ParametersUnitObserved + WFDEIWFDEIERA-ICORDEX-SMHI
c– −5.51 −6.09 −6.64 −5.99 
c– 0.44 0.25 0.14 −0.20 
c– −0.05 −0.05 −0.05 −0.07 
tx °C −1.68 −1.38 −0.36 0.94 
wind scale m/s 2.01 4.70 4.84 3.40 
FDR days 11.59 12.03 9.94 10.75 
SDR days 39.98 32.66 29.04 31.55 
p_corr_factor – 1.34 1.09 0.41 0.4 

Evaluation of discharge simulation using different forcing datasets

Figure 7 compares daily observed and simulated discharges using the different forcing datasets. Calibration and validation results in terms of error statistics are shown in Table 6. Simulations from Observed + WFDEI forcing give the best results in terms of error statistics during both calibration (NSE = 0.90, NSEsqrt = 0.91, KGE = 0.94, Gbench = 0.43, and Dv = −0.91%) and validation (NSE = 0.90, NSEsqrt = 0.93, KGE = 0.92, Gbench = 0.56, and Dv = −1.29%) periods. Dv in both calibration and validation was also best for Observed + WFDEI and found to be less than −2%. The second best performance was achieved for WFDEI (Table 6). For NSE, NSEsqrt, and KGE, the ERA-I and WFDEI datasets were similar, but the volume differences during calibration were comparatively higher for ERA-I than WFDEI. Comparatively poorer model performance and higher volume differences (Table 6) were found for the model calibrated with the CORDEX-SMHI dataset. Goodness of fit with respect to the benchmark (Gbench) for ERA-I and CORDEX-SMHI during calibration periods was found negative, indicating poorer performance than the benchmark series. However, it should be noted the peak discharge was best simulated with CORDEX-SMHI data (Figure 7).

Table 6

Performance statistics for different forcing datasets

Error parametersObserved + WFDEI
WFDEI
ERA-I
CORDEX-SMHI
CalibrationValidationCalibrationValidationCalibrationValidationCalibrationValidation
NSE 0.90 0.90 0.83 0.80 0.82 0.83 0.60 0.48 
KGE 0.94 0.92 0.86 0.81 0.82 0.81 0.73 0.52 
Dv −0.91 −1.29 −1.72 −5.68 −4.01 −5.53 −19.9 −42.4 
NSEsqrt 0.91 0.93 0.88 0.84 0.88 0.87 0.70 0.57 
Gbench 0.43 0.56 0.01 0.07 −0.03 0.19 −1.34 −1.47 
Error parametersObserved + WFDEI
WFDEI
ERA-I
CORDEX-SMHI
CalibrationValidationCalibrationValidationCalibrationValidationCalibrationValidation
NSE 0.90 0.90 0.83 0.80 0.82 0.83 0.60 0.48 
KGE 0.94 0.92 0.86 0.81 0.82 0.81 0.73 0.52 
Dv −0.91 −1.29 −1.72 −5.68 −4.01 −5.53 −19.9 −42.4 
NSEsqrt 0.91 0.93 0.88 0.84 0.88 0.87 0.70 0.57 
Gbench 0.43 0.56 0.01 0.07 −0.03 0.19 −1.34 −1.47 
Figure 7

Calibrated and simulated daily discharge for the Narayani river catchment for the period 2000–2009.

Figure 7

Calibrated and simulated daily discharge for the Narayani river catchment for the period 2000–2009.

In Figure 8, the year 2004 was plotted to show in more details the representation of the annual cycle by the different forcing datasets. Overall, all datasets were able to simulate the cycle well; however, the CORDEX-SMHI simulation deviates somewhat from the rest and shows a higher peak flow early in the wet season.

Figure 8

Daily observed and simulated discharges from different forcing models for the year 2004.

Figure 8

Daily observed and simulated discharges from different forcing models for the year 2004.

Quantile–Quantile (QQ) plots for the different forcing datasets are shown in Figure 9. The QQ-plot for the calibration period is shown in Figure 9(a), for the validation period in Figure 9(b). To highlight low flows, Figures 9(c) and 9(d) show the QQ-plot on log-scale. The QQ-plot is a graphical technique to determine if two datasets come from the same population with a common distribution (Renard et al. 2010). Departures from the 1:1 reference line indicate discrepancies between the simulated and the observed discharge. Figure 9(a) shows that all simulations slightly overestimate discharge (with respect to the observed discharge) up to a value of around 4,000 m3/s. Between 4,000 and 8,000 m3/s, the CORDEX-SMHI simulations were higher than observed, whereas other simulations were lower than observed. After 8,000 m3/s, CORDEX-SMHI simulations were able to capture peaks where others simulations were lower than observed. Although the Observed + WFDEI forcing dataset resulted in discharge simulations with the best performance, it was not able to capture the highest observed discharge (simulated discharge points after 95% quantile). The simulations from CORDEX-SMHI were able to capture relatively high peak-discharge events (above the 99% quantile); however, they have large deviations from the observations between 90% and 95% quantiles. For the low flows (i.e. 1–10% quantile), the simulation from the WFDEI forcing captures observation better than the Observed + WFDEI forcing (Figure 9(c)). However, low-flow simulation from ERA-I until 5% quantile was better than the rest.

Figure 9

QQ-plots of simulated versus observed discharge for the Narayani river catchment during model (a) calibration and (b) validation. Similarly, the QQ-plot for log-transformed observed and simulated discharges during model (c) calibration and (d) validation. The vertical lines denote the 1%, 5%, 10%, 90%, 95%, and 99% quantiles of the observed daily discharge.

Figure 9

QQ-plots of simulated versus observed discharge for the Narayani river catchment during model (a) calibration and (b) validation. Similarly, the QQ-plot for log-transformed observed and simulated discharges during model (c) calibration and (d) validation. The vertical lines denote the 1%, 5%, 10%, 90%, 95%, and 99% quantiles of the observed daily discharge.

Similarly, in the validation period (Figure 9(b)), simulated discharge and observed daily discharge for values less than 90% quantile from all forcing datasets fit well to the observations. Above 95% quantile, most forcing datasets fail to capture these highest discharges. Interestingly, the CORDEX-SMHI, which generally provides lower performance (Table 6), notable in the lower flow range, manages to capture the highest flows best. However, these simulations were mostly higher than the reference line, indicating a general overestimation during the full range of simulations. Similar to the calibration periods, low-flow simulation from WFDEI was also found closest to the reference (observations) during the validation periods (Figure 9(d)).

Water balance analysis

Figure 10 compares the catchment average (2000–2009) water balance components for the four different forcing datasets (the exact values are given in Table 7). Annual averages over the 10-year simulation period, 25.1%, 24.8%, 30.7%, and 17.8% of precipitation, were lost to evapotranspiration for Observed + WFDEI, WFDEI, ERA-I, and CORDEX-SMHI datasets, respectively. Among all forcing datasets, Observed + WFDEI shows lower runoff than the others. Higher runoff from CORDEX-SMHI forcing was observed, and it is associated with higher average precipitation during the monsoon period (Figure 3). The lowest (72.8 mm) total glacier melt during the simulation period was observed for the WFDEI dataset, while the highest (297.3 mm) was observed from ERA-I. The change in storage (Equation (5)) shows that 73.1, 32.7, 68.7, and 29.2 mm of water were surplus (a positive change in storage) in the catchment for Observed + WFDEI, CORDEX-SMHI, ERA-I, and WFDEI forcing data, respectively. The largest glacier melt contribution (19%) to total runoff was observed for ERA-I forcing, while the lowest was observed for WFDEI (4%). Different water balance components and their associated changes in storage are shown in Table 7.

Table 7

Water balance components for different forcing (numbers in mm/year)

Water balance componentsObserved + WFDEIWFDEIERA-ICORDEX-SMHI
Precipitation 1,771.1 1,965.7 1,901.4 2,234.1 
Evapotranspiration 449.6 487.4 584.6 399.5 
Glacier melt 242.2 72.8 297.3 132.5 
Runoff 1,490.7 1,521.9 1,545.4 1,934.4 
Change in storage (Equation (5)) 73.0 29.2 68.7 32.7 
Water balance componentsObserved + WFDEIWFDEIERA-ICORDEX-SMHI
Precipitation 1,771.1 1,965.7 1,901.4 2,234.1 
Evapotranspiration 449.6 487.4 584.6 399.5 
Glacier melt 242.2 72.8 297.3 132.5 
Runoff 1,490.7 1,521.9 1,545.4 1,934.4 
Change in storage (Equation (5)) 73.0 29.2 68.7 32.7 
Figure 10

Annual average (2000–2009) water balance components (Evap, actual evapotranspiration; Prec, precipitation calculated by using p_corr_factor) for different forcing datasets in the Narayani catchment.

Figure 10

Annual average (2000–2009) water balance components (Evap, actual evapotranspiration; Prec, precipitation calculated by using p_corr_factor) for different forcing datasets in the Narayani catchment.

DISCUSSION

Discussion on model parameters

As with any hydrologic model, the model used in this study is sensitive to the parameters; particularly the Kirchner coefficients, precipitation correction factor, and the threshold temperature (Teweldebrhan et al. 2018). While calibrating the model, different values of the parameters were obtained for different forcing datasets (Table 5). An examination of the model parameters revealed that the calibrated values of the precipitation correction factor, wind scale, threshold temperature, and slow albedo decay rates differ a lot between the forcing datasets. Larson & Peck (1974) and Lakew et al. (2017) pointed out that the model simulations are highly sensitive to the precipitation correction factor, which helps adjusting systematic bias embedded in the precipitation forcing. The precipitation correction factor may account for both under catch and lack of representative stations (Lakew et al. 2017) and is directly related to the bias in precipitation data. Prior studies by Engeland et al. (2016) and Teweldebrhan et al. (2018) suggested that the Kirchner parameters are highly sensitive since they have a large influence on a given simulation. We found that the variance for Kirchner response parameters was less than 0.05. An explanation is that these parameters are not dependent on the forcing datasets but rather on the physical catchment characteristics and thereby the concentration time of the catchment. Kirchner response parameters are dependent on the physical catchment characteristics and thereby the concentration time of the catchment and are therefore less sensitive to the forcing.

Potential factors controlling hydrological model efficiency during model calibration and validation

Different forcing datasets lead to different levels of performance efficiency in terms of discharge simulations (Table 6). Explaining these differences is not easy, as each forcing dataset has different characteristics such as the spatial resolution and the assimilation methods. Figures 35 present the catchment average monthly mean forcing variables for the different forcing datasets. From the figures, we can see that in general, the spatio-temporal distribution of precipitation from WFDEI was quite similar to the interpolated observations. The extent of higher precipitation than observation during pre-monsoon and monsoon seasons was larger for ERA-I and CORDEX-SMHI. Overestimation of precipitation during the monsoon season identified for CORDEX-SMHI in the Himalayan region (Ghimire et al. 2018) and for the ERA-I over the Indus catchment (Dahri et al. 2016). Moulin et al. (2009) and Zhu et al. (2014) indicate that a poor representativeness of precipitation over the catchment is a major source of uncertainty in discharge simulation.

Despite some of the biases in the forcing data, the performance of the model during calibration and validation showed that all datasets performed reasonably well. The good agreement between simulated and observed discharges during the model calibration and validation periods (Figure 7) is to a large extent due to the application of the precipitation correction factor. The multiplicative correction factor to the precipitation removes a substantial proportion bias in both precipitation and simulated discharges. It should be clarified here that such an adjustment of bias in discharge (through the input forcing correction) is an ‘ad hoc’ procedure to overcome the quality of forcing data. The ‘p_corr_factor’ shown in Table 5 indicates that the CORDEX-SMHI was reduced by a factor of 0.4, and ERA-I precipitation was reduced by a factor of 0.41 during the calibration period to minimize bias in discharge estimation. These numbers indicate that there were significant biases in these two forcing datasets and suggest that some kind of bias correction should be considered in any application. We also notice that the precipitation correction factor for WFDEI was near one. This is reasonable and encouraging, as this dataset is already bias-corrected (using observed precipitation). The precipitation correction factor for observed datasets was 1.34, indicating that the areal precipitation is underestimated. This underestimation is likely caused by a combination of undercatch of the precipitation gauges and the lack of representative observation stations. The undercatch depends on the precipitation gauge, wind, and precipitation phase (Mekonnen et al. 2015; Zhao et al. 2015), and the non-representativeness is caused by the limited number of observation stations at high altitudes. Better performance in discharge simulation by the application of a correction factor to precipitation was also shown by previous studies Larson & Peck (1974), Lawrence et al. (2009), and Lakew et al. (2017) and is a standard practice in the hydrologic analysis. The improvement of simulation from all forcing datasets versus prior results presented by Bhattarai et al. (2018) using the HBV model (Bergström 1995) in the same catchment is likely attributed to the distributed nature of Shyft.

Though we achieved high NSE, NSEsqrt, and KGE using WFDEI, ERA-I, and Observed + WFDEI forcing datasets with reasonable Dv (Table 6), the simulations were still not able to capture the highest discharges as seen in the QQ-plots in Figure 9. As compared to high discharge, low discharge events were well captured by all forcing datasets. Local simulated peak discharge during the monsoon season may be attributed to less representativeness of the input precipitation data. Particularly, the lack of ability to capture intense localized precipitation events could be the main reason for an inability to simulate high discharge concurrent with the observations. The CORDEX-SMHI simulation generated relatively higher discharge peaks as compared with the other datasets but has overall higher positive biases during the pre-monsoon season (Figure 3). The CORDEX-SMHI forcing datasets could provide an opportunity for the peak discharge analysis in the Himalayan catchment.

Discussion on the water balance analysis

The catchment water balance component analysis (Table 7) revealed notable differences among the forcing datasets. The discrepancies can, to some degree, be attributed to the uncertainty inherent in each forcing dataset. As mentioned by Nkiaka et al. (2017), different precipitation in each forcing dataset can strongly influence the optimized parameters that control the rates and threshold of hydrological processes in the catchment. Although measured evapotranspiration data are not available for comparison, the percentage of evapotranspiration estimates from WFDEI (24.79%) was found to be similar to the previous study by Sakai et al. (2004), Giertz et al. (2006), Sintondji et al. (2014), and Ragettli et al. (2015) in similar Himalayan catchments. The lowest evapotranspiration (399.5 mm/year) was found for the CORDEX-SMHI forcing dataset, which can be associated with high RH, and lower average shortwave radiation (Figure 3) than other forcing data. Lower evaporation (449.6 mm/year) and lower precipitation (1,771.1 mm/year) from the Observed + WFDEI than WFDEI may partly result from precipitation under-catch in the catchment.

Higher positive storage changes were found for the ERA-I and Observed + WFDEI forcing datasets as compared to CORDEX-SMHI and WFDEI. As compared to Observed + WFDEI forcing dataset, a lower storage change for ERA-I was observed and results from higher evapotranspiration in ERA-I. A change in storage was found smaller than the glacier melt for all forcing datasets. A smaller change in storage is due to the fact that it originates from snow that accumulates over the year and contributes to the glacier mass balance. Smaller changes in storage also indicate that the glacier mass balance was negative over the simulation period. Previous studies by Kulkarni (1992), Cogley et al. (2010), Bolch et al. (2012), Wagnon et al. (2013), and Gurung et al. (2016) also showed that the Himalayan glaciers are experiencing a higher retreat rate in recent decades. The highest glacier melt contribution was observed for ERA-I (19%), while the lowest glacier melt contribution is observed for WFDEI (4%). A study by Gupta et al. (2019) on the Marshyangdi River catchment (with 24% glacier cover area) shows that glacier melt contributes to 11.8% of the total discharge. A relatively lower percentage of snow and glacier melt contribution to the total discharge from the Modi River catchment (with 12% glacier cover area) is presented by Kayastha & Shrestha (2019). Similarly, a study by Nepal (2016) on the Dudhkoshi catchment (with 15% glacier cover area) shows that glacier melt contribution is 17% of the total discharge. A similar percentage of glacier melt contribution to the total discharge for Observed + WFDEI (16%) was observed for the Narayani river catchment, although the glacier cover area was only 8%. As suggested by Bonekamp et al. (2019), differences in glacier melt contribution to total discharge result from differences in meteorological forcing, as we also observed in our study.

Uncertainty in the model simulation and observation

Hydrological projections are subject to considerable uncertainty (Dobler et al. 2012) and are easily affected by various factors, including local and climatic conditions, optimized parameters (Shen et al. 2012), and the quality of forcing data (Teweldebrhan et al. 2018). Different factors also have a varying degree of impact on the discharge simulation. The challenge using a model with many model parameters is that we might get a good model fit but a less robust model for prediction and forecasting due to over-fitting. However, the uncertainty in the simulations will almost certainly be higher due to the increased uncertainty in the parameter values. In our study, the discharge was highly seasonal and particularly precipitation dependent. Inclusion of the precipitation correction factor as a calibrating parameter improves the predictability of the model at the expense of increases to the uncertainty in water balance components like evapotranspiration and snowmelt. However, only a single value of the precipitation correction factor (p_corr_factor) was applied. More complex precipitation correction could include seasonality and impact of orographic effects, though we feel this would result in greater uncertainty. Nonetheless, uncertainties remain.

First, uncertainties in the observed data result from the uneven spatial distribution and few monitoring meteorological stations, which are mainly located in lower elevation regions (Figure 2). Furthermore, many observations are manual. This poses a unique source of uncertainty in measurements, not necessarily greater or less from automated. Furthermore, there are numerous sources of errors associated in the establishment of the river stage and parameterization of rating curves. During model calibration, we assumed that observed discharge is correct. Some of deviations between model simulations and observations might actually be explained by errors in the observations and not in the model or the forcings, to assess the uncertainty in streamflow observations with required detailed knowledge about local river profiles and data used to establish the rating curve and it is outside the scope of this paper. Second, uncertainties might originate from the estimation of AE and other water balance components, which again is influenced by the precipitation estimates and other forcing data used for driving the model. And thirdly, there are the uncertainties in the model formulation itself. The model in this study assumed that the discharge in the river depends on the amount of water stored in the catchment, and we did not consider the impact of industrial and households' water consumption or any other regulation. Despite this, the simulation performance achieved by the model presented here is quite good, though the results are limited to the one catchment in Nepal. A further examination should evaluate whether the selected forcing datasets could be applied to simulated discharge at other mountain regions with longer time periods.

CONCLUSIONS

We aimed to identify the quality of discharge simulation for the Narayani catchment of Nepal based on different forcing datasets. The forcing datasets WFDEI, CORDEX-SMHI, ERA-I, and ground-based observation combined with WFDEI (i.e. Observed + WFDEI) were used to estimate the Naryani catchment discharge on a daily basis for the period of January 2000–December 2009. Because of the uniform coverage and data consistency, global forcing datasets from global and regional climate models were considered an important supplement to station data. In this study, the distributed Shyft hydrologic simulation platform was selected as a discharge simulation tool. The forcing data were interpolated using IDW for generating daily observed precipitation, while a gradient-based Kriging method was used for generating temperature fields matching the model resolution. In this study, Shyft was calibrated for the period of 2000–2004 and validation was done from 2005 to 2009.

Our analysis showed that large differences exist between different forcing datasets particularly in the amount of precipitation. Precipitation from ERA-I and CORDEX-SMHI was unrealistic leading to poor model performance. To improve the existing bias in precipitation particularly over the high mountain regions, further validation and algorithm improvements are required. With the application of a precipitation correction factor (p_corr_factor), a significant amount of bias in precipitation could be mitigated. Still, the performance of the bias-corrected forcing data, i.e. WFDEI, was better than the climate model datasets. Comparing model performance during model calibration and validation periods, relatively higher NSE, NSEsqrt, and KGE with lower Dv were found for the Observed + WFDEI and WFDEI forcing datasets. The water balance analysis shows that higher evapotranspiration with large glacier melt is observed for the ERA-I forcing dataset, while the average evapotranspiration calculated from the WFDEI (24.79%) is similar to previous studies.

Therefore, based on the different results from the Shyft, we conclude that, in the data-poor Himalayan catchment, the WFDEI forcing dataset may be the best choice for water resource planning and hydropower inflow calculations. Discharge simulations resulting from the WFDEI forcing data were particularly promising for hydropower estimation and water resource assessment in data-scarce or ungauged regions. CORDEX-SMHI data captured higher peaks and may be suitable for the peak discharge analysis. However, to use CORDEX-SMHI and ERA-I in ungauged catchments for the water balance analysis, bias correction is required. Further analysis by implementing Shyft in the different Himalayan catchments for different forcing datasets is a recommended further step to assess the regional extensibility of the current results.

FUNDING

This work is a contribution to the Strategic Research Initiative ‘Land Atmosphere Interaction in Cold Environments’ (LATICE) of the University of Oslo and partially supported through the Norwegian Research Council's INDOR program under the Hydrologic sensitivity to Cryosphere-Aerosol interaction in Mountain Processes (HyCAMP) project (NRF No. 222195)

ACKNOWLEDGEMENT

We are thankful to the Department of Hydrology and Meteorology, Government of Nepal for providing observed hydrological and meteorological datasets. We also gratefully acknowledge Special Issue Editor Dr. Kolbjørn Engeland and anonymous reviewer for their comments on the research and manuscript.

SUPPLEMENTARY MATERIAL

The Supplementary Material for this paper is available online at https://dx.doi.org/10.2166/nh.2020.079.

REFERENCES

REFERENCES
Adhikari
D.
2006
Hydropower development in Nepal
.
NRB Economic Review
18
,
70
94
.
Bakken
T. H.
,
Killingtveit
Å.
,
Engeland
K.
,
Alfredsen
K.
&
Harby
A.
2013
Water consumption from hydropower plants – review of published estimates and an assessment of the concept
.
Hydrology and Earth System Sciences
17
(
10
),
3983
4000
.
Bergström
S.
1995
The HBV model
. In:
Computer Models of Watershed Hydrology
(
Singh
V. P.
, ed.).
Water Resources Publications
,
Highlands Ranch, CO
,
USA
.
Berrisford
P.
,
Dee
D. P.
,
Poli
P.
,
Brugge
R.
,
Fielding
K.
,
Fuentes
M.
,
Kbackslashr allberg
P. W.
,
Kobayashi
S.
,
Uppala
S.
&
Simmons
A.
2011
The ERA-Interim Archive Version 2.0
.
ECMWF
,
Shinfield Park
,
Reading
, p.
23
.
Bharti
V.
&
Singh
C.
2015
Evaluation of error in TRMM 3B42V7 precipitation estimates over the Himalayan region
.
Journal of Geophysical Research: Atmospheres
120
(
24
),
12458
12473
.
Bhattarai
S.
,
Zhou
Y.
,
Shakya
N. M.
&
Zhao
C.
2018
Hydrological modelling and climate change impact assessment using HBV light model: a case study of Narayani River Basin, Nepal
.
Nature Environment and Pollution Technology
17
(
3
),
691
702
.
Bhattarai
B. C.
,
Burkhart
J. F.
,
Stordal
F.
&
Xu
C. Y.
2019
Aerosol optical depth over the Nepalese cryosphere derived from an empirical model
.
Frontiers in Earth Science
7
,
178
.
Bolch
T.
,
Kulkarni
A.
,
Kääb
A.
,
Huggel
C.
,
Paul
F.
,
Cogley
J. G.
,
Frey
H.
,
Kargel
J. S.
,
Fujita
K.
,
Scheel
M.
,
Bajracharya
S.
&
Stoffel
M.
2012
The state and fate of Himalayan glaciers
.
Science
336
,
310
314
.
Bonekamp
P. N.
,
de Kok
R. J.
,
Collier
E.
&
Immerzeel
W. W.
2019
Contrasting meteorological drivers of the glacier mass balance between the Karakoram and central Himalaya
.
Frontiers in Earth Science
7
,
107
.
Burkhart
J. F.
,
Helset
S.
,
Abdella
Y. S.
&
Lappegard
G.
2016
Operational Research: Evaluating Multimodel Implementations for 24/7 Runtime Environments. American Geophysical Union, Fall General Assembly 2016, Abstract ID H51F-1541
.
Channan
S.
,
Collins
K.
&
Emanuel
W. R.
2014
Global Mosaics of the Standard MODIS Land Cover Type Data
, Vol.
30
.
University of Maryland and the Pacific Northwest National Laboratory
,
College Park, MD
,
USA
.
Chowdhury
M. D. R.
&
Ward
N.
2004
Hydro-meteorological variability in the greater Ganges–Brahmaputra–Meghna basins
.
International Journal of Climatology
24
(
12
),
1495
1508
.
Cogley
J. G.
,
Kargel
J. S.
,
Kaser
G.
,
van der Veen
C. J.
,
Paul
F.
,
Cogley
J. G.
,
Frey
H.
,
Kargel
J. S.
,
Fujita
K.
,
Scheel
M.
,
Bajracharya
S.
&
Stoffel
M.
2010
Tracking the source of glacier misinformation
.
Science
327
(
5965
),
522
.
Conrad
O.
,
Bechtel
B.
,
Bock
M.
,
Dietrich
H.
,
Fischer
E.
,
Gerlitz
L.
,
Wehberg
J.
,
Wichmann
V.
&
Böhner
J.
2015
System for automated geoscientific analyses (SAGA) v. 2.1.4
.
Geoscientific Model Development Discussions
8
(
2
),
2271
2312
.
Dahal
R. K.
&
Hasegawa
S.
2008
Representative rainfall thresholds for landslides in the Nepal Himalaya
.
Geomorphology
100
(
3–4
),
429
443
.
Dahri
Z. H.
,
Ludwig
F.
,
Moors
E.
,
Ahmad
B.
,
Khan
A.
&
Kabat
P.
2016
An appraisal of precipitation distribution in the high-altitude catchments of the Indus basin
.
The Science of the Total Environment
548–549
,
289
306
.
Daly
C.
,
Neilson
R. P.
&
Phillips
D. L.
1994
A statistical-topographic model for mapping climatological precipitation over mountainous terrain
.
Journal of Applied Meteorology
33
(
2
),
140
158
.
Dee
D. P.
,
Uppala
S. M.
,
Simmons
A. J.
,
Berrisford
P.
,
Poli
P.
,
Kobayashi
S.
,
Andrae
U.
,
Balmaseda
M. A.
,
Balsamo
G.
,
Bauer
P.
,
Bechtold
P.
,
Beljaars
A. C. M.
,
van de Berg
L.
,
Bidlot
J.
,
Bormann
N.
,
Delsol
C.
,
Dragani
R.
,
Fuentes
M.
,
Geer
A. J.
,
Haimberger
L.
,
Healy
S. B.
,
Hersbach
H.
,
Hólm
E. V.
,
Isaksen
L.
,
Kållberg
P.
,
Köhler
M.
,
Matricardi
M.
,
McNally
A. P.
,
Monge-Sanz
B. M.
,
Morcrette
J.-J.
,
Park
B.-K.
,
Peubey
C.
,
de Rosnay
P.
,
Tavolato
C.
,
Thépaut
J.-N.
&
Vitart
F.
2011
The ERA-Interim reanalysis: configuration and performance of the data assimilation system
.
Quarterly Journal of the Royal Meteorological Society
137
(
656
),
553
597
.
Dehecq
A.
,
Gourmelen
N.
,
Gardner
A. S.
,
Brun
F.
,
Goldberg
D.
,
Nienow
P. W.
,
Berthier
E.
,
Vincent
C.
,
Wagnon
P.
&
Trouvé
E.
2019
Twenty-first century glacier slowdown driven by mass loss in High Mountain Asia
.
Nature Geoscience
12
(
1
),
22
27
.
DeWalle
D. R.
&
Rango
A.
2008
Principles of Snow Hydrology
.
Cambridge University Press
,
Cambridge
.
Dobler
C.
,
Hagemann
S.
,
Wilby
R. L.
&
Stätter
J.
2012
Quantifying different sources of uncertainty in hydrological projections in an Alpine watershed
.
Hydrology and Earth System Sciences
16
,
4343
4360
.
Duan
Q.
,
Sorooshian
S.
&
Gupta
V.
1992
Effective and efficient global optimization for conceptual rainfall-runoff models
.
Water Resources Research
28
(
4
),
1015
1031
.
Edenhofer
O.
,
Pichs-Madruga
R.
,
Sokona
Y.
,
Seyboth
K.
,
Kadner
S.
,
Zwickel
T.
,
Eickemeier
P.
,
Hansen
G.
,
Schlömer
S.
,
von Stechow
C.
&
Matschoss
P.
2011
Renewable Energy Sources and Climate Change Mitigation: Special Report of the Intergovernmental Panel on Climate Change
.
Cambridge University Press
,
Cambridge
.
Engeland
K.
,
Steinsland
I.
,
Johansen
S. S.
,
Petersen-Øverleir
A.
&
Kolberg
S.
2016
Effects of uncertainties in hydrological modelling. A case study of a mountainous catchment in Southern Norway
.
Journal of Hydrology
536
,
147
160
.
Essou
G. R. C.
,
Sabarly
F.
,
Lucas-Picher
P.
,
Brissette
F.
&
Poulin
A.
2016
Can precipitation and temperature from meteorological reanalyses be used for hydrological modeling?
Journal of Hydrometeorology
17
(
7
),
1929
1950
.
Giertz
S.
,
Jaeger
A.
,
Schopp
M.
,
Diekkrügerdiekkr¨diekkrüger
B.
&
Jaeger
A.
2006
An interdisciplinary scenario analysis to assess the water availability and water consumption in the Upper Ouémé catchment in Benin
.
Advances in Geosciences
9
(
2006
),
3
13
.
Giorgi
F.
,
Jones
C.
&
Asrar
G. R.
2009
Addressing climate information needs at the regional level: the CORDEX framework
.
World Meteorological Organization (WMO) Bulletin
58
(
3
),
175
195
.
Gupta
H. V.
,
Kling
H.
,
Yilmaz
K. K.
,
Martinez
G. F.
&
Kling
H.
2009
Decomposition of the mean squared error and NSE performance criteria: implications for improving hydrological modelling
.
Journal of Hydrology
377
,
80
91
.
Gupta
A.
,
Kayastha
R. B.
,
Ramanathan
A. L.
&
Dimri
A. P.
2019
Comparison of hydrological regime of glacierized Marshyangdi and Tamor river basins of Nepal
.
Environmental Earth Sciences
78
(
14
),
427
.
Gurung
S.
,
Bhattrai
B. C.
,
Kayastha
R. B.
,
Stumm
D.
,
Joshi
S.
&
Mool
P. K.
2016
Study of annual mass balance (2011–2013) of Rikha Samba Glacier, Hidden Valley, Mustang, Nepal
.
Sciences in Cold and Arid Regions
8
(
4
),
0311
0318
.
Hegdahl
T. J.
,
Tallaksen
L. M.
,
Engeland
K.
,
Burkhart
J. F.
&
Xu
C.-Y.
2016
Discharge sensitivity to snowmelt parameterization: a case study for Upper Beas basin in Himachal Pradesh
.
Hydrology Research
47
(
4
),
683
700
.
Hock
R.
2003
Temperature index melt modelling in mountain areas
.
Journal of Hydrology
282
(
1–4
),
104
115
.
Iizumi
T.
,
Okada
M.
&
Yokozawza
M.
2014
A meteorological forcing data set for global crop modeling: development, evaluation, and intercomparison
.
Journal of Geophysical Research: Atmospheres
119
(
2
),
363
384
.
Immerzeel
W. W.
,
van Beek
L. P. H.
&
Bierkens
M. F. P.
2010
Climate change will affect the Asian water towers
.
Science
328
(
5984
),
1382
1385
.
Jain
K.
,
Jain
S. K.
,
Jain
N.
&
Xu
C.-Y.
2017
Hydrologic modeling of a Himalayan mountain basin by using the SWAT model
.
Journal Hydrology and Earth System Sciences
.
https://doi.org/10.5194/hess-2017-100
Kansakar
S. R.
,
Hannah
D. M.
,
Gerrard
J.
&
Rees
G.
2004
Spatial pattern in the precipitation regime of Nepal
.
International Journal of Climatology
24
,
1645
1659
.
Karki
R.
2010
Status of automatic weather stations in Nepal and comparison of air temperature and precipitation data between automatic weather station and manual observation
. In
IOM-104
. .
Kauffeldt
A.
,
Wetterhall
F.
,
Pappenberger
F.
,
Salamon
P.
&
Thielen
J.
2016
Technical review of large-scale hydrological models for implementation in operational flood forecasting schemes on continental level
.
Environmental Modelling & Software
75
,
68
76
.
Kayastha
R. B.
&
Shrestha
A.
2019
Snow and ice melt contribution in the daily discharge of Langtang and Modi Rivers, Nepal
. In:
Environmental Change in the Himalayan Region
(S. Anup & T. Pankaj Thapa, eds).
Springer
,
Cham
, pp.
1
21
.
Khanal
N. R.
,
Mool
P. K.
,
Shrestha
A. B.
,
Rasul
G.
,
Ghimire
P. K.
,
Shrestha
R. B.
&
Joshi
S. P.
2015
A comprehensive approach and methods for glacial lake outburst flood risk assessment, with examples from Nepal and the transboundary area
.
International Journal of Water Resources Development
31
(
2
),
219
237
.
Kirchner
J. W.
2009
Catchments as simple dynamical systems: catchment characterization, rainfall-runoff modeling, and doing hydrology backward
.
Water Resources Research
45
(
2
).
https://doi.org/10.1029/2008WR006912
Kripalani
R. H.
&
Sontakke
N. A.
1996
Rainfall variability over Bangladesh and Nepal: comparison and connections with features over India
.
International Journal of Climatology
16
,
689
703
.
Kulkarni
A. V.
1992
Mass balance of Himalayan glaciers using AAR and ELA methods
.
Journal of Glaciology
38
(
128
),
101
104
.
Kumar
D.
2017
River Ganges – historical, cultural and socioeconomic attributes
.
Aquatic Ecosystem Health & Management
20
(
1–2
),
8
20
.
Lambert
A.
1972
Catchment models based on ISO-functions
.
Journal of the Institution of Water Engineers
26
,
413
422
.
Larson
L. W.
&
Peck
E. L.
1974
Accuracy of precipitation measurements for hydrologic modeling
.
Water Resources Research
10
(
4
),
857
863
.
Lawrence
D.
,
Haddeland
I.
&
Langsholt
E.
2009
Calibration of HBV hydrological models using PEST parameter estimation
.
NVE Report
.
Available from: www.nve.no
.
Li
H.
,
Xu
C.-Y.
,
Beldring
S.
,
Tallaksen
L. M.
&
Jain
S. K.
2016
Water resources under climate change in Himalayan basins
.
Water Resources Management
30
(
2
),
843
859
.
doi:10.1007/s11269-015-1194-5
.
Li
H.
,
Haugen
J. E.
&
Xu
C.-Y.
2018
Precipitation pattern in the Western Himalayas revealed by four datasets
.
Hydrology and Earth System Sciences
22
,
5097
5110
.
Li
L.
,
Shen
M. X.
,
Hou
Y. K.
,
Xu
C.-Y.
,
Lutz
A. F.
,
Chen
J.
,
Jain
S. K.
,
Li
J. J.
&
Chen
H.
2019
Twenty-first-century glacio-hydrological changes in the Himalayan headwater Beas River basin
.
Hydrology and Earth System Sciences
23
,
1483
1503
.
Lombraña
J. U.
2017
Evaluation of Snow Simulations in SHyFT
.
Norwegian University of Science and Technology (NTNU)
.
Martinec
J.
,
Rango
A.
,
Roberts
R.
&
Baumgartner
M. F.
1998
Snowmelt Runoff Model (SRM) User's Manual
.
University of Berne, Department of Geography
,
Berne
.
Matt
F. N.
,
Burkhart
J. F.
&
Pietikäinen
J.-P.
2018
Modelling hydrologic impacts of light absorbing aerosol deposition on snow at the catchment scale
.
Hydrology and Earth System Sciences
22
,
179
201
.
Mekonnen
G. B.
,
Matula
S.
,
Doležal
F.
&
Fišák
J.
2015
Adjustment to rainfall measurement undercatch with a tipping-bucket rain gauge using ground-level manual gauges
.
Meteorology and Atmospheric Physics
127
(
3
),
241
256
.
Ménégoz
M.
,
Gallée
H.
&
Jacobi
H. W.
2013
Precipitation and snow cover in the Himalaya: from reanalysis to regional climate simulations
.
Hydrology and Earth System Sciences
17
,
3921
3936
.
Moulin
L.
,
Gaume
E.
&
Obled
C.
2009
Uncertainties on mean areal precipitation: assessment and impact on streamflow simulations
.
Hydrology and Earth System Sciences
1
3,
99
114
.
Nayava
J. L.
1974
Heavy monsoon rainfall in Nepal
.
Weather
29
(
12
),
443
450
.
Nayava
J. L.
1980
Rainfall in Nepal
.
Himalayan Review
12
,
1
18
.
Nkiaka
E.
,
Nawaz
N.
,
Lovett
J.
,
Nkiaka
E.
,
Nawaz
N. R.
&
Lovett
J. C.
2017
Evaluating global reanalysis datasets as input for hydrological modelling in the Sudano-Sahel region
.
Hydrology
4
(
1
),
13
.
Oliver
M. A.
&
Webster
R.
1990
Kriging: a method of interpolation for geographical information systems
.
International Journal of Geographical Information System
4
(
3
),
313
332
.
Omani
N.
,
Srinivasan
R.
,
Karthikeyan
R.
,
Smith
P.
,
Omani
N.
,
Srinivasan
R.
,
Karthikeyan
R.
&
Smith
P. K.
2017
Hydrological modeling of highly glacierized basins (Andes, Alps, and Central Asia)
.
Water
9
(
2
),
111
.
Pellicciotti
F.
,
Buergi
C.
,
Immerzeel
W. W.
,
Konz
M.
&
Shrestha
A. B.
2012
Challenges and uncertainties in hydrological modeling of remote Hindu Kush–Karakoram–Himalayan (HKH) Basins: suggestions for calibration strategies
.
Mountain Research and Development
32
(
1
),
39
50
.
Peña-Arancibia
J. L.
,
Zhang
Y.
,
Pagendam
D. E.
,
Viney
N. R.
,
Lerat
J.
,
van Dijk
A. I. J. M.
,
Vaze
J.
&
Frost
A. J.
2015
Streamflow rating uncertainty: characterisation and impacts on model calibration and performance
.
Environmental Modelling & Software
63
,
32
44
.
Pradhananga
N. S.
,
Kayastha
R. B.
,
Bhattarai
B. C.
,
Adhikari
T. R.
,
Pradhan
S. C.
,
Devkota
L. P.
,
Shrestha
A. B.
&
Mool
P. K.
2014
Estimation of discharge from Langtang River basin, Rasuwa, Nepal, using a glacio-hydrological model
.
Annals of Glaciology
55
(
66
),
223
230
.
Priestley
C. H. B.
&
Taylor
R. J.
1972
On the assessment of surface heat flux and evaporation using large-scale parameters
.
Monthly Weather Review
100
(
2
),
81
92
.
QGIS Development Team
2016
QGIS Geographic Information System
.
Open Source Geospatial Foundation
.
Retrieved from: http://qgis.osgeo.org
.
Ragettli
S.
,
Pellicciotti
F.
,
Immerzeel
W. W.
,
Miles
E. S.
,
Petersen
L.
,
Heynen
M.
,
Shea
J. M.
,
Stumm
D.
,
Joshi
S.
&
Shrestha
A.
2015
Unraveling the hydrology of a Himalayan catchment through integration of high resolution in situ data and remote sensing with an advanced simulation model
.
Advances in Water Resources
78
,
94
111
.
Raimonet
M.
,
Oudin
L.
,
Thieu
V.
,
Silvestre
M.
,
Vautard
R.
,
Rabouille
C.
&
Moigne
P. L. E.
2017
Evaluation of gridded meteorological datasets for hydrological modeling
.
Journal of Hydrometeorology
18
(
11
),
3027
3041
.
Renard
B.
,
Kavetski
D.
,
Kuczera
G.
,
Thyer
M.
&
Franks
S. W.
2010
Understanding predictive uncertainty in hydrologic modeling: the challenge of identifying input and structural errors
.
Water Resources Research
46
,
5
.
Rochester
R. E. L.
2010
Uncertainty in Hydrological Modelling: A Case Study in the Tern Catchment
.
UCL
,
Shropshire
,
UK
. .
Sakai
A.
,
Fujita
K.
&
Kubota
J.
2004
Evaporation and percolation effect on melting at debris-covered Lirung Glacier, Nepal Himalayas, 1996
.
Bulletin of Glacier Research
21
,
9
15
.
Sanjay
J.
,
Krishnan
R.
,
Shrestha
A. B.
,
Rajbhandari
R.
&
Ren
G. Y.
2017
Downscaled climate change projections for the Hindu Kush Himalayan region using CORDEX South Asia regional climate models
.
Advances in Climate Change Research
8
(
3
),
185
198
.
Sapkota
A.
2016
Regional Modelling of the Narayani Basin in Nepal
.
Norwegian University of Science and Technology, Trondheim
.
Seibert
J.
2001
On the need for benchmarks in hydrological modelling
.
Hydrological Processes
15
(
6
),
1063
1064
.
Seiller
G.
,
Anctil
F.
&
Perrin
C.
2012
Multimodel evaluation of twenty lumped hydrological models under contrasted climate conditions
.
Hydrology and Earth System Sciences
16
(
4
),
1171
1189
.
Shepard
D.
1968
A Two-Dimensional Interpolation Function for Irregularly-Spaced Data
.
Proceedings of the 1968 ACM National Conference, New York, 27–29 August 1968
, pp.
517
524
.
Shrestha
A. B.
,
Wake
C. P.
,
Mayewski
P. A.
,
Dibb
J. E.
,
Shrestha
A. B.
,
Wake
C. P.
,
Mayewski
P. A.
&
Dibb
J. E.
1999
Maximum temperature trends in the Himalaya and its vicinity: an analysis based on temperature records from Nepal for the period 1971–94
.
Journal of Climate
12
(
9
),
2775
2786
.
Sintondji
L. O.
,
Zokpodo
B.
,
Ahouansou
D. M.
,
Vissin
W. E.
&
Agbossou
K. E.
2014
Modelling the water balance of Ouémé catchment at the Savé outlet in Benin: contribution to the sustainable water resource management
.
International Journal of Agriscience
4
,
74
88
.
Teweldebrhan
A. T.
,
Burkhart
J. F.
&
Schuler
T. V.
2018
Parameter uncertainty analysis for an operational hydrological model using residual-based and limits of acceptability approaches
.
Hydrology and Earth System Sciences
22
,
5021
5039
.
Wagnon
P.
,
Vincent
C.
,
Arnaud
Y.
,
Berthier
E.
,
Vuillermoz
E.
,
Gruber
S.
,
Ménégoz
M.
,
Gilbert
A.
,
Dumont
M.
,
Shea
J. M.
,
Stumm
D.
&
Pokhrel
B. K.
2013
The cryosphere seasonal and annual mass balances of Mera and Pokalde glaciers (Nepal Himalaya) since 2007
.
The Cryosphere
7
,
1769
1786
.
Weedon
G. P.
,
Gomes
S.
,
Viterbo
P.
,
Shuttleworth
W. J.
,
Blyth
E.
,
Österle
H.
,
Adam
J. C.
,
Bellouin
N.
,
Boucher
O.
&
Best
M.
2011
Creation of the WATCH forcing data and its use to assess global and regional reference crop evaporation over land during the twentieth century
.
Journal of Hydrometeorology
12
(
5
),
823
848
.
Weedon
G. P.
,
Balsamo
G.
,
Bellouin
N.
,
Gomes
S.
,
Best
M. J.
&
Viterbo
P.
2014
The WFDEI meteorological forcing data set: WATCH forcing data methodology applied to ERA-Interim reanalysis data
.
Water Resources Research
50
(
9
),
7505
7514
.
Xu
H.
,
Xu
C. Y.
,
Sælthun
N. R.
,
Zhou
B.
&
Xu
Y.
2015
Evaluation of reanalysis and satellite-based precipitation datasets in driving hydrological models in a humid region of Southern China
.
Stochastic Environmental Research and Risk Assessment
29
(
8
),
2003
2020
.
Zhang
Y. G.
,
Hao
Z. C.
,
Xu
C.-Y.
&
Lai
X. D.
2019
Response of melt water and rainfall runoff to climate change and their roles in controlling streamflow changes of the two upstream basins over the Tibetan Plateau
.
Hydrology Research
.
doi: 10.2166/nh.2019.075
.

Supplementary data