Discharge sensitivity to snowmelt parameterization: a case study for Upper Beas basin in Himachal Pradesh, India

The large Himalayan rivers, Indus, Ganges, Brahmaputra, Yangtze and Yellow river, are the main water supply for some of the most densely populated areas of the world. The glacier- and snowmelt contribution to discharge in these regions is important and influences the availability of water for domestic, agricultural and industrial use, as well as for hydropower production (Winiger et al. 2005). The contribution of glacier- and snowmelt to annual discharge varies widely throughout the Himalayas. For the smaller Upper Indus basins, studies report that snowmelt contributes up to 60% of annual discharge (Jeelani et al. 2012; Li et al. 2014), whereas other studies indicate glacier contribution alone in the order of 30% (e.g. Bookhagen & Burbank 2010). Due to the importance of glacier- and snowmelt in discharge generation, the Intergovernmental Panel on Climate Change fourth assessment report (IPCC 2007) expressed great concern for the water availability in the Himalayan region, based on the climate change and glacier retreat prognoses. The report erroneously referred that most of the Himalayan glaciers could disappear by 2035 under current climatic conditions (Cruz et al. 2007). The prediction was later strongly modified suggesting a much longer time horizon for the Himalayan glaciers (IPCC 2014). However, it created a focus on the region in which there is limited knowledge and available data. Overall, the glaciers of Himalaya are retreating though some glaciers show a different behavior of mass gain, referred to as the Karakoram anomaly (Hewitt 2005; Bhambri et al. 2011; Bolch et al. 2012; Kääb et al. 2012).

Hydrological models are simplifications of the physical processes that take place in the natural environment. Different methods for snowmelt calculations exist and the choice of model complexity is determined by the objective of the study, the available data and the nature of the environment to which the model is applied (Haan et al. 1982). The temperature-index method uses temperature as the sole parameter controlling snow melt, and it is widely used due to its simplicity and good performance (WMO 1986; Hock 2003). The latter is because a close correlation between temperature and shortwave radiation exists in most regions. Sicart et al. (2008) showed, however, that for a high altitude and low latitude glacier in Bolivia, the energy budget was controlled by the net shortwave radiation and had low correlation with temperature. This indicates that in some regions, glacier- and snowmelt cannot accurately be described by temperature alone, and improved glacier- and snowmelt modelling can be achieved by taking into account the shortwave radiation. A combined temperature-radiation-index method is suggested by Pellicciotti et al. (2012) for modelling snowmelt in the Hunza River basin in Karakoram Himalaya. The challenge using the temperature-radiation-index or a complete energy balance method is the extended need for input data. Both methods need shortwave radiation; in addition, the energy balance method needs longwave radiation, wind, and humidity data. These climate observations are not readily available, specifically not for the remote areas of the Himalayas. An alternative is the use of global data sets or remotely sensed data from which certain input data may be derived.

In this study, a gridded hydrological model is set up for a high mountainous basin in the Indian Himalayas, the Upper Beas basin. The aim of this study is two-fold: (i) to evaluate the sensitivity of simulated discharge to increased complexity in the snowmelt routine; and (ii) to assess the water balance components contribution to seasonal and annual water balance. The model applies local as well as global data sets, and three different snowmelt routines are evaluated.

The glacier coverage varies among the sub-basins, from 1.5% in Tirthan to 28.3% in Parvati. Table 1 provides some key characteristics of the five sub-basins. Input data consist of daily observations of precipitation, temperature, and humidity from a network of six climate stations (Figure 1). The stations are located at elevations ranging from 904 masl at Pandoh Dam, to 1,971 masl at Manali. Pandoh Dam is the only station measuring potential evaporation, whereas Banjar and Sainj only measure precipitation. In Table 2 the annual mean of the observed variables (station values) in the Upper Beas basin is presented for the calibration and validation period separately. Gridded daily wind and shortwave incoming radiation were obtained from the WATCH Forcing Data (WFD) set, a biased corrected version of the ERA-40 reanalysis data set with a spatial resolution of 0.5 ° (Weedon et al. 2011). In WFD no bias correction is applied to wind, whereas the shortwave incoming radiation is corrected with respect to the average cloud cover and effects of changing atmospheric aerosol loading given in the CRU TS2.1 data set (New et al. 1999, 2000; Mitchell & Jones 2005).

