Reconciling high-altitude precipitation in the upper Indus basin with glacier mass balances and runoff

Mountain ranges in Asia are important water suppliers, especially if downstream climates are arid, water demands are high and glaciers are abundant. In such basins, the hydrological cycle depends heavily on high-altitude precipitation. Yet direct observations of high-altitude precipitation are lacking and satellite derived products are of insufficient resolution and quality to capture spatial variation and magnitude of mountain precipitation. Here we use glacier mass balances to inversely infer the high-altitude precipitation in the upper Indus basin and show that the amount of precipitation required to sustain the observed mass balances of large glacier systems is far beyond what is observed at valley stations or estimated by gridded precipitation products. An independent validation with observed river flow confirms that the water balance can indeed only be closed when the highaltitude precipitation on average is more than twice as high and in extreme cases up to a factor of 10 higher than previously thought. We conclude that these findings alter the present understanding of high-altitude hydrology and will have an important bearing on climate change impact studies, planning and design of hydropower plants and irrigation reservoirs as well as the regional geopolitical situation in general.


Introduction
Of all Asian basins that find their headwaters in the greater Himalayas, the Indus basin depends most strongly on highaltitude water resources (Immerzeel et al., 2010;Lutz et al., 2014) The largest glacier systems outside the polar regions are found in this area and the seasonal snow cover is the most extensive of all Asian basins (Immerzeel et al., 2009).In combination with a semi-arid downstream climate, a high demand for water as a result of the largest irrigation scheme in the world and a large and quickly growing population, the importance of the upper Indus basin (UIB) is evident (Immerzeel and Bierkens, 2012).
The hydrology of the UIB (4.37 × 10 5 km 2 ) is, however, poorly understood.The quantification of the water balance in space and time is a major challenge due the lack of measurements and the inaccessibility of the terrain.The magnitude and distribution of high-altitude precipitation, which is the driver of the hydrological cycle, is one of its largest unknowns (Hewitt, 2005(Hewitt, , 2007;;Immerzeel et al., 2013;Mishra, 2015;Ragettli and Pellicciotti, 2012;Winiger et al., 2005).Annual precipitation patterns in the UIB result from the intricate interplay between synoptic scale circulation and valley scale topography-atmosphere interaction resulting in orographic precipitation and funnelling of air movement (Barros et al., 2004;Hewitt, 2013).At the synoptic scale, annual precipitation originates from two sources: the south-eastern monsoon during the summer and moisture transported by the westerly jet stream over central Asia (Mölg et al., 2013;Scherler et al., 2011) during winter.The relative contribution of westerly disturbances to the total annual precipitation increases from south-east to north-west, and the anomalous behaviour of Karakoram glaciers is commonly attributed to changes in winter precipitation (Scherler et al., 2011;Yao et al., 2012).(Lehner et al., 2008), basin hypsometry and three gridded precipitation products.(a) shows the digital elevation model, the location of the major glacier systems (area > 5 km 2 ), the available stations in the Global Summary of the Day (GSOD) of the World Meteorological Organization (WMO) and the hydrological stations used for validation.Panel B shows box plots of the elevation distribution of the basin, the large glacier systems, the GSOD meteorological stations and the average elevation of the catchment area of each hydrological station.(c) to (e) show the average gridded annual precipitation between 1998 and 2007 for the APHRODITE (Yatagai et al., 2012), TRMM (Huffman et al., 2007) and ERA-Interim (Dee et al., 2011) data sets.
At smaller scales the complex interaction between the valley topography and the atmosphere dictates the spatial distribution of precipitation (Bookhagen and Burbank, 2006;Immerzeel et al., 2014b).Valley bottoms, where stations are located, are generally dry and precipitation increases up to a certain maximum altitude (HMAX) above which all moisture has been orographically forced out of the air and precipitation decreases again.In westerly dominated rainfall regimes HMAX is generally higher, which is likely related to the higher tropospheric altitude of the westerly airflow (Harper, 2003;Hewitt, 2005Hewitt, , 2007;;Scherler et al., 2011;Winiger et al., 2005).
Gridded precipitation products are the de facto standard in hydrological assessments, and they are either based on observations (e.g.APHRODITES; Yatagai et al., 2012), remote sensing (e.g. the Tropical Rainfall Monitoring Mission; Huffman et al., 2007) or re-analysis (e.g.ERA-Interim; Dee et al., 2011) (Fig. 1c-e).In most cases the station data strongly influence the distribution and magnitude of the precipitation in those data products; however, the vast majority of the UIB is located at elevations far beyond the average station elevation (Fig. 1a-b).The few stations that are at elevations above 2000 m are located in dry valleys and we hypothesise that the high-altitude precipitation is considerably underestimated (Fig. 1c-e).Moreover, remotesensing-based products, such as the Tropical Rainfall Measuring Mission (TRMM), are insufficiently capable of capturing snowfall (Bookhagen and Burbank, 2006;Huffman et al., 2007) and the spatial resolution (25-75 km 2 ) of most rainfall products (and the underlying models) is insufficient to capture topography-atmosphere interaction at the valley scale (Fig. 1c-e).Thus, there is a pressing need to improve the quantification of high-altitude precipitation, preferably at large spatial extents and at high resolution.
A possible way to correct mountain precipitation is to inversely close the water balance.Previous studies in Sweden and Switzerland have shown that it is possible to derive vertical precipitation gradients using observed runoff in a physically realistic manner (Valéry et al., 2009(Valéry et al., , 2010)).Earlier work at the small scale in high mountain Asia suggested that the glacier mass balance may be used to reconstruct precipitation in its catchment area (Harper, 2003;Immerzeel et al., 2012a).Figure 1a and b show that UIB glaciers are located at high elevations that are not represented by station data.Therefore, the mass balances of the glaciers may contain important information on high-altitude accumulation in an area that is inaccessible and ungauged, but very important from a hydrological point of view.In this study we further elaborate this approach by inversely modelling average annual precipitation from the mass balance of 550 large (> 5 km 2 ) glacier systems located throughout the UIB.We perform a rigorous uncertainty analysis and we validate our findings using independent observations of river runoff.

Methods
We estimate high-altitude precipitation by using a glacier mass balance model to simulate observed glacier mass balances.We use a gridded data set from valley bottom stations as a basis for our precipitation estimate and we compute a vertical precipitation gradient (PG; % m −1 ) until observed mass balances match the simulated mass balance.We repeat this process for the 550 major glacier systems in the UIB, and the resulting PGs are then spatially interpolated to generate a spatial field that represents the altitude dependence of precipitation.We use this field to update the APHRODITE precipitation and generate a corrected precipitation field that is able to reproduce the observed glacier mass balance.We validate the findings independently with a water balance approach.Estimated (annual) runoff, based on the corrected precipitation, actual evapotranspiration based on four gridded products and the observed glacier mass balance, is compared with an extensive set of UIB runoff observations.We also analyse the physical realism of our simulations by deriving a Turc-Budyko plot using precipitation, measured runoff and potential evapotranspiration.A rigorous uncertainty analysis is also conducted on the six most critical model parameters including potential effects of spatial correlation.

Glacier mass balance and outlines
Glacier mass balance trends based on NASA's Ice, Cloud and land Elevation Satellite (ICESat) (Kääb et al., 2012a) are recomputed for the period 2003 until 2008 for the three major mountain ranges in the UIB: the Karakoram, the Hindu Kush and the Himalaya (Fig. 1).For each zone the mass balance is computed including a regional uncertainty estimate (Kääb et al., 2012a).From the zonal uncertainty (σ z ) we estimate the standard deviation between glaciers within a zone (σ g ) as where n is the number of glaciers within a zone.The σ g values used in the uncertainty analysis are shown in Table 1.
The glacier boundaries are based on the glacier inventory of the International Centre for Integrated Mountain Development (Bajracharya and Shrestha, 2011).

Precipitation and temperature
The daily APHRODITE precipitation (Yatagai et al., 2012) and air temperature data sets (Yasutomi et al., 2011) from 2003 until 2007 are used as reference data sets to ensure maximum temporal overlap with the ICESat-based glacier mass balance data set (Kääb et al., 2012a).The precipitation data set is resampled from the nominal resolution of 25 km 2 to a resolution of 1 km 2 using the nearest neighbour algorithm.The air temperature data set is then bias corrected using monthly linear regressions with independent station data to account for altitudinal and seasonal variations in air temperature lapse rates (Fig. 3).

Runoff and evapotranspiration
We use runoff data, potential evapotranspiration (ETp) and actual evapotranspiration (ETa) data for the validation of our results.For runoff we compiled all available published data, which we complemented with data made available by the Pakistan Meteorological Department and the Pakistan Water and Power Development Authority.
Evapotranspiration is notoriously difficult to monitor and there are few direct measurements of ETa in the upper Indus.In earlier UIB studies, ET was estimated using empirical formulae based on air temperature but was only applied to the Siachen glacier (Bhutiyani, 1999;Reggiani and Rientjes, 2014).We take into account the uncertainty in ET in our streamflow estimates and develop a blended product based on re-analysis data sets, a global hydrological model and an energy balance model.Four gridded ETa and three gridded ETp products were resampled to a 1 km 2 resolution at which we perform all our analyses: -ERA-Interim re-analysis (Dee et al., 2011): ERA-Interim uses the HTESSEL land-surface scheme (Dee et al., 2011) to compute ETa.For transpiration a distinction is made between high and low vegetation in the HTESSEL scheme and these are parameterised from the Global Land Cover Characteristics database at a nominal resolution of 1 km 2 .
-Modern-Era Retrospective Analysis for Research and Applications (MERRA) re-analysis (Rienecker et al., 2011): the MERRA re-analysis product of NASA applies the state-of-the-art GEOS-5 data assimilation system that includes many modern observation systems in a climate framework.MERRA uses the GEOS-5 catchment land surface model (Koster et al., 2000) to compute actual ET.For the MERRA product ETp is not available.

Model description
We use the PCRaster spatial-temporal modelling environment (Karssenberg et al., 2001) to model the mass balance of the major glaciers in each zone and subsequently estimate precipitation gradients required to sustain the observed mass balance.The model operates at a daily time step from 2003 to 2007 and a spatial resolution of 1 km 2 .For each time step the total accumulation and total melt are aggregated over the entire glacier surface.Only glaciers with a surface area above 5 km 2 are included in the analysis (Karakoram is 232 glaciers, Hindu Kush is 119, Himalaya is 204 glaciers), as the ICESat measurements do not reflect smaller glaciers.

Nov
Elevation (m asl) ∆T(OBS − APHRO) ( ° C) q q q q q q q q q q q q q q 1000 3000 5000 −4 −2 0 the precipitation is negatively lapsed from its maximum at HMAX with the same PG according to for HREF < H ≤ HMAX, and for H > HMAX.
HREF and HMAX values are derived from literature (Table 1) and uncertainty is taken into account in the uncertainty analysis.HMAX varies per zone and lies at a lower elevation in the Himalayas than in the other two zones (Table 1).We spatially interpolate HMAX from the average zonal values to cover the entire UIB.The melt is modelled over the glacier area using the positive degree day method (Hock, 2005), with different degree day factors (DDFs) for debris-covered (DDFd) and debrisfree (DDFdf) glaciers derived from literature (Table 1).To account for uncertainty in DDFs, the DDFd and DDFdf are taken into account separately in the uncertainty analysis.
At temperatures below the critical temperature of 2 • C ( Immerzeel et al., 2013;Singh and Bengtsson, 2004), precipitation falls in the form of snow and contributes to the accumulation.Avalanche nourishment of glaciers is a key contributor for UIB glaciers (Hewitt, 2005(Hewitt, , 2011) ) and to take this process into account, we extend the glacier area with steep areas directly adjacent to the glacier with a slope over an average threshold slope (TS) of 0.2 m m −1 .This average threshold slope is derived by analysing the slopes of all glacier pixels in the basin (Fig. 4).To account for uncertainty in TS, this parameter is taken into account in the uncertainty analysis.
For each glacier system, we execute transient model runs from 2003 to 2007 and we compute the average annual mass balance from the total accumulation and melt over this period.We make two different runs for each glacier system with two different PGs (0.3 and 0.6 % m −1 ) and we use the simulated mass balances of these two runs and the observed mass balances based on ICESat to optimise the PG per glacier, such that the simulated mass balance matches the observed.
To interpolate the glacier-specific PG values to PG spatial fields over the entire domain we use geostatistical conditional simulation (Goovaerts, 1997) Reconciling high-altitude precipitation in the upper Indus basin q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 2500 3000 3500 4000 4500 5000 5500 6000 6500 7000 centroid.The semi-variogram has the following parameters: nugget = 0, the range = 120 km, sill is the variance of PGs.

Uncertainty analysis
A rigorous uncertainty analysis is performed to take into account the uncertainty in parameter values and uncertainty in regional patterns.To account for parameter uncertainty, we perform a 10 000 member Monte Carlo simulation on the parameters given in Table 1.For each run we randomly sample the parameter space based on the average (µ) and the standard deviation (σ ), which are all based on literature values.For the positively valued parameters, we use a log-Gaussian distribution and a Gaussian distribution in case parameter values can be negative.We take into account uncertainty in the following key parameters (HREF, HMAX, DDFd, DDFdf, TS) for the PG as well as uncertainty in the mass balance against which the PG is optimised (mass balance, MB).We randomly vary the five parameters (HREF, HMAX, DDFd, DDFdf, TS) 10 000 times and calculate the PG for each glacier for each random parameter set drawn, thus resulting in 10 000 PG sets for each glacier considered.
For each of the 10 000 PG sets, we then use conditional simulation (see above) to arrive at 10 000 equally probable spatial PG fields, taking account of parameter's uncertainty, mass-balance uncertainty and the interpolation error.Note that for each of the 10 000 sets, the variogram is scaled with the variance of the PGs associated with the specific parameter combination drawn.Finally, based on the results of the 10 000 simulations we derive the average-corrected precipitation field and the associated uncertainty in the estimates Using the 10 000 combinations of parameters and associated PGs, we ran a multi-variate linear regression analy-sis to determine relative contribution of each parameter to the spread in the PG to understand which parameter has the largest influence on the PG.
It is possible that certain parameters used in the model are spatially correlated.To account for uncertainty in this spatial correlation and the presence of spatial patterns in the parameters, we perform a sensitivity analysis where we consider three cases: -Fully correlated: we assume the parameters are spatially fully correlated within a zone, e.g. for each of the 10 000 simulations a parameter has the same value within a zone.
-Uncorrelated: we assume the parameters are spatially uncorrelated and within each zone each glacier system is assigned a random value.
-Intermediate case: we use geostatistical unconditional simulation (Goovaerts, 1997) with a standardised semivariogram (nugget = 0, sill is the variance of parameter, range = 120 km) to simulate parameter values for each glacier system.

Validation
We estimate the average annual runoff (Q) for sub-basins in the UIB from where P cor is the average corrected precipitation, ET the average annual evapotranspiration based on the four products described above and MB the glacier mass balance expressed over the sub-basin area in mm yr −1 .We then compare the estimated runoff values to the observed time series (Table 2).
For the three zones (Himalaya, Karakoram and Hindu Kush) we also perform a water balance analysis to verify whether the use of the corrected precipitation product results in a more realistic closure of the water balance.Finally, we test the physical realism of the corrected precipitation product using a non-dimensional Turc-Budyko plot as described in Valéry et al. (2010).This plot is based on two assumptions: (i) the mean annual runoff should not exceed the mean annual precipitation and (ii) the mean annual runoff should be larger than or equal to the difference between precipitation and potential evapotranspiration.By plotting P / ETp versus Q / P on a catchment basis, it is tested whether the use of corrected precipitation results in more physically realistic values.

Corrected precipitation
The average annual precipitation based on 10 000 conditionally simulated fields reveals a striking pattern of high-altitude precipitation.The amount of precipitation required to sustain the large glacier systems is much higher than either the station observations or the gridded precipitation products imply.For the entire UIB the uncorrected average annual precipitation (Yatagai et al., 2012) for 2003-2007 is 437 mm yr −1 (191 km 3 yr −1 ), an underestimation of more than 200 % compared with our corrected precipitation estimate of 913 ± 323 mm yr −1 (399 ± 141 km 3 yr −1 ; Fig. 5).
The greatest corrected annual precipitation totals in the UIB (1271 mm yr −1 ) are observed in the elevation belt between 3750 to 4250 m (compared to 403 mm yr −1 for the uncorrected case).In absolute terms the main water-producing region is located in the elevation belt between 4250 and 4750 m where approximately 78 km 3 of rain and snow precipitates annually.
In the most extreme case, precipitation is underestimated by a factor 5 to 10 in the region where the Pamir, Karakorum and Hindu Kush ranges intersect (Fig. 5).Our inverse modelling shows that the large glacier systems in the region can only be sustained if snowfall in their accumulation areas totals around 2000 mm yr −1 (Hewitt, 2007).This is in sharp contrast to precipitation amounts between 200 and 300 mm yr −1 that are reported by the gridded precipitation products (Fig. 1).Our results match well with the few studies on high-altitude precipitation that are available.Annual accumulation values between 1000 and 3000 mm have been reported for accumulation pits above 4000 m in the Karakoram (The Batura Glacier Investigation Group, 1979;Wake, 1989;Winiger et al., 2005).Our results show that the highest precipitation amounts are found along the monsoon-influenced southern Himalayan arc with values up to 3000 mm yr −1 , while north of the Himalayan range the precipitation decreases quickly towards a vast dry area in the north-eastern part of the UIB (Shyok sub-basin).In the north-western part of the UIB, westerly storm systems are expected to generate considerable amounts of precipitation at high altitude.
Our results reveal a strong relation between elevation and precipitation with a median PG for the entire upper Indus of 0.0989 % m −1 , but with large regional differences.Median precipitation gradients in the Hindu Kush and Karakoram ranges (0.260 and 0.119 % m −1 , respectively) are significantly larger than those observed in the Himalayan range, e.g 0.044 % m −1 (Fig. 6).In the Hindu Kush, for example, for every 1000 m elevation rise, precipitation increases by 260 % with respect to APHRODITE, which is based on valley floor precipitation.Higher HMAX in the Hindu Kush and the Karakoram (e.g.5500 m versus 4500 m in the Himalayas; Hewitt, 2007;Immerzeel et al., 2014a;Putkonen, 2004;Seko, 1987;Winiger et al., 2005) suggests that westerly airflow indeed has a higher tropospheric altitude and that the interplay between elevation and precipitation is stronger for this type of precipitation.Further research should thus focus on the use of high-resolution cloud-resolving weather models (Collier et al., 2014;Mölg et al., 2013) for this region to further resolve seasonal topography-precipitation interaction at both synoptic and valley scales.The estimated precipitation is considerably higher than what was reported in previous studies.Several studies have used TRMM products to quantify UIB precipitation (Bookhagen and Burbank, 2010;Immerzeel et al., 2009Immerzeel et al., , 2010) ) and they show average annual precipitation values around 300 mm.It was also noted that the water balance was not closing and average annual river runoff at Tarbela exceeded the TRMM precipitation (Immerzeel et al., 2009).Two possible reasons have been suggested to explain this gap: (i) the high-altitude precipitation is underestimated and (ii) the glaciers are in a significant negative imbalance (Immerzeel et al., 2009).Since the ICESat study and several other geodetic mass balance studies (Gardelle et al., 2013;Kääb et al., 2012b) it has become clear that the glaciers in this region are not experiencing a significant ice loss and that this cannot be the explanation for the missing water in the water balance.This supports our conclusion that it is the high-altitude precipitation that has been underestimated.A study based on long-term observations of Tarbela inflow also confirm our results (Reggiani and Rientjes, 2014).In this study the total UIB precipitation above Tarbela is estimated to be 675 ± 100 mm yr −1 and the difference remaining between our results may stem from the fact that the UIB we consider is twice the size of the area above the Tarbela, the different approach used to estimate actual ET, the different period considered and their assumption that ice storage has not changed between 1961 and 2009.

