Water Balance Assessment under Di ﬀ erent Glacier Coverage Scenarios in the Hunza Basin

: The potential impact of glacier recession on river discharge from the Hunza river basin was estimated as an indicator for downstream changes in the Indus river system. The J2000 model was used to analyze the water balance in the basin and simulate the contribution of snow and ice melt to total discharge at present and under three scenarios of glacier recession. Precipitation was corrected using virtual weather stations created at a higher elevation and a precipitation gradient. Snowmelt from the whole basin contributed, on average, 45% of the total river discharge during the modeling period and 47% of the ice melt from the glacier area. Total ice melt declined by 55%, 81%, and 96% under scenarios of glacier recession to 4000, 4500, and 5000 masl, respectively. The contribution of ice melt to river discharge decreased to 29%, 14%, and 4% under the three scenarios, while total discharge from the Hunza river decreased by 28%, 40%, and 46%. The results suggest that glacier recession in the Hunza river basin could have serious implications for downstream water availability. Understanding melt contribution in the basin based on ongoing and projected future climatic change can play a crucial role in future water resource management.


Introduction
The Hindu Kush Himalaya (HKH) has some of the largest glaciers in the world [1]. They act as natural water storage-reservoirs that store precipitation in the winter in the form of snow and release it in the summer as meltwater [2]. A major part of the flow in the Indus river system in the western Himalaya comes from snow and glacier fed river catchments in the Karakoram region [3]. Modeling studies suggest that the contribution of snow and glacier melt in the Indus river is as high as 50-80% [4][5][6][7].
Climate change is expected to affect various components in the hydrological cycle [8]. Recent changes identified at high elevations in the Karakoram include shifts in seasonal temperatures, snowfall, and snow cover [9]. Temperature change is expected to adversely affect the glacier and ice reserves of the Himalayan region with a potential shift in the equilibrium line altitude (ELA), recession of glacier termini to higher elevations, and reduction in glacier area and ice volume [10]. Khattak, et al. [11] reported an increase in winter maximum temperatures between 1976 and 2005 in the upper, the middle, and the lower parts of the Indus basin of 1.79, 1.66, and 1.20 • C, respectively. Nepal and Shrestha [12] also reported an overall gradual increase in temperature in the Indus basin but with some differences in the reported seasonal trends. Bolch, et al. [13] project that HKH glaciers will lose more than one-third of their volume, even if warming is kept to 1.5 • C. In the first instance, rising temperatures are likely to lead to an increased melt rate and runoff [14]. Archer [15] suggests that a 1 • C rise in mean summer temperature would result in a 16% increase in summer runoff into the Hunza and the Shyok rivers The Pakistan Water and Power Development Authority (WAPDA) has installed three climate stations with precipitation gauges at different elevations, as shown in Figure 1 and Table 2. Additional virtual weather stations were created for this study, and they are described below in Section 2.3.2. The Hunza basin has an arid to semi-arid climate characterized by cold winters and hot summers at lower altitudes, with a wide variation between temperature extremes. The high mountains limit the intrusion of the monsoon, whose influence weakens northwestward [36]. The westerly circulation brings weak intensity precipitation in both summers and winters, which is the primary source of precipitation, contributing about two-thirds of high-altitude snowfall in the Karakoram region [7]. The basin has close to 80% snow cover in winters and 30% in summers [3]. About 90% of the total glaciated area in the Karakoram range lies at 5000-6000 masl, and most of this area is in the accumulation zone [37]. Maximum precipitation occurs at an elevation of around 5500 masl [37,38].
The distribution of elevation range in a catchment area is crucial in hydrological modeling. Table  1 shows the hypsometry of the Hunza basin. Close to 80% of the basin area lies in the elevation range of 3500-5500 masl, thus any changes that occur in this elevation range could have a significant impact on the basin's hydrological dynamics. Similarly, more than 50% of the glacier area lies in the elevation range of 4500-5500 masl, and this elevation range is particularly important for glacio-hydrological processes. The distribution of elevation range in a catchment area is crucial in hydrological modeling. Table 1 shows the hypsometry of the Hunza basin. Close to 80% of the basin area lies in the elevation range of 3500-5500 masl, thus any changes that occur in this elevation range could have a significant impact on the basin's hydrological dynamics. Similarly, more than 50% of the glacier area lies in the elevation range of 4500-5500 masl, and this elevation range is particularly important for glacio-hydrological processes.

