Effect of baseline snowpack assumptions in the HySIM model in predicting future hydrological behaviour of a Himalayan catchment

In mountainous and glacial catchments, the magnitude and seasonality of river discharge are greatly influenced by snowpack development and behaviour (Ye et al. 2003). Hydrological modelling in such mountainous regions is, however, significantly constrained with model biases and input data uncertainties regarding snow and ice reserves (Fontaine et al. 2002); thus affecting the reliability of snow/glacier melt-runoff predictions and water resources and flood risk management. For example, snow depth is an important parameter for hydrological simulation (Bell et al. 2016), especially for peak floods caused by snowmelt (Bergeron et al. 2014). Accurate observation and measurement of snowpack and glaciers at the catchment or river basin scale, including data on water-equivalent storage, temporal depth–area relationships, albedo, etc., is challenging in such poorly instrumented, remote and often inaccessible environments. Given the importance of these properties for physically based, spatially distributed hydrological modelling, modellers may rely on the increasing array of remote sensing data sets that are becoming available, such as the NOAA Geostationary Orbiting Environmental Satellite (GOES) or MODerate resolution Imaging Spectroradiometer (MODIS) on the NASA Terra satellites (Schmugge et al. 2002), along with snow water equivalent (SWE) retrieval methods and algorithms (Dai et al. 2012).

However, practitioners still use various indirect observations, assumptions, expert knowledge and inverse modelling skills for conceptual modelling in ‘melt-dominated’ catchments. These include assumptions of linear relationships between terrain elevation and snow depth, estimation of SWE through evaluation of snow-covered area (Liston 1999), use of passive microwave data (Dong et al. 2005) and incorporation of interpolation of ground-based measurements with remote sensing data (Kelly et al. 2003; Kelly 2009), but none of these approaches are free from uncertainty in actual SWE storage.

Despite these challenges in defining snow boundary conditions, conceptual hydrological models with snowmelt components, using either energy budget (e.g., Anderson 1976), temperature-index or degree-day methods for snow accumulation and melting (e.g., Martinec et al. 1983) are widely used in snow-dominated catchments. Examples include HBV (Akhtar et al. 2008), SRM (Kult et al. 2012), GERM (Farinotti et al. 2012), GSM-SOCONT (Schaefli & Huss 2011), J-2000 (Nepal et al. 2014), OEZ (Hagg et al. 2007), SNOWMOD (Singh et al. 2006), UBC Watershed (Loukas et al. 2002), WATFLOOD (Comeau et al. 2009) and HySIM (Manley and Water Resource Associates Ltd 2006). Glacier mass balance is a major climate indicator in glacio-hydrological modelling using index-based approaches in conceptual models (e.g., TOPKAPI-ETH, SRM-model, etc.). Previous work in the Alps in Europe has evaluated the influence of snow and ice reserves on simulated hydrological behaviour, demonstrating the importance of prior knowledge on ice-thickness distribution, total ice volume amount and spatial distribution of winter snow and glacier area change in constraining snow and ice melt for runoff projections (e.g., Bavay et al. 2013; Huss et al. 2014), compared to assuming unlimited snow/glacier melt. In this study we have selected the HySIM model due to its: (i) capability to simulate surface water resources in both data-limited and snow-dominated environments; (ii) ability to simulate daily surface runoff, percolation to groundwater and river flow; (iii) in-built multiple parameter optimization with four objective functions; (iv) inclusion of advanced hydraulic and hydrologic parameters in the model; and (v) graphics facility for rapid visualization of inputs and outputs.

It is common within climate change studies on glacier-dominated catchments for glacier extent to be assumed to be in steady state (Hagg et al. 2013). However, Hagg et al. (2007) consider that treating glaciers as static ice bodies is a major limitation in climate change impact studies. Unsteady state and seasonal snow model setups may lead to enhanced hydrological seasonality due to larger direct runoff in wet seasons and lower glacier or snowmelt in dry seasons (Juen et al. 2007).

