www.the-cryosphere.net/5/349/2011/ doi:10.5194/tc-5-349-2011 © Author(s) 2011. CC Attribution 3.0 License. The Cryosphere

Abstract. Mass loss of Himalayan glaciers has wide-ranging consequences such as changing runoff distribution, sea level rise and an increasing risk of glacial lake outburst floods (GLOFs). The assessment of the regional and global impact of glacier changes in the Himalaya is, however, hampered by a lack of mass balance data for most of the range. Multi-temporal digital terrain models (DTMs) allow glacier mass balance to be calculated. Here, we present a time series of mass changes for ten glaciers covering an area of about 50 km2 south and west of Mt. Everest, Nepal, using stereo Corona spy imagery (years 1962 and 1970), aerial images and recent high resolution satellite data (Cartosat-1). This is the longest time series of mass changes in the Himalaya. We reveal that the glaciers have been significantly losing mass since at least 1970, despite thick debris cover. The specific mass loss for 1970–2007 is 0.32 ± 0.08 m w.e. a−1, however, not higher than the global average. Comparisons of the recent DTMs with earlier time periods indicate an accelerated mass loss. This is, however, hardly statistically significant due to high uncertainty, especially of the lower resolution ASTER DTM. The characteristics of surface lowering can be explained by spatial variations of glacier velocity, the thickness of the debris-cover, and ice melt due to exposed ice cliffs and ponds.


Introduction
Recent debate on whether Himalayan glaciers are shrinking faster than those in other parts of the world (Cogley et al., 2010) highlighted the lack of knowledge about the glaciers in this region.The best measure of glacier health is mass balance, which can be directly linked to climate and compared Correspondence to: T. Bolch (tobias.bolch@geo.uzh.ch) to other regions.The topographic setting and the glacier hypsometry modify the way the climatic signal is translated into glacier mass balance, but this can calculated using a DTM (Paul and Haeberli, 2008).However, only a few in-situ mass balance measurements have been made on Himalayan glaciers, and existing data series are short (Kulkarni, 1992;Fujita et al., 2001;Wagnon et al., 2007;Dobhal et al., 2008).Comparisons of digital terrain models for different years can complement field measurements, and allow regional mass balance to be estimated (Bamber and Rivera, 2007).However, to date it has only been applied to some glaciers in the western Himalaya for 1999 to 2004 (Berthier et al., 2007) and to four glaciers at Mt. Everest for 1962 to 2002 (Bolch et al., 2008b).Broader and more detailed knowledge of glacier mass balance is also needed to decrease the high uncertainty about the importance of Himalayan glaciers for water resources (e.g.Immerzeel et al., 2010) and sea level rise (e.g.Braithwaite and Raper, 2002).Finally, improved knowledge of glacier recession is needed to better estimate hazards of glacier lake outburst floods (GLOFs; Richardson and Reynolds, 2000).
The aims of this study are threefold.First, we aim to evaluate the results of the pilot study by Bolch et al. (2008b) using independent data sets.This study revealed surface lowering of 0.33 ± 0.24 m a −1 by analysing a 1962 Corona and an ASTER DTM generated based on 2001, 2002 and 2003 data but had high uncertainties.The second aim is to present mass balance estimates for a larger sample of glaciers around Mt. Everest, including Imja Glacier which is of high interest due to the proglacial lake which formed in the 1960s and has rapidly grown since (Bolch et al., 2008a;Fujita et al., 2009).In addition, the mass balance of the entire Khumbu Glacier will be presented for the first time.Thirdly, we wanted to produce the first time-series of glacier volume changes around Mt. Everest to show the suitability of different optical imagery to derive mass balance variability over time and to discuss the possible causes of the surface changes.The tongues of nine of the ten studied glaciers are heavily covered by supraglacial debris (Fig. 1), with a debris area of about 36.3 km 2 and an average debris thickness increasing downglacier (Nakawo et al., 1999;Hambrey et al., 2008).The glaciers are mainly nourished by snow and ice avalanches which accumulate cones below the steep headwalls.Only Khumbu Glacier has an extensive accumulation area (called "Western Cwm").Glacier equilibrium line altitudes (ELAs) are roughly estimated to be situated above 5600 m (Asahi, 2001).Ice velocities typically decrease downglacier from the ELA and most of the studied glaciers have extensive stagnant ice in their lower reaches (Bolch et al., 2008a;Scherler et al., 2008;Quincey et al., 2009).Between 1962 and 2005 the overall glacier area loss in the study area was ∼5% with an increasing debriscovered area (∼2.5%) but almost stable terminus positions (Bolch et al., 2008b).

