Integrating a glacier retreat model into a hydrological model – Case studies of three glacierised catchments in Norway and Himalayan region

Glaciers are crucial in many countries where meltwater from glaciers is an important source of water for drinking water supply, irrigation, hydropower generation and the ecological system. Glaciers are also important indicators of climate change. They have been significantly altered due to the global warming and have subsequently affected the regional hydrological regime. However, few models are able to parameterise the dynamics of the glacier system and consequent runoff processes in glacier fed basins with desirable performance measures. To narrow this gap, we have developed an integrated approach by coupling a hydrological model (HBV) and a glacier retreat model (Dh-parameterisation) and tested this approach in three basins with different glacier coverage and subject to different climate and hydrologic regimes. Results show that the coupled model is able to give satisfactory estimations of runoff and glacier mass balance in the Nigardsbreen basin where the measured data are available to verify the results. In addition, the model can provide maps of snowpack distribution and estimate runoff components from glaciers. 2015 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http:// creativecommons.org/licenses/by/4.0/).


Introduction
Glaciers are essential in the water system. They store about 75% of the fresh water on the earth. Approximately 99.5% of the ice volume store in ice sheets and the remaining 0.5% in mountain glaciers (Khadka et al., 2014). Although the mountain glaciers are relatively small, they are critically important to the humanity (Beniston, 2003;Arora et al., 2014). Firstly, mountain glaciers are relevant to humankind because of their proximity to populated areas (Singh et al., 2006). During 14th to 19th century, a period known as the Little Ice Age, advancing glaciers periodically extended down into valleys and destroyed communities, crops, and livelihoods (Carey, 2007). Glacier recession also induces a series of natural hazards (Carey, 2005). Secondly, mountain glaciers indeed substantially contribute to the development of society and economy (Beniston, 2003). For example, glaciers only cover approximately 1% area of mainland Norway. However, 15% of its electricity is generated by runoff from these glacierised basins (Andreassen et al., 2005). Thirdly, glacier meltwater is an important water resource (Burlando et al., 2002) and a supplement to streamflow under drought conditions (Marshall, 2014).
Furthermore, glaciers are considered as one of the most sensitive indicators of climate change (Burlando et al., 2002). That is because runoff from a glacierised basin is more dependent on the available energy than from a non-glacierised basin. Therefore, a modification of the prevalent climate, particularly of air temperature, can considerably affect the hydrologic regime (Burlando et al., 2002;Radić and Hock, 2014).
Glacier responses to climate change are of concern for both scientific research and public communities (Hagg et al., 2007;Barnett et al., 2005). In many regions of the earth, glaciers are retreating and seasonal snowfalls are diminishing as shown by observations and modelling (Bolch et al., 2012;Barnett et al., 2005). Zemp et al. (2008) reported that the global average annual mass was À0.58 m water equivalent (w.e.) year À1 for the decade [1996][1997][1998][1999][2000][2001][2002][2003][2004][2005]. Assuming a minor change in precipitation, ice-and snow-covered areas are predicted to decrease at accelerating rates due to increased melting of snowpack and ice as well as due to reduced accumulation of snow and decreased surface albedo . The expected climate change and its serious implications demand interdisciplinary research on different topics including glaciers.
The response of glaciers to global warming is generally stated as: due to the release of water from glacial storage, runoff initially increases especially during late spring and summer, and after complete removal of the glacier ice, runoff stabilizes and drops below the previous level (Huss et al., 2008). Theoretically, the total amount of runoff cannot be changed by the disappearance of glaciers compared to stable glacier conditions. However, changes in seasonal distribution of runoff have significant impacts on water availability, ecosystems and human well-being (Jain et al., 2011).
The inhomogeneity among mountain glaciers, however, makes this generalisation not applicable for all glaciers (Jain et al., 2011). This inhomogeneity is significant among different climate regimes, even among the glaciers in a close neighbourhood (Horton et al., 2006). The rate of glacier retreat depends on both large-scale and time-invariant factors, and small-scale and time-dependent factors (Burlando et al., 2002;Khadka et al., 2014). Geographical, topographical and climatic conditions act at large scales and relatively stable, whereas the size, location of glaciers and study period of research are at a small spatial scale and time dependent (Bolch, 2007;Marshall, 2014).
The processes of melt water from glaciers draining to the catchment outlet are very complex. Glaciologists and hydrologists have made many efforts to describe and model them by employing glaciological methods and hydrological methods or by integrated knowledge. The methods used can be categorised into three types: glaciological, hydrological and interdisciplinary approaches.
From a glaciological point of view, it is essential to simulate response of glaciers to climate change and contribution to streamflow by glacier wastage. This approach mainly uses models to relate meteorological measurements to accumulation and ablation rates on the glacier surface. These mass balance models vary in complexity, including conceptual models and or more complicated models based on the energy balance (Gottlieb, 1980;Arnold et al., 1996;Klok and Oerlemans, 2002;Hock, 2005) and models for ice dynamics (Greuell, 1992;Hubbard et al., 1998;Jouvet et al., 2009). They have been validated by in-situ glacier measurements, such as surface mass balance or ice thickness change. Then glacier mass loss can be compared with discharge downstream of the glacier (Kotlarski et al., 2010). However, glacier melt water is only a raw volume input into the hydrological system, similar to rainfall. The discharge is further modified by, e.g. evaporation and groundwater storage. Therefore, the direct comparison leads to an overestimation of glacier contribution to discharge (Kaser et al., 2010), unless contributions from non-glacierised areas are also considered.
Many diverse approaches have been used in developing hydrological models for a better representation of snow processes and ice melting (Bergström, 1976;Shrestha et al., 2012;Wrede et al., 2013). The models mainly use the meteorological data, satellite data of snow cover and climate scenarios to study effects of climate change on glaciers (Khadka et al., 2014;Jasper et al., 2004). Among the models, the Snowmelt Runoff Model (SRM) (Khadka et al., 2014;Immerzeel et al., 2009) and the HBV model (Mayr et al., 2013;Hagg et al., 2006;Akhtar et al., 2008) are widely used. However, these studies did not take into account the change in the glacier area to the imposed climatic changes. Although these studies provide valuable insights into the possible range of the future options, they suffer from a large uncertainty about the plausibility of the future evolution of glaciers (Immerzeel et al., 2012). This limits the time scale over which the results could apply.
Several studies using the HBV model considered the retreat of glaciers (Hagg et al., 2006;Akhtar et al., 2008) but, since the change of glacier extent were arbitrarily assigned, they can be considered as rather simple and non-dynamic experiments not allowing to account for the transient evolution.
The coupled models with both descriptions of glacier and hydrology are quite new. Several of such models have been reported to link glacier dynamics and hydrological processes. Horton et al. (2006) studied the climate change impacts on the runoff regimes in the Swiss Alps by a conceptual reservoir-based hydrological model. The glacier extent was updated by a conceptual model based on the accumulation area ratio (AAR) method, which assumes that the accumulation area of a glacier occupies a fixed proportion of the total glacier area. However, the AAR method is not able to reproduce the transient response of glaciers to a changing climate as it assumes the glacier to be permanently in steady-state (Huss et al., 2008). Huss et al. (2008) ran a distributed model to simulate the daily surface mass balance on three highly glacierised catchments in Switzerland. The model included the water balance calculation and a parameterisation of glacier retreat. The model was validated with monthly runoff and decadal ice volume change. Stahl et al. (2008) proposed and applied a methodology for estimating changes in streamflow associated with the coupled effects of climatic change and associated glacier response, with a specific focus on transient responses. The authors combined the HBV model with a model of the glacier area evolution based on volume-area scaling. However, as shown by Lüthi (2009) and Radić and Hock (2014), volume-area scaling cannot describe a critical time lag between the area and the volume response to the prevailing climate. Immerzeel et al. (2012) suggested a combined approach consisting of a simple precipitation-runoff model and an ice dynamics model including basal sliding for a glacierised Himalayan catchment in Nepal. The model was first calibrated using the recent location of the glacier terminus and then by using discharge observations. Naz et al. (2014) published a physically-based, spatially distributed hydrological model coupled to a shallow ice dynamics model. A common drawback of the currently available physically-based models is a high demand of input data and computational power, which hampers their applicability to large-scale catchments.
This paper aims at presenting a coupled model for hydrology and glacier retreat suitable for large catchments with low data demand. The combined model explicitly simulates glacier evolution and major hydrological processes at a high spatial resolution.