Uncertainty
We estimated the uncertainty in the modelled precipitation field with the standard deviation (σ ) of the 10 000 realisations (Fig. 5c).The signal-to-noise ratio is satisfactory in the entire domain, e.g. the σ is always considerably smaller than the average precipitation with an average coefficient of variation of 0.35.The largest absolute uncertainty is found along the Himalayan arc and this coincides with the precipitation pattern found here.Strikingly, the region where the underestimation of precipitation is largest, at the intersection of the three mountain ranges in the northern UIB, is also an area where the uncertainty is small even though precipitation gradients are large.By running a multiple regression analysis after optimising the PGs, we quantify the contribution of each parameter to the total uncertainty.The largest source of uncertainty in our estimate of UIB high-altitude precipitation stems from the MB estimates, followed by the DDFdf, DDFd, TS, HREF and HMAX, although regional differences are considerable (Fig. 7).The MB constrains the precipitation gradients and thereby exerts a strong control over the corrected precipitation fields, in particular because the intra-zonal variation in MB is relatively large (Table 1).Thus, improved spatial monitoring techniques of the MB are expected to greatly improve precipitation estimates.
Figure 8 shows the result of uncertainty analysis associated with the spatial correlation of the parameters, which reveals that the impact on the average-corrected precipitation is limited.Locally there are minor differences in the corrected precipitation amounts, but overall the magnitude and spatial patterns are similar.However, there are considerable differences in the uncertainty.The lowest uncertainty is found for the fully uncorrelated case, the fully correlated case has the highest uncertainty whereas the intermediate case is in between both.For the fully correlated case all glacier systems have the same parameter set for the specific realisation and this results in a larger final uncertainty.In the uncorrelated case each glacier system has a different randomly sampled parameter set and this reduces the overall uncertainty as it spatially attenuates the variation in precipitation gradients.