The purpose of this study is to explore the implications of model-setup assumptions regarding the behaviour of historical snowpack reserves for water resources simulation in the Himalayan Beas river basin of north India and to see how the model assumptions on snow and glaciers make an impact on hydrological model simulations under changing climate. The two research questions are: (i) How do differing assumptions on baseline snowpack behaviour affect model parameter values and performance? (ii) How do the baseline snowpack assumptions affect simulated future climate change impacts?

Hydrological modellers often ignore inadequacies in the input data or data assumptions which can affect baseline/future model performance, instead focusing more on the quantification of model parameters. Therefore, in this study, the consequences of such baseline set-up assumptions for simulated future hydrological response to changed temperature and precipitation under two future snowpack scenarios were assessed. The paper provides valuable guidance for conceptual hydrological modelling in snow-dominated regions, considering the effects of data inadequacies in snowpack characterization and its implications in climate change-related hydrological impact studies. However, it must be noted that the purpose of this paper is not to predict the time scale of glacial retreat or loss, recognizing the difficulty and uncertainty in estimating rates of loss of glacial area or volume in the Himalayas and Tibetan Plateau (Yao et al. 2012).

This study focuses on the perennial River Beas in the north-western Himalayan region of India. The river basin upstream of Pong Dam, bounded by latitude 31°28′–32°26′N and longitude 75°56′–77°48′E, has a drainage area of 12,560 km² (Figure 1) which varies in elevation from 245 to 6,617 m above sea level. The major land cover classes include forest, glacier and bare rock, with about 65% of the area covered with snow during winter (Singh & Bengtsson 2003). The Beas basin receives around 70% of the annual rainfall during the summer monsoon between June and September, and the basin is characterized by mean minimum and mean maximum winter temperatures of −1.6 °C and 7.7 °C, respectively (Singh & Ganju 2008) although there are significant spatial differences due to the elevation range. Daily gauged inflows to the Pong reservoir are available from January 1998 to December 2008 (11 years) used for this study.

This study has used Tropical Rainfall Measuring Mission (TRMM) 3B42V7 daily gridded precipitation data for 1998 to 2008 to represent the precipitation variability of the catchment. The TRMM data, which are widely used in environmental and hydrological research, have a spatial resolution of 0.25° × 0.25° and covers the latitudinal band of 50° N-S. TRMM data sets have been used in studies of many catchments in the Himalayan region (Bookhagen & Burbank 2010). Moreover, Xue et al. (2013) demonstrated that 3B42V7 provides better basin-scale agreement with observed (2001–2010) monthly and daily rain gauge data and improved rainfall intensity distribution than the earlier 3B42V6.

Daily reference evapotranspiration (ET₀) was calculated using the FAO Penman–Monteith method (Allen et al. 1998) from the 0.312 degree (∼38 km) gridded meteorological variables (daily maximum temperature, daily minimum temperature, daily wind velocity, daily average relative humidity and daily average solar radiation) from the NCEP Climate Forecast System Reanalysis (CFSR) data for 1998–2008. The standard normal probability density functions and standard normal cumulative density functions of ET₀ for the three sub-basins (Figure 2) demonstrate that the CFSR data adequately represent the influence of elevation.

HySIM is a continuous, daily, conceptual rainfall–runoff model that has been extensively used in mountainous catchments including in climate change studies (Wilby 2005). The HySIM model contains two sub-routines for simulating the river basin hydrology and channel hydraulics. The hydrology in each sub-basin is simulated using seven storages representing snow, vegetation, soil layers, unsaturated and saturated zones, while the channel hydraulic sub-routine uses kinematic routing flows within sub-basins. Details of the model parameters are given in Pilling & Jones (1999). The empirical degree-day approach is used to calculate daily snowmelt or accumulation. HySIM uses precipitation, ET₀ and a temperature-based snowmelt model to simulate streamflow.

The Beas river basin was sub-divided into three sub-catchments on the basis of the river network and elevation (Figure 1) to appropriately represent the catchment structure/river network in such a conceptual model. These represent the glacier and permanent snow-covered areas of the basin (Upper), the seasonal snow cover areas (Middle) and the remaining lower elevation parts of the catchment that receive rainfall-only (Lower). The upper, middle and lower basins areas are 5,720, 3,440 and 3,350 km², respectively. Areal averaging was used to transform the precipitation and evapotranspiration data from the different resolution grids. Soil parameters were initialized based on spatially weighted values from the Harmonized World Soil Database (FAO/IIASA/ISRIC/ISSCAS/JRC 2012) although model values for the soil hydraulic parameters were calibrated.