Methodology
The general framework is to couple a distributed HBV model with a mass-conserving glacier retreat model. The HBV model calculates the accumulation and ablation of snow and glacier ice for every grid at a daily time step and the glacier retreat model describes glacier surface elevation changes and updates glacier area in response to the total amount of glacier mass change calculated by the HBV model.

Hydrological model
The HBV model concept was initially developed for runoff simulation in the Scandinavian countries in the 1970s (Bergström, 1976). So far, the HBV model has been applied in more than 80 countries and is used as a standard tool for flood forecasting and for simulating inflow to hydropower reservoirs in Norway as well as many other European areas. The scientific basis and applicability of the model have been verified at many sites in Europe (Uhlenbrook et al., 1999;Götzinger and Bárdossy, 2005; Hailegeorgis and Alfredsen, in press) and Asia (Akhtar et al., 2008;Li et al., 2014b). The algorithms have been described in details (e.g., Bergström, 1976;Lindström et al., 1997;Beldring et al., 2003); therefore, we present here only the equations of special significance to our research.
The model version used in this study was first published by Beldring et al. (2003). It is a distributed model considering the spatial distribution of meteorological data at a daily time step. The main inputs are precipitation, temperature and surface elevation. The station precipitation is firstly corrected for under catch (Eq. (1)). Then the corrected precipitation and air temperature series are interpolated for every grid using the inverse distance weighting (IDW) method considering elevation effects (Eq. (2)). Finally water balance is computed by the HBV algorithms for every grid.
where P c is corrected precipitation (mm day À1 ) and P o is precipitation measured at the weather station (mm day À1 ). K r and K s are free parameters respectively for rainfall and snowfall correction.
where P g and T g are precipitation (mm day À1 ) and temperature ( C) at each grid. W i is the weight of station i calculated by the IDW method. H g and H i are the elevations in meters above mean sea level (m amsl) of grid and station i. n is the number of stations used for interpolation. c t ( C per 100 m) and c p are free parameters that describe elevation effects on temperature and precipitation, respectively. Snow accumulates when the air temperature is below the threshold temperature of snow accumulation (T acc ; C) and there is precipitation. Snow melting is calculated by a degree-day method based on temperature according to Eq. (3) (Lindström et al., 1997). Glacier melting starts when there is no snow coverage and the air temperature is higher than the melting threshold temperature (T t ). The melting rate of ice (Melt ice ) is based on the same method, but another melting factor (MF ice ) as shown in Eq. (4). At the end of melting season (i.e. 31st August) existing snow on the glacier is transformed into ice, and then the glacier gains mass.
where T is air temperature ( C) and T t is the threshold temperature of snow melting; Melt snow and Melt ice are melting rate of snow and glacier (m day À1 ) obtained with the degree-day factors MF snow and MF ice , respectively (both in m day À1°CÀ1 ).

