Climate Change Impacts on the Upper Indus Hydrology: Sources, Shifts and Extremes

The Indus basin heavily depends on its upstream mountainous part for the downstream supply of water while downstream demands are high. Since downstream demands will likely continue to increase, accurate hydrological projections for the future supply are important. We use an ensemble of statistically downscaled CMIP5 General Circulation Model outputs for RCP4.5 and RCP8.5 to force a cryospheric-hydrological model and generate transient hydrological projections for the entire 21st century for the upper Indus basin. Three methodological advances are introduced: (i) A new precipitation dataset that corrects for the underestimation of high-altitude precipitation is used. (ii) The model is calibrated using data on river runoff, snow cover and geodetic glacier mass balance. (iii) An advanced statistical downscaling technique is used that accounts for changes in precipitation extremes. The analysis of the results focuses on changes in sources of runoff, seasonality and hydrological extremes. We conclude that the future of the upper Indus basin’s water availability is highly uncertain in the long run, mainly due to the large spread in the future precipitation projections. Despite large uncertainties in the future climate and long-term water availability, basin-wide patterns and trends of seasonal shifts in water availability are consistent across climate change scenarios. Most prominent is the attenuation of the annual hydrograph and shift from summer peak flow towards the other seasons for most ensemble members. In addition there are distinct spatial patterns in the response that relate to monsoon influence and the importance of meltwater. Analysis of future hydrological extremes reveals that increases in intensity and frequency of extreme discharges are very likely for most of the upper Indus basin and most ensemble members.