Data and methodology
We used 1970 Corona KH-4B (declassified US spy imagery) data, 1984 aerial photographs (camera: Wild RC 10) (Altherr and Grün, 1990) and 2007 Cartosat-1 (Indian Remote Sensing Satellite, IRS P5) images (Table 1).In addition, we used previously generated 1962 Corona and 2002 ASTER DTMs (Bolch et al., 2008b).We did not consider the DTM data from the Shuttle Radar Topography Mission (SRTM) due to large data gaps especially at the clean ice area of Khumbu Glacier and the coarser spatial resolution (∼90 m) in comparison to the ASTER DTM (30 m).In addition, using only stereo optical data results in a methodologically consistent data set.The SRTM DTM data in contrast are based on C-band adar beams that penetrate into snow and can cause higher uncertainties in the snow-covered accumulation areas.We applied the Remote Sensing Software Package  Graz (RSG, http://dib.joanneum.at/rsg/)6.13 for processing Corona, PCI Geomatica OrthoEngine 10.2 for Cartosat, and Leica Photogrammetry Suite (LPS) 9.1 for the aerial images.We used 14 non-differential GPS points acquired in 2006 and 2008 and ∼200 points from the National Geographic 1:50k topographic map (Altherr and Grün, 1990) as ground control points (GCPs).The GPS points were mostly measured along the main trekking routes within a height range between 3900 m and 5600 m.Their horizontal accuracy is about 7.9 m and their vertical accuracy in comparison to topographic map height points is about 20.6 m.The GCPs based on the map are almost randomly distributed in accordance with the un-systematic distribution of significant peaks and other geomorphological forms marked as height points.In addition, these points represent different slope angles and aspects well (Pieczonka et al., 2011).The RMSE z and RSME x,y of the map were computed based on the GPS points to be 20.6 m and 17.8 m, respectively.This matches the results achieved by Altherr and Grün (1990).In addition, we used automatically selected tie points (TPs) to improve the sensor models.The overall quality of the generated raw DTMs appears promising as the glacier tongues are almost fully represented (Figs. 2, 3).Data gaps occur mainly due to snow cover and cast shadow.In order to measure glacier elevation changes as precisely as possible it is recommended to adjust the DTMs relative to each other (Nuth and Kääb, 2011).Tilts which occurred especially in the Corona DTM were corrected using trend surfaces calculated based on stable non-glacierised areas of the DTMs (Fig. 2, Pieczonka et al., 2011).We observed slight horizontal shifts of the generated DTMs although we used the same GCPs for all images whenever possible.In order to avoid biases introduced thereby and to improve the z-accuracy, we choose the Cartosat-1 DTM as the master reference as it has a high spatial resolution and showed the lowest mean elevation difference (5.9 m) and standard deviation (18.3 m) relative to the SRTM3 DTM.We co-registered the other DTMs to it by minimizing the standard deviation of the elevation differences (Berthier et al., 2007).The applied shifts varied between 5 and 30 m. Altitudinal differences which exceeded ±100 m (usually around data gaps and near DTM edges) were omitted assuming that these values represent outliers.Similar assumptions were made by Berthier et al. (2010).We resampled all DTMs bilinearly to the pixel size of the coarsest DTM (30 m) in order to reduce the effect of different resolutions.
The relative uncertainties of the DTMs prior to and after the adjustments were calculated based on the ice free terrain relative to the 2007 master DTM.The mean difference between the final adjusted DTMs was in the range −0.1 to −1.8 m while the RMSE z was 7.8 to 19.8 m (Table 1).To address the uncertainty of the elevation differences of the glacierised areas we calculated statistical parameters for the differences of ice covered and the non-ice covered areas separately (Table 2).The standard deviation of the non glacier area (STDV noglac ) or the RSME z can be used as a first esti-mate of the uncertainty, but would probably overestimate it as it does not account for the reduction of the error due to averaging over larger regions (Berthier et al., 2007(Berthier et al., , 2010)).Another suitable measure for the uncertainty is the standard error of the mean (SE) defined as while n is the number of the included pixels.We choose a decorrelation length of 600 m for the ASTER DTMs with an effective spatial resolution of 30 m and a length of 400 m for all other higher resolution DTMs to minimize the effect of auto-correlation.These numbers are slightly more conservative than the average value of 500 m utilized by Berthier at el. (2010) for DTMs with coarser resolution (mostly 40 m).Koblet et al. ( 2010) suggested a decorrelation length of 100 m for DTMs based on aerial images of a spatial resolution of 5 m.We used the standard error (SE) and the mean elevation difference (MED) of the non glacier area (Table 2) as an estimate of the uncertainty according to the law of error propagation: Volume change was calculated for each glacier assuming that the density profile remains unchanged and that only ice is lost or gained (Paterson, 1994;Zemp et al., 2010).To convert volume changes into mass change, we assumed an ice density of 900 kg m −3 and assigned an additional uncertainty of 7% due to lack of ground truth (Zemp et al., 2010).We interpolated small data voids (<10 pixel) within the ice covered  areas using a spline algorithm.We did not fill the larger data gaps e.g. on steep slopes.The glacier tongues, the avalanche cones, and Western Cwm are represented in the DTMs of 1970, 2002, and 2007 (Fig. 3, Tables 3, 4) which allow estimation of the mass balance for the entire glacier.Only Changri Nup, Duwo, and the debris-free Chukhung Glacier have larger data gaps.Detailed investigations on Khumbu Glacier are limited to the tongue below ∼5700 m (mainly the ablation area) due to the small coverage of the aerial images.
3 Volume changes and mass losses

Periods 1970-2007 and 2002-2007 for the whole study area
Between 1970 and 2007 significant surface lowering occurred on all investigated glaciers (Fig. 3, Table 2).The greatest lowering was on Imja/Lhotse Shar Glacier.Except for this glacier, which displays surface lowering throughout the terminus, most glaciers show maximum lowering in their mid ablation zones, with a negligible change near their termini.Overall ice loss is estimated to be >0.6 km 3 with an average surface lowering of 0.36 ± 0.07 m a −1 or a specific www.the-cryosphere.net/5/349/2011/The Cryosphere, 5, 349-358, 2011  3).The surface lowering for the debris-covered parts only is −0.39±0.07m a −1 , clearly showing that significant thickness loss occurs despite thick debris cover.Most glaciers also experienced surface lowering between 2002 and 2007 (Fig. 3, Table 3).The specific mass balance of all 10 glaciers has possibly doubled compared to 1970-2007 (−0.79 ± 0.52 m w.e. a −1 ).However, the uncertainty is high.