In this study, seven parameters were calibrated: rooting depth (mm) [RD], permeability – horizon boundary (mm hr⁻¹) [PHB], permeability – base lower horizon (mm hr⁻¹) [PBLH], interflow rate within the upper soil layer (mm hr⁻¹) [IU], interflow rate within the lower soil layer (mm hr⁻¹) [IL], snow temperature threshold (°C) [ST] and snowmelt rate for each degree of temperature above the threshold (mm day⁻¹°C⁻¹) [SM]. The permanent Himalayan snow and ice cover within the upper sub-catchment in the model was initialized with an ice/snow depth of 25 m informed by past research (Kulkarni et al. 2005; Kulkarni & Karyakarte 2014). HySim uses the commonly applied empirical degree-day approach to simulate snowmelt and accumulation (Pilling & Jones 1999). In this model, when the mean air temperature (T) falls below snow threshold (ST) parameter, any precipitation will be assumed as snow and added to snow storage. Similarly, daily potential snowmelt (M_s, mm/day) will be released from snow storage when the mean air temperature (T) exceeded base temperature (T_b, assumed as 0 °C), as indicated by snowmelt rate parameter (SM, this is similar to degree-day factor) based on equation [M_s = SM × (T − T_b)].

The methodology adopted in this study is shown in Figure 3.

1. Two separate HySIM models were built using different assumptions as to the baseline behaviour of the glacier and permanent snowpack in the Upper sub-basin.

   • a.

     Model Assumption 1: constrained snowmelt [glacier and permanent snowpack are in steady state]

     In this model, the annual snowmelt in the Upper sub-basin is limited to the seasonal snow accumulation so that there is no net loss of snow and ice reserves over the baseline period.

   • b.

     Model Assumption 2: unconstrained snowmelt [glacier and permanent snowpack are in an unsteady state]

     In this model, the snowmelt is only energy-limited so that there can be net loss (or gain) of snow and ice reserves over the baseline period. The model was initialized with a snowpack depth of 25 m in the Upper sub-catchment to ensure that annual snowmelt is not limited by snowpack availability.

2. Both models had a 2-year warm-up period (1998–1999) and then were independently calibrated, by modifying the selected parameter values, against observed discharge data (2000–2004, inclusive) and validated (2005–2008, inclusive). The automated optimization procedure in HySIM, using the reduced error of estimate technique, was used to generate one hundred parameter sets for each model within the parameter uncertainty range. The hydrological performances were then evaluated by comparing the simulated and measured daily discharge using the Nash–Sutcliffe efficiency criterion (NSE, Nash & Sutcliffe 1970) and the percent bias (PBIAS) goodness-of-fit measures, as recommended by Moriasi et al. (2007).

3. After identifying the two parameter sets which provided the best calibration/validation performance for the two contrasting models, the parameter values in both models were fixed for all future simulations.

4. A climate ‘scenario neutral framework’ was set up that assesses the catchment response to a plausible range of future climate changes, while avoiding the application of time varying GCM/RCM scenarios simulated under particular assumptions of social/economic/environmental policies (Prudhomme et al. 2010; Remesan & Holman 2015). Ranges of future annual temperature and annual precipitation changes were informed by the regional summary results from 25–39 GCMs given in Christensen et al. (2013). Six temperature change factors between ΔT = 0 °C and ΔT = +5 °C (in steps of 1 °C) and seven precipitation change factors from ΔP = −10% to ΔP = +20% (in steps of 5%) were used. Each absolute temperature change factor was added to the historical NCEP data to provide modified temperature and subsequently ET₀ time series, assuming all other weather variables were unchanged. The relative changes in precipitation were applied to the TRMM historical time series.