The J2000 Model
The J2000 hydrological model is a process-oriented hydrological model for hydrological simulations of mesoscale and macroscale catchments [39,40] implemented in the Jena Adaptable Modeling System (JAMS). JAMS is a software framework for component-based development and application of environmental models [41,42]. The simulation of different hydrological processes is carried out in encapsulated process modules, which are, to a great extent, independent of each other. The flexibility of this model allows changing, substituting, or adding of individual modules or processes without having to restructure the model from the start. A glacier module was integrated into the J2000 modeling system to simulate glacier runoff from a glacierized area. Figure 2 shows the system of different modules within the J2000 model.  [40] and Nepal et al. [33]) . Note: ET = evapotranspiration, P = precipitation, T = air temperature, W = wind speed, RH = relative humidity, SH = sunshine hours, LPS = large pore storage, MPS = middle pore storage, DPS = depression storage, RO = runoff.
The model takes into account all the important hydrological processes and components, such as interception, evapotranspiration, snow and glacier melt, soil water, groundwater, and routing. Details of the modeling application, model input data, and calibration parameters are provided in [33].  The model takes into account all the important hydrological processes and components, such as interception, evapotranspiration, snow and glacier melt, soil water, groundwater, and routing. Details of the modeling application, model input data, and calibration parameters are provided in [33]. Table 2 shows the hydro-meteorological data available for the Hunza basin in the installed stations and a gauging station at Dainyor bridge. Wind speed (m/s) data were obtained from the Cultural Area of Khunjerab station installed by the University of Bonn close to Ziarat station. Sunshine hours data were obtained from the Meteorological Department of Pakistan. All data had a daily temporal resolution. Note: P = precipitation (mm), T min, max = minimum and maximum temperature ( • C), RH = relative humidity (%).

Precipitation Correction Using Lapse Rate
Precipitation uncertainty is one of the main challenges for hydrological simulation in the Hunza basin. The meteorological stations are all valley-based and lie below the elevation of maximum precipitation in the Karakoram region at 5000-6000 masl [43][44][45], thus the precipitation recorded at the stations is likely to represent under-catch [27]. To address this, three virtual weather stations (VWSs) were created at elevations around the maximum precipitation zone: VWS 1 at 5351 masl, VWS2 at 5393 masl, and VWS3 at 5364 masl (for location, see Figure 1). Precipitation at the VWSs was calculated using a vertical precipitation gradient (PG) of 0.04% m −1 , as suggested in various studies [46][47][48][49][50], using the values from a linked observation station at lower elevation as the base. The precipitation at VWS1 was calculated using precipitation at Naltar station as a reference, at VWS2 using precipitation at Khunjerab station as a reference, and at VWS3 using precipitation at Ziarat station as reference. Table 3 shows details of the stations and average annual precipitation measured or calculated at each. The average annual precipitation in the Hunza basin during the modeling period was calculated to be 731 mm (using data from both the installed stations and the virtual weather stations). This is close to the estimate of 795 mm for the Hunza river catchment made by Khan and Koch [51] using a new approach for interpolation and regionalization of observed precipitation data across the upper Indus basin with adjustment for orographic effects and changes in glacier storage. All six precipitation stations were used as precipitation input to the J2000 model, i.e., three physical stations and three virtual stations.