The topographic data for the model was built from a digital elevation model based on the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) using the Global Digital Elevation Model version 2 (GDEM2). ASTER GDEM2 has an overall accuracy of about 17 m at a 95% confidence level (METI/NASA 2011). The original GDEM2 resolution of 30 m was resampled to a 1 km resolution. Glacier covered area (GCA) is derived from the glacier inventory of the Upper Beas basin (Berthier et al. 2007; Frey et al. 2012) downloaded from the Global Land and Ice Measurements from Space database (www.glims.org). Land use and land cover information for the basin have been prepared based on the Indian Remote Sensing (IRS) satellite IRS-P6 data (NRSC/ISRO 2008).

ENKI (https://bitbucket.org/enkiopensource/enki) is a modelling framework that allows the running of a customized hydrological precipitation-runoff model, built from a library of subroutines, for a region that is specified by geographical information system (GIS) data (grids and point networks) as described in Table 3. The model was run on daily time steps for a grid with a spatial resolution of 1 × 1 km². The model parameterizations were applied globally to each grid cell. Precipitation was interpolated using an inverse distance weighting routine with a fixed elevation gradient (PrecGrad, Table 4), and the temperature was interpolated using 3D kriging and a temperature lapse rate (PriTgrad, PriSDtgrad,Table 4). Humidity and shortwave radiation were interpolated using an inverse distance weighting routine. The actual evapotranspiration was calculated by adjusting mean monthly values of potential evaporation (measured at Pandoh Dam) with respect to anomalies in temperature, wind, moist air during rainfall, vegetation height, and presence of snow-cover. The actual evapotranspiration increases with temperature and wind, and is reduced during rainfall events, for grid-cells with low vegetation, and for snow covered area. The snow accumulation and snowmelt routine was based on the snow depletion curve (SDC) model (Kolberg & Gottschalk 2006; Kolberg et al. 2006). The SDC is represented by a gamma distribution, specified by a coefficient of variation and a possible initial bare ground fraction. Daily snow covered area is calculated based on daily updated SDC throughout the melting season.

Snowmelt was calculated by three separate routines: (i) the complete energy balance method (GamSnow routine), (ii) the temperature index method (GamDDF routine) and (iii) the combined radiation and temperature index method (GamSRF routine). All parameter values are given in Table 4.

The parameters FastDecayRate (days) and SlowDecayRate (days) can be interpreted as the time it takes for the albedo to decrease 95% of the value defined by the difference between the maximum and minimum albedo. After snowfall larger than a predefined depth (ResetSnowDepth), the albedo is reset to maximum albedo (MaxAlbedo). Glacier albedo is set to a constant.

Incoming and outgoing longwave radiations, L_in and L_out, are calculated following the Stephan–Boltzmann law. L_in is based on air temperature T_a, whereas snow surface temperature T_ss, calculated as T_ss = 1.16·T_a–2.09, with 0 °C as maximum value, is used to estimate the outgoing longwave radiation. The latent and sensible heat fluxes, H_L and H_SE, are calculated using a bulk-transfer approach that depends on wind speed, temperature and air humidity. The two parameters defining the wind profile, intercept (windconstant) and slope (windscale), are determined by calibration. H_L and H_SE are assumed always positive, and set to zero for temperatures below 0 °C. The subsurface flux E_G is calculated assuming a linear temperature profile in the snow surface layer, where temperature at the bottom of the layer is set to 0 °C. Snowmelt in mm is derived by multiplying the available energy ΔE with the latent heat of fusion for water.

Glacier melt occurs when glaciers are exposed due to a reduction in snow cover, i.e. when the snow covered area is less than the GCA within a grid-cell. The glacier melt routine has unique parameters for α, DDF and SRF, all determined by calibration.

The soil moisture and response routines adopted are schematically similar to the HBV model (Bergstrom 1976). An overview of all parameters included their calibration range and optimized values (for those that were calibrated), are given in Table 4. Simulated discharge is calculated by accumulating runoff from all grid cells within a basin for each time step, without any delay. In the following, the different versions of the model are named by the snowmelt routine used, i.e. GamSnow, GamSRF and GamDDF models.

To separately assess the effect of the snowmelt routine on the simulated discharge, a two-step calibration procedure was performed. First, the parameters listed in Table 4 and those related to the GamSnow routine (Table 5) were calibrated. Subsequently, the calibrated parameter values related to the soil and response routines were held constant (optimized values specified in Table 4), whereas parameters related to snowmelt calculations (GamDDF and GamSRF) models were calibrated (Table 5). Four observed discharge series were included in the calibration of the Upper Beas basin: Thalout, representing the whole Upper Beas basin, and the sub-basins Tirthan, Sainj and Bhuntar.

The three models were calibrated using observed daily discharge for the period 01.08.1991 to 31.07.1996 and validated from 01.08.1996 to 31.07.2001. A warm-up period was set from 01.01.1989 to 31.07.1991.

The model performances in terms of NSE, Pbias and Pearson's r are shown in Table 6 for the Upper Beas basin and the three sub-basins Tirthan, Sainj and Bhuntar. Discharge data for Bhuntar was used for the calibration, whereas its sub-basins Manali and Parvati were used in the more detailed analyses, resulting in a total of four sub-basins (Tirthan, Sainj, Parvati and Manali). Acceptable NSE values (∼0.7) are found for the Upper Beas basin for the calibration period, whereas lower values are seen for the validation period (∼0.55). Notably low NSE was obtained for Bhuntar and to some degree also for Sainj for the validation period, whereas Tirthan showed consistent values for both the calibration and validation period, and all sub-basins (∼0.60–0.64). The volume deviation, Pbias, varies from negative values for Sainj (–15.7%) and the GamDDF model, to positive values (8.2%) for Upper Beas and the GamSRF model during the calibration period. Only positive values are seen for the validation period, ranging from 4.8 to 55%, where the highest deviation is found for Bhuntar for the GamSRF model. Pearson's r is in the range 0.8–0.9 for all models and sub-basins, both for the calibration and the validation period.

Overall, the performance criteria shows that all models are able to satisfactorily simulate daily discharge for the calibration period. The NSE values for the calibration period are similar to results obtained for the same data set in other studies using different hydrological models (Li et al. 2013, 2014). This indicates that the calibration result is limited by the quality of the hydro-metrological data and to a lesser degree by the choice of model. A reduced model performance for the validation period is reflected in the values for NSE and Pbias, and less in the correlation coefficient (daily values). The three models evaluated, representing different snowmelt routines, mainly deviate in their volume prediction, whereas there are smaller differences in the timing of flow events, reflected in relatively high Pbias and high correlation coefficient, respectively.

At the seasonal scale, the largest volume deviation in model prediction was found for July and August for all sub-basins, reflecting the high melt contribution during these months. For the highest elevated basins (Sainj and Parvati), the largest discharge in the monsoon period in Sainj was obtained by GamSnow, whereas GamDDF gave the highest and GamSnow the lowest simulated discharge in Parvati. This suggests a lack of consistency among models in simulating snowmelt discharge for the highest elevated basins during the monsoon months.

None of the models was able to predict the negative trend in annual discharge observed for Sainj, Parvati and Manali. The negative Pbias for annual values for these sub-basins could implicate error in the input data, other explanations are glacier dynamics or interpolation of precipitation. However, all models acceptably simulated Tirthan, which experienced no trend in discharge.

Highest precipitation was estimated for the northern sub-basins Manali and Parvati, whereas lowest precipitation was obtained in Tirthan, at the lowest altitude. The precipitation dependency with altitude is described by a power-law relationship in the interpolation routine. Negi (2002) suggests that a shadow effect occurs for the high mountains east in the Beas basin, where typically the valleys lying northeast of the basin, in Lahul and Spiti, are mountain desert areas with annual precipitation around 400 mm. It is likely that a reduction in precipitation would occur above some elevation. Overestimation of precipitation at the highest elevations will, due to all year round low temperature, be retained as snow storage, not affecting simulated discharge when temperature is stable. However, the overestimation of precipitation and the development of too large snow storage will potentially become a source of error in a warmer, future climate. In the high-elevated Parvati sub-basin, precipitation has been declining since 1998, as has observed discharge. A similar reduction is not seen in simulated discharge, which is assumed to be sustained by questionably large snow storage developed during the first years of simulation. A smaller fraction of precipitation is stored as snow in the model for Manali sub-basin. Due to its lower mean altitude the discharge is more directly connected to precipitation. Manali does however have the largest volume deviation between simulated and observed discharge. Precipitation and observed discharge in Manali would indicate a volume increase in the reservoirs of perennial snow and glacier growth. This implies that the volume deviation between observed and simulated discharge in Parvati and Manali could arise from erroneous precipitation interpolations. For Tirthan and Sainj, the simulated discharge is larger than can be sustained by precipitation alone. The basins have no snow storage from previous years, suggesting that excess water originates from glacier melt. This is supported by the negative storage change calculated for the two southern sub-basins for most years. Sainj show the same pattern as Parvati, as the simulated and observed discharge has different behavior in the later part of the period. Observed discharge follows the reduction in precipitation, whereas this is not the case for simulated discharge which is kept high for both sub-basins, probably by simulating melt at higher elevations.

Evaporation is an overall minor component in the annual water balance budget for all sub-basins. The highest evaporation is found in the valleys at lower altitude regions where annual mean temperatures are higher. Hence, Tirthan has the highest average evaporation (215 mm/year), followed by Manali (133 mm/year), Sainj (102 mm/year) and Parvati (44 mm/year).

The observed data show a decrease in discharge (Figure 4) whereas there is an increase in both precipitation and temperature (Table 2). This indicates a non-stationarity in precipitation and discharge. Li et al. (2015) highlight non-stationarity in the precipitation to discharge ratio analyzing the annual water balance for the 1997–2005 period for the Bhuntar sub-basin. Their results show that the linear relationship between precipitation and discharge during the early period (1997–2001) is different from the relationship during the later period (2002–2005). This could imply deficiencies in the data quality, but also it raises the question as to how the two dominant weather systems, the summer Monsoon and Western Disturbance, will affect the annual precipitation/discharge ratio. This could be an additional explanation to the deviation we see in the data for the different years.

The observed and measured reduction in GCA in the Beas basin during 1991–2001 (Kulkarni et al. 2005) is not accounted for in the model, which assumes a static GCA with an infinite volume of ice. The Parbati glacier is located within the Parvati sub-basin, and is one of the glaciers with the largest areal reduction in the Western Himalaya. The glacier snout has retreated to about 6,000 m in the period 1962–1990, with a somewhat slower retreat from 1992 to 2001 (Kulkarni et al. 2005, 2011). A reduction in the GCA and increased elevation of the snowline support an increase in discharge due to additional melt. If, however, the rate of reduction is slowing down or the retraction has passed a critical point, further retraction will reduce the melt volume (Jansson et al. 2003). Reduced glacier melt may explain the trend in observed discharge for the study period seen for the northern sub-basins Manali and Parvati, which cannot be explained by a reduction in precipitation and, further, has not been captured by the model.

The pre-monsoon simulated discharge in March and April was overestimated by all models, indicating a too early start of model snowmelt compared to observed discharge (Figure 7). Precipitation in the pre-monsoon period is typically low, so discharge is mainly fed by snowmelt and hence, is dependent on simulating snowmelt satisfactory. Thus, an overestimation of pre-monsoon discharge suggests the presence of snow at low elevations or a too early start of the melt season in the model.

Important factors contributing to the uncertainty in our modelling study include: (i) spatial interpolation based on a sparse network of station date; (ii) downscaling coarse resolution WDF data in challenging terrain; and (iii) model structure and parameterization.

Generally, the quality of input data will have a large influence on model performance. Data used in this study is largely drawn from manual measurements and hence, is subject to human errors. Precipitation measurements might under-estimate the real precipitation, and losses due to extreme precipitation or snow events under windy conditions are particularly difficult to account for (Wolff et al. 2015). The Upper Beas basin is characterized by large elevation differences and lack of observations in the highest elevated areas. The input temperature and precipitation is based on extrapolation of their respective laps-rates. This implies that the interpolation of climate variables is highly uncertain in the region and depends on good observation networks, both in terms of spatial representation and quality of observations. However, the lack of observation sites at high altitude is considered the most critical factor in determining uncertainty in interpolated input data, precipitation in particular, as suggested by overestimation of precipitation for the highest elevations in this basin. In addition, the marked seasonality in hydroclimatology suggests that the model simulations could gain from seasonal specific parameters related to the interpolation of temperature and precipitation. Thayyen et al. (2005) found a non-linear temperature lapse-rate for a valley in the Din Gad basin in Central Himalaya. They demonstrated that the valley lapse-rate was smaller compared to the alpine lapse-rates at higher elevations during the monsoon months, and that the alpine zone showed, in addition, a higher annual variation than the lower valley zone. This could give additional information to the large inter-annual differences found in the precipitation to discharge ratio in this study and by Li et al. (2015). However, these findings are not necessarily representative for our study region as large regional variations can be expected since the wet and dry adiabatic lapse-rates depend on atmospheric moisture content and hence, on precipitation and temperature in the area (Aguado & Burt 2010).

The high snow accumulation simulated for Paravati is questionable when confronted with the areal retreat observed for glaciers in this sub-basin (Kulkarni et al. 2005, 2011). Rather, it may be a result of an excessive precipitation gradient in the model overestimating precipitation at high altitudes. Additional information from remotely sensed data of snow covered areas might help determining more realistic values. Alternative data sources like Tropical Rainfall Measuring Mission (Huffman et al. 2007) and WFD (Weedon et al. 2011) were considered. Previous studies (Li et al. 2012, 2013) have used these data to force hydrological models, the general conclusion is that the model performance is still not comparably good to that of using observed data. Different downscaling methods for regional climate models (RCM) for climate projections studies was tested on the Beas basin by Li et al. (2016). The bias correction significantly reduced the difference from RCM to in situ observations, and similar methods should be investigated for the downscaling of global data sets. The use of wind and radiation from WFD for energy balance modelling in this region with a complex topography is challenging, and local topography was not accounted for when interpolating these inputs. Validation of downscaled wind and radiation by surface observations is not possible, since local data is missing.

The Beas River with tributaries is prone to large annual variations and extensive flooding, and both low and high flows can be inaccurate. There is an extensive sediment transport in the rivers (Jain et al. 2007), implying unstable river-profiles, and by this it is assumed difficult to establish stable rating curves.

In GamSnow a simplified algorithm for the calculation of snow surface temperature assumes turbulent fluxes always to be positive, an empirical relationship fitted for Norwegian conditions. Energy balance modelling at the Chhota Shigri glacier just north of the Beas basin (Pithan 2011), found a strong negative latent heat contribution in the pre-monsoon period as also supported by findings for a low latitude glacier in Bolivia (Sicart et al. 2008). Hence, these results are in conflict with the snow surface temperature calculations adopted in our model setup. Wind is highly variable and dependent on local roughness in addition to the overall topography. This will influence the calculations of turbulent fluxes, and hence the surface temperature. Shortwave radiation will also be influenced by the topographic parameters, slope and aspect, which modifies the exchange surface for radiation, in addition to steep terrain leading to shading (DeWalle & Rango 2008). These effects were ignored in the interpolation of climate input data and may result in reduced model performance for GamSnow and GamSRF model. The importance for this model setup is however questionable, since all three models experience problems with the water balance, independently of the method used for snowmelt calculations.

Glacier dynamics are not implemented in the models, and the static GCAs are hence an infinite water source. For the southern, low altitude, sub-basins Sainj and Tirthan glacier melt contributed to simulated discharge in most years, indicated by a consistent negative change in storage. The glacier melt contribution was, especially for Tirthan, high relative to the comparable small GCA. This underlines the problem with a relatively coarse resolution in steep terrain. In case only a small part of the grid cell is glaciated, its temperature, which is derived based on the mean altitude for the grid, may overestimate the melting. Methods to reduce erroneous glacier melt, for example introduce model limitations for the maximum glacier melt, from each grid and time step based on the glacier-covered area should be investigated.

Water from glacier- and snowmelt is important for both early spring flow and peak flow during the summer monsoon in the high altitude river basins of the Himalaya. Three different parameterizations of the snowmelt routine with varying complexity have been evaluated for the Upper Beas basin including four sub-basins. The models were run on a daily time step and a spatial resolution of 1 × 1 km² using a combination of local and global data sets. The following conclusions are drawn from the study:

• The three snowmelt parameterization routines performed equally well with respect to simulated discharge for the calibrated period with NSE as objective function. Variability in input data, input data processing and interpolation, likely have a larger influence on model performance than the choice of snowmelt routine.

• None of the models was able to predict the negative trend in discharge observed for the sub-basins Sainj, Parvati and Manali for the study period. This is likely caused by (i) static glacier parameterization, (ii) precipitation interpolation routine overestimating high altitude precipitation and thus, snow storage or a combination. It is noted that Tirthan, experiencing no trend in observed discharge, was better predicted.

• The water balance components showed minor seasonal differences among the models. However, discharge is overestimated in the pre-monsoon period for all models and sub-basins, suggesting a too early melt onset in the model.

Increased model complexity did not enhance the model performance for the Upper Beas basin, implicating that downscaling the WFD gives too rough estimates for wind and radiation to be representative for challenging regions with complex topography. Validation of downscaled radiation and wind by surface observations is difficult since local data is missing. Further studies evaluating optimum downscaling strategies for this region should be encouraged. Implementation of a dynamic glacier module and different precipitation interpolation routines should be further investigated for improved simulation of glacier- and snowmelt. Seasonality in model parameters and a better differentiation in sub-basins parameter values should be considered in future calibration schemes.

We would like to thank Dr Sharad K. Jainj for valuable insight and local knowledge. This study was jointly supported by the Research Council of Norway projects JOINTINDNOR (203867), and INDNOR (222195): Hydrologic sensitivity to Cryosphere-Aerosol interaction in Mountain Processes (HyCAMP).