5. The two calibrated/validated HySIM models were then initialized with two future contrasting snow/glacier scenarios:

   • a.

     Medium-term scenario (M): Permanent snow/ice reserves persist in the Upper sub-basin, with future snowmelt unconstrained by snow and ice reserves;

   • b.

     Long-term scenario (L): Insignificant permanent snow/ice reserves remain in the Upper sub-basin so that future annual snowmelt is limited to the seasonal accumulation of snow.

We selected these two scenarios to represent two realistic but distinct behavioural systems for the Himalayas – a nearer-term scenario (which we have termed ‘Medium’), in which the permanent snow/ice reserves persist in the Upper sub-basin and a longer-term future scenario in which there are insignificant permanent snow/ice reserves remaining in the Upper sub-basin. Given the major uncertainties in the future temporal responses of the Himalayan glaciers to climate change, we have deliberately not ascribed time horizons to our two scenarios. Where appropriate, we have used the notations of 1 M, 2 M, 1 L and 2 L to denote the different baseline model assumptions and future scenarios.

• 6.

  The four model builds (each combination of calibrated/validated Model 1 and Model 2 with the medium- and long-term future snowpack scenarios) were then each run for the 42 combinations of changed temperature (6) and precipitation (7), with IRS produced of the changes in daily low flows [Q90] and high flows [Q10] under these plausible future changes in temperature and precipitation, compared to those from 2000 to 2008 (i.e., with zero temperature/precipitation change). The daily flow that is equalled or exceeded 10 and 90% of the time (annual high flow (Q10) and annual low flow (Q90)), respectively, are commonly used in designing hydropower projects (IITR 2011).

This section presents and discusses the effect of the conceptual model structure assumptions of historical snowpack behaviour on the calibrated parameter values and baseline hydrological performance of the two models in the Beas river basin. Then the influence of the resultant baseline parameterization on simulated climate change impacts in the medium- and long-term snowpack scenarios are presented and discussed. It is seen that the glaciers in the Himalayas (excluding the Karakorum) are believed to be more sensitive to future warming than those in the Tibetan Plateau due to seasonal change of snow accumulation (Yao et al. 2012). This is consistent with the estimated current relative rates of loss in the glaciers of Himachal Pradesh (Kulkarni et al. 2007; Yao et al. 2012).

The time series of observed and modelled flow for both constrained and unconstrained models are given in Figure 4, which demonstrates that both models can describe the observed Beas river basin hydrology with the contrasting snowpack assumptions. The calibration and validation results of the simulated daily discharge in Table 1 show that both models give almost identical NSE values for both periods, which fall within the ‘Very good’ model performance range of Henriksen et al. (2003) for the calibration period and on or within 0.01 of the boundary between ‘Good’ and ‘Very good’ model performance for the validation period. There is greater divergence with the results for PBIAS, which are ‘Excellent’ for Model 1 [constrained snowmelt] for both periods but ‘Very good’ (calibration) and on the ‘Good’ boundary (validation) for Model 2 [unconstrained snowmelt]. However, the Q-Q plot for 2000–2008 in Figure 5 shows that both models under-predict the extreme high flow and low flow events in the Beas river, although this is likely to partially represent the limitations of the TRMM precipitation data. The probability of exceedence curves of observed flow, unconstrained modelled flow and constrained modelled flow are given in Figure 6, which demonstrates that both models reproduce the observed distribution of daily flows reasonably well.

The slight performance disparity between Model 1 and 2 is related to the model assumptions. The differences in the calibrated parameter values between the two models (Table 2) indicate that the calibration process has sought to compensate for the hydrological consequences of the differing snowpack behaviour assumptions between the models, given that input data uncertainties (e.g., in the TRMM daily precipitation data) and model structural deficits are constant across the models. However, hydrological modellers often assume that modelling performance issues are substantively caused by inadequacies in the quantification of model parameters themselves, rather than inadequacies in the input data or data assumptions; but this is not always the case (Mukhopadhyay & Dutta 2010; Remesan & Holman 2015). In such cases, catchment models may have unrecognized limitations in simulating streamflow under scenarios of climate or land-use change that are outside of the calibration conditions, due to what might be thought of as compromised parameter values.

A common assumption in hydrological modelling is that a successfully calibrated and validated hydrological model contains process representation that is suitable for climate change impact studies. This ignores the potential expansion of modelling uncertainty associated with the varying assumptions on initial boundary conditions and resultant parameter values evident in Table 2.

