Glaciers as a Proxy to Quantify the Spatial Distribution of Precipitation in the Hunza Basin

Abstract Accurate quantification of the spatial distribution of precipitation in mountain regions is crucial for assessments of water resources and for the understanding of high-altitude hydrology, yet it is one of the largest unknowns due to the lack of high-altitude observations. The Hunza basin in Pakistan contains very large glacier systems, which, given the melt, cannot persist unless precipitation (snow input) is much higher than what is observed at the meteorological stations, mostly located in mountain valleys. Several studies, therefore, suggest strong positive vertical precipitation lapse rates; in the present study, we quantify this lapse rate by using glaciers as a proxy. We assume a neutral mass balance for the glaciers for the period from 2001 to 2003, and we inversely model the precipitation lapse by balancing the total accumulation in the catchment area and the ablation over the glacier area for the 50 largest glacier systems in the Hunza basin in the Karakoram. Our results reveal a vertical precipitation lapse rate that equals 0.21 ± 0.12% m−1, with a maximum precipitation at an elevation of 5500 masl. We showed that the total annual basin precipitation (828 mm) is 260% higher than what is estimated based on interpolated observations (319 mm); this has major consequences for hydrological modeling and water resource assessments in general. Our results were validated by using previously published studies on individual glaciers as well as the water balance of the Hunza basin. The approach is more widely applicable in mountain ranges where precipitation measurements at high altitude are lacking.


Introduction
The Indus basin in Pakistan contains the largest continuous irrigation system in the world, and millions of people are dependent on the food produced in this area (Immerzeel et al 2010;Kreutzmann 2011). The Karakoram mountain range plays a crucial role in sustaining the Indus basin's water resources. Some of the world's largest glacier systems outside the polar regions are found in this area, and glacial and snow melt are important hydrological processes. A key uncertainty in this region, as in many mountain ranges in the world, is the spatial distribution of precipitation: a crucial factor when trying to understand high-altitude hydrology. The relation between precipitation and the Asian mountain ranges remains poorly defined due to remoteness and the lack of reliable rainfall networks (Bookhagen and Burbank 2006). Over short horizontal distances, precipitation can vary enormously due to orographic effects, but precipitation gauge networks are virtually nonexistent at high altitudes, and, even if they exist, they are mostly located in the valley bottoms and have difficulty in capturing snowfall. Snow-accumulation measurements by using snow pillows, snow courses, pits, and cores from accumulation zones are also scarce and usually confined to short observation periods.
With the advent of remote sensing-derived estimates of precipitation, such as the Tropical Rainfall Monitoring Mission (TRMM) (Huffman et al 2007), the quantification of spatial patterns has improved, but strong biases with rain gauges remain. For the entire Indus basin, TRMM data were bias corrected by using station data, which showed a significant overall improvement (Cheema and Bastiaanssen 2012). However, the coarse spatial resolution (,25 km) is insufficient in most hydrological basin-scale studies in mountainous terrain where orographic precipitation plays a major role. Techniques have been developed to downscale TRMM further to a 1-km resolution based on the response of vegetation (Immerzeel et al 2009b), but this is unsuitable in the Karakoram given the sparseness of vegetation. Therefore, the applicability of TRMM and other satellite-derived data sets in mountainous areas remains unreliable, owing to the spatial resolution and the poor capacity to detect snowfall when using satellite sensors.