Hydrological Response Units (HRU)
The land cover data [52] were derived from the 300 m global land cover map from the European Space Agency GlobCover initiative. The fifteen GlobCover land cover classes found in the study area were reclassified to five classes that have similar effects on hydrological dynamics at a resolution of 500 m. The permanent snow and ice layers in GlobCover were combined and designated as bare land, because snow accumulation is calculated within the J2000 model. A glacier layer derived by ICIMOD [53] for the period 2005 ± 3 was overlaid on the land cover data.
Harmonized World Soil Database (HWSD) version 1.2 data [54] were used to determine the soil parameters in the model. Four major different soil types were found in the study area. Due to the lack of a geological dataset, geological information for the study area was derived from the basic geological characteristic of the region that controls maximum percolation rates and ground water storage. Three regional classes of geological information from the soil data and the literature were used. In the model, there is no infiltration or percolation to the soil in glacier areas, thus these are regarded as no soil areas.
The Shuttle Radar Topographic Mission (SRTM) digital elevation model (DEM) at 90 m resolution from the CGIAR Consortium for Spatial Information [55] was used to define the topography and delineate watershed boundaries and hydrological response units (HRUs). The preparation and the quantification of model parameter files were done in a GIS environment. For preparation of the soil parameter file, soil texture information in percentages was provided as input data to the software "HYDRUS 1D" to understand the pedo-transfer function of soils under three different pressure scenarios (0 mbar, 60 mbar, and 15,000 mbar). The parameter files and their values were static in the modeling application. All the input data (i.e., soil, land-cover, and DEM) were prepared in raster format at the same resolution of 500 m; this mainly controls the number of HRUs to be formed without losing the heterogeneity of a catchment. These data were used to delineate HRUs using the web-based tool (http:/intecral.uni-jena.de/hruweb).

Calibration with Snow Cover
The J2000 simulated snow cover was compared with Moderate Resolution Imaging Spectroradiometer (MODIS) snow cover data [56] and the comparison used to calibrate the snow and glacier related parameters. The snow cover comparison was carried out from March 2000 to November 2010. MODIS snow products have been used by others to understand snow extent in the Indus basin [3,4,14]. Both the coefficient of determination (R 2 ) and visual inspection were used to compare the modeled snow extent with MODIS snow extent.

Calibration and Validation with Discharge
The ice melt and rainfall-runoff parameters were then calibrated by comparing the simulated discharge with the observed discharge measured at Dainyor bridge. Simulated discharge data from the model for the period 2000-2004 were used for calibration and for 2008-2010 for validation. The model performance was evaluated using both efficiency criteria and visual inspection related to systematic (over and under prediction) and dynamic (timing, rising/falling limb, and base flow) behaviors of the model. A combination of four different efficiency criteria-Kling-Gupta efficiency (KGE), Nash-Sutcliffe efficiency (ENS), logarithm of Nash-Sutcliffe (LNS), and coefficient of determination (R 2 )-was used, as a single criterion often cannot provide a complete picture of model performance [57]. KGE is a decomposition of ENS that facilitates analysis of the relative importance of its different components in the context of hydrological modeling. ENS uses the observed mean as a baseline, which can lead to

Discharge under Different Glacier Recession Scenarios
The Karakoram glaciers are in a transition phase from positive to negative mass balance [60]. There is a slightly negative mass balance in glaciers up to 5000 masl elevation [61]. Although most of the terrain is high enough to have temperatures sufficiently far below freezing to maintain glacier mass, climate models project that the warmer climate will lead to a reduction in annual snowfall of 20-40% in the upper Indus [62]. Reference [11] reported that the winter maximum temperature in the upper Indus basin increased at an average rate of 1.79 • C per 39 years between 1987 and 2005. If this type of temperature increase continues, it could lead to increased glacier melt and recession of glaciers, converting the glacier area to bare land.
This study considered three future scenarios for the upward movement of glacier termini: shift to 4000 masl (Scenario 1), 4500 masl (Scenario 2), and 5000 masl (Scenario 3). The impact on glacier area under each of these scenarios is shown in Table 4 and Figure 3.

Scenario Conditions
Change from Baseline Scenario 1 Glacier termini recede to 4000 masl Glacier area decreases by 10% Scenario 2 Glacier termini recede to 4500 masl Glacier area decreases by 22% Scenario 3 Glacier termini recede to 5000 masl Glacier area decreases by 41%

Snow Cover Simulation
The J2000 model's ability to simulate seasonal snow cover was assessed by comparing the simulated snow cover area (calculated using the corrected precipitation values) with MODIS snow