Impact response surfaces have therefore been constructed to understand how the differing parameterization resulting from the contrasting baseline snowpack assumptions translates into uncertainty in the future hydrological response of the River Beas to a plausible range of changes in future annual precipitation and temperature change to 2100 for the region (Christensen et al. 2013) under medium-term (continued presence of permanent snowpack) and long-term (loss of permanent snowpack) conditions.

Both models simulate increasing average annual discharge with increasing temperature and precipitation, with the temperature increase of +5 °C (through its influence on snowmelt) being more able to offset the simulated annual precipitation decrease of −10% (Table 3). The change in average annual discharge is generally larger with Model 2 (unconstrained snowmelt), with the difference between the two models ranging from 0.6 to 1.5 × 10⁹ m³/yr (compared to an observed average annual discharge of 7.8 × 10⁹ m³/yr), while the percentage change difference between the models is up to 45% of the simulated baseline [with ΔT = 0 °C, ΔP = 0%]. Although it has not been considered in this study, it is noted that shrinkage in glacier area and melt duration would impact glacier melt runoff as it can connect to a function of contributing glacier area and melt rate (Bliss et al. 2014).

Figure 7 shows that there is broad similarity in both the Q10 impact response surfaces when the two models are applied to medium-term conditions (i.e., 1 M and 2 M settings) in which future snowmelt is only energy-limited due to a continuing permanent snowpack. In general, Q10 increases in both models with increasing temperature as there is increased snowmelt of the permanent snow/ice reserves over a longer period, although there is a tipping point in Model 1 (constrained snowmelt) at around a temperature increase of 1 °C, above which, increases in snowmelt outweigh the effect of increased actual evapotranspiration. Q10 also tends to increase with increasing precipitation with a similar sensitivity in both models, although the effect decreases with increasing temperature from a range of around 40% of the baseline Q10 for ΔT = +0 °C to 30% for ΔT = +5 °C. Q10 for the unconstrained snowmelt model (Model 2 and 2 M setting) is more sensitive to changing temperature and precipitation than the constrained snowmelt model (Model 1 and 1 M setting), exhibiting generally larger changes (up to an increase of 79% of the baseline Q10, compared to 40% with Model 1) and a larger range of responses across the scenario-neutral climate space (of 90% of the baseline Q10, compared to 58%).

However, there are very significant differences in the responses in Q90 between the two models. In Model 2 (unconstrained snowmelt, 2 M setting), Q90 increases with increasing temperature and precipitation across the scenario-space by up to 133% of the baseline Q90 due to the increasing duration and magnitude of snow/ice melt, although there is a minor inflection line around ΔT = +1 °C. With small temperature increases (<1 °C) and a decrease in precipitation, Q90 decreases by up to −12% arising from decreased recharge due to increased actual ET not being offset by increased snow/glacier melt or precipitation.

In contrast, warming leads to general decreases in Q90 in Model 1 (constrained snowmelt, 1 M setting) across most of the scenario space by up to −73% of the baseline Q90, despite the increase in average annual discharge (Table 3). This arises due to the effect of the calibrated parameterization on the interplay between a number of hydrological processes within the sub-catchments (Table 2). First, the snowmelt rate in the Upper catchment is lower in Model 1 compared to Model 2 (1.0 compared to 1.2 mm/°C/d) leading to smaller increases in dry season snowmelt from this sub-catchment. Second, seasonal snow accumulation in the Middle catchment melts sooner and quicker (due to a lower snowmelt threshold and higher snowmelt rate) exhausting this seasonal reserve before the dry season, but the potential for increased recharge into the groundwater store to support dry season flows is limited by the smaller lower boundary permeability which leads to increased lateral flow. Finally, the Middle catchment in Model 1 has a greater calibrated rooting depth of 2.2 m (compared to 1.2 m in Model 2), so that more water can be extracted from the soil moisture reserves to meet evapotranspirative demand leading to reduced recharge and therefore baseflow.