Introduction
The water resources supplied by the upper Indus basin (UIB) are essential to millions of people and future changes in both demand and supply may have large impacts [1]. The UIB provides water for the world's largest continuous irrigation scheme through several large reservoirs (e.g. the Tarbela and Mangla dams, Fig 1), which depend for more than 50% of their annual inflow on snow and glacier melt water [2][3][4][5]. In combination with variable precipitation patterns, the intra-annual variation in streamflow is high [6] and so is the supply to the downstream areas. Water demands are high, primarily because of water consumption by irrigated agriculture [7], and hydropower generation [8]. At the same time the downstream part of the basin is characterized by very dry conditions [9], making it largely dependent on water supply from the upstream areas. The downstream demands exceed the supply and on an annual basis groundwater resources are depleted by an estimated 31 km 3 [10], which makes the Indus basin aquifer the most overstressed aquifer in the world [11,12]. The uncertainty in the future mountain water resources combined with the Indus basin's large dependency on these upstream water resources, makes the Indus basin a climate change hotspot [1,13]. This vulnerability is enhanced by an expected large regional population growth in the coming decades [14], associated with increases in water and energy demand in the basin. Extreme weather events are likely to become more frequent in the region in the future [15,16], which poses serious threats for a region which is already facing severe flooding events [17] and other natural hazards.
The UIB has a complex climate. Several studies investigated historical trends in precipitation and temperature in the UIB. Trend analysis on precipitation for 17 stations throughout the UIB showed statistically significant increasing trends in precipitation for several stations in annual, summer and winter precipitation between 1961 and 1990 [19]. Air temperature trends between 1961 and 2000 were also assessed and it was found that (i) the diurnal temperature range is increasing consistently in all seasons, (ii) winter mean and maximum temperatures show significant increases and (iii) mean and minimum summer temperatures show a consistent declining trend [20]. These findings were confirmed also for a more recent period  for roughly the same stations [21]. Trend analysis on the ERA40 reanalysis dataset for the Baltoro region in the Karakoram (Fig 1) showed negative summer temperature trends from 1958 until 1990 and a positive trend from 1991 to 2001 [22]. The authors also found an increasing trend in annual precipitation from 1970 to 1990 and a decreasing trend during the 1990s. Trend analysis on several gridded precipitation products did not confirm these findings [23]. Studies on the winter westerly disturbances, being the major source of winter precipitation, indicate strong intra-seasonal variability and a trend of enhanced frequency and strength of these disturbances in the Karakoram and western Himalaya between 1979 and 2010, leading to increased heavy winter precipitation [24,25].
There is great debate on the response of glaciers in the UIB to climate change during the last decade. The glaciers in the Himalayan range are seemingly losing mass at rates similar to other mountainous regions in the world, however the glaciers in the Karakoram and Pamir mountain ranges have neutral mass balances on average and are characterized by a large number of surging glaciers [26][27][28][29][30][31]. This so called Karakoram anomaly has not been explained, but a possible reason could be a combination of a decrease in summer temperatures and an increase in precipitation. However this is still speculative and requires further study and understanding of atmospheric processes leading to high-altitude precipitation. This hypothesis is supported by an increasing trend in snow cover that was found in the Hunza basin based on MODIS snow cover analysis [32,33] and that the water balance of the UIB can be closed without large negative glacier mass balances [6]. On the other hand, decreasing trends in snow cover for the most westerly-influenced subbasins, including Hunza, and increasing trends for the more monsooninfluenced subbasins were found [34]. A trend analysis of snow cover in the monsoon-dominated Sutlej basin indicated a trend of snow cover reduction between 2000 and 2009 [35]. Other research concludes that the Karakoram is protected from reductions in annual snowfall under climatic warming because the seasonal cycle is dominated by non-monsoonal winter precipitation [36].
Rising temperatures in basins strongly dependent on glacier melt are likely to result in an increase in stream flow in the near future and a decline in the far future. This is caused by the fact that the total amount of glacier melt is a trade-off between increasing melt rates on one hand and reduced glacier volumes on the other hand. The moment when the trend in glacier melt changes from positive to negative is highly variable [37,38]. Analysis of a 1961 to 2009 record of reservoir inflow at Tarbela, which is the largest reservoir on the main stem of the Indus river (Fig 1), indicates a declining trend, although statistically insignificant [6]. Further upstream trend analysis on streamflow records at different locations identified stable or declining trends in runoff too [32,39,40]. These studies indicate that large parts of the UIB are (not yet) experiencing accelerated melt, which could indeed be partly attributed to the Karakoram anomaly. However, contrary to these findings, a recent study in the Shigar river basin reports rising river flows [41]. However, the authors do not relate this to the existence of the Karakoram anomaly. Instead, they argue that an increase in runoff is possible under neutral glacier mass balance conditions as a result of increasing temperature and precipitation, i.e. the mass turnover of the glacier is increasing, yet the mass balance remains neutral.
Climate simulations are used to generate projections of future climate change in the UIB. Analysis of precipitation change signals in a large number of General Circulation Model (GCM) runs indicates that an increase in summer precipitation and on average no significant change in winter precipitation are likely [42]. However, the spread in the precipitation changes from the GCM ensemble is large, because the complex UIB climate is difficult to simulate [43]. Analysis with regional climate models (RCM) reveals consistent warming until the end of the century with greater warming in the upper Indus than in the lower Indus. Precipitation projections show a non-uniform change with increases projected for the upper parts and decreases for the lower parts [44,45]. However care needs to be taken in using RCMs directly in impact studies. A recent study that analyzed the uncertainty of the CORDEX South Asia regional climate models showed that the RCMs exhibit large uncertainties in both temperature and precipitation, that they exhibit a large cold bias and that they are unable to reproduce observed warming trends [46]. Empirical-statistical downscaling, which may be better suited under such complex conditions, is another approach to generate forcing for climate change impact models, where climate model output is statistically corrected using transfer functions with local observations during a historical period. Empirical-statistical downscaling of GCMs in the UIB based on an ensemble of selected GCMs showed a modest increase in precipitation and a consistent warming, which is stronger in the upper parts of the basin [3,37]. The application of a stochastic weather generator to downscale RCM data in the northern UIB lead to a projection of yearround increasing precipitation, with increased intensity during the wettest months and yearround uniformly increasing temperatures [47].
Hydrological impact studies have been conducted for the UIB at various spatial scales and key assumptions in those studies relate to (i) the reference climate dataset being used, (ii) the future climate forcing and downscaling method, (iii) the type and complexity of the hydrological model, (iv) the treatment of glacier evolution in the future and (v) the calibration and validation strategy. Hydrological projections based on different approaches indicate likely increases in flow at least during the first half of the 21 st century for particular subbasins [37,38] and at the basin scale [3,44,45,48]. Projections of changes in hydrological extremes in the UIB are very limited [49], but are at the same time very much desired [3,38,50].
In this study we systematically assess the present day hydrology of the UIB and the impacts of climate change using a new fully distributed cryospheric-hydrological model at a high spatial resolution (1 km 2 ) that includes all relevant components of the high altitude water balance [51]. We introduce several novel components which may advance our understanding of the complex impact of climate change on the UIB hydrology: • A new historical precipitation dataset [52] that corrects for the underestimation of high altitude precipitation is used.
• The model is calibrated on river runoff at several locations, as well as MODIS based snow cover estimates and geodetic glacier mass balance data.
• An advanced statistical downscaling technique for climate change scenarios until 2100 is used that accounts for changes in precipitation extremes.
• The analysis is focusing on changes in sources of runoff, changes in seasonality and changes in hydrological extremes.
Pakistan, India and China makes the UIB a transboundary river basin in a geopolitically complex region. It is among the most glaciated areas on Earth, with *22.000 km 2 of glacier surface area [53]. The climate of the UIB is complex and is the result of an intricate interaction between monsoon circulation, westerlies and the topography [46,[54][55][56][57]. The interaction between topography and precipitation manifests itself at various scales ranging from a synoptic scale of several hundreds of kilometers to an orographic meso-scale of less than 30 kilometers [58]. Along the Himalayan arc the monsoon influence is largest but this influence decreases in north-western direction where mid-latitude westerlies become increasingly important, e.g. at the junction of the Karakoram, Pamir and Hindu-Kush mountain ranges (Fig 1). Precipitation from the westerlies is highest in winter when low-pressure systems reach the western margin of the greater Himalaya. This supply of moisture reaches higher elevations than the summer monsoon, which might be related to the higher tropospheric extent of the westerly airflow [55].

Cryospheric-hydrological model
We use the high resolution, raster-based, fully distributed Spatial Processes in Hydrology (SPHY) cryospheric-hydrological model (open source, version 2.0) [51], which was applied in a river basin-scale study on climate change impacts for water availability in five major Asian river basins before [3]. The model runs at 1 km 2 spatial resolution with a daily time step. The actual runoff which is calculated for each grid cell consists of four possible contributing factors: rainfall-runoff, snow melt, glacier melt, and baseflow. For each grid cell the total runoff generated per time step (Q TOT ) is calculated: where Q GM is runoff from glacier melt, Q SM is runoff from snow melt, Q RR is rainfall-runoff and Q BF is baseflow. To determine the contribution of each of the four components to the total runoff within a grid cell, a subgrid parameterization is used in which for each cell the fractional ice cover (G F ), ranging from 0 (no ice cover) to 1 (complete ice cover), is determined. Glacier melt is simulated using a degree-day modelling approach [59]. A differentiation in debris-covered and debris-free glaciers is made based on thresholds for elevation and terrain slope [60], and different degree day factors are used for both glacier types (Table 1). For the remaining fraction of the grid cell, the model maintains a dynamic snow and soil water storage. The glacier fraction per grid cell is adapted dynamically in time. A variable groundwater storage is maintained for the entire grid cell. A part of the glacier melt generated in the glacierized cell fraction is treated as surface runoff and the remaining part is treated as groundwater recharge. Runoff from snow melt consists of the snow melt released from the snow storage, which is simulated using a degree-day modelling approach. Besides accumulation and melt, refreezing of snow melt and rain water within the snow storage is included in the model. Gravitational snow transport between grid cells is simulated with the SnowSlide routine [61]. Snow sublimation is estimated using an elevation-dependent potential sublimation function. We assume that the majority of sublimation comes from snowblown sublimation with highest wind speeds prevailing at higher elevations, and therefore potential sublimation is assumed to increase linearly with elevation above 3000 m a.s.l. by a calibrated factor. The actual sublimation is the potential sublimation limited by the snow storage present in the grid cell. Rainfall-runoff consists of the surface runoff from rainfall and lateral flow released from the soil water storage. Soil water processes are simulated for a topsoil and subsoil layer and processes simulated include evapotranspiration, infiltration, percolation, capillary rise, surface runoff and lateral flow. Baseflow is released from the groundwater storage. Each of these four runoff types is routed downstream using a digital elevation model (DEM) and a routing recession function. A detailed description of the model has been published before [51].

Datasets
Meteorological observations from stations are sparse in the mountains, in particular in the upper Indus. Data mostly originates from valley stations which are not representative for high altitude precipitation. The few stations that are located at higher elevations are typically subject to undercatch in case of snow [62]. Therefore, meteorological datasets consistently seem to underestimate precipitation in the UIB [23,63,64]. As forcing for the SPHY-model we therefore use a corrected precipitation dataset [52], which uses the observation-based APHRODITE [65] dataset as a basis. In the corrected dataset, the raw APHRODITE precipitation data have been corrected by using the glacier mass balance of the major glacier systems as a proxy to estimate high altitude precipitation. Details of the methodology and dataset have been published before [52]. The correction factors that were found in the referred study [52] for 2003-2007 are applied to the daily APHRODITE data for 1971-2000 to generate a 30-year reference climate dataset at 1 km 2 spatial resolution and daily timestep. By using this dataset we aim to overcome the fundamental problem of underestimated precipitation in distributed modeling of highmountain hydrology. Because the 2003-2007 period for which the correction factors were derived [52] does not overlap with the 30-year reference period, and we cannot establish the correction factors for earlier periods due to the lack of IceSAT data [66], we validate the correction factors for the 1971-2000 period by comparing the corrected precipitation amounts to observed discharge. We compare average annual precipitation amounts to average annual discharge amounts during periods falling entirely within 1971-2000 (Table 2). Fig 2 shows that the corrected precipitation amounts are in most cases higher than the observed discharge amounts, whereas the uncorrected precipitation amounts are almost all lower than the observed discharge amounts. Because of the systematic underestimation in high altitude precipitation, we conclude that the corrected precipitation dataset is appropriate to be used as historical precipitation forcing in our study. As digital elevation model (DEM) we use the 15 arc-second void-filled and hydrologically conditioned HydroSHEDS DEM [67], which is based on the SRTM DEM [18]. This DEM is resampled to 1 km 2 spatial resolution. Glacier outlines are extracted from the Randolph Glacier Inventory [68], and they are recalculated to a fractional glacier cover per 1 km 2 grid cell. Land use characteristics are extracted from the MERIS Globcover product [69], and quantitative soil properties are derived from the Harmonized World Soil Database (HWSD, [70]) using pedotransfer functions [71].
MODIS snow cover data [72] and geodetic glacier mass balance data [73] are used for model calibration. IceSat derived glacier mass balance data [27] is used for calibration of a basin-scale parameterization of glacier changes [74]. Discharge observations provided by the Pakistan Water and Power Development Authority (WAPDA) are used for model calibration and validation.

Calibration and Validation
The model is calibrated using a systematic three-step approach to overcome equifinality problems [75,76]. First, parameters related to glacier melt are calibrated using geodetic mass balance data for the Hunza basin (Figs 1 and 3a). The geodetic mass balance data indicates differences in glacier surface elevation, from differencing SRTM [18] and ASTER [77] DEMs. The SRTM DEM was acquired in February 2000, but due to radar penetration it underestimates glacier elevations and is likely to be more representative of the elevation of glaciers at the end of the 1999 melt season [78]. The ASTER DEMs were collected in late September and early October 2008. The elevation differences are transformed to average annual glacier mass balances (m w.e. yr -1 ) for 30 individual glaciers by using an average ice density of 850 kg m -3 [79]. The 30 individual glaciers only include glaciers with a surface area covering at least 5 km 2 (5 model grid cells) to avoid scale problems, and fractional glacier cover for the individual model grid cells are extracted from an updated version of the ICIMOD glacier inventory, which includes distinction of debris-free and debris-covered ice surfaces for the Hunza basin (courtesy of S.R. Bajracharya). Using the model temperature and precipitation forcing, the glacier mass balances for the individual glaciers are simulated for October 1999 to September 2007. Accumulation is calculated as solid precipitation falling on the grid cells with fractional glacier cover and the Table 2. Runoff stations used for validation of the corrected precipitation dataset. Catchment areas are delineated based on the SRTM DEM [18]. P uncor is uncorrected APHRODITE [65]. P cor is corrected APHRODITE using published correction factors, which were derived for 2003-2007 [52]. Numbers in parentheses behind station names are for reference in Figs 1 and 2.

Station Lat Lon Area (km 2 ) Observed Q (m 3 s -1 ) ObservedQ (mm yr -1 ) P uncor (mm yr -1 ) P cor (mm yr -1 ) Period
Tarbela inflow a adjacent grid cells with a slope steeper than 0.2 towards the glacier surface [52]. The simulated data does not coincide completely with the geodetic mass balance data because the forcing data is only available until 2007. The model parameters related to glacier melt (DDFci, DDFdc, see Table 1) are then optimized for agreement between the simulated and observed glacier mass balances and different melt parameters are used for debris-covered glaciers and debris-free glaciers. Parameters are optimized by running the model with different combinations of manually sampled parameter values from the parameter ranges listed in Table 1, and the combination of parameters that yields the best agreement between the simulated and observed glacier mass balances is selected.
Second, parameters related to snow storage and melt (DDFs, SnowSC, Sm, ShdMin, SubPot, see Table 1) are calibrated independently by comparing SPHY simulated snow cover and MODIS remotely sensed snow cover [72]. Remotely sensed snow cover has proven to be useful to improve model calibration in areas with high snow cover [76,80]. The same processed MOD10C2 dataset is used as has been used in another study [81]. From the beginning of 2000 until halfway 2008 the snow cover imagery is averaged for 46 different periods of 8 days (5 days for the last period) to generate 46 different average snow cover maps. That means period 1 is the average snow cover for 1-8 January for 2000 until 2008, period 2 is the average snow cover for 9-16 January for 2000 until 2008, etc. Because the MODIS snow cover product is  Table 2 available at 0.05°x 0.05°spatial resolution, SPHY model snow cover output is averaged over 8 day periods, resampled and projected from 1 km 2 spatial resolution to the same time periods, resolution and geographic projection as the MODIS product. Parameters related to snow melt are optimised to minimize the difference between SPHY simulated snow cover and MODIS observed snow cover. Parameter values are optimized by running different manually sampled parameter values from the parameter ranges listed in Table 1, and the combination of parameters that yields the smallest difference between simulated and observed snow cover is selected.
Third, after calibration of the model parameters related to glacier melt and snow melt, remaining parameters related to baseflow and routing (αGW, kx, see Table 1) are calibrated to observed discharge at two locations in the UIB. The selection of locations is primarily dictated by data availability and data access. Secondly the selection is made such that it is a representative subset of the UIB, with different climatic and hydrological characteristics. Calibration is performed at a daily time step for the same periods for which stream flow records are available. Parameter optimization is done using PEST parameter estimation software (freeware, version 6.05) [82]. The calibrated parameters are assumed to be spatially uniform, i.e. one set of parameters is calibrated and assumed to be applicable to the entire UIB.
The calibrated SPHY model is independently validated to observed discharge at two locations that are not used in model calibration.

GCM downscaling
We select two ensembles containing four GCM runs from the CMIP5 database [83]: one ensemble for the medium stabilization scenario RCP4.5 and one ensemble for the very high radiative forcing scenario RCP8.5. We did not include the mitigation scenario leading to a very low radiative forcing level (RCP2.6). It is unlikely that this RCP can be met, since it requires an immediate drastic decline of emissions followed by ongoing carbon sequestration in the second half of the 21 st century, whereas the future emissions expected to come from existing capital are large [84,85]. As we aim to present robust, realistic projections in our study we choose not to include RCP2.6 in the climate model ensemble. By selecting RCP4.5 and RCP8.5 we cover the entire range of radiative forcing resulting from RCP4.5, RCP6 and RCP8.5. To include a wide range of possible futures and because for our study area there is no particular GCM performing best, and no GCM is able to simulate all aspects of the precipitation dynamics in the region satisfactory [42,86,87], we choose to use the entire range of projections available. For both ensembles we therefore select four GCM runs covering the entire spectrum of projected changes in temperature and precipitation, as projected by all the CMIP5 GCM runs with output available for that RCP (Table 3). We select the models closest to the 10 th and 90 th percentile values of the projections, to avoid the inclusion of outliers, similar as in other studies [3,37,88].
The selected GCM runs are statistically downscaled by applying the Advanced Delta Change (ADC) method [89]. ADC has the advantage over the classical delta change method [90,91] that it is not based on changes in the mean, but changes in the entire precipitation distribution, including extreme precipitation, which is a prerequisite for the assessment of changes in hydrological extremes. This is achieved by applying a non-linear transformation to five-day sums of precipitation data. Five-day sums are considered because extreme discharge events usually depend on multiple days of extreme precipitation. The transformation parameters are determined from the GCM control and future runs. Because of the large difference in resolution between the historical dataset (1 km 2 ) and the GCM data (*1.0 to 2.5°), both datasets are interpolated to a common grid of 25 km 2 resolution. Because ADC focuses on increasing detail in the high end of the precipitation distributions, two different equations are used for the transformation of the observed 5-day precipitation sums, based on the 90% quantile. This quantile is determined per calendar month over the entire reference period for every 25 km 2 grid cell. The two transformation equations are: Where P Ã represents the transformed 5-day sums, P the reference climate dataset 5-day sums, P 90 the 90% quantile and a and b are the transformation coefficients. The superscripts O , C and F denote whether the variable represents respectively the reference climate time series ( O ), the GCM control series ( C ) or the GCM future series ( F ). For 5-day precipitation sums that exceed the P 90 of their month an excess value (E) is determined: E = P − P 90 . The mean future excess (E F ) and mean control excess (E C ) in Eq 3 are determined per calendar month over the entire future or control period for every 25 km 2 grid cell: The linear scaling of the transformed precipitation with the ratio of future and control excess in Eq 3 expresses a change in the slope of the extreme value plot of the five-day maximum precipitation amounts [89].
The transformation coefficients a and b are derived from the 60% and 90% quantiles by: Bias correction factors g 1 and g 2 account for systematic differences in P 60 and P 90 in the reference climate time series and GCM control series, and are determined by: and To reduce sampling variability in the transformation coefficients, the P 60 and P 90 are smoothed temporally by using a weighted mean with weights of 0.25, 0.5 and 0.25 on respectively the previous, current and next month. The mean excesses are smoothed temporally in a similar manner. A detailed description of the ADC-method has been published before [89].
Because the variability in precipitation within a common grid cell in the UIB is much larger than in the Rhine basin, for which the ADC-method was originally developed, the a and b parameters are additionally capped to avoid unrealistically high transformed daily precipitation values, which can occur due to the non-linear transformation of the precipitation value. This is done by constraining the a parameter and associated b parameter as follows: This constraining is based on the distribution of a and b parameter values observed in the transformation in the Rhine basin [92].
The transformation parameters are determined and five-day sums are transformed for each future period spanning 10 years, by using a moving 30-year window from the GCM future series centered around the 10 year future period under consideration. For example, in the calculation of the transformation parameters for 2061-2070, the GCM future series for 2051-2080 is used. For the last future ten-year period (2091-2100) the GCM future series for 2071-2100 is used, similar as for 2081-2090, because most GCM runs do not go beyond 2100. After transformation a change factor can be determined for each five-day sum, which can be subsequently used to transform the individual days that belong to that specific sum. The change factor R is determined as: To generate a baseline daily time series spanning 100 years from 2001 to 2100 that can be subjected to the change factors, a series of 100 years of daily precipitation is randomly selected from the 1971-2000 reference climate dataset (P O d ). The change factor (R) is used to transform the individual daily precipitation sums to future daily precipitation (P F d ): Due to the non-linear transformation of precipitation, the mean climate change signal in the bias-corrected downscaled data is modified from the mean climate change signal in the raw GCM data. Such modification of the mean climate change signal is often observed in statistical bias-correction and downscaling methods and may be considered as an undesired deficiency of a bias-correction and downscaling method [93], although this is a current topic of discussion [94][95][96]. We choose to correct for this effect and therefore the transformed daily precipitation values are scaled for each future ten year period at the grid cell level at monthly scale to the ratio of future and reference precipitation sum according to the raw GCM data as follows: With corP F d being the final transformed daily precipitation value, P F being the future precipitation sum in the GCM future run, P C being the precipitation sum in the GCM control run, P O being the precipitation sum in the reference dataset and P F d being the initially transformed precipitation Eq (11).
The temperature transformation, in contrast to that of precipitation, is linear in nature and has the form [89]: where T Ã represents the transformed temperature; T the temperature in the reference climate dataset; T O , T C and T F the monthly mean of respectively the reference, GCM control and GCM future temperature; σ C and σ C the standard deviations of the daily GCM control and GCM future temperature calculated per calendar month. The temperature transformation is applied to daily temperature values directly. The same series of 100 years of randomly selected years from the reference period as for precipitation is used for the transformation of air temperature data. The transformation is applied to mean, maximum and minimum air temperature separately (i.e.: T in Eq 13 can be replaced by Tmean, Tmax or Tmin). Each of the downscaled GCM scenarios is used to force the hydrological model with transient runs from 1 January 2001 until 31 December 2100.

Future glacier changes
Future glacier changes are simulated at large scale for the UIB divided in three sub-regions: one sub-region for each of the mountain ranges Hindu Kush, Karakoram and Himalaya (Fig 1). For each of these three regions a regionalized glacier mass balance model is used to estimate changes in the regional glacier extent as a function of the glacier size distribution in the subregions and the downscaled future climate data [74]. This glacier mass balance model is specifically developed for implementation in large-scale hydrological models, where the spatial resolution does not allow for the simulation of individual glaciers and data scarcity is an issue. The model is initially forced with the climatic forcing for the reference period and calibrated to subregion-averaged glacier mass balance data derived from IceSAT data [27], before it is used to calculate sub-region-scale glacier changes for each of the downscaled GCM ensemble members from 2001 until 2100. The Randolph Glacier Inventory [68] is assumed to be representative for the state of the glacier extent at the start of the future simulation in 2001.

Calibration and validation
After calibration of the degree-day factors of debris-covered and debris-free glaciers ( Table 1) the area-weighted mean glacier mass balance (+0.11 m we yr -1 ) matches very well with the observed area-weighted mean glacier mass balance (+0.12 m we yr -1 ) (Fig 3b). The interquartile range is also similar. However, the total spread within the sample of 30 individual glaciers is larger in the simulation than in the observations. The larger spread in the simulation stems most probably from the quite coarse model resolution at 1 km 2 . The calibrated values for the degree day factors for debris-free and debris-covered glaciers (Table 1) fall well within the range of values derived in field experiments in the greater Hindu-Kush-Himalayan region [98]. Given the large scale and the fact that we use a fixed parameter set for all glaciers we conclude that the calibrated parameters can be considered representative for the UIB.
Averaged over the UIB, the calibrated SPHY model simulates snow cover reasonably well (Fig 3c and 3d). The largest overestimates occur in the Karakoram range and the Himalayan range in the most southeastern part of the UIB. The largest underestimates occur in the Hindu Kush and mountain ranges to the south of the Karakoram. At the basin scale, there is also a slight overestimation of snow cover during most parts of the year (Fig 3d). Overestimates may well be related to the fact that snow redistribution by wind from one grid cell to another is not included in the SPHY model. Another explanation could be related to the simple approach used to estimate sublimation, whereas sublimation can potentially be an important component of the high-altitude water balance in the HKH region [99]. Studies in other areas revealed that blowing snow sublimation plays a larger role than ground sublimation from the snow pack, and that sublimation losses can be in the order of tens of percents of the total snow accumulation, and up to *90% on very windy ridges [100][101][102].
Calibration to observed discharge shows that averaged over the two locations, the Nash-Sutcliffe efficiency [103] calculated at a daily time step equals 0.81, whereas Pearson's correlation coefficient equals 0.92 (Table 4). For the location at Tarbela, covering a large part of the Indus basin, there is a positive bias of 9.7% in the simulation. The bias is largest during the months with high contribution of snow melt to the discharge (Fig 4a), and is thus likely related to the overestimate of snow cover on the part of the Tibetan Plateau that is part of this catchment (Fig 3d), which in turn relates to the high precipitation forcing in spring. For the Jhelum basin upstream of Mangla reservoir, with a large contribution of snow melt to the stream flow, the model simulates the seasonal patterns in stream flow well (Fig 4b). For the Hunza basin upstream of Dainyor bridge, which harbours the highest and most scarcely monitored part of the UIB and is used for model validation, simulated stream flow is slightly underestimated during the peak season in July and August, and overestimated during September and October (Fig  4c). During these months the stream flow is dominated by glacier melt, which is driven by air temperature. This suggests that the APHRODITE temperature fields may lack some accuracy for this area where observations are very scarce. Besides the model slightly underestimates snow cover in the northern part of the Hunza basin (Fig 3c), which may contribute to the underestimate of the flow peak. For the Chenab basin, located to the southeast of the Jhelum basin, the model underestimates the flow during July and August, leading to a large negative bias (Fig 4d). In this case the bias is most likely related to a shortage of precipitation in the forcing data, being 1222 mm yr -1 for the validation period whereas the observed discharge is only slightly lower (1100 mm yr -1 ).
The annual water balance for 2003-2007, largely coinciding with the period covered by Ice-SAT derived glacier mass balances for three sub-zones in the UIB (Fig 1), is plausible for the Indus upstream of Besham Qila with precipitation input being 664 mm yr -1 , the negative glacier mass balance contributing 25 mm yr -1 , and evapotranspiration, sublimation and discharge being 174 mm yr -1 , 139 mm yr -1 and 367 mm yr -1 , respectively on the other side of the water balance. The gap in the water balance of 9 mm yr -1 is negligible and can be attributed to changes in storages in the soil, snow cover and groundwater. Given the complexity of high mountain hydrology, the scale of the application, the use of one parameter set for the entire basin, and large uncertainties in the meteorological model forcing, we conclude that the model performance is satisfactory for our purpose to estimate the impacts of climate change for the UIB's future hydrology.

Present day hydrology
The stream flow compositions during the reference period have a large spatial variation in the UIB (Fig 5). Strong south to north and east to west gradients are visible in the intensity of the rainfall-runoff generation, consistent with the intensity of the monsoon that comes in from the southeast during the monsoon season. In the monsoon-dominated Sutlej basin the contribution of rainfall to the total flow at the outlet is 74%, whereas this is 33% for the Indus at Tarbela. Snow melt has highest importance in the water coming from the Hindu Kush mountains in the Kabul basin, which receive large amounts of solid precipitation from westerly disturbances during the winter months. Glacier melt contribution is highest in the most glaciated Karakoram subbasins, like Hunza (85%) and Shigar (43%), and the upstream reaches of Kunar. This makes the Indus river the most melt-water dependent river leaving the UIB (55% glacier-and snow melt at Tarbela). The lower-latitude Satluj, Beas, Ravi, Chenab and Jhelum rivers are dominated by input from rainfall, most of which falls during the monsoon season. The Jhelum river also has a substantial snow melt component (32%). Our estimates of stream flow composition match reasonably well with what others found based on a conceptual model [4,5].
The results from this study show a slight shift in runoff composition towards higher contribution of snow melt and rainfall and lower contribution of glacier melt compared to our earlier findings [3]. This is because the current study is focused on the UIB only whereas the earlier study comprised the upstream basins of five Asian rivers, and most importantly because in the current study we use precipitation forcing that is corrected for the underestimate of high altitude precipitation, whereas this was not available at the time of the 2014 study. In the 2014 results, the shortage of precipitation input is compensated by higher glacier melt rates when calibrated only to observed stream flow, a common problem in the simulation of mountain hydrology [76]. Associated to the large differences in the stream flow composition between the tributaries are also differences in the intra-annual distribution of river discharge. Although the peak of glacier melt largely coincides with the peak in monsoonal rains, the snow melt peak occurs during spring. These contrasts in hydrological regimes and stream flow composition of the different tributaries feeding the downstream basin may lead to different responses to future climate change.

Future climate
The downscaled GCM ensembles for RCP4.5 and RCP8.5 show that the future climate in the UIB is highly uncertain. Both ensembles indicate strong warming (Fig 6a and 6b), with significantly stronger warming for the parts of the basin with the highest elevation. The difference in warming can be up to *1°C (RCP4.5) and *2°C (RCP8.5) between the lowest and the highest areas in the UIB. This is well in line with presently observed elevation-dependent warming [104,105]. Comparing the average warming in the UIB (+2.1 to +8.0°C between 1971-2000 and 2071-2100) to the global average (+1.8 to +4.4 between the same periods for the same RCPs [106]), also demonstrates that the UIB is likely to warm stronger than other parts in the world. The uncertainty in warming is largest in the eastern and northern parts of the UIB. Seasonal differences in the temperature projections are limited. For both RCPs, strongest temperature increases are projected for January and June (Fig 6e and 6f). These projected temperature changes are well in line with what was found in other studies for the Indus basin [44,45], although the different scenarios and climate models used in those studies make a direct comparison difficult. The precipitation projections are highly uncertain. The RCP4.5 mean projection shows clear contrasting trends of precipitation increase in the southeastern part and precipitation decrease in the northwestern part of the UIB (Fig 6c). This contrast is also observed for the RCP8.5 mean projection (Fig 6d), although the area with a projected increase in precipitation is larger. Besides, the precipitation increase is much higher for RCP8.5 than for RCP4.5, and the magnitude of precipitation decrease is smaller. The range of the precipitation projections however is very large. For both ensembles, for each geographical location there are both ensemble members that predict an increase and a decrease in precipitation. The ensemble range is much larger for RCP8.5 compared to RCP4.5 and can be up to 100% for the most downstream parts. The mean projected precipitation trend in the southeastern parts of the UIB suggest that monsoon intensity increases, and that the monsoon protrudes further to the northwest with increasing temperatures. Averaged over the UIB, significant seasonal patterns can be observed in the precipitation projections (Fig 6e and 6f). In general, although subject to a large uncertainty, the mean projection in both ensembles is indicating precipitation decrease during February-May and increase in October-January, with the increase being strongest in October. Especially the months October through January have very large uncertainties in RCP8.5. Despite this large uncertainty, an increase in precipitation is likely. Shifts in precipitation patterns originating from westerly disturbances are more difficult to interpret. The increasing precipitation during winter months combined with decreases during the early spring months could suggest that the westerly disturbances set in earlier, however the spatial pattern reveals mostly a precipitation decrease (on annual scale) in those areas where westerly disturbances are the main contributor to total precipitation. The trends, and large uncertainties of precipitation change we find in our ensembles are similar to what was found in an analysis of 32 CMIP5 GCMs over the Hindu-Kush-Karakoram-Himalaya region [42], and once more demonstrate the need for improvement of climate simulations in this region, to lower the uncertainty in the future's climate.

Future glacier extent
The large uncertainty in the climate change scenarios translates in the projected changes in glacier extent (Table 5). Even though the wet scenarios project large increases in precipitation, glacier area decreases considerably during the 21st century throughout the basin, since the precipitation increases cannot compensate for the ample rises in temperature. Our projections are in the same order as projections made in recent other studies at large scale [107,108].

Future hydrology
The uncertainty in UIB's future climate evidently also reflects in the projections of the future hydrology. Nevertheless, several remarkably consistent patterns of projected hydrological changes can be observed across the range of scenarios.
Stream flow composition. The contribution of glacier melt is projected to decrease by the end of the century across all scenarios (Fig 7a, 7d, 7g and 7j). For RCP8.5 the decrease is strongest for the wet, warm scenario and smallest for the dry, cold scenario. The changes in snow melt contribution also show a consistent signal across scenarios, but with high spatial variation (Fig 7b, 7e, 7h and 7k). The strongest decreases are projected for the Hindu Kush mountain range consistent with the high warming rates. In the Karakoram and in the Zanskar subbasin, the contribution of snow melt increases in favor of glacier melt, since the glacier area is reduced but seasonal snow still provides a considerable amount of melt water. Although the strongest precipitation increases are projected for the winter months (Fig 6e and 6f), all year increases in temperature lead to a shift in the precipitation regime to more precipitation falling as rain instead of snow. For the ensemble means, averaged over the UIB the portion of the precipitation falling as rain changes from 58% during 1971-2000 to 66% during 2071-2100 for RCP4.5 and 75% for RCP8.5, consistent with earlier projections of changes in UIB snowfall [109]. Remarkably, despite this shift in precipitation regime, snow melt contribution to total runoff increases or stays equal in most parts of the UIB except for the Kabul basin (Fig 6b, 6e and 6h). This can be explained by the combined effect of increased evapotranspiration due to higher temperatures and increased water availability in the soil and a reduction of sublimation due to decreases in snow cover. For the western part of the Kabul basin the strongest increases in temperature are projected (Fig 6a and 6b), leading to a reduction in snow melt contribution and increase in rainfall-runoff contribution across scenarios. The RCP8.5 wet & warm scenario leads to largest increases in rainfall-runoff contribution (Fig 7l) and for this scenario the contribution of snow melt is mostly reduced (Fig 7k).
Water availability and intra-annual shifts. Changes in stream flow composition are also related to hydrological changes in different times of the year (Fig 8). For most catchments (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13) in the near future (2021-2050) for RCP4.5, flows show little changes during the high flow season and increase during autumn and especially spring. This is most likely due to an increase in autumn and winter precipitation (Fig 6e and 6f) and earlier onset of snow-and glacier melt. Despite the projected decrease in annual precipitation in most of the main Indus branch's basin (Fig 6c), annual-averaged water availability is unchanged for the locations in the Indus river (6)(7)(8)(9)(10)(11)(12), and increases slightly for the upstream subbasins of Hunza, Shigar and Shyok. For Hunza and Shigar this is most likely related to increased glacier melt, and for the Shyok basin upstream of Yogo it is a combination of precipitation increases for the mean of the scenarios (Fig 6c) and increased glacier melt. The lower altitude subbasins (14)(15)(16)(17) show a different pattern of seasonal shifts, with strong decreases in flow during June and July and often also for the spring months. Autumn and winter flows increase slightly, and annual-averaged water availability decreases slightly for these sub-basins. These basins have large rainfall-runoff and snow melt components, and decreases in precipitation during spring and the monsoon season combined with higher evapotranspiration rates, most likely cause runoff to decrease during those months, whereas precipitation increases during the winter months, cause increasing runoff during winter (Fig 6e). For the end of the century (2071-2100), the mean projection for RCP4.5 shows similar changes in intra-annual water distribution as for the near future, but much more pronounced. As glacier areas have reduced significantly by then, the amount of glacier melt water decreases substantially, causing reductions in discharge during the summer months. In addition, flows in the high flow season decline further by reduced precipitation during the monsoon season, and total water availability decreases for the entire UIB due to reduced precipitation in combination with increased evapotranspiration. Flows in spring tend to increase more strongly due to earlier onset of snow and glacier melt during these months. Only for the Satluj river, being the most rain-dominated river in the UIB, increases in water availability are projected for the far future according to the RCP4.5 ensemble mean since precipitation is projected to increase for this part of the UIB (Fig 6c).
In terms of total water availability, the RCP8.5 ensemble shows quite contrasting projections with increases in annual water availability in the near and far future. The patterns in the shifts in the ensemble mean projection however are consistent with RCP4.5, implying a transition to a more attenuated hydrograph. That flows increase during all parts of the year, including the high flow season, is most likely because precipitation is projected to increase for all seasons except spring and glacier melt rates in the near future increase stronger compared to RCP4.5 due to stronger temperature increase. Earlier onset of melt in spring causes runoff to increase during spring despite reduced precipitation input during this season. The projection for the far future shows that despite strong precipitation increases, the glacier-melt dominated Chitral, Hunza, Gilgit and Shigar subbasins experience reductions in flow during the high flow season, since the glacier extent has decreased strongly by then ( Table 5). The similar contrasting shifts between the high-altitude and lower altitude subbasins as for RCP4.5 can be observed. Besides, the contrast in the precipitation projections between the Kabul subbasin and the remaining part of the UIB (Fig 6d) are also visible in the projections of changes in total water availability. The remarkable strong year-round increase in flows in the near future as well as far future for the Shyok subbasin, can most likely be explained by the fact that projected precipitation increases are strongest in this subbasin (Fig 6d). Similarly strong year-round increases in flow  Upper Indus Hydrology and Climate Change for the rain-dominated Satluj river can be explained by strong precipitation increases in this subbasin.
Our results of intra-annual changes are in line with the projections made for the Shigar catchment [38]. There, the initial increase of summer flows is projected halfway through the century followed by a decline at the end of the century, that is accompanied by increasing flows in spring. A previous study projects increasing flow in the UIB until the end of the 21 st century, with more rapid increase during the first half of the century [45]. The authors assessed changes for RCP4.5 and RCP8.5, and used one downscaled GCM and one RCM for their projections. They project stronger increase in winter flows compared to summer flows, consistent with our results. Similar results were found using the previous generation IPCC scenarios A2 and B2 for one RCM [110]. Accurate comparisons to the cited studies is however hampered by the use of different scenarios and climate models.
The patterns are consistent for both RCPs, but the uncertainty is large: for the combined RCP4.5 and RCP8.5 ensembles total water supply from the UIB in 2071-2100 changes by -15% to +60% with respect to 1971-2000. Large uncertainties in hydrological projections have also been found earlier for the Shigar catchment [37,38], and at larger scale [111]. Striking is the particularly large uncertainty observed for the Gilgit subbasin in both RCPs (Fig 8), which is most likely caused by a particularly large uncertainty in monsoon and autumn precipitation and summer air temperatures in both RCPs for this subbasin.
Hydrological extremes. Large changes in extreme discharges can be expected for most parts of the UIB (Fig 9). For most rivers, the highest water levels occur during the coinciding melting and monsoon season, and therefore the changes in return levels are largely determined by the projected climate changes during those months. However, the peak flows are also significantly determined by the meltwater from snow and glacier melt stemming from winter precipitation, forming a basic flow level during the melting and monsoon season which is exacerbated by runoff originating from extreme precipitation events.
Remarkably, the return levels for extreme discharges in the very upstream Hunza river with its highly glaciated basin increase for all scenarios, including the scenarios projecting overall dryer conditions. For the Hunza river with a large contribution of glacier melt, the decreases in glacier extent play a large role, lowering the continuous flow from glacier melt during the melting season. Nevertheless, even for the far future, when the contribution of glacier melt and the total flow has significantly decreased, the extremes in discharge are clearly increasing, due to increases in extreme precipitation, across scenarios. Earlier work in the Shigar subbasin to the east of Hunza also indicates that hydrological extremes may considerably increase until the end of the century [38,49].
The return levels for extreme discharges at Tarbela, where the main Indus branch leaves the UIB, increase for most scenarios as well, except the RCP8.5 dry & warm scenario, because precipitation events are projected to be more intense across the climate model ensemble. At Tarbela, the most extreme changes in return levels are projected for the RCP8.5 wet & warm scenario, with 100-years return level increasing by more than 100% between 1971-2000 and 2071-2100.
For the Shyok basin, return levels clearly increase most for the RCP8.5 wet & warm scenario, which also project the largest precipitation increases. Interestingly, the RCP4.5 wet & warm scenario projects stronger return level increases for the near future compared to the far future, despite increasing precipitation intensity. Since the Shyok river has a large glacier melt contribution, this is related to the lower continuous flow from glacier melt during the melting season in the far future.
At the outlet of the rainfall-runoff dominated Satluj basin the range of projected changes in return levels is largest. As this is the most rainfall-runoff dominated river, the discharge of this river is also most sensitive to changes in extreme precipitation. In addition, precipitation projections have large spread for this part of the UIB (Fig 6c and 6d), which may also imply a large spread in the projected precipitation extremes.
The model runs forced with the RCP8.5 MIROC5 and CSIRO-Mk3 GCMs clearly stand out from the other model runs for Satluj, Shyok, and the Indus at Tarbela, projecting the wettest future and strongest increases in precipitation intensity. Since the uncertainty in future climate is larger for the RCP8.5 ensemble compared to the RCP4.5 ensemble, it is not surprising that the range of the projected return levels is also larger for the RCP8.5 ensemble.

Uncertainty
This study sheds light on the propagation of uncertainty in the future climate for the future hydrology. We emphasize that the future climate in the upper Indus basin is highly uncertain as none of the current state-of-the art GCMs and RCMs satisfactory simulates the monsoon and westerly dynamics in the region [46,86,112], making the reliability of future scenarios questionable. We stress the importance of improvement in the representation of the complex climate in High Mountain Asia in order to be able to narrow down the uncertainty in future projections.
The ensemble of selected climate models determines the outcome of a climate change impact study to a large extent. In this study, climate models to be included in the climate model ensemble have been selected based on the envelope of projections of mean air temperature and precipitation totals. Other selection approaches, based on climate model skill [e.g. 113], changes in multiple climate properties [e.g. 114,115], or a combination of envelope and skill [e.g. 15] will lead to different climate model ensembles and thus to different projections of climate change and its impacts.
There are many different statistical downscaling approaches and choosing the most appropriate method is challenging, especially for areas with complex climate and terrain like the UIB. A study comparing different empirical-statistical downscaling methods for precipitation in the Austrian Alps concludes that methods which apply a non-linear transformation, like ADC does, are among the best performing approaches over terrain with complex topography [116]. A successful application of a method [117], on which ADC is based, has been reported for mountainous parts of the Rhone basin [118]. Although the ADC method was not specifically developed for application in the Hindu Kush-Himalayan region and we could not validate its performance this does provide confidence in its applicability. As in any study projecting future changes in extremes, the uncertainty of the downscaled projections increases with the return period of the events, because these only occupy a very small part of the precipitation intensity distribution. In addition, the ADC method focuses on proper transformation of multi-day precipitation events, assuming these to be the main driver of flooding events. Although we believe that multi-day precipitation events are also the main driver of flooding events in the UIB, future studies should also pay particular attention to changes in multi-day events of high air temperature during the melting season, as these are probably also important drivers of flooding events in the upstream parts of the UIB. In that respect it would be recommendable to test the application of the Quantile Mapping [e.g. 119] downscaling methodology to climate change projections for the UIB, which correct each quantile of the precipitation and air temperature distributions separately and have demonstrated good performance over terrain with complex topography [37,116].
Besides the uncertainty within climate model ensembles, climate models themselves, and empirical-statistical downscaling, additional uncertainties are introduced in the hydrological model forcing and other data, parameters, and representation of physical processes.
Precipitation products show a large range of precipitation amounts for High Mountain Asia [23]. Although we use climate model forcing that is corrected for the underestimate of highaltitude precipitation [52], these data can still have large biases. For example, the UIB-averaged corrected precipitation is estimated to be 913±323 mm yr -1 between 2003 and 2007. Another study assessed precipitation in the UIB based on station data and precipitation estimates for the accumulation zones of major glaciers [64]. The study found similar precipitation totals as reported in the dataset which we used for our study [52], providing further confidence in the used precipitation input. Further narrowing down of the uncertainty in historical precipitation data is a prerequisite for better estimates of future climate change impacts [120]. Similarly, air temperature datasets show a large spread over the UIB. An analysis of the sensitivity of glacier melt water amounts to the large variation in baseline air temperature data shows that variation can lead to an estimated glacier melt contribution to total flow varying from 4 to 78% for the baseline climate, and even larger variations in future projections [121]. Important data used in this study that also introduce uncertainty are the subregion-averaged glacier mass balance data derived from IceSAT data [27], since they are used for the calibration of the large scale glacier change parameterization.
Uncertainties are introduced by using a uniform set of calibrated model parameters for the entire domain. The values of most of the calibrated parameters vary in space and time in reality. However, due to the lack of data which can be used for calibration that covers the entire UIB, we are limited to calibrate the parameters for different subareas and limited temporal periods of our model domain and extrapolate them to larger areas and time periods. This is a common problem in model calibration for data-scarce areas like the UIB [76]. The increasing availability of geodetic glacier mass balance data can help to calibrate spatial variation in degree-day factors for glacier melt, and this can potentially also be done for degree-day factors of snow using MODIS snow cover data. The calibration of parameters which are calibrated using observed discharge are limited to the subbasins which have discharge data available. Different sets of parameters could have been calibrated for each of the gauged subbasins, but then difficulties would arise in assigning parameter values to the ungauged part of the basin. This could in the future potentially be (partly) overcome by categorizing subcatchments by different characteristics, such as climatic differences, degree of glaciation and catchment size. Parameters could then be calibrated for gauged catchments and transferred to ungauged catchments. The potential of such approaches [122] in the UIB region are to be explored in the future.
Parameters themselves have their own uncertainties, which are ideally all taken into account. A study compared three sources of model uncertainty (model parameters, climate projections, natural climate variability) for future projections for a hydrological modeling study in the Hunza subbasin [123]. The study showed that, for heavily glacierized basins, the uncertainty stemming from parameter uncertainty often exceeds the uncertainty stemming from uncertainty in the future climate and natural climate variability. In the cited study an ensemble of three GCMs was used. When a larger ensemble of climate models with a wide range of projections is used, like in our study, the uncertainty stemming from uncertainty in the future climate is probably larger. By calibrating model parameters in a three-step approach using geodetic glacier mass balance data, snow cover data, and observed discharge, the parameter uncertainty is reduced, because the model parameters are constrained to the processes they are affecting, reducing equifinality problems [76].
Uncertainties in the mapping of glacier extent have implications for the simulated contribution of glacier melt to the total flow. For the Karakoram range, the total glacier area according to three different glacier inventories varies from 21193 to 26018 km 2 (i.e. the highest estimate is *23% higher than the lowest estimate) [124], and thus the inventory used will have large consequences for the simulated amounts of glacier melt during the reference period and may also lead to different simulated hydrological responses to climate change. The glacier inventory also determines the initial glacier volume, which is used as a starting point for the simulation of glacier change projections, which is based on volume-area scaling in our approach [74,125]. Besides, the method which is used to estimate initial ice volume from glacier outlines can also vary. A comparison of six different methods to calculate ice volume showed ice volume in the Karakoram ranging from 1683 to 2827 km2 (i.e. the highest estimate being *68% larger than the lowest estimate) [126]. More uncertainties are involved in the glacier change projections as discussed in earlier published work [74]. Currently no basin-wide map with distinction of debris-free and debris-covered glaciers is available for the UIB, and thus the differentiation of both glacier surface types is based on assumptions of elevation and slope constraints controlling the glacier surface type. A map with distinction of both glacier surface types would solve this key issue. Another key issue is the limited understanding of the role of sublimation in the high mountain water balance [99,101].

Conclusions
In this study we use a distributed hydrological model which we force with the latest suite of climate models using an advanced statistical downscaling technique. This study stands out from previous work as for the first time shifts in seasonal water availability are assessed in combination with changes in hydrological extremes at basin scale for the upper Indus basin.
Assessing future hydrological changes in the upper Indus basin is complicated by large uncertainties in the historical and future climate, uncertainties in glacier extent and glacier mass balance, and uncertainties in hydrological model processes and parameters.
From our results we can conclude that the upper Indus basin faces a very uncertain future in terms of water availability in the long run. Projections of changes in water availability from the upper Indus basin at the end of the 21 st century range from -15% to +60% with respect to 1971-2000. This uncertainty mainly stems from the large spread in the projections of precipitation change throughout the 21 st century. Therefore, formulating adequate adaptation measures which take into account the uncertain future is of vital importance, thus requiring hydrological projections to be made based on an ensemble of climate models representing all possible futures.
Despite the large uncertainties in future climate and water availability, basin-wide patterns and trends of intra-annual shifts in water availability are consistent across climate change scenarios. These trends mainly consist of minor increases in summer flows combined with increased flows during other seasons in the near future (2021-2050) and decreases in summer flows combined with stronger increasing flows during the other seasons in the far future (2071-2100). Furthermore, increases in intensity and frequency of extreme discharges are found for most of the UIB and for most scenarios and models considered, implying increases in flooding events during the 21 st century.
Population growth in combination with increasing standards of living and associated increases in energy and food production will continue to expand the downstream water and energy demand [127,128]. This implies a growing dependency on the uncertain future water resources, which calls for sound basin-wide adaptation strategies to be developed across sectors that take into account the changing demand and supply in the Indus basin as well as uncertainties therein.
financial support of core research at ICIMOD. This work is partly carried out by the Himalayan Adaptation, Water and Resilience (HI-AWARE) consortium under the Collaborative Adaptation Research Initiative in Africa and Asia (CARIAA) with financial support from the United Kingdom's Department for International Development (DFID) and the International Development Research Centre (IDRC), Ottawa, Canada. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement 676819). In addition, this work was partially supported by core funds of ICIMOD contributed by the governments of Afghanistan, Australia, Austria, Bangladesh, Bhutan, China, India, Myanmar, Nepal, Norway, Pakistan, Switzerland, and the United Kingdom. DFID partly funds ICIMOD's Indus Basin Programme, under which this study was undertaken. These funds partly supported the contributions of AFL, WWI and ABS. DFID and IDRC funds the HI-AWARE consortium, of which ICIMOD and FutureWater are consortium members. These funds partly supported the contributions of AFL, WWI and ABS.
ERC partly funded the contribution of WWI through grant agreement 676819. ICIMOD core funds partly funded the contribution of ABS. AFL and WWI are employed by Future-Water. FutureWater provided support in the form of salaries for authors AFL and WWI, but did not have any additional role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. The specific roles of these authors are articulated in the 'author contributions' section.
We have the following interests: AFL and WWI are employed by FutureWater. There are no patents, products in development or marketed products to declare. This does not alter our adherence to all the PLOS ONE policies on sharing data and materials. We are grateful to T. Bolch for providing the geodetic mass balance data and S. Bajracharya for providing detailed glacier outlines for the Hunza basin. We acknowledge the World Climate Research Program's Working Group on Coupled Modeling, which is responsible for CMIP5, and we thank the climate modelling groups for producing and making available their model output. We thank the Pakistan Water and Power Development Authority and the Pakistan Meteorological Department for making available discharge and meteorological station data. We thank two anonymous reviewers for their constructive comments that helped to improve the manuscript significantly.