Detailed multi-temporal investigations on Khumbu Glacier
The surface of the ablation area of Khumbu Glacier lowered in all investigated time periods (Table 3).DTM differencing (Figs. 3 and 4) and longitudinal profiles, in particular for 1970-2007, show almost no surface lowering in the clean ice zone below Khumbu Icefall (Figs. 5 and 6, section A); an increasing downwasting in the debris-covered part, with the highest lowering between 2 and 8 km from the terminus (B,  For 1970For -1984 only lowering between ∼1.5 and 5.5 km of the terminus is significant (Fig. 5).Between 1970 and 2007, average surface lowering rate in the ablation area was −0.38 ± 0.07 m a −1 .The rate for 1984-2002 is higher than for 1962-1970and 1970-1984, and comparison of the recent DTMs (2002-2007) ) suggests further accelerated surface lowering (Table 4).However, these differences are statistically not significant.Comparing the periods 1970-1984 and 1984-2007, however, shows an almost statistically significant increase in the surface lowering rate (0.18 ± 0.29 m a −1 in comparison to 0.56±0.13m a −1 ).The accumulation zone of Khumbu Glacier has possibly also thinned during the investigated time, while there might be a slight thickening in recent time (2002)(2003)(2004)(2005)(2006)(2007).

Discussion
Stereo capability, acquisition in the 1960s and 1970s and relatively high spatial resolution make Corona KH-4, KH-4A and B imagery a valuable source for geodetic mass balance estimations.While the earlier KH-4 data were already found to be suitable for this task (Bolch et al., 2008b) we have shown that the later higher resolution KH-4 results in lower uncertainties.The accuracy and even distribution of the ground control points used for the rectification of the imagery are crucial for the resultant accuracy of the DTM.Hence, it could be expected that an even better accuracy than ours could be achieved if more precise GCPs are available.
Other reconnaissance images such as Hexagon KH-9 from the 1970s and early 1980s are also suitable for this task if coverage with no clouds and little snow cover are available (Surazakov and Aizen, 2010).Surazakov and Aizen (2010) have shown a RSME z for mountainous terrain of ∼20 m, which is the same range as our results using KH-4B data despite lower resolution of the KH-9 data.Reseau marks on the KH-9 images facilitate minimizing the distortion.Corona and Hexagon data are also of high value for mapping glaciers and investigating area changes (Bhambri et al., 2011;Bolch et al., 2010;Narama et al., 2010).The generation of mass balance time series estimation of different data sets with different resolution requires careful co-registration and adjustment.This is especially true for the 1984 DTM based on aerial imagery, highlighted by a clear positive-negative trend on hillslopes with opposing aspects for the periods 1970-1984 and 1984-2002 (visible on the lower right of the images on Fig. 4).However, a shifting to adjust this area would worsen the accuracy of the steep slope in the northern part of the DTM.Hence, the distortion of the 1984 DTM is more complex and could not be fully adjusted relative to the others.Although inaccuracies remain on steep slopes with all DTMs, most of the glacier area is not affected by these biases.The inclusion of the lower resolution 2002 ASTER DTM provides some insight in the recent volume loss of the glaciers.However, the volume change of the glaciers is too small for a significant signal to be detected.In addition, the large scatter in the volume change amongst the investigated glaciers needs further investigation.
The quality of the DTMs, however, is supported by local details, such as the area on Khumbu Glacier arrowed on Fig. 4. A large thinning occurred in this area between 1970 and 1984, which can be attributed to the growth of a lake basin visible on the 1984 aerial photograph (Supplement, Fig. 4).Melting and calving around the margins of lakes is well known to produce locally high ablation rates (Sakai et al., 2000;Benn et al., 2001).Thickening occurred in this region between 1984 and 2002, attributable to drainage of the lake reducing ablation and ice inflow from upglacier.
The utilized imagery also allows length changes of the glaciers to be examined in detail.Bajracharya and Mool (2009) argued that Khumbu Glacier has undergone terminus retreat, based on comparison of recent imagery with old topographic maps.Our images, however, indicate that the terminus region of Khumbu Glacier has undergone very little change in recent decades (Fig. 5b, c and Supplement Fig. 1).This highlights the fact that debris-covered glaciers may have stable terminus positions even though they are in negative balance, and that volume changes are more reliable indicators of glacier health than area changes (cf.Scherler et al., 2011).
The calculated average 1970-2007 thickness changes for the whole study area of −0.36 ± 0.07 m a −1 , based on independent data sets, is very close to the value of −0.33 ± 0.24 m a −1 presented by Bolch et al. (2008b) for 1962-2002, but with smaller uncertainty.The wider coverage of this study including Imja Glacier, and the accumulation area of Khumbu Glacier, as well as lower uncertainties and the multi-temporal coverage, allow greater insight into decadal glacier changes and the influence of debris cover.
The longitudinal profile of glacier thinning of Khumbu Glacier shows similar characteristics to those presented by Nakawo et al. (1999) based on estimated ice flow and thermal properties derived from Landsat data.Based on field measurements, Kadota et al. (2000) found a slight surface lowering (∼5 m) about 1 km upglacier of the terminus and higher surface loss (more than 10 m in places) up to the pinnacle zone close to Everest base camp for the period 1978-1995.This is in good agreement with our data.These observations increase confidence in the observed patterns of downwasting, despite the remaining uncertainties.
The pattern of surface lowering on Khumbu Glacier can be explained in terms of ice dynamics and surface melt rates.Sustained high rates of ice delivery below the icefall largely offset melt in the upper ablation area, where debris cover is thin or absent (Figs. 5 and 6a, section A).Further downglacier, thin debris cover increases the ice melt (section B), in line with field measurements of increased surface lowering (Takeuchi et al., 2000).Thinning rates remain high downglacier despite an increasing debris thickness, due to very low glacier velocities and ablation associated with supraglacial lakes and exposed ice cliffs (Sakai et al., 2000(Sakai et al., , 2002) ) (section C, Supplement, Fig. 2).Almost no thinning was observed within 1 km of the terminus (section D), which may reflect either a thick, complete debris cover or indicate that ice loss is already complete.The possible slight surface lowering in the accumulation area of the glacier from 1970 to 2007 might be due to less snowfall.This is consistent with an ice core record at the East Rongbuk Glacier north of Mt.Everest that indicates decreasing snow accumulation for 1970-2001 (Kaspari et al., 2008).
The pattern of downwasting has resulted in the development of a concave glacier profile on Khumbu Glacier, with very low surface gradients in the lower ablation zone (Fig. 5b).This trend indicates that a glacial lake could develop about 1.5 to 3 km upstream of the terminus, as was predicted in simulations using a 1D-coupled mass balance and flow model by Naito et al. (2000).The formation of a large lake on Khumbu Glacier would have major consequences for outburst flood hazards downstream.
The greatest amount of mass loss occurred on Imja Glacier, which can be at least partly attributed to the proglacial Imja Lake, which enhances ice losses by calving.This lake grew significantly since its formation in the late 1960s, reaching ∼0.9 km 2 in 2008 (Supplement, Fig. 3, Fujita et al., 2009;Bolch et al., 2008a).The comparatively thin debris cover of Imja Glacier, apparent in exposed ice cliffs, is likely to be another reason for the higher mass loss.At Imja Glacier, a slight thinning is also observed on the relict ice-cored moraine situated below Imja Lake.This is in line with recent field measurements in this area which revealed a lowering of about 1 m a −1 (Fujita et al., 2009).The mass loss of the smaller glaciers, such as Amphu Laptse Glacier, amounts to half of that of Imja Glacier but is still significant.
Khumbu Glacier, for which data from the year 1984 was available, showed a higher surface lowering for the period 1984-2007 in comparison to 1962-1984.The recent trend of more negative mass balances since 2002, however, needs further investigation, as it is not statistically significant.Accelerated thinning could reflect decreasing velocity (Quincey et al. 2009), higher air temperatures (Prasad et al., 2009), decreasing snow accumulation (Kaspari et al., 2008), or a combination of these factors.The influence of black carbon (BC) as summarized by Ramanathan and Carmichael (2008) cannot be excluded but is negligible for the ablation zones with debris-cover as BC does not lead into changes in the albedo there.The catchment topography plays an important role for the glacier flow regimes in Khumbu Himalaya (Quincey et al., 2009).For example, Kangshung Glacier with a large high altitude accumulation area flowing from Mt. Everest towards the east, shows flow activity across the entire snout while all other avalanche fed glaciers and glaciers with a lower accumulation area such as Khumbu Glacier contain large stagnant parts.Our result shows comparatively little scatter of the mass balance amongst the investigated, mainly avalanche fed glaciers in the same study region.Mass balance estimate of Kangshung Glacier would reveal how the catchment topography and the higher flow influence the mass balance.Unfortunately, only the upper part of Kangshung Glacier is covered by our multi-temporal DTMs making a mass balance estimate impossible.
The total mass balance of Khumbu Glacier for 1970-2007, at −0.27 ± 0.08 m w.e. a −1 is lower than that observed on several other Himalayan glaciers including Chhota Shigri Glacier (−0.98 m w.e. a −1 , 2002-2006 Wagnon et al.,  2007 and -1.02 to -1.12 m w.e. a −1 , 1999-2004 Berthier et  al., 2007) or the small debris-free Glacier AX010 (−0.6 to −0.8 m w.e. a −1 , 1978-1999 Fujita et al., 2001), but similar to Dokriani Glacier (−0.32 m w.e. a −1 , 1992-2000 Dobhal  et al., 2008).However, the different observation times and glacier sizes have to be considered, and Khumbu Glacier also has possibly a more negative mass balance in recent years.The tendency towards increased mass loss has also been observed worldwide and for the few other Himalayan glaciers with mass balance estimates (Cogley, 2011).The mass loss of the investigated glaciers is similar to the average mass loss of the 30 reference glaciers worldwide for 1976-2005 (−0.32 m w.e. a −1 ) (Zemp et al., 2009).