Across the entire range of temperature and precipitation changes, the long-term future average annual flows are lower than in the medium term for both models (data not shown). This arises as the volume of snowmelt is limited to the seasonal winter precipitation in the Upper and Middle catchments (as there is no diminishing permanent snowpack as in the medium-term future) with the result that the increased temperature increases the evapotranspiration and thus decreases rainfall–runoff and recharge. Both models (i.e., 1 L and 2 L settings) show an increasing average annual discharge with increasing precipitation (compared to [ΔT = 0 °C, ΔP = 0%]), the magnitude of which decreases with increasing temperature (in contrast to the medium-term response) (Table 3). There is a greater range of changes in average annual discharge across the scenario-neutral space with Model 2 (unconstrained snowmelt), with the difference between the two models in the magnitude of the simulated change from the baseline [ΔT = 0 °C, ΔP = 0%] ranging between around 650–1,000 × 10⁶ m³/yr, while the percentage change difference between the models is only up to 15% of the baseline. Both of these impact uncertainties between the two models were lower than in the medium-term future.

There are much more similar hydrological responses in low and high flows for the long-term future across the scenario-neutral space between the two models (Figure 8), than for the medium-term impact response surfaces (Figure 7). Q10 and Q90 both decrease with increasing temperature and decreasing precipitation, as rainfall and seasonal snow accumulation decrease and evapotranspiration increases. As the snowmelt in both the Upper and Middle catchments are limited to the seasonal snow accumulation, the largest decreases in Q10 of −50% (Model 1, 1 L setting) and −62% (Model 2, 2 L setting) are associated with the lowest precipitation (ΔP = −10%) and highest temperature increase (ΔT = +5 °C). Similarly, the largest increases in Q10 of +41% (Model 1) and +51% (Model 2) are associated with the greatest precipitation (ΔP = +20%) and lowest temperature increase (ΔT = +0 °C). Under the highest temperature scenario (ΔT = +5 °C), changes in Q10 range from the most optimistic changes of −2% (Model 1) to −9% (Model 2) for ΔP = +20%, to the most pessimistic changes of −50% (Model 1) to −62% (Model 2) for ΔP = −10%.

As would be expected with the loss of the permanent snow/ice cover in the Upper sub-basin, Q90 decreases with increasing temperature, although Q90 can increase under certain conditions of increasing precipitation. For the most pessimistic combination [ΔT = +5 °C, ΔP = −10%], Q90 decreases by between −68% (Model 1) and −74% (Model 2). This decrease arises from the earlier loss of the seasonal snowmelt contribution (as the duration of the melt period shortens) so that river discharge becomes reliant on the groundwater baseflow. However, the increased evapotranspiration associated with warming leads to decreased recharge in the lower and middle sub-basins and thus to further reduced baseflow. In contrast, the hottest and wettest scenario [ΔT = +5 °C, ΔP = +20%] leads to increases in Q90 of between 42% (Model 1) and 13% (Model 2) as the increased precipitation leads to increased recharge during the monsoon period when soils are at field capacity, outweighing the effects of the increased evapotranspiration during that period.

The long-term uncertainty introduced by the two baseline snowpack behaviour assumptions across the scenario-neutral climate space ranges between −14 to +13% for Q10 and 5 to −31% for Q90. This is about 10% less than the medium-term uncertainty in Q10 (of 36% of the baseline Q10) and reflects the limitation imposed on the hydrological response of the two models by the loss of the permanent snow/glacier reserve and the influence of the limited seasonal snow accumulation.