Snow Cover Simulation
The J2000 model's ability to simulate seasonal snow cover was assessed by comparing the simulated snow cover area (calculated using the corrected precipitation values) with MODIS snow cover values (Figure 4). In J2000, the maximum snow cover area of the 8-day interval matching the MODIS time period was chosen. The snow cover area simulated by the model was in good agreement with the area inferred from the MODIS snow cover data (R 2 = 0.65; positive bias of 6%). The MODIS data for the Hunza basin show a maximum average snow cover area of 81% (11,213 km 2 ) in March and a minimum average snow cover area of 42% (5835 km 2 ) in August for 2000-2010. The simulated values from the J2000 model for maximum snow cover area (11,904 km 2 ) in March were similar to the observed values; those for the minimum snow cover area (4128 km 2 ) in August were lower than the observed values, but there was good agreement in the period of snow melt. In some years (2002,2005,2008), the snow accumulation was underestimated by the J2000 model. Overall, the snow cover area was simulated well by the J2000 model, but snow cover variability was less well represented.

Hydrograph Analysis during Calibration and Validation
A split sample procedure was used for calibration and validation of the simulated discharge values against observed discharge [63]. The period for which data was available was not very long, thus a larger part (2000)(2001)(2002)(2003)(2004) was used for calibration to ensure calibration was meaningful and the remainder [2008-2010 (no data available for 2005-2007)] was used for validation. Figure 5 shows the simulated and the observed daily discharge values for the calibration and the validation periods together with the corrected precipitation, and Figure 6 shows the average monthly simulated and observed discharge values and a scatter plot of the daily observed and simulated discharge values. Table 5   The MODIS data for the Hunza basin show a maximum average snow cover area of 81% (11,213 km 2 ) in March and a minimum average snow cover area of 42% (5835 km 2 ) in August for 2000-2010. The simulated values from the J2000 model for maximum snow cover area (11,904 km 2 ) in March were similar to the observed values; those for the minimum snow cover area (4128 km 2 ) in August were lower than the observed values, but there was good agreement in the period of snow melt. In some years (2002,2005,2008), the snow accumulation was underestimated by the J2000 model. Overall, the snow cover area was simulated well by the J2000 model, but snow cover variability was less well represented.

Hydrograph Analysis during Calibration and Validation
A split sample procedure was used for calibration and validation of the simulated discharge values against observed discharge [63]. The period for which data was available was not very long, thus a larger part (2000)(2001)(2002)(2003)(2004) was used for calibration to ensure calibration was meaningful and the remainder [2008-2010 (no data available for 2005-2007)] was used for validation. Figure 5 shows the simulated and the observed daily discharge values for the calibration and the validation periods together with the corrected precipitation, and Figure 6 shows the average monthly simulated and observed discharge values and a scatter plot of the daily observed and simulated discharge values. Table 5 shows the values of the different efficiency criteria derived from comparison of the simulated daily discharge values with observed values over the calibration and the validation periods using the observed precipitation values and the corrected precipitation values. For the water balance assessment (discussed in the later sections), we took the data from the period of 2000-2004 and 2008-2010 (a total of eight years).    The graphical and the statistical evaluation showed that the J2000 model was able to reproduce the overall hydrological dynamics fairly well. The base flow was well simulated during both the calibration and the validation periods. The rising and the falling limbs were also well simulated by the model. However, there were some high peaks shown in the simulated discharge in April 2002 and 2003 that were not identified in the observed discharge. Both visual inspection and the efficiency criteria values indicated that the simulated discharge values (especially base flow and discharge peaks) were considerably closer to the observed values when using the corrected precipitation in both the calibration and the validation periods.

Contribution of Ice Melt to Total Discharge
In the model, glacier melt runoff (melt from the glacier area) is the sum of snowmelt runoff (seasonal snowfall), ice melt runoff, and rain runoff (from rainfall over the glacier area). Glacier ice melt begins in a glacier HRU after the seasonal snowfall has melted and the seasonal snow storage on that HRU is zero. Figure 7 shows the simulated monthly average contribution of the ice melt component to total simulated discharge from the basin. The monthly values of proportional contribution to basin total discharge for the individual components are shown in Table 6. The graphical and the statistical evaluation showed that the J2000 model was able to reproduce the overall hydrological dynamics fairly well. The base flow was well simulated during both the calibration and the validation periods. The rising and the falling limbs were also well simulated by the model. However, there were some high peaks shown in the simulated discharge in April 2002 and 2003 that were not identified in the observed discharge. Both visual inspection and the efficiency criteria values indicated that the simulated discharge values (especially base flow and discharge peaks) were considerably closer to the observed values when using the corrected precipitation in both the calibration and the validation periods.

