Brief Communication : Glaciers in the Hunza Catchment ( Karakoram ) are in balance since the 1970 s

Previous geodetic estimates of mass changes in the Karakoram revealed balanced budgets or a possible slight mass gain since the year ~2000. Indications for longer-term stability exist but no mass budget analyses are available before 1


Introduction
Glacier melt water is of high importance for the run-off of the Indus River (Immerzeel et al., 2010) but the exact glacier share is not known.This is partly due to the lack of knowledge about precipitation, snow cover, and snow water equivalent, but also about glacier mass balance, their characteristics and their responses to climate change.Karakoram glaciers, which occupy a large portion of the glacierized area of the Indus basin, have recently shown unusual behaviour: on average no significant area changes but frequent advances and surge activities have been observed during the last decades (Bhambri et al., 2013;Bolch et al., 2012;Hewitt, 2011;Copland et al., 2011).Geodetic mass estimates revealed balanced glacier mass budgets or even slight mass gain since ~2000 (Rankl and Braun, 2016;Gardelle et al., 2013;Kääb et al., 2015).However, no mass budget analyses are available for the period before 2000.Herreid et al. (2015) found no significant change in debriscoverage of the glaciers in the Hispar and Shimshal sub-regions of the Hunza River basin for the period 1977 until 2014 and concluded that this might be due to balanced glacier budgets during this period.Temperature measurements that are available since 1961 in the Karakoram show, in contrast to many other regions of the globe, a consistent decline in summer and an increase during winter (Fowler and Archer, 2006).Hence, these measurements would support the assumption that glaciers would be in balanced or slightly positive conditions over the last several decades of the 20 th century.
Declassified stereo satellite images from the 1960s and 1970s such as Corona KH-4 and Hexagon KH-9 have been proven to be suitable to generate digital terrain models and assess glacier mass changes since the 1960s (Bolch et al., 2008;Pieczonka et al., 2013).Hence, the aim of this study is to revisit existing information and extend the time series back to some of the The Cryosphere Discuss., doi:10.5194/tc-2016Discuss., doi:10.5194/tc- -197, 2016 Manuscript under review for journal The Cryosphere Published: 9 September 2016 c Author(s) 2016.CC-BY 3.0 License.earliest available satellite imagery.We focus on the Hunza catchment in the Central Karakoram (Figure 1) where high heterogeneity of glacier behaviour was found in previous studies (e.g.Bolch et al., 2012;Quincey and Luckman, 2014).
Moreover, suitable Hexagon KH-9 data from the 1970s and recent stereo data from ~2010 such as ASTER and Cartosat-1 data were available.The Hunza River is a tributary to the Gilgit River, which flows into the upper Indus.The area of the basin is about 13,715 km² and approximately 25% of the basin is covered by glaciers.These glaciers constitute more than 15% of the glacierized area of the entire Karakoram.Frequent surges reported for several glaciers in this basin (Quincey and Luckman, 2014;Copland et al., 2011;Rankl et al., 2014) complicatethe analysis of mass budgets as often only a certain part of a surge is captured.

Data
The SRTM digital terrain model (DTM) version 4, with a spatial resolution of 1 arc second (~30m, SRTM1) was utilized as reference dataset.The SRTM1 DTM was acquired by the use of two C-Band radar antennas (operating in interferometric mode) during 11 -22 February 2000 and is frequently used for glaciological investigations.It can be assumed that the represented ice surface of the ablation region is close to the surface at the end of the 1999 ablation period, assuming a full penetration of the radar beam into snow (Paul and Haeberli, 2008).However, deeper radar penetration can be expected in the accumulation region (Berthier et al., 2006).Data voids which are mainly restricted to some accumulation areas were filled with ASTER GDEM2 data.
The Cryosphere Discuss., doi:10.5194/tc-2016Discuss., doi:10.5194/tc- -197, 2016 Manuscript under review for journal The Cryosphere Published: 9 September 2016 c Author(s) 2016.CC-BY 3.0 License.This DTM is a merge of several ASTER scenes covering the period ~2000 -2010 (Slater et al., 2011).We used ten on demand generated ASTER DTMs (product AST14DEM) from five different acquisition dates of the year 2008-2010.A small missing stripe was filled with a DTM from 2013 scenes (Table S1, Fig. S1).ASTER scenes were visually checked and from the most promising available (no clouds, minimum snow cover) the respective DTMs were ordered and used for DTM differencing.
To estimate rates of elevation change, an acquisition year was assigned for each glacier.We used the mean year for glaciers covered by two scenes acquired in different years.Two high-resolution Cartosat-1 stereo scenes captured on 11 July 2010 (Table 1) were used to compare the results gained with the lower resolution ASTER DTM.Cartosat-1 (IRS-P5) was launched by Indian Space Research Organisation (ISRO) in May 2005.The satellite has two high resolution (2.5 m) panchromatic sensors recording stereo images along the track (Titarov, 2008).The major advantage of this dataset is besides the high spatial resolution and also the 12 bit pixel depth.On the other hand, the spatial coverage is relatively small (25 x 25 km) and our two stereo pairs cover only one large glacier (Khurdopin Glacier, glacier nr.6 in Figure 1) in full.Declassified Hexagon KH-9 imagery allowed us to extend the time series back to 1973 (Table 1).The KH9-Hexagon mission was part of the US keyhole reconnaissance satellite program whose images were declassified in 2002.Imagery from this program have already been applied to investigate glacier mass changes (e.g.Pieczonka et al., 2013).
The ICIMOD glacier inventory (Bajracharya and Shrestha, 2011), also available through the GLIMS data base (www.glims.org),was used as a baseline data set and manually adjusted based on the utilized optical imagery for DTM generation and Landsat ETM+ scenes of the years 2000 and 2001.