The modelling results presented have demonstrated that the long-term transition to a Himalayan river whose hydrological response is dominated by rain and the shorter-duration melt of seasonal snow will lead to significant temporal changes in the water balance and flow dynamics. Here, although both of the models show this gross difference between the two future periods, there are important differences in simulated future hydrological response that arise as a consequence of their different baseline parameterization (Table 2). Figure 9 shows how the differences in calibrated parameter values associated with the differing assumptions of baseline snowpack behaviour impact on the change in the flow duration curves for four selected future climates. Also shown is the uncertainty range for the two models for the historical period. The selected scenarios are [ΔP = −10%, ΔT = 0 °C] (highest reduction in precipitation and no change in temperature), [ΔP = −10%, ΔT = 5 °C] (highest reduction in precipitation and highest increasein temperature), [ΔP = +20%, ΔT = 0 °C] (highest increase in precipitation and no change in temperature), [ΔP = +20%, ΔT = 5 °C] (highest increase in precipitation and temperature from baseline climate). These show that precipitation changes in the absence of temperature increases tend to lead to the similar percentage changes in discharge throughout the flow duration curve, whereas temperature increases cause the largest increase in discharge around the 25th to 50th exceedance probability reflecting increases in pre- and post-monsoon snowmelt. However, the largest differences between the responses of the two models to the same change in climate, arising from the assumptions regarding the historical snowpack behaviour, occur in the medium-term future between about the 60th and 100th exceedance probability reflecting the different contributions of snowmelt to dry season flows.

The assumptions regarding the historical snowpack behaviour lead to an uncertainty between the two baseline models across the historical flow regime that is less than 30% of the observed baseline flow (with the exception of the 97th–100th exceedance probability). When the future discharge uncertainty across the range of flow exceedance between the two models for each of the two future periods is compared with this, it is apparent that the future uncertainty only exceeds the historical uncertainty for the medium term under the warming scenario space, indicating that the magnitude of the baseline model uncertainty is not conserved but is magnified under such futures. In contrast, in the longer term, there is no general expansion of the uncertainty due to the behavioural constraints imposed by the lack of a permanent snow/glacier. The results show that, given the current uncertainty in spatio-temporal dynamics of the glaciers and permanent ice and snowfields, the choices that a modeller makes in the baseline model build can lead to differences in the simulated magnitude of future changes in average annual discharge and hydrological response that are significantly larger than the baseline uncertainty.

The results show an average loss of −0.37 m/yr SWE (unconstrained Model 1) and average gain of +0.005 m/yr SWE (constrained Model 2) during 2000–2008 (negative indicates mass loss and positive values indicate mass gain). The total simulated mass balance loss during the 2000–2008 period is −3,323 mm with uncertainty bounds of [−131 to −691 in mm/year SWE corresponding to the year 2005 and year 2008, respectively] for the unconstrained model. Similarly, the constrained model (with no permanent snow assumption) has shown a total mass change of +46 mm (i.e., seasonal accumulation) during 2000–2008 period with uncertainty bounds of [−37 mm/year to +50 mm/year SWE for the year 2002 and year 2004, respectively]. The losses are comparable with other studies in the Himalayan region, e.g., −0.89 m/yr in Shaune Garang basin for 2001–2008 (Kumar et al. 2016), −0.44 ± 0.09 m/yr for Lahaul and Spiti region for 1999–2011 and the regional review of Pritchard (2017). The streamflow contribution from snow/ice melt is shown in Figure 10 for both constrained and unconstrained models. One can clearly see from this figure that unlike the unconstrained model with permanent snow assumption (i.e., Model 2), the snow/ice contribution from the constrained model is relatively less during low flow periods but compensated with other water balance components. In the case of the Model 2 simulation, the total snow/glacier melt is more than in the Model 1 simulation because of our conceptual assumptions. The water balances for the constrained and unconstrained model simulation are represented in Figure 11, indicating how surface runoff and low AET are compensating for the effect of the lack of a permanent snowpack in the constrained model. For better understanding of major water balance gains and losses in Beas basin with different model assumptions see Figure 12. It is clearly seen that maximum snowpack melting is occurring in the hot monsoon season (June, July and August months). Variations in soil moisture changes (in the upper and lower layer) and differences in monthly groundwater gains/losses are clearly observed in these figures with compensation effect on AET on both constrained and unconstrained models.