Contribution of Ice Melt to Total Discharge
In the model, glacier melt runoff (melt from the glacier area) is the sum of snowmelt runoff (seasonal snowfall), ice melt runoff, and rain runoff (from rainfall over the glacier area). Glacier ice melt begins in a glacier HRU after the seasonal snowfall has melted and the seasonal snow storage on that HRU is zero. Figure 7 shows the simulated monthly average contribution of the ice melt component to total simulated discharge from the basin. The monthly values of proportional contribution to basin total discharge for the individual components are shown in Table 6.  Figure 7. Monthly average contribution of ice melt to total discharge during the modeling period. Table 6. Average monthly melt contribution from glacier area to total discharge from the Hunza basin. Melt from the glacier area contributed 60% (on average) to the total discharge from the basin during the modeling period; 47% of the total from ice melt, 12% from seasonal snowmelt, and 1% from rain runoff. The maximum average contribution of glacier melt to total discharge was in August (73%); there was no contribution in December, January, or February. The seasonal contribution of ice melt was particularly significant at 52% of the average monthly discharge in the summer period (June to October).

Month Melt Runoff from Glacier Area (mm) Snowmelt (mm) Ice Melt (mm) Rain Runoff (mm)
Snowmelt from outside of the glacier area comprised 39% of total runoff in the basin, but a part of this infiltrated to soil and evaporated (about 6%), thus the total contribution of snowmelt from outside the glacier area to discharge was 33%. The contribution of all snowmelt to total discharge was 45% (33% snowmelt from outside the glacier area and 12% snowmelt from the glacier area). Table 7 shows the contribution of snowmelt, glacier ice melt, and total snow and glacier melt to discharge identified in various modeling studies in and close to the Hunza basin. Glacier ice melt is given separately to meet the different definitions used in the studies. For example, Refs. [14,29] do not include snowfall or rainfall on the glacier surface and only consider glacier ice melt, whereas the present study used seasonal snowfall and rainfall as well as glacier ice melt to assess the contribution to discharge from the glacier area (contribution to glacier ice melt alone was 47%, and both snow and ice melt was 59%). Our results suggest that snow and glacier melt from the whole basin contribute 92% of total discharge, which is close to the values given by [4,14,29], the studies that provide the most direct comparison in terms of area and basin size. The approach used by [4] was similar to that used in the J2000, with seasonal snowfall on the glacier taken into consideration, although the value calculated for glacier ice melt in our study was slightly higher. The values calculated for glacier ice melt by [14,29] were higher again than those calculated using our approach because of the different methods used to realize glacier melt. Note: * the percentage contribution of ice melt runoff from the glacier area to river discharge (excludes snowmelt and rain runoff from the glacier area).

Water Balance Analysis
The results of the overall water balance analysis for the Hunza basin are shown in Figure 8. The average annual input to total water during the modeling period was 731 mm (70%) from precipitation and 315 mm (30%) from glacier ice melt, giving a total of 1046 mm. Actual evapotranspiration (actEt) returned 162 mm (16%) of total input to the atmosphere, 672 mm (64%) left the basin as total river discharge, and 20% remained in different forms of storage, such as snow, soil, and groundwater.
average annual input to total water during the modeling period was 731 mm (70%) from precipitation and 315 mm (30%) from glacier ice melt, giving a total of 1046 mm. Actual evapotranspiration (actEt) returned 162 mm (16%) of total input to the atmosphere, 672 mm (64%) left the basin as total river discharge, and 20% remained in different forms of storage, such as snow, soil, and groundwater.  Figure 9 shows the results of the analysis of the impact of glacier recession on ice melt runoff as well as the resultant total discharge calculated using the J2000 model. The analysis focused on the impact on ice melt, as this is likely to be the component most affected. The snowfall and rainfall components are not affected by the change in glacier area (the only component changed in the scenario), as they will simply runoff from bare land instead of the glacier surface. In a warming world,  Figure 9 shows the results of the analysis of the impact of glacier recession on ice melt runoff as well as the resultant total discharge calculated using the J2000 model. The analysis focused on the impact on ice melt, as this is likely to be the component most affected. The snowfall and rainfall components are not affected by the change in glacier area (the only component changed in the scenario), as they will simply runoff from bare land instead of the glacier surface. In a warming world, the change in temperature will also affect snowfall distribution and other hydrological processes, but in this study, we only looked at the glacier recession scenario.