DTM generation, postprocessing, and uncertainty
All KH-9 DTMs were generated with Erdas Imagine 2014 Photogrammetry Suite using the frame camera model with a focal length of 30.5 cm.Image pre-processing, comprising the elimination of internal distortions and reseau grid removal, has been performed following Pieczonka et al. (2013).GCPs have been collected from Landsat 7 ETM+ imagery with SRTM1 as a vertical reference (Pieczonka et al., 2013).Fiducials were measured manually considering the principal point in the image centre.All stereo images have been processed with a RMS of <~1.5 pixels.The final Hexagon KH-9 DTMs cover the entire Hunza basin.A small gap of about 20 pixels between the two generated DTMs exists.
The Cartosat-1 stereo pairs have also been processed using PCI Orthoengine 2014 with 32 and 35 GCPs.To improve the quality of the DTMs, image enhancement techniques were applied prior to DTM generation in order to overcome low image contrast and temporal differences in image acquisition.The root mean squared error (RSME) varied between 0.3 and 3.9 pixels (Table S2).The spatial resolution of all DTMs was chosen as 30 m.
In order to obtain reliable results on glacier surface elevation changes, the DTMs must be properly co-registered (Nuth and Kääb, 2011).As we observed tilts when differencing the original DTMs, we first minimized elevation differences between the different DTMs with respect to the SRTM1 master DTM based on spatial trend corrections.We considered only elevation differences (Δh) between ±150m over non-glacierized terrain with slopes less than 15° (Bolch et al., 2008;The Cryosphere Discuss., doi:10.5194/tc-2016-197, 2016 Manuscript under review for journal The Cryosphere Published: 9 September 2016 c Author(s) 2016.CC-BY 3.0 License.Pieczonka et al., 2013).Subsequently, all DTMs were further co-registered following the approach by Nuth and Kääb (2011).The final displacement between all DTMs and SRTM1 were less than or equal to one pixel (≤30 m) on average.Data voids and mismatches that result in incorrect elevation values can occur in areas with low image contrast such as cast shadows and bright snow.Mismatches due to snow in the accumulation regions led often to unrealistic low elevation values using KH-9 data that would subsequently lead to unrealistic surface lowering values in parts of the accumulation region (Pieczonka and Bolch, 2015).However, thickness change distributions for glaciers with negative mass budgets typically have a minimum lowering at the glacier head with increasing values towards the glacier front following a non-linear trend (Huss et al., 2010).This pattern is different for surging glaciers that often exhibit high positive Δh values at comparatively low elevations, strong surface lowering around the ELA, and then decreasing Δh values towards the upper reaches (Gardelle et al. 2013, Rankl andBraun, 2016).Elevation change patterns can also be affected by thick debris cover where the highest lowering usually does not occur close to the terminus but upstream (Bolch et al., 2008).
As both debris-covered glaciers and surge-type glaciers are common in the investigated region we could not apply a general threshold to remove Δh outliers, but used the general assumption that lower elevations show stronger Δh variability than higher elevations.This should still be true for surging glaciers and those with balanced conditions.The related calculations followed the approach by Pieczonka and Bolch (2015) and used a sigmoid function allowing a larger range of Δh values in the middle part of the ablation region to preserve the signal of surging glaciers and a narrower range (-20≤Δh≤20 m) at the glacier head.We filled all data gaps (including the gap between the two KH9 DTMs) by means of ordinary kriging in order to get the weighted moving average based on neighbouring pixel values.
The penetration of the radar beams into firn and snow has to be considered in case of the comparison of DTMs generated from microwave data such as the SRTM1.However, the value can only be estimated as it depends on several unknown parameters (e.g.snow depth and characteristics) and is therefore one major source of uncertainty (Kääb et al., 2015;Gardelle et al., 2013).We applied the correction suggested by Kääb et al. (2012), who analysed the beam penetration of the C-band SRTM data in a similar region of the Karakoram and found a penetration of 2.4 ±1.4 m.The conversion of volume to mass changes needs to consider the combined ice and snow density.As both are unknown, we used a density of 850 ±60 kg/m³ as a reasonable and widely used assumption for a longer time period (Huss, 2013).
There is no best method to estimate the uncertainty (e) of the DEM differencing when no precise and well distributed GCPs are available.Here, we followed a widely applied approach and calculated the uncertainty following Gardelle et al., (2013): where EΔhi is the standard deviation of the mean elevation change of the non-glacierized terrain of each altitude band and Neff is the effective number of observations.The latter is calculated using the total number of observations Ntot, the pixel size R (30 m), and d is the distance of spatial autocorrelation of the elevation change maps (1025 m) The Cryosphere Discuss., doi:10.5194/tc-2016Discuss., doi:10.5194/tc- -197, 2016 Manuscript under review for journal The Cryosphere Published: 9 September 2016 c Author(s) 2016.CC-BY 3.0 License.
The overall uncertainty of the DEM difference is the average of EΔh weighted by the glacier hypsometry.A further uncertainty to be considered is the uncertainty of the mapped area of the glaciers E a .We assume an uncertainty of 5% which is towards the upper bound of published estimates of the uncertainty of mapped glaciers based on similar satellite data (e.g.Paul et al. 2013).The final uncertainty is calculated considering also the uncertainty of the radar penetration (E p , ±1.4 m) and of the volume to mass conversion (E m , ±7% of the elevation change).
We did not apply a seasonality correction as most of the images were acquired close to the end of the ablation period but assume the effect is well within the considered uncertainty.