Given the focus of our study on the model response to the differing baseline snowpack assumptions alone, it must also be acknowledged that there are many other sources of uncertainty which affect the impact of climate change on the hydrological response of glaciated mountainous river basins. First, the simulated changes in flows resulting from the effect of increased temperatures on snow/ice melt in our two extreme cases do not take into account the change in glacier area which may lead to over-estimation of the snowmelt-temperature response, although our results are consistent with the increases in annual average streamflow of 10% and 18% under warming scenarios of 2 °C and 4 °C (with no change in precipitation) simulated by Nepal et al. (2014) in the Himalayan Dudh Kosiriver basin. Second, although there are continuing debates regarding the temperature index approach to calculating snowmelt and accumulation (e.g., Pellicciotti et al. 2012; Hock 2003; Mackintosh et al. 2017), it continues to be a widely used method, particularly in data-poor basins. Third, although more sophisticated methods exist, the change factor method used to perturb the historical climatology is still widely used in impact analysis studies (e.g., Anandhi et al. 2011; Fatichi et al. 2015). Finally, the limited duration comparison of 9 years of baseline simulation against 9 years of scenario analysis, due to TRMM V7 and flow data availability, is acknowledged.

The paucity of data and the uncertainty in understanding the snow/ice resources and their spatio-temporal dynamics in poorly instrumented mountainous regions remains a significant challenge for conceptual hydrological modelling in snow and glacier-dominated catchments. This often necessitates model assumptions to be made regarding initial snow/ice boundary conditions. This study evaluated how contrasting baseline permanent snowpack assumptions in a semi-lumped conceptual model affected simulated hydrological response to changes in future climate. Two conceptual model settings in which the glacier and permanent snowpack are either in steady state (Model 1, in which annual snowmelt is constrained by the seasonal snow accumulation so that there is no net loss over the baseline period) or are not in equilibrium (Model 2, so that snowmelt is only energy-limited and there can be a net loss or gain over the baseline period) were calibrated and validated in the mountainous Beas river catchment in northern India. The findings from the study are as follows:

1. Despite the contrasting assumptions about historical snowpack response to baseline climate variability, both models had generally very good and similar baseline hydrological performance. The compensation effect of parameter fitting was demonstrated to obscure the effect of uncertain boundary conditions on process and model behaviour.

2. The effect of differing parameterization on baseline hydrological process behaviour is hidden, so that models can have unrecognized limitations in simulating streamflow under scenarios that are outside of the baseline conditions. We have demonstrated that the potential inadequacies of such resultant parameters values only become evident when the models are subjected to changing climatic conditions.

3. Across a scenario-neutral climate change space of changing annual temperature of 1–5 °C and changing annual precipitation of −10 to +20%, the model parameterization under the two baseline snowpack behaviour assumptions introduced medium-term discharge uncertainty across the scenario-neutral climate space of between 36% (ΔQ10 of −3 to −39%) for Q10 and 163% (ΔQ90 of −141% to +22%) for Q90. These uncertainty ranges obtained from the difference between Model 1 and Model 2 during 9-year scenario analysis are significantly greater than the modelled baseline uncertainty of around 30%.

4. Under longer-term changes associated with the loss of the permanent snow/glacier, the uncertainty across the scenario-neutral climate space amounts was between 28% (ΔQ10 of −15% to +13%) for Q10 and 37% (ΔQ90 of −31 to +6%) for Q90.

5. Although the uncertainty arising from the baseline parameterization reduces with the transition from a melt-dominated river to a rain-dominated river, associated with the loss of area and volume of permanent snow/icepack resources, the significant temporal changes in flow dynamics and magnitudes (Nepal 2016) will have significant impacts on simulated future dry season water resources for irrigation, hydropower and livelihoods.

6. In this particular case study, although the NSE values are similar for both assumptions, we acknowledge the differences in the PBIAS values (especially in the validation phase). Nevertheless, our study suggests that, where there is considerable uncertainty in historical snowpack reserves and dynamics, an ensemble of hydrological model-builds calibrated to the different assumptions should be used to inform the understanding of the resultant effect of parameter biases on climate change impact studies.

We thank Dr Sanjay Jain (National Institute of Hydrology, Roorkee) and the Bhakra Beas Management Board for supplying discharge data and Water Resource Associates for the use of the HySIM software. This work was supported by the Natural Environment Research Council (NERC) – United Kingdom (grant number NE/I022329/1 and NE/N015541/1), as part of a joint NERC – Indian Ministry of Earth Sciences research programme. All data supporting this study are openly available at https://doi.org/10.17862/cranfield.rd.6969836. No conflict of interest is declared.