Conclusions
This study presents the longest time series of geodetically derived mass-balance and glacier volume change estimates obtained to date in the Himalaya.Geodetic mass-balance estimates based on early stereo Corona and recent satellite data are suitable for tracking glacier changes through time, thus filling major gaps in glaciological knowledge of the Himalaya and other mountain regions.However, careful relative adjustments of the DTMs are necessary to obtain suitable accuracies of DTMs based on different data sources with different resolutions.Mass balance information is urgently needed to improve knowledge of the response of Himalayan glaciers to climate change and to allow prediction of future glacier change and its influence on water resources, river runoff, sea level rise, and glacial hazards.
Glaciers south of Mt.Everest have continuously lost mass from 1970 until 2007, with a possibly increasing rate in recent years.All glaciers lost mass despite partly thick debriscover.The highest loss was observed at Imja Glacier which terminates into a lake.The specific mass balance of the investigated glaciers of −0.32 ± 0.08 m w.e. a −1 is not higher than the global average.

Fig. 1 .
Fig. 1.Study area; location, names and debris-covered portion of the glaciers in the study area, coverage of the utilized satellite data.Background: SRTM3 CGIAR, Vers.4, study area: ASTER DTM; glacier outlines based on Bolch et al. (2008b).

Fig. 5 .
Fig. 5. (A): Khumbu Glacier based on the 2007 Cartosat-1 image including the profile and the different sections of surface lowering (see Fig. 6A).(B, C): Terminus of Khumbu Glacier based on the 1970 Corona KH-4B and on the 2007 Cartosat-1 image.

Fig. 6 .
Fig. 6. (A): Profiles of the DTM differences of Khumbu Glacier.(B): Longitudinal profiles of the surface elevation of Khumbu Glacier 1970 and 2007.See Figs. 3 and 5 for the location of the profiles.

Table 1 .
Utilized imagery and derived DTM characteristics based on the ice free area relative to the 2007 master DTM.

Table 2 .
Statistics of the DTM differences for the investigated periods.

Table 3 .
Fujita et al. (2009)and mass balance 1970Glacier volume loss and mass balance  -2007Glacier volume loss and mass balance  , and 2002Glacier volume loss and mass balance  -2007..For the year 2007, the estimated volume of Imja lake (37.8 × 10 6 m 3 , calculated based on the area extent 0.91 km 2 and the average depth of 41.4 m) was added to the elevation difference of1970 -2007 .For 2002 -2007we added the volume difference between 2007 and 2002 (2.0 × 10 6 m 3 ).Data is based onFujita et al. (2009).
1We assumed an ELA of 5700 m based on Ashai (2001) and interpretation of the satellite images. 2

Table 4 .
Volume loss of the ablation area of Khumbu Glacier1962-2007.