Impact of Different Scenarios on Ice Melt Runoff and River Discharge
The simulated monthly average ice melt runoff during the modeling period was used as the baseline for the glacier shrinkage scenarios. Total ice melt declined by 55%, 81%, and 96%, respectively, under the scenarios of recession to 4000, 4500, and 5000 masl; the average annual contribution to total discharge declined from 47% at the baseline to close to zero (29%, 14%, and 4%, respectively), and the total average annual river discharge declined by 28%, 40%, and 46%, respectively. The reduction was mainly in summer discharge. Discharge actually increased very slightly in some months during the accumulation period (December, January, February, and March) because, as glaciers recede, the snow falls on bare land, infiltrates after melt, and contributes to river discharge through interflow and groundwater.

Uncertainties and Limitations
Climate data quality plays a very crucial role in hydrological modeling. The climate stations in the Hunza basin are valley-based and do not sufficiently represent the spatial distribution of precipitation in the mountain range. The precipitation in this study was corrected using a precipitation gradient, but these were calculated values, and some uncertainty remains. The validation of the J2000 model simulated snow cover area using MODIS snow cover data helped to constrain the snow and glacier related parameters in the model, which reduced the parametric uncertainties. However, the MODIS snow cover data also contained some uncertainty, as they were limited to 500 m resolution over an 8-day period. Finally, both accuracy and length of discharge data are essential to calibrate and validate hydrological models, but limited data availability meant that the typical length of the modeling period in this study was only eight years. the change in temperature will also affect snowfall distribution and other hydrological processes, but in this study, we only looked at the glacier recession scenario. The simulated monthly average ice melt runoff during the modeling period was used as the baseline for the glacier shrinkage scenarios. Total ice melt declined by 55, 81, and 96%, respectively, under the scenarios of recession to 4000, 4500, and 5000 masl; the average annual contribution to total discharge declined from 47% at the baseline to close to zero (29%, 14%, and 4%, respectively), and the total average annual river discharge declined by 28%, 40%, and 46%, respectively. The reduction was mainly in summer discharge. Discharge actually increased very slightly in some months during the accumulation period (December, January, February, and March) because, as glaciers recede, the snow falls on bare land, infiltrates after melt, and contributes to river discharge through interflow and groundwater.

Conclusions
The key findings of the analysis of the water balance in the Hunza river basin using the J2000 hydrological model and precipitation values corrected using virtual weather stations were as follows: • The snow cover area of the basin simulated by the model was in good agreement (R 2 = 0.65) with the snow cover area calculated from the 8-day interval MODIS data. The snow cover validation helped to constrain the snow and glacier related parameters, which gave additional confidence to the water balance simulations compared to validation of discharge data only. • The J2000 model can be used successfully for snow and glacier melt simulation in the western Himalaya. Snowmelt from the whole Hunza basin contributed, on average, 45% of the total river discharge during the modeling period, and glacier ice melt contribution was 47%.

•
The contribution of ice melt to river discharge decreased to 29%, 14%, and 4% under the three scenarios, respectively, while total discharge from the Hunza river decreased by 28%, 40%, and 46%.
The drastic reduction in ice melt contribution to river flow with glacier recession suggests that glacier storage in the Hunza basin is crucial for sustaining river flow. The Hunza river is one of the main contributors of snow and glacier melt to the Indus river. If storage capacity is lost, the flow in the whole Indus river could be affected adversely, which could have a marked impact on all the downstream systems that depend on the river water.
Funding: This research received no external funding.