Validation
The corrected precipitation is validated independently by a comparison to published average annual runoff data of 27 stations (Table 2).Runoff estimates based on the corrected precipitation agree well with the average observed annual runoff (Fig. 9, top panel).It is interesting to note that the higher catchments (r = 0.98, red outline) show a better cor- relation with the observed runoff than the lower catchments (r = 0.76, black outline).The runoff estimated for the uncorrected APHRODITE is consistently lower than the observed runoff, and in some occasions even negative.Runoff estimates were also made based on the ERA-Interim and TRMM precipitation products.The TRMM results yield a similar underestimation as the uncorrected APHRODITES product, but the runoff estimates based on the ERA-Interim precipitation agrees reasonably well with the observations.However, the coarse resolution (∼ 70 km 2 ) (Fig. 1) is problematic and cannot be used to reproduce the mass balance (Fig. 11).Averaged over large catchments the precipitation may be applied for hydrological modelling, but at smaller scales there are likely very large biases.As a result, the observed glacier mass balances cannot be reproduced when the ERA-Interim data set is used.Although the ERA-Interim data set may not be used to reproduce the glacier mass balances, it can be used to verify the atmospheric convergence as the product is based on a data assimilation scheme and the ECMWF IFS forecast model that includes fully coupled components for atmosphere, land surface and ocean waves, including closure of the atmospheric water balance.The total precipitation sum from 1998 to 2007 of the ERA-Interim data set over the entire UIB is 947 mm, which is very close to our corrected precipitation sum of 913 mm.This indicates that the westerlies and monsoon circulation transport sufficient moisture into the region to account for the precipitation we estimate.The source of precipitation in the upper Indus is a mixture of the Arabian Sea (westerlies), Bay of Bengal (Monsoon) and potentially also intra-basin moisture recycling (Tuinenburg et al., 2012;  Wei et al., 2013); however, further research with atmospheric models is required to quantify these contributions further.The zonal water balance analysis (Fig. 9, bottom panels) reveals that the water balance is much more realistic when the corrected precipitation is used.Although the uncertainties are considerable, our analysis shows that the Himalaya and Hindu Kush zones are about twice as wet as the Karakoram zone.For all three zones the glacier mass imbalance only plays a marginal role in the overall water balance and about 60 % of the total precipitation runs off while 40% is lost through evapotranspiration.Notable the values for Corg, which represents the water balance gap in case the uncorrected precipitation is used, are approximately 500 mm yr −1 in all three zones.Our validation does not take into account groundwater fluxes and we have assumed that over the observed period from 2003 to 2007 there is no net loss or gain of groundwater in the upper Indus basin.We do acknowledge that groundwater may play an important role in the hydrology.A study in the Nepal Himalaya shows that fractured basement aquifers fill during the monsoon and they purge in the post-monsoon thus causing a natural delay in runoff (Andermann et al., 2012).However, this does not imply sig-nificant net gains or losses over multiple year periods, which is what we consider.A second component that we have not considered and that may play a role in the high-altitude water cycle is sublimation.There are some indications that wind redistribution and sublimation may play a considerable role in the high-altitude water balance (Wagnon et al., 2013).However, our PGs are constrained on the observed mass balance; hence, our precipitation can be considered as a net precipitation and sublimation losses are thus accounted for.
In Fig. 10 the Turc-Budyko plot is shown to confirm the physical realism of our results.Those dots located in the hatched part of the graph are physically less realistic.For the uncorrected case almost all dots (open dots) are above the Q / P = 1 line.For the corrected case the Q / P values are much more plausible; however, there many catchments that are located slightly to the right side of the theoretical Budyko line, meaning that the Q is smaller than the difference between P and ETp.Possible deviation can potentially be explained by uncertainties in observed flows and ETp estimates; the fact that in glacierised catchments the theoretical Budyko curve may be different because of a glacier imbalance can be an additional water balance term that is unac- q q q q q q q q q q q q q Himalaya Flux (mm/y) The box plots with a red outline have an average elevation higher than 4000 m. and the box plots with a black outline have an elevation lower than 4000 m.Bottom panels: water balance components of each zone (P cor is corrected precipitation, P org is uncorrected APHRODITES precipitation, ET is actual evapotranspiration, Mass is glacier mass balance, Q cor is estimated runoff and C org is water balance gap in case the P org is used).
counted for, a too short time period that is used to construct the water balance and, finally, that some of the discharge observations do not align in time with the rest of the water balance terms.Overall we conclude though that the use of the corrected precipitation results in physically more realistic results, where the water balance could be closed and no significant amount of precipitation input is missing.
Figure 11 shows how the average simulated zonal glacier mass balance using the corrected, the APHRODITES, the ERA-Interim and the TRMM precipitation data sets.It shows that none of the precipitation products can reproduce the observed mass balance.Mostly the mass balances are underestimated, which is consistent with the underestimation of the precipitation.The ERA-Interim data set overestimates the mass balance in the Himalaya, but underestimates the mass balances in the other two zones as a result of the coarse resolution.