Systems knowledge
Mountain Research and Development (MRD) An international, peer-reviewed open access journal published by the International Mountain Society (IMS) www. mrd-journal.org In many hydrological studies, simulation models are used. These models are forced by precipitation data observed at precipitation gauges located in valleys (Singh and Bengtsson 2004;Rees and Collins 2006;Singh et al 2006;Immerzeel et al 2011). Precipitation measurements are subject to large systematic errors that can add up to 30% and even more in those cases in which a significant part of precipitation falls in the form of snow and undercatch prevails (Rubel and Hantel 1999;Cheema and Bastiaanssen 2012). In addition, valley stations in mountain regions are generally not representative for basin precipitation because of strong vertical precipitation lapse rates. Several regional studies show that this is particularly valid for the Karakoram. During an expedition in 1974-1975, The Baltura Investigation Group already observed 1000-1300 mm of annual snow accumulation in the accumulation zone of the Batura glacier and reported a factor 10 difference between valley precipitation and what was observed in the accumulation zone (The Batura Glacier Investigation Group 1979). Goudie et al (1984) also identified this difference and estimated that total precipitation in the glacier source areas in the Hunza basin are in the order of 1000-2500 mm/y. Winiger et al (2005) combined in situ measurements of snow depth and water equivalent, remotely sensed data on snow cover, and runoff characteristics to estimate the vertical spatiotemporal distribution of precipitation. They concluded that precipitation at 5000-m elevation in the Western Karakoram can be a factor 6 higher than what is observed at valley stations.
In other studies, it is claimed that, in the Karakoram, there is a 5-to 10-fold increase in precipitation between the glacier termini around 2500 m and the accumulation zones above 4800 masl, and maximum precipitation occurs between 5000 and 6000 m (Hewitt 2005(Hewitt , 2007. In the Baltoro catchment, located to the northwest of the Hunza basin, the average annual precipitation at 3200 m was estimated at 300 mm, increasing to up to 2500 mm at 8000 m . Several other studies deployed different proxies, such as geochemical analysis and isotopes, to quantify precipitation distribution with elevation in the Karakoram. They all consistently report strong vertical lapse rates and a peak in precipitation between 5000 and 6000 mm (Wake 1989;Young and Schmok 1989;Young and Hewitt 1990;Hewitt 2011).
Besides this direct evidence, the water balance also provides an important proxy for the existence of strong vertical precipitation lapse rates. In a recent study, the hydrological regime of the Hunza basin in the Karakoram was assessed by using a combination of satellite-derived snow cover, streamflow observations, and observed meteorological data. The observed streamflow at Dainyor Bridge was compared with observed precipitation at Gilgit; it was found that, over the period 1980-2005, the average streamflow equals 750 mm/y, whereas, the average observed precipitation at Gilgit is only 100 mm/y (Tahir et al 2011). This is not a recent phenomenon; Ferguson (1984) also concluded that observed discharges in the Hunza cannot be explained by precipitation observed in Gilgit. Although much less extreme, this is further supported by a reported mismatch between precipitation and observed runoff for the entire Upper Indus Basin (Immerzeel et al 2009a).
We concluded that satellite-derived precipitation estimates and uncorrected precipitation gauge observations are unsuitable to force glaciohydrological models in water resource assessment studies in mountainous regions at the small to intermediate spatial scales. There is a strong need for robust correction methods of precipitation time series that do not depend on expensive field campaigns in accessible terrain.
In this article, we propose such a method by using the glacial mass balance as a proxy to estimate vertical lapse rates in precipitation in the Hunza basin in the Karakoram. We optimized vertical precipitation lapse rates until ablation and accumulation in the catchment areas of the 50 largest glacier systems in the basin corresponded to average observed glacial mass balances. This lapse rate can be used as a first-order approximation of the required correction of precipitation time series. This approach can be applied in similar glacierized regions and can contribute to an improved understanding of hydrological processes due to a more reliable forcing. It is expected that an approximation of the vertical precipitation lapse rate will reduce equifinality problems pertinent to many hydrological modeling studies in which different parameter sets result in similar model outcomes (Pellicciotti et al 2012).

Study area
The Hunza basin is located in the upstream part of the Indus basin ( Figure 1) and has a total catchment area of 14,234 km 2 ; 27.6% of the area is glacierized (3930 km 2 ). Large parts of the ablation areas are generally debris covered, with a total area of 1239 km 2 (31.5% of the total glacier area). The Hunza basin contains several large glacier systems such as the Hispar glacier (521.2 km 2 ) and the Batura glacier (336.7 km 2 ). The elevation ranges from 1394 masl to 7849 masl, and average daily temperature ranges from 29.0uC in January to 14.6uC in July, based on observations at Ziarat (3669 masl). The snow cover in winter is approximately 80%; this decreases to 30% in summer, based on a 2000-2004 Moderate Resolution Imaging Spectroradiometer (MODIS)-derived data set (Tahir et al 2011).
The precipitation regime of the Hunza basin is complex, with 2 principal sources of precipitation. The dominant source in terms of total amounts of precipitation delivered is the monsoon, bringing storm systems from the south from late July through September. Monsoon rains result in heavy precipitation in the greater Himalayan region, but sometimes the monsoon is strong enough to penetrate the central Karakoram. The second source of precipitation originates from depressions that come from the west. This is the major moisture source during the winter months. A secondary summer accumulation season at high elevation in the Karakoram and Hindu Kush usually derives from the same source. These westerly sources provide the dominant nourishment for the glacier systems of the Karakoram (Young and Hewitt 1990). These patterns are confirmed by the observed average monthly precipitation from 2001 to 2003 at 3 meteorological stations in the Hunza basin: Naltar, Khunjerab, and Ziarat (Figures 1, 2). Naltar is located in the south of the basin and has the highest precipitation and the strongest influence of the monsoon (June, July, August [JJA]). The JJA variability is also the highest in Naltar, underlining the strong variability in monsoon strength. Precipitation in Ziarat and Khunjerab is more evenly distributed over the seasons, with less variability, which confirms reliable precipitation that results from westerly disturbances.
The average precipitation observed at the 3 different stations from 2001-2003 in the Hunza basin was 331 mm/y and the average streamflow at Dainyor bridge during the same period was 310 m 3 /s, which equals 688 mm/y. The average streamflow is thus more than twice as high as the observed precipitation, which illustrates the strong undercatch at the valley stations. The streamflow generated in the Hunza eventually flows into the Tarbela

Methodology
Whereas, elsewhere in Asia, glaciers have shown, similar to global trends, predominantly negative mass balances over the past decade, the glaciers in the Karakoram seem to be stable or even slightly increasing in mass. This has been attributed to increases in winter precipitation and lowering of summer temperatures (Archer and Fowler 2004;Hewitt 2005;Fowler and Archer 2006;Tahir et al 2011). Our hypothesis is that, under the assumption of a neutral glacial mass balance over the years 2001-2003, we can estimate the vertical precipitation lapse rate by balancing the total accumulation in the catchment areas of the glaciers with the total ablation on the glacier surface. All analyses performed to test this hypothesis were conducted in the PCRaster dynamic modeling language (Karssenberg et al 2001).
The 50 largest glacier systems were extracted from a glacier inventory (Mool et al 2005); for each of the glacier systems, the total catchment area was determined by using the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Map (GDEM), which was resampled and projected from the nominal 30-m resolution to a 100-m resolution and to an equal area projection. All calculations were performed on this 100-m resolution raster.
To account for horizontal differences in precipitation and temperature, we used daily data from 3 different meteorological stations in the Hunza basin: Khunjerab, Ziarat, and Naltar ( Figure 1). First, we interpolated daily precipitation by using a spline interpolator, which resulted in a smooth precipitation field (P obs (x,y)). We then used the same interpolator to construct a reference elevation field (h r (x,y)) by interpolating the station elevations (Table 1). The precipitation lapse rates are relative to h r (x,y). We then assumed that precipitation increases linearly to an elevation of 5500 masl (Young and Hewitt 1990) and decreases with the same rate at higher elevations. For elevations below 5500 masl, we derived the spatial precipitation field through: For elevations higher than 5500 masl, we derived the spatial precipitation field through: where P(x,y) is the precipitation in millimeters at location x, y, P obs is the observed precipitation at Ziarat, h (x,y) is the elevation (masl), and c is the precipitation lapse rate (%/m). Daily temperatures at the station locations were first converted to sea level temperatures by using a temperature lapse rate of 20.0068uC/m, and then interpolated by using the same spline interpolator as for precipitation. We then added to this temperature field the product of the digital elevation model and the temperature lapse rate. We then estimated the total accumulation (A) for each glacier system by aggregating daily snowfall over its catchment area. Subsequently, we determined glacial melt by using a simple degree day factor (DDF) approach with different DDFs for debriscovered (DDFDC) and clean-ice (DDFCI) glaciers (Table 2), and we aggregated total melt (M) over each glacier's surface. By using a linear forecast function, we estimated c by matching total accumulation with total ablation over all glacier systems considered: where i is an index for each glacier system under consideration. The uncertainty in the optimized precipitation lapse due to uncertainty in DDFDC, DDFCI, and temperature lapse rate (TL) was taken into account by sampling 20 parameter sets, by assuming that the parameter values were normally distributed (Table 2). For each parameter set, c was determined, and we estimated the error in c by calculating the standard deviation of the 20 simulations. Average DDFs and associated standard deviations for clean-ice and debris-covered glaciers under similar climatic conditions were derived from the literature (The Batura Glacier Investigation Group 1979;Mihalcea et al 2006;Zhang et al 2006;Konz et al 2007;Immerzeel et al 2011). The average TL and its standard deviation were based on the range between the dry and saturated adiabatic lapse rates (Table 2).

Results and discussion
In Figure 3, the distribution of average annual mass balances is shown for all 50 glacier systems. When the precipitation lapse rate is ignored, there was a strong negative overall mass balance (21.48 m/y). After optimization, the overall mass balance averaged over all glacier systems is 0.00 m/y. There is substantial variation between the different glacier systems, as melt rates, aspect, nourishment from avalanches, debris cover, and topographical differences together define the behavior of each system (Table 3). However, by averaging over all glacier systems, it is possible to derive a representative precipitation lapse rate and the optimized c is 0.21%/m. This lapse rate is positive up to an elevation of 5500 m and negative at higher elevations. The average annual precipitation map with the south to north gradient clear visible is shown in Figure 4, which corresponds to the higher precipitation observed at Naltar than at Khunjerab and Zialtar (Table 1). The Batura wall forms a precipitation divide; north of the Batura wall precipitation is much lower than on the southern side (Winiger et al 2005). The map also shows the positive relation between elevation and precipitation up to elevation to 5500 m and a reverse signal at higher elevations.
The precipitation distribution per elevation band is shown in Figure 5. Averaged over the entire basin, the elevation peaks at an elevation around 6000 masl at 1174 mm/y. It is interesting to note that this is higher than 5500 masl because of the elevation distribution in the basin. The drier northern part has a relatively large area around the 5500 masl elevation band, and this reduces the average. In the Hunza valley, precipitation is only around 180 mm/y. Locally, in particular, in the vicinity of the Batura wall, precipitation can be much higher and values up to 3000 mm/y are estimated at the south side of the Batura wall. The overall basin precipitation after taking into account the vertical precipitation lapse rate is 828 mm/y. If the data of the 3 stations are interpolated without accounting for elevation effects, then the total basin precipitation equals 319 mm/y. These findings are in good agreement with previously published studies.
There is large variation in the melt rates, so we used a simple DDF approach with different DDFs for debriscover and clean-ice glaciers. For individual glacier systems, the actual melt rates may vary considerably from the assumed melt rates because a DDF approach is a simple approximation of a full energy balance approach, variation in surface conditions, glacier aspects, and thickness of debris cover; this variation is shown in Figure 3. These factors all result in random errors; therefore, the integration of values over the 50 major glacier systems does provide a plausible estimate of the average precipitation lapse rate. To further account for both this uncertainty in DDFs and in temperature lapse rate, we derived different 20 parameters sets, which were used to estimate the standard deviation in the precipitation lapse rate. This analysis shows that the standard deviation in c equals 0.12%/m.
In this study, we derived precipitation fields by using a single vertical precipitation lapse rate; this means that precipitation will linearly increase with elevation. Several articles indicate that, in the central Himalaya, monsoon precipitation peaks at an elevation around 3000-4000 m (Burbank et al 2003;Putkonen 2004;Bookhagen and Burbank 2006), but, for the Karakoram, this may be different, and precipitation may increase with elevation, up to much higher elevations (Young and Hewitt 1990;   Hewitt 2007). This may be because monsoon precipitation is mostly convective, whereas, in the Karakoram, stratiform precipitation is also abundant. Although our method is applicable in all glacierized mountain basins, the elevation of maximum precipitation is unique to each climate zone and cannot be generalized.
Our results are validated by comparing the interpolated observed precipitation and the lapsed precipitation fields with the observed river flow at Dainyor bridge ( Figure 6). The monthly river flow regime shows a peak in July, whereas precipitation peaks in February-March and in August. This shows that the Hunza river flow regime is dominated by meltwater.  The average river flow over the 2001-2003 period equals 688 mm/y, but the interpolated observed precipitation is only 319 mm. In other words, the total river flow is 215% of interpolated observed precipitation. If we assume that the gap between observed river flow and observed precipitation is closed by excess glacial melt, then this would mean a negative specific mass balance for all glaciers in the Hunza of 1.34 m/y, which is unrealistic given the neutral or slightly increasing mass balance that has been reported in this period. The lapsed precipitation, however, equals 828 mm/y, and, based on the observed river flow, this would result in a runoff coefficient of 83%, which is plausible for a mountainous basin with steep slopes and limited evapotranspiration.

Conclusions
In this study, we estimated vertical precipitation lapse rates in the Hunza basin in the Karakoram. By assuming a neutral glacier mass balance, we inversely modeled the required precipitation to balance total accumulation and total ablation for 50 major glacier systems, which allows us to draw the following key conclusions. Using glacier mass balances as proxy to derive vertical precipitation lapse rates can be a very powerful technique; it is more widely applicable in regions where high-altitude precipitation measurements are scarce. This method is straightforward, does not require extensive field campaigns, and can be used to correct precipitation time series for hydrological simulation models and water resource assessment studies. There is a very strong positive relation between precipitation and elevation in the Hunza basin in Pakistan. We estimated that the vertical precipitation lapse rate is equal 0.21 6 0.12% m 21 under the assumption that precipitation increases with elevation up to 5500 masl and the lapse rate reverses at higher altitudes. In this case, this implies that the total basin precipitation is 260% higher than what is observed at the valley stations (319 mm/y), which FIGURE 5 Precipitation distribution per elevation zone for the entire Hunza river basin. obviously has major consequences for hydrological modeling and water resource assessments in general. Our results were validated by using previous published studies on individual glaciers and were based on the water balance at the outlet of the Hunza basin. The uncertainty in the estimated lapse rate is considerable, and the coefficient of variation is 57%. The main sources of uncertainty are the melt rates of clean-ice and debris-covered glaciers and the temperature lapse rates.