Glacier volume and mass changes
Results for the period 1999 -~2009 show heterogeneous glacier behaviour, with several surging glaciers in the study region (Figure 2).A northern tributary of Hispar Glacier thickened by approximately 150 m at the confluence of the glaciers.
Khurdopin and its neighbouring glaciers in south-eastern Shimshal Valley both show significant thickening and thinning within their tongues which is typical of a surge occurring during the study period (Figure 2).In contrast, the large debriscovered Batura Glacier west of Shimshal Valley showed surface lowering throughout the tongue leading to an overall volume loss (Figure 2, 3).For the entire study region we found no significant mass changes (Table 1) which is in line with previous results for the period after 1999 in Central Karakoram (e.g.Gardelle et al., 2013).In addition, we add information for glaciers east of their study region (e.g.Batura Glacier) and cover the entire Hispar Glacier, one of the largest glaciers (length 50 km) in the Karakoram.For this region we cannot confirm positive mass budgets but our study indicates a slight mass loss similar to Rankl and Braun (2016) and Kääb et al. (2015).However, our uncertainties are larger due to the utilization of the lower resolution, but freely available ASTER DTM.In addition, we want to emphasize that the analysed glaciers and period slightly differ which could be a reason for the (not significant) differences to the existing studies.
The Cryosphere Discuss., doi:10.5194/tc-2016Discuss., doi:10.5194/tc- -197, 2016 Manuscript under review for journal The Cryosphere Published: 9 September 2016 c Author(s) 2016.CC-BY 3.0 License.Our extended time series show that glaciers in Hunza Valley experienced no significant mass changes but a heterogeneous behaviour also for the period 1973 to 1999 (Table 1).Hence, central Karakoram glaciers were, on average, in balance for at least the last 40 years.Although longer-term balanced budgets could be assumed based on existing information about (insignificant) area changes since the 1970s (Bhambri et al., 2013), on average similar debris coverage since 1977 (Herreid et al., 2015), we confirm for the first time balanced-budgets using elevation differencing.
Most glaciers experienced similar mass budgets for both investigation periods.However, it seems that some glaciers had more negative budgets in the second period than before 2000.This is especially true for the debris-covered Batura Glacier whose tongue showed significant lowering during 1999-2009.Over the entire study period there is also no significant difference in the mass budgets of surge-type and non-surge-type glaciers, a result also found by Gardelle et al. (2013).Surgetype glaciers, however, showed different surge stages in the two periods and more negative values in the recent period.Kurdophin Glacier, for example, experienced a significant thickening near the snout and a significant lowering around the ELA, both combined resulting in an about zero mass budget (Figure 3).Different tributaries of Hispar Glacier also show surge events (Figure 2).
A further source of uncertainty are the data voids in the original SRTM data as about 10% of the glacier area in our study regions are affected and previous studies showed that there can be significant deviations to reality in the data used to fill the voids (e.g.Kääb et al. 2012).In addition the time stamp of the data is often unknown.We compared therefore the results of

Figure 1 :
Figure 1: Overview map of the study region

Figure 2 :
Figure 2: Elevation difference between the Hexagon KH-9 and the SRTM1 DTMs (left) and the SRTM1 and ASTER DTMs (right).Black dots indicate surge-type glaciers, the numbers selected larger glaciers: 1: Batura, 2: Pasu, 3: Barpu, 4: Hispar, 5: Yazghil, 6: Khurdopin, 7: Vijerab the void filled and the non-void filled version.The resultant value for both periods differs only by about 0.02 m a -1 , which is well within the estimated uncertainty.The minor differences are, besides the possibility that the ASTER GDEM data used for void filling are quite reliable for the study region, due to the fact that most of the void are located in the accumulation regions where we restricted the maximum possible deviation.Only one larger glacier is also affected in the ablation region.This is Shishpar Glacier a south exposed glacier located south of Batura Glacier.The void-filled data allowed to detect the surge activity between 1973 and 1999 with an estimated mass budget of +0.04 ± 0.19 m w.e. a -1 .

Figure 3 :
Figure 3: Longitudinal profiles of surface elevation changes for selected glaciers for the two periods.