Conclusions
In this study we inversely model high-altitude precipitation in the upper Indus basins from glacier mass balance trends derived by remote sensing.Although there are significant uncertainties, our results unambiguously show that highaltitude precipitation in this region is underestimated and that the large glaciers here can only be sustained if high-altitude accumulation is much higher than most commonly used gridded data products.
Our results have an important bearing on water resources management studies in the region.The observed gap between precipitation and streamflow (Immerzeel et al., 2009) (with q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Aphrodite Corrected Dif period Partial overlap Same period Figure 10.Non-dimensional graphical representation of catchments using their mean runoff, Q, precipitation, P , and potential evapotranspiration, PE.The grey line in the empty centre area represents the theoretical Budyko relationship in the non-dimensional graph.The size of the dots is scaled to the catchment area.streamflow being larger) cannot be attributed to the observed glacier mass balance (Kääb et al., 2012a), but is most likely the result of an underestimation of precipitation, as also follows from this study.With no apparent decreasing trends in precipitation (Archer and Fowler, 2004), the observed negative trends in streamflow in the glacierised parts of the UIB should thus be primarily attributed to decreased glacier and snowmelt (Sharif et al., 2013) and increased glacier storage (Gardelle et al., 2012).In a recent study the notion of negative trends in UIB runoff was contested and based on a recent analysis  it was concluded that runoff of Karakoram rivers is increasing (Mukhopadhyay and Khan, 2014b).The study suggests that increase glacier melt during summer is the underlying reason, which in combination with positive precipitation trends in summer does not contradict the neutral glacier mass balances in the region.From all of these studies it becomes apparent that precipitation is the key to understanding behaviour of glacier and hydrology at large in the UIB.The precipitation we estimate in this study differs considerably, in magnitude and spatial distribution, from data sets that are commonly used in design of reservoirs for hydropower and irrigation and as such it may have a significant impact and improve such planning processes.
The water resources of the Indus River play an important geopolitical role in the region, and our results could contribute to the provision of independent estimates of UIB precipitation.We conclude that the water resources in the UIB are even more important and abundant than previously thought.Most precipitation at high altitude is now stored in the glaciers, but when global warming persists and the runoff regime becomes more rain dominated, the downstream impacts of our findings will become more evident.