Glacier retreat model
The Dh-parameterisation for modelling the changes in surface elevation and extent for retreating glaciers owes its origin to varying thinning rates over a glacier (Huss et al., 2010). For a certain mass change, surface elevation changes are smallest in the accumulation area and the largest near the terminus of mountain glaciers. This has been confirmed by measurements (Arendt et al., 2002;Bauder et al., 2007) and numerical modelling (Jóhannesson et al., 1989). The Dh-parameterisation describes the surface elevation changes as a function of elevation and the ice volume change, as Eq. (5).
where h is the glacier surface elevation (m amsl); h max and h min are the maximum and minimum elevations (both in m amsl); h r is the normalised elevation; Dh is the normalised surface elevation change. c; a; b and c are parameters and their values can be derived from glacier surface maps of different years, or calibration (Huss et al., 2010). Integration of the dimensionless Dh function (Eq. (5)) over the entire glacier, taking into account the area distribution and the ice density q ice (900 kg m À3 ) must equal the total glacier mass change B a (kg) calculated by the HBV model in a given time interval: where Dh i is the normalised elevation change of grid i; A i is its area (m 2 ); n is the number of glacier grids. Further, f s is a factor that scales the magnitude of the dimensionless ice thickness change. It is chosen for each time interval such that Eq. (6) is satisfied. Then the glacier surface is changed as: where h 0 is the elevation (m amsl) in the previous time step and h 1 is the updated elevation (m amsl). In the last step, the change in glacier extent is determined by the updated glacier surface elevation. The glacier disappears for grids where the thickness is not greater than zero. The required inputs are the initial ice thickness and surface elevation. The algorithm runs for individual glaciers, but all glaciers in a basin share the same values of the parameters. In this research, the glacier flow sheds are constructed from the surface elevation according to the D8 (deterministic eight-node) algorithms (O'Callaghan and Mark, 1984) and the grids having the same terminus point belong to one glacier. Glacier area has a slower climate response than runoff (Lüthi, 2009). Therefore, glacier extent is updated at a yearly scale, at the end of every melting season, i.e. 31st August.
The Dh-parameterisation avoids high demand of field data and computation resources such as an ice dynamics model (Jouvet et al., 2009) and can provide similar estimates as a 3D finite element ice flow model as shown by Huss et al. (2010).

Calibration and criteria
The combined model is calibrated by a model independent calibration package (PEST). PEST is available for free download (Doherty, 2005) and has been widely used in environmental and hydrological modelling (Doherty and Johnston, 2003;Immerzeel et al., 2012;Li et al., 2014a). The objective function is to obtain the minimum of a weighted least square sum of the discrepancies between simulated and observed series, i.e. daily discharge (m 3 s À1 ) and if available, glacier annual mass balance (m year À1 ) where the weight of an annual mass balance value is assumed 10,000 times of a discharge value. When glacier data are not available, the sum of the discharge series is additionally used to control the volume error.
The criteria used for evaluating precision of streamflow simulation are the Nash-Sutcliffe efficiency (NSE) (Nash and Sutcliffe, 1970) and relative mean error (RME), while the Pearson correlation coefficient (COR) is used to evaluate mass balance simulation. The criteria are shown by Eq. (8).
where O i is the observed series; S i is the simulated series; n is the length of series; O and S are the mean values of observed series and simulated series.

Study sites and data
Data of three basins are used to test the coupled model: (1) the Nigardsbreen basin in Norway (Fig. 1), (2) the Chamkhar Chhu basin in Bhutan (Fig. 2) and (3) the Beas basin in India (Fig. 2). The area of the basins ranges from 65 to 3201 km 2 and the glacier-covered fractions are between 15% and 73% ( Table 1). Hydrographs of all basins show a similar discharge pattern characterised by a peak occurring between May and August (Fig. 3).

Nigardsbreen basin
The Nigardsbreen basin is a glacierised watershed in the high mountains in western Norway (Fig. 1). It has a small area of 65 km 2 , but with a large range of elevation. The highest point is 1957 m amsl and the lowest in only 285 m amsl. The climate is humid, as it is influenced by the moist currents from the ocean. In addition, the climate is locally modified by the presence of the glacier. The mean annual air temperature is À0.47 C and the mean precipitation reaches 3736 mm year À1 , with a large amount falling in winter as snow (Andreassen et al., 2005). Streamflow is largely determined by melting of snow and ice in the warm period of the year.
Nigardsbreen is one of the largest outlet glaciers from Jostedalsbreen (the largest ice cap in mainland Europe). The glacier is exposed towards the southeast and extends from approximately 315 to 1957 m amsl. Approximately 73% of the basin area is covered by ice (Andreassen et al., 2012). Nigardsbreen is well-known in the scientific community as it is one of the most studied alpine glaciers in Northern Europe (Oerlemans, 1986(Oerlemans, , 1997(Oerlemans, , 2007Engelhardt et al., 2012).
The climate data are spatially interpolated maps of precipitation and temperature produced by the Norwegian Meteorological Institute (met.no) using 24-h mean temperature and accumulated precipitation measured at meteorological stations (NVE, 2015). These datasets are in a grid format at 1 km resolution and have been evaluated and used in many studies covering mainland Norway, such as hydrological modelling (Beldring et al., 2003;Li et al., 2014a), glacier mass balance estimation (Engelhardt et al., 2012(Engelhardt et al., , 2013, permafrost evolution (Gisnaeas et al., 2013) and snow depth estimation (Vormoor and Skaugen, 2013). The areal mean value of the Nigardsbreen basin is assigned to a virtual station located at the middle of the basin.
The daily discharge is obtained from the national hydrological database, which is collected and managed by NVE. The series are transformed from the in-situ measured water depth by the Bayesian Rating Curve Fitting method (Petersen-Øverleir et al., 2009) and its quality is ensured by the NVE hydrological quality control system. The mass changes of the glacier are derived from the direct glaciological method or stakes-and-pits method, which is traditionally used by glaciologists to measure glacier mass balance (Østrem and Brugman, 1991;Hagg et al., 2004). Stakes are used to measure the changes in thickness and pits are selected sites for measuring the density at different depth. The measurements are interpolated over the glacier to obtain the total amount of glacier mass change. According to the experience of NVE, the measurements are suspected to be subject to a slight positive bias for the Norwegian coastal glaciers and the reason is not clear yet (Andreassen et al., 2011). However, no better dataset is available at present. The mass balance data have been extensively used for glacier modelling purposes (Oerlemans, 1986(Oerlemans, , 1997(Oerlemans, , 2007Engelhardt et al., 2012). Maps of bedrock and glacier surface elevation are provided by the Norwegian Water Resources and Energy Directorate (NVE). The profiles were measured by the Radio-Echo Sounding (RES) method during spring and early summer in the years 1981, 1984and 1985(Saeetrang and Wold, 1986. The bedrock map is obtained through interpolation of the measured profiles and the open valley. The glacier surface elevation at a spatial resolution of 25 m are derived from complete or partial aerial photos taken within the  (Joshi, 2007(Joshi, , 2008(Joshi, , 2011 The Beas basin at the Bhuntar gauging station. (c) The Chamkhar Chhu basin at the Kurjey gauging station. The range of DEM in (a) is assigned to give a better presentation rather than the minimum and maximum for the displaying area. Other signs are same as in Fig. 1.

Table 1
A short summary of the basins. Note: A is area of the basin in km 2 . ME is median elevation in m amsl. GF is glacier fraction in %. P is annual precipitation in mm year À1 . T is annual mean temperature in°C.
period 1984-2009. The ice thickness is calculated as the difference between the glacier surface and bedrock for each grid. The thickness map is further aggregated at a resolution of 100 m, which is the spatial resolution of the HBV model for the Nigardsbreen basin. The starting date of the HBV model simulations is 1st September, 1989, which is roughly the middle of the period of the aerial photos (Table 2).

Chamkhar Chhu basin
The Chamkhar Chhu basin is located in central Bhutan (Fig. 2) and is the source of one of the major national rivers. The basin has three branches originating from the glaciers of the Gangkar Punsum region and the glaciers of the Monla Karchung La region. The river flows south-easterly and finally joins the Brahmaputra River in India. The basin area above the Kurjey gauging station is 1353 km 2 . The elevation ranges from 6653 m amsl in the upper northern glacial region to 2643 m amsl in the southern low-land. The northern part above 4000 m amsl is mainly occupied by glaciers (Fig. 2(c)) whereas the southern low-land is covered by forests.
The climate is strongly influenced by elevation and monsoon; therefore it varies from the southeast to the northwest. Based on observations in the period from 1998 to 2008, the mean precipitation is 1786 mm year À1 and the mean annual air temperature is 1.75 C. The monsoon normally starts in June and lasts until early September. It brings significant amounts of rainfall and warm weather. Subsequently, river flow rises due to the rainfall and melting of snow and ice. As the monsoon proceeds or retreats, there are four clear seasons, spring (March to May), summer (June to August), autumn (September to November), and winter (December to the following February).
The climate data are measured by in-situ meteorological stations, as shown in Fig. 2. Data of seven stations are used; however none of them lies inside the basin. The mean of daily maximum and minimum temperatures is taken by the HBV model as input of daily air temperature. The discharge series of the Kurjey station are obtained from the national authorities, which are in charge of collecting hydrological data. The measurements are available for the period 1998-2008 and have been used for evaluating climate change impacts on hydropower development (Beldring, 2011).
Glacier ice thickness distribution for the Chamkhar Chhu basin is a part of the global dataset produced by Huss and Farinotti (2012) using a method based on glacier mass turnover and principles of ice-flow mechanics (Farinotti et al., 2009). Required input data are a digital elevation model and glacier outlines. For each individual glacier ice thickness distribution is determined for about the year 2000 depending on the date of the utilised glacier inventory data (Pfeffer et al., 2014). The starting date for modelling is mainly determined by the observation period of the meteorological data and is set to 1st September, 1993.
The elevation and land use data are obtained from the Department of Hydromet Services, Bhutan. The original data are at a spatial resolution of 25 m and are rescaled at 1 km. The land use data are reclassified into three broad classes, high biomass (including broadleaf forest, coniferous forest and scrub), low biomass (including erosion, pasture, rock) and human affected (agriculture and urban).

Beas basin
The Beas River is an important branch of the Indus River system in northern India (Fig. 2). It originates at the southern side of the Rohtang Pass in the Himalaya. The river is 470 km long and has a drainage area of 12,916 km 2 (Gupta et al., 1982). This area is extremely rich in hydroelectricity resources. In total, there are 11 hydroelectric plants projects and three of them are ongoing or have been finished (SANDRP, 2015). To avoid the effects of flow regulation and to have a large study area, the Bhuntar gauging station is selected. The station lies downstream of the confluence with the eastern branch, the Parbati River. Only the Malana Hydel Scheme, with a capacity of 86 MW, was running during the study period of 1997-2005. The area above the station is 3202 km 2 . The elevation decreases from 6288 m amsl in the northern mountains to 1055 m amsl. The area above 4500 m amsl is occupied by permanent snow and glaciers (Fig. 2(b)).
The climate is a result of a combined effect of elevation and monsoon. The northern part is much colder and drier than the low valleys. Influenced by the monsoon, there are four seasons, winter (January to March), pre-monsoon (April to June), monsoon (July to September) and post-monsoon (October to December) . The monsoon strength is a major indicator of magnitude of precipitation and temperature. The air currents  of the monsoon originate in the Bay of Bengal and they are weaker after striking east Himalaya and a long westward travel than in the Chamkhar Chhu basin . Therefore, the precipitation decreases with height. Based on the observations during the period from 1997 to 2005, the mean precipitation is 1116 mm year À1 and the mean annual air temperature is À1.04 C. Three meteorological stations measure precipitation, and daily maximum and minimum temperature. Two of the stations are located in the selected area, respectively at the north valley and the Bhuntar station. The mean of daily maximum and minimum temperature is utilised by the HBV model as temperature input. The daily discharge series only cover the period 1997-2005 and the data quality is largely unknown. Quality control is visually conducted and points of suspicious errors are not included in model calibration and validation.
The ice thickness data are also a part of the dataset provided by Huss and Farinotti (2012). The DEMs are downloaded from the Hydro1k global datasets (EROS, 1996).

Nigardsbreen basin
The model is run at a daily time step and at a spatial resolution of 100 m. The model is calibrated for the period from 1991 to 2002 using discharge and annual mass balance and is validated for the period from 2003 to 2012. The numerical values of NSE and RME are tabulated in Table 3 which shows a high model efficiency. Fig. 4 visually shows the good performance of the model in reproducing the historical daily flow series. The annual mass balance simulations also fit well with the observations (Fig. 5). The cumulative mass change in the period 1991-2012 is +7255 mm w.e. by observations and +7231 mm w.e. according to the simulation.
As shown in Fig. 6, summer is the season with the least snow cover, with approximately 84% of the total basin area. For autumn,  Table 3.  the model produces a thinner snowpack than for summer, but a larger snow-covered area. The thinner pattern is because on 31st August, the model removes snow on glaciers and accumulation and ablation are relatively small for the autumn months.
However, the summer pattern is a result of accumulation and ablation effects for each hydrological year. The spatial distribution of snow is in agreement with the seNorge snow modelling product (NVE, 2015).

Chamkhar Chhu basin
The model is run at a daily time step at a spatial resolution of 1 km. The calibration period is from 1998 to 2004 and the validation period is from 2005 to 2008 against discharge. The graphs of observed and simulated discharge are shown in Fig. 7. The NSE values for the calibration and validation are not less than 0.85 (Table 3), which indicates a very good simulation of runoff.
As shown in Fig. 8, the glaciers continued to retreat and mass balance was negative, roughly À0.51 m w.e. year À1 for the period 1998-2008. The maps of seasonal mean snowpack are shown in Fig. 9. They show similar model response as the Nigardsbreen basin.

Beas basin
The model is run at a daily time step on a spatial resolution of 1 km. The model is calibrated against discharge in the period from 1997 to 2002 and validation is performed for the period 2003-2005. The model efficiency is the lowest among the three basins (Table 3). As shown in Fig. 10, the model cannot capture the peaks in the observed hydrograph and significantly overestimates the low flow. The low quality of the data is likely one of the main reasons for this, as the extreme peaks are more than three times of normal high flow and most of them are only sustained for one day.
To further explore the reasons of the relatively low model efficiency, we plot the relationship between annual precipitation and runoff in Fig. 11. The year 2002 can be considered as a turning point and the relationships of precipitation and runoff are substantially modified since this year. The reasons cannot be conclusively stated by examining the annual temperature and precipitation. The shift is possibly due to some sub-daily events or changes in glacier energy balance. For the Himachal Pradesh (Western Himalaya, India) region, Azam et al. (2012) also reported a particularly negative mass balance since the year 2002 compared to the data collected during the period 1987-1989. Fig. 12 shows the simulated snow pattern of the Beas basin. In addition to the similar model responses of the Nigardsbreen and  Chamkhar Chhu basins, the eastern part in the summer receives more snow than in the spring. This confirms that significant snow accumulation occurs during the summer months. The snow pattern of the Beas basin is the most complex among the three basins.

Model performance
The model evaluation is based on discharge and available mass balance measurements. The results given in Table 3 show that the coupled model is able to accurately simulate both daily discharge and annual glacier mass balance in the Nigardsbreen basin. Though there are no mass balance measurements in the Chamkhar Chhu basin, the numerical criteria of NSE are higher than or equal 0.85 for the calibration and validation modes, which is considered as very good simulation. The model efficiency of the Beas basin is much lower than for the other two basins, but results are considered as acceptable, considering the availability and quality of the data.
Mapping the spatial distribution of internal variables, such as snowpack, is an advantage of distributed modelling. In the three basins, the snow storage increases from the lower elevations to the higher elevations. As stated before, the model transforms the snow remaining on glaciers to ice by 31st August, which means   that snow on glaciers is zero at the end of every summer. However, if glaciers are misrepresented by the model, the simulated snow could be misleading and should be interpreted carefully. The misrepresentation can be introduced by errors in the initial ice thickness data or generated during simulating glacier area.

Glacier runoff
The glacierised areas are considered to be very valuable for water resources and hydropower management (Barnett et al., 2005). Most glaciers are located at high elevation, which results in a high hydropower potential (hydraulic head). Additionally, the evaporation rate is small and the specific runoff is high in these areas.
Glacier runoff is defined as the runoff from the glacierised area, including all melt of snow and glacier ice/firn, and rain water in the most general sense Bliss et al., 2014). Fig. 13 shows the monthly mean runoff of the three basins. The ranking of glacier contribution to runoff in the descending order is the Nigardsbreen basin (92.5%), the Chamkhar basin (48.1%) and the Beas basin (27.5%). The main reason for the highest contribution for the Nigardsbreen basin is its largest glacierisation. Additionally, the glacier contribution also depends on climate, especially the spatial distribution of precipitation within a watershed. The more precipitation falls on the glacierised areas, the higher the glacier's contributions to runoff. As shown in Fig. 14, the precipitation decreases from upstream to downstream in the Nigardsbreen basin and in the Chamkhar Chhu basin. However, in the Beas basin, more precipitation falls in the valleys compared to the high mountains. As a result of these spatial distribution patterns of precipitation, the ratio of the glacier's contribution is higher in the Chamkhar Chhu basin than in the Beas basin, even though the latter has a higher glacierisation.
There is no doubt that runoff of these basins is strongly defined by both glaciers and precipitation. Glacier melting is dominant in the Nigardsbreen basin, whereas precipitation is the main contributor in the Chamkhar Chhu and Beas basins. For a given increase in temperature, the three basins are expected to respond differently in terms of runoff and glacier existence. The sensitivity of a basin to climate change and effects on runoff can be studied by running the model with climate projections and this is an on-going research.

Uncertainties
Uncertainties can inherit from data, models and their parameters; small uncertainties can accumulate into considerable uncertainties in the target outputs through the modelling processes (Seiller and Anctil, 2014;Bastola et al., 2011).
The quality of climate data is very important in hydrological modelling. In the studies of the Chamkhar Chhu and Beas basins, the meteorological stations are sparsely distributed and mainly located at low elevation places. This distribution leads to a low representativeness of the spatial precipitation and temperature. Though in the interpolation of the station measurements, elevation effects have been considered, the exposure to the sun also influences snow and ice melt. For very high mountain glaciers, air temperature is seldom above the melting point and due to thin atmospheres, it is a humble indicator of energy (Hock, 2003;Sicart et al., 2008). The accumulated ice flows to lower elevations due to ice dynamics. The accumulation and melting should be carefully interpreted in their total amount rather than the spatial details.
The discharge series are essential to calibrate and validate hydrological models. Therefore, their accuracy and length should be adequate. However, the available discharge series of the Himalayan basins typically have a length of only ten years. In addition, the accuracy is affected by the measurement methodology and quality is not ensured, particularly for the Beas basin.
Initial conditions provide a starting point for a model system. The initial conditions of the hydrological models are determined by model simulations, usually referred to as model ''spin-up''. However, this method cannot be used to construct the glacier state in current climate. A glacier forms over many years and recent warming has induced melting for most glaciers globally. Considerable uncertainties exist in the initial state of glaciers in the three basins.
For the Dh-parameterisation, a limitation is that this model is only valid for retreating glaciers, which is generally true for most mountain glaciers at present and for the near future. At present no scheme is implemented that allows describing the expansion of glaciers following positive mass change; they just increase the total glacier volume. The four parameters of the Dh-parameterisation are only calibrated according to discharge in the Chamkhar Chhu and Beas basins. How well they are identified needs further evaluations. For the HBV parameters, the optimised values are given in Table 4, which shows that the parameters values are in general within their empirical ranges except few parameters for the Beas basin.

Conclusions
Glaciers strongly influence hydrology in glacierised areas. However, these effects have not been well understood and well simulated quantitatively at present due to imperfect tools and data scarcity. In this study, we integrated a simple approach for simulating glacier geometry change, the Dh-parameterisation into the HBV model, to study the dynamics of hydrological processes in glacierised basins. The combined model requires easily accessible inputs data (precipitation, temperature and initial ice thickness) and provides the dynamics of glacier system and consequent runoff processes. The model can be used for estimating runoff, snow and glacier mass balance in past and future.
The coupled model is tested in three high mountain basins with different climates and data quantity and quality. Among the basins, the Nigardsbreen basin has the longest time series of discharge data and observed annual glacier mass balance. Comparisons of model simulations with the observations show that the model yields very high model efficiency in simulation of discharge and annual mass balance for the Nigardsbreen basin. The least model efficiency is found for the Beas basin possibly due to uncertainties in input data and the changed precipitation-runoff relationship since the year 2002. Additionally, the model can provide maps of snowpack distribution and estimate runoff components from glaciers.