Figure 1 .
Figure1.Overview of the UIB(Lehner et al., 2008), basin hypsometry and three gridded precipitation products.(a) shows the digital elevation model, the location of the major glacier systems (area > 5 km 2 ), the available stations in the Global Summary of the Day (GSOD) of the World Meteorological Organization (WMO) and the hydrological stations used for validation.Panel B shows box plots of the elevation distribution of the basin, the large glacier systems, the GSOD meteorological stations and the average elevation of the catchment area of each hydrological station.(c) to (e) show the average gridded annual precipitation between 1998 and 2007 for the APHRODITE(Yatagai et al., 2012), TRMM(Huffman et al., 2007) and ERA-Interim(Dee et al., 2011) data sets.

Figure 4 .
Figure 4. Box plots of slopes of glacierised areas per elevation bin.

Figure 5 .
Figure 5. Corrected precipitation and estimated uncertainty for the UIB for the case with intermediate spatial correlation between model parameters.(a) shows the average modelled precipitation field based on 10 000 simulations for the period 2003-2007, (b) shows the ratio of corrected precipitation to the uncorrected APHRODITE precipitation for the same period, (c) shows the standard deviation of the 10 000 simulations and (d) shows the average precipitation gradient.

Figure 6 .
Figure 6.Box plots of precipitation gradients for the entire UIB and for the three regions separately.

Figure 7 .
Figure7.Normalised weights of multiple regression of the precipitation gradients the predictors slope (slope threshold for avalanching to contribute to accumulation), HREF (base elevation from which lapsing starts), HMAX (elevation with peak precipitation), DDFd (degree day factor for debris-covered glaciers), DDFdf (degree day factor for debris-free glaciers) and the MB (mass balance of the glacier).

Figure 8 .
Figure 8. Impact of spatial correlation of parameters on the corrected precipitation field and associated uncertainty.The top panels show the corrected precipitation field (a) and uncertainty (b) for the fully uncorrelated case.The middle panels (d, e) for the fully correlated case and the bottom panels (e, f) for the intermediate case.

Figure 9 .
Figure 9. Validation of the precipitation correction using observed discharge (Table2).Top panel: the box plots are based on the runoff estimate based on 10 000 corrected precipitation fields (grey: stations for which the observed record does not coincide with the 2003-2007 period; yellow: stations for which the 2003-2007 period is part of the observational record; green: stations for which the observations are based precisely on the 2003-2007 period).The black dots and red diamonds (estimated runoff below 50 m 3 s −1 ) show the estimated runoff based on the uncorrected precipitation.The box plots with a red outline have an average elevation higher than 4000 m. and the box plots with a black outline have an elevation lower than 4000 m.Bottom panels: water balance components of each zone (P cor is corrected precipitation, P org is uncorrected APHRODITES precipitation, ET is actual evapotranspiration, Mass is glacier mass balance, Q cor is estimated runoff and C org is water balance gap in case the P org is used).

Figure 11 .
Figure 11.Reconstructed mass balances based on the corrected, APHRODITE, ERA-Interim and TRMM data sets.The black horizontal dotted line shows the observed mass balance for each zone.

Table 2 .
Khattak et al. (2011)for validation.Catchment areas are delineated based on SRTM DEM (Shuttle Radar Topographic Missiondigital elevation model).a is calculated based on discharge provided by the Pakistan Water and Power Development (WAPDA), b is based on Mukhopadhyay and Khan (2014a), c is based on Sharif et al. (2013), d is based on Archer (2003) and e is based onKhattak et al. (2011).