Heterogeneity in glacier response in the upper Shyok valley , northeast

Glaciers in the Karakoram show long-term irregular behaviour with comparatively frequent and sudden advances. A glacier inventory of the upper Shyok valley situated in northeast Karakoram has been generated for the year 2002 using Landsat ETM + and SRTM3 DEM as baseline data for the investigations and subsequent change analysis. The upper Shyok valley contained 2123 glaciers (larger than 0.02 km 2 in size) with an area of 2977.9± 95.3 km2 in 2002. Out of these, 18 glaciers with an area of 1004.1 ± 32.1 km2 showed surge-type behaviour. Change analysis based on Hexagon KH-9 (years 1973 and 1974) and Landsat TM/ETM+ (years 1989, 2002 and 2011) images had to be restricted to a subset of 136 glaciers (covering an area of 1609.7 ± 51.5 km2 in 2002) due to adverse snow conditions. The area of the investigated glaciers, including the 18 surge-type glaciers identified, showed no significant changes during all studied periods. However, the analysis provides a hint that the overall glacier area slightly decreased until about 1989 (area 1973: 1613.6 ± 43.6 km2; area 1989: 1602.0 ± 33.6 km2) followed by an increase (area 2002: 1609.7± 51.5; area 2011: 1615.8 ± 35.5 km2). Although the overall change in area is insignificant, advances in glacier tongues since the end of the 1980s are clearly visible. Detailed estimations of length changes for individual glaciers since the 1970s and for Central Rimo Glacier since the 1930s confirm the irregular retreat and advance.


Introduction
Meltwater discharge from Himalayan glaciers can play a significant role in the livelihood of people living in the downstream areas.Qualitatively, glacier meltwater in the Western Himalaya and Karakoram region is less influenced by the summer monsoon compared to the central and eastern parts of the mountain chain (Immerzeel et al., 2010;Bolch et al., 2012).Recent studies revealed that most parts of northwestern Himalaya showed less glacier shrinkage than the eastern parts of the mountain range during the last decades (Bhambri and Bolch, 2009;Bolch et al., 2012;Kääb et al., 2012).Conversely, glaciers in the western and central Karakoram region showed long-term irregular behaviour with frequent advances and possible slight mass gain since the 2000s (Hewitt, 2011;Copland et al., 2011;Bolch et al., 2012;Gardelle et al., 2012Gardelle et al., , 2013;;Kääb et al., 2012;Minora et al., 2013).Moreover, individual glacier advances in eastern Karakoram have also been reported in the Shyok valley during the last decade (Raina and Srivastva, 2008).These individual advances and mass gain episodes could be attributed to surging (Barrand and Murray, 2006;Hewitt, 2011;Copland et al., 2011;Quincey et al., 2011), winter temperature decrease (Shekhar et al., 2010) and increased solid precipitation in the accumulation areas (Fowler and Archer, 2006).Copland et al. (2011) reported an increase of glacier surge activities after 1990 in the western and central Karakoram region.There are few published studies on surging phenomenon on individual glaciers such as Rimo, Chong Kumdan, Kichik Kumdan and Aktash glaciers in the eastern

R. Bhambri et al.: Heterogeneity in glacier response in the upper Shyok valley
Karakoram (Raina and Srivastva, 2008;Tangri et al., 2013).Periodic advances of Chong Kumdan Glacier blocked the flow of Shyok River on several occasions, resulting in the formation of an ice dam.Several hazardous situations in downstream areas due to sudden collapse of ice dams have been reported (Mason, 1930;Hewitt, 1982;Raina and Srivastava, 2008;Hewitt and Liu, 2010).These glacier lake outburst floods (GLOFs) in upper Shyok valley, eastern Karakoram have been documented since as early as 1533 (Hewitt and Liu, 2010).Hence, it is essential to frequently monitor Karakoram glaciers not only to understand the surging process but also to mitigate hazardous situations arising due to GLOFs.Conventional field surveys are laborious and can be dangerous in high mountainous regions, whereas multispectral and multi-temporal satellite data offer great potential to monitor these glaciers at regular intervals (e.g.yearly or decadal).Topographic mapping which was carried out based on ground survey for eastern Karakoram, including Rimo Glacier (Dainelli, 1932), allows us to extend the investigation back in time.Topographic information on glaciers (e.g.minimum and maximum elevation) are largely unknown and longer-term direct climate measurements are non-existent for this eastern part of the Karakoram.In addition, to date, this region is not yet covered by the Global Land Ice Measurements from Space (GLIMS) database (www.glims.org,Raup et al., 2007).
Therefore, the main objective of this study is to (1) generate a complete and up-to-date glacier inventory for the upper Shyok valley and to provide information on the general glacier characteristics for the GLIMS data base and the Randolph Glacier Inventory (RGI; Arendt et al., 2012); (2) analyse glacier changes in upper Shyok valley during the last 40 yr and beyond; and (3) identify surge-type glaciers in upper Shyok valley.

Study area
The Shyok valley is a part of the Karakoram and covers the eastern part of the upper Indus basin (Mason, 1938).Shyok River is a major tributary of Indus River, which originates from the snout (∼ 4950 m a.s.l.) of Rimo Glacier and meets Indus River at Keris (2275 m a.s.l.) near the Skardu (Kaul, 1998).The entire study area includes glaciers of the Rimo Muztagh, Saser Muztagh, Depsang plains and the Chang Chenmo valley (Fig. 1).The Rimo, Kumdan, Saser and Kunzang group of ridges form a divide between the Shyok and Nubra catchments.Our study area covers ∼ 14 200 km 2 and an elevation range of ∼ 3600-7600 m a.s.l.The geology, geomorphology and sub-climatic regimes differ between the Karakoram Batholith terrain of the Sasir-Rimo Mustagh, and the eastern tributary basins between the Pangong Suture Zone and Karakoram Pass (Searle, 1991).The study area is situated under winter snow regime (Thayyen and Gergan, 2010) but is probably influenced by monsoon activity dur-ing the summer (Wake, 1989).Glaciers of eastern Karakoram receive solid precipitation from westerly air masses originating from Mediterranean Sea and/or Atlantic Ocean (Shrivastava, 2006).Maximum snowfall occurs between November and April (55 % of total snowfall) and slightly less during the May to October (45 % of total snowfall) (Shrivastava, 2006).On the basis of cosmogenic radionuclide (CRN) dating, Dortch et al. (2010) have identified three glacial stages (Deshkit 1, Deshkit 2 and Deshkit 3 stages) in the Nubra and Shyok valley confluence, which are controlled by the mid-latitude westerlies and oscillation in the Northern Hemisphere ice sheets and oceans.

Data sources
Multispectral orthorectified Landsat TM/ETM+ (spatial resolution 30 m), panchromatic Hexagon KH-9 (spatial resolution ∼ 7 m) and panchromatic OrbView-3 (spatial resolution 1 m) satellite data from 1973 to 2011 were acquired from United States Geological Survey (USGS, http://earthexplorer.usgs.gov/)(Table 1).The void-filled SRTM3 DEM from the Consortium for Spatial Information -Consultative Group for International Agriculture Research (CSI -CGIAR), version 4 (http://srtm.csi.cgiar.org/)(90 m spatial resolution) was used as a reference DEM for semiautomatic delineation of drainage basins and extraction of topographic glacier parameters.The ETM+ scene (year 2002) was selected as the reference image due to very little seasonal snow cover during the ablation period (Table 1).In addition, this year is closest to the year 2000, for which a global glacier inventory is recommended by Paul et al. (2009) for generation of a baseline data set to facilitate global glaciological applications.Since marginal areas of the studied valleys are not covered by the reference Landsat scene, an additional scene of September 2001 was used (Table 1).All the utilized Landsat TM/ETM images were available orthorectified in the processing level L1T, except the 1989 TM scene which was only processed in L1G.This scene had a slight horizontal shift of ∼ 30 m compared with the 2002 reference image.Hence, we co-registered the 1989 scene to the reference image using the projective transformation algorithm of Erdas Imagine.This technique allows the inclusion of DEM information from the SRTM3 DEM to enhance rectification as compared to the simple polynomial algorithm.The image geometry of Hexagon KH-9 (1973, 1974) is complex (Surazakov and Aizen, 2010;Pieczonka et al., 2013) and a proper co-registration is crucial.The Hexagon KH-9 image data is large in size (∼ 2 GB) and the KH-9 scenes were divided in 8 parts to optimize the processing speed.Each part of Hexagon KH-9 was co-registered based on ∼ 50 ground control points (GCPs) derived from the 2002 Landsat ETM+ imagery by spline adjustment using ESRI ArcGIS 9.3.(cf.Bhambri et al., 2011b).The average horizontal shift between reference 2002 Landsat ETM+ image and corresponding subsets of Hexagon (1973Hexagon ( , 1974) ) images was found to be ∼ 20 m.The topographic map of eastern Karakoram for the year 1930 (Dainelli, 1932; Fig. 2) has been used and coregistered to the 2002 reference image by spline adjustment.The glacier map is too rough for detailed change assessment but is suitable to estimate the length change of Central Rimo Glacier.

Glacier mapping, inventory and changes
Debris-free glaciers were mapped by the well-established semi-automated TM3/TM5 (visible band red [RED]/short wave infra-red band [SWIR]) band ratio approach followed by a 3 by 3 median filter to eliminate isolated pixels (Paul and Kääb, 2005;Racoviteanu et al., 2009;Bolch et al., 2010a;Bhambri et al., 2011a).The resultant glacier polygons were visually checked to eliminate misclassified pro-glacial lakes, seasonal snow cover, rocky outcrops, and shadow areas and were subsequently manually improved.The precise identification of the few exiting debris-covered glacier termini was rather difficult using Landsat TM/ETM+ imagery alone.Therefore, high-resolution panchromatic OrbView 3 images were used as an additional source of information.Several indicators such as signs of movement (identified with overlays of multi-temporal images) and streams at the end of the ter-minus helped to determine the most likely position of the glacier termini.
The contiguous ice masses were separated into their drainage basins using an automated approach for hydrological divides as proposed by Bolch et al. (2010a) and subsequently visually checked and improved where required.It is assumed that the ice divides did not change over the study period.This approach also minimizes the chance of errors that may occur due to different delineation of the upper glacier boundary, which could be due to varying snow conditions.The minimum size of mapped glaciers to be included in the inventory was set to 0.02 km 2 according to Frey et al. (2012).The Hexagon images and Landsat TM scenes for the years 1989 and 2011 were partially suitable for glacier mapping due to fresh snow cover.Therefore only 136 glaciers, ranging in size from 0.2 to 270 km 2 , could be compared for spatial area variation.For this purpose, the 2002 outlines were manually adjusted using the 1973 Hexagon and the 1989 and 2011 TM images.The adjustments were restricted to the lower part of the glaciers as no visible changes could be identified in the accumulation region, and snow cover hampered correct identification at some high-elevation areas.
Four selected glaciers (Central Rimo and Kumdan glaciers: Chong Kumdan, Kichik Kumdan and Aktash, Fig. 1) were studied in detail.The selection of glaciers is based on previous studies (Mason, 1930;Hewitt, 1982;Raina and Srivastava, 2008;Hewitt and Liu, 2010) that have documented surging activities during the last two centuries.This has enabled us to identify the surge cycles, taking into account the Kumdan glaciers which blocked the Shyok River several times and created hazardous situations for downstream areas.In addition, Kichik Kumdan and Aktash glaciers were analysed on a yearly scale based on Landsat satellite images from 1998 to 2011.
For the calculation of length changes along the termini of Kichik Kumdan, Aktash and Central Remo glaciers, stripes with 50 m spacing were drawn parallel to the main flow direction of the glacier.Length change was calculated as the average length from the intersection of the stripes with the glacier outlines (cf.Koblet et al., 2010).Several glaciers showed surge-type behaviour and were visually identified on multi-temporal satellite images (called "surge-type glaciers" hereafter).Indications were typical folded moraines and heavily crevassed surfaces at the glacier front and/or rapid frontal advances in a short period.The influence of topographic parameters (minimum, maximum, median, mean, aspect and range elevation) on surge-type and non-surge-type glacier area change is also examined, as these parameters could be important factors to explain glacier area change (Paul et al., 2009).
Our study includes various data sources at different spatial and temporal resolutions.Hence, error assessment is crucial to ascertain the accuracy and significance of the results.The potential uncertainty of the multi-temporal analysis generally results from location and mapping errors.Glacier mapping uncertainty mainly depends on the spatial and radiometric resolution of the used image and the conditions at the time of the acquisition (e.g.seasonal snow).In perfect conditions some studies achieved less than half-pixel accuracy (e.g.Bolch et al., 2010b).We used 2002 ETM+ scene with a 15 m panchromatic band as the basis of the inventory and estimated the mapping uncertainty for each glacier based on a buffer around the glacier margins as suggested by Granshaw and Fountain (2006) and Bolch et al. (2010b).The buffer size was set to 7.5 m for the ETM+ and we assume that the uncertainty is similar for the adjusted outlines of the TM scenes.We chose a slightly higher buffer size of 10m for the Hexagon images taking the co-registration inaccuracy into account.The average mapping uncertainty was 3.2 % for the 2002 ETM+ image, 2.1 % and 2.2 % for the 1989 TM and 2011 TM images, respectively and 2.7 % for the 1973/1974 Hexagon images.It is to be noted that there are fewer (136) and, on average, larger glaciers for the area change analysis, causing the lower uncertainty than for the 2002 data.For evaluation of these more theoretical measures of uncertainty we also performed a comparison between glacier outlines derived from the multispectral bands with 30 m resolution of a 2006 Landsat ETM+ image and those derived from a highresolution OrbView-3 image of the same year for a sample of four selected glaciers as suggested by Paul et al. (2013).The resultant uncertainty was ±2.7 % which is within the range of the above mentioned uncertainty and uncertainties reported by previous studies (Bolch et al., 2010a, b;Paul et al., 2013).
The area change uncertainty is estimated according to standard error propagation to be the root sum square of the uncertainty for the outlines mapped from different sources.With this approach, both the uncertainty of the upper glacier area, where the area changes are usually minor but could not be adjusted at some instance due to adverse snow con-ditions, and the mapping uncertainty of outlines of the individual years could be addressed.The location uncertainty for the common outlines is insignificant and well within the uncertainty calculated with the buffer method as we used the 2002 outlines as the basis and adjusted the glacier extents of the other years manually.
The uncertainty of the length change investigations is estimated by the following formula as suggested by Hall et al. (2003).
where, x1 is the pixel resolution of image 1, x2 the pixel resolution of image 2 and E reg represents the registration error.
Here, the registration error uncertainty for Hexagon scenes is considered to be 20 m and for the remaining Landsat scenes to be 5 m as these scenes almost perfectly matched with each other based on visual checks along the glacier terminus.
The uncertainty for the 1930 map is significantly higher and depends mainly on the scale and the mapping accuracy of the glacier outline.Unfortunately, there is no information about the mapping accuracy given in Dainelli (1932).Dainelli visited the glacier and used the existing Survey of India topographic map as a basis for his 1 : 750 000 scale map.For the present study, it is assumed that the uncertainty of the position of the glacier tongue is within the line thickness, which is ∼ 0.4 mm and represents 480 m on the ground for the utilized scale.The co-registration accuracy is ∼ 100 m and hence the overall uncertainty of the 1930-1973/1974 period is estimated to be ∼ 600 m.

Glacier characteristics
The 2002 inventory of the upper Shyok basin includes 2123 glaciers larger than 0.02 km 2 , with a total area of 2977.9 ± 95.3 km 2 .Out of these, glaciers smaller than 1.0 km 2 comprise 79 % (1683) of the total number, but encompass only 13 % of the total area (Fig. 3).This conforms to the pattern in other regions of Himalaya (Raina and Srivastava, 2008;Bajracharya and Shrestha, 2011;Frey et al., 2012), Tibet, Western Canada or the Alps (e.g.Kääb et al., 2002;Bolch et al., 2010a, b).The mean size of glaciers in the study area is 1.4 km 2 which is slightly higher than the glacierized basins of Ganga (1.1 km 2 ) and Brahmaputra (1.2 km 2 ) (Bajracharya and Shrestha, 2011).Eighteen (18) surge-type glaciers have been identified with a total area of 1004.1 ± 32.1 km 2 covering ∼ 34 % of the glacierized area of the upper Shyok valley.Hence, the average size of the surgetype glaciers is significantly larger (∼ 56 km 2 ) than the average size of the non-surge-type glaciers (∼ 1 km 2 ).Five surgetype glaciers (total area: 64.3 ± 2.0 km 2 ) are located east of the Shyok main stream and the remaining thirteen glaciers in Rimo and Saser Muztagh ranges (Fig. 1).The mean glacier elevation ranges from 5275 to 6456 m a.s.l., with an average elevation of 5830 m a.s.l. in the study area.This elevation is higher than in western Himalayan basins like Jhelum (4432 m a.s.l.) and Ravi (4761 m a.s.l.) (Frey et al., 2012) and Garhwal Himalaya such as upper Bhagirathi (5629 m a.s.l.) and Saraswati/Alaknanda (5490 m a.s.l.) (Bhambri et al., 2011b).It exhibits to follow the pattern of the central and western Himalayan regions, where the mean glacier elevation is increasing from southwest to northeast (Frey et al., 2012).This could be due to the increase in precipitation from the Indus-Ganges lowlands to the Himalayan mountains and a decrease towards the Trans-Himalaya (Bookhagen and Burbank, 2006;Bolch et al., 2012).The glaciers were divided into five size classes namely, < 1, 1-5, 5-10, 10-50 and > 50 km 2 .The size class > 50 km 2 contains only 6 glaciers covering an area of 833.0 ± 26.6 km 2 .All these large glaciers (e.g.Central and South Rimo, North and South Shukpa Kunchang, Chong Kumdan and Kichik Kumdan glaciers) are located in the western part of study area (Rimo and Saser Muztagh ranges).
Here, the glaciers have a higher average elevation range of 565 m as compared to glaciers of the eastern part (Depsang Plains and Chang Chenmo region: 393 m).
Thirty-five (35) glaciers have been identified that are partly covered with debris in their ablation zones.However, in most cases the debris is restricted to medial moraines.The overall area covered by debris was found to be 37.0 ± 1.1 km 2 (∼ 1.2 % of the whole ice cover) in entire study area.Most of the glaciers with debris cover (91 % of the total debriscovered ice area) are located in the western part (Rimo and Saser Muztagh) of the upper Shyok valley and the remaining are to the east of main Shyok stream.Similar results were also reported by Frey et al. (2012) in western Himalaya where the debris-covered ice area decreases from southwest to northeast.

Glacier area and length changes
On average, the glacier area changes are within the uncertainty range.However, some individual glaciers show significant changes (e.g.Kumdan group of glaciers).Most of glacier area changes occur below the equilibrium line altitude (ELA) and are clearly visible at the frontal part of terminus (Tables 2, 3, Figs. 2, 4, 6).Between 1973 and 1989, the overall glacier area probably decreased slightly (−0.7 ± 3.4 %) and also the overall area of the surge-type glaciers showed indications of a shrinkage in a similar range (−0.8 ± 3.4 %).
In the period from 1989 to 2002, all glaciers exhibited on average signs of area increase (+0.5 ± 3.8 %).This increase was mainly due to the surge-type glaciers as the non-surgetype glaciers were more or less stable (−0.1 ± 3.8 %).During the last decade (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011), also the non-surge-type glaciers showed indications of a slight area increase.Overall, ∼ 75 % of the total area change of +0.4 ± 3.9 % can be attributed to surging glaciers during 2002-2011.This demonstrates that the total absolute glacier area changes in the study area were dominated by the surge-type glaciers.

Glacier inventory and changes
Overall, 2123 glaciers with a total area of 2977.9 ± 95.3 km 2 in the upper Shyok valley were identified for the year 2002.However, Geological Survey of India (GSI) reported only 1263 glaciers with an area of 2608 km 2 in the same region (Sangewar and Shukla, 2009).To investigate the reason for the difference, a small-scale map of the GSI glacier inventory was scanned and overlaid on the outlines based on the 2002 satellite image.This highlights that the difference in glacier numbers can be attributed to the inclusion of smaller glaciers with a minimum threshold of 0.02 km 2 in the present glacier inventory, but the omission of smaller ice bodies in the case of the GSI inventory.The GSI glacier inventory for the Shyok valley is based on 1 : 50 000 scale Survey of India (SOI) topographic maps (which are mostly from the 1960s), aided by satellite images and aerial photographs with limited field checks (Raina and Srivastava, 2008;Sangewar and Shukla, 2009).However, these SOI topographic maps do not show all ice bodies due to scale limitations, and some glacier outlines might have been misinterpreted from aerial photographs due to seasonal snow and debris cover (Raina and Srivastava, 2008;Bhambri and Bolch, 2009).In addition, contiguous ice masses may be counted as a single entity or could be subdivided into multiple glaciers depending on the purpose of the inventory.
The mean glacier elevation in the present study area is located at an elevation of about 5830 m a.s.l., whereas the mean glacier elevation of the entire Karakoram has a clearly lower elevation of about 5330 m a.s.l.(Bolch et al., 2012).Similarly, the lowest glacier terminus in the study area is located at an elevation of about 4600 m a.s.l.while the lowest glacier termini of the entire Karakoram, which is located in the Hunza valley, reaches down to about 2300 m a.s.l.(Hewitt, 2011).This shows that glaciers are situated at higher elevations in the eastern Karakoram, which is located in the rain shadow of both the monsoon and the westerly precipitation as compared to western Karakoram (cf.Bolch et al., 2012).
All the glaciers in the study area are almost debris free (2088 count; debris-free area: 2941 km 2 ) and only about 1 % of the glacier area (35 glaciers; area: 36.9 km 2 ) is covered with debris.This is likely due to local lithology (Karakorum plutonic complex) of rocks in the Shyok basin (Searle, 1991) and the fact that most of the headwalls are covered with ice and are broader and not so steep as compared to regions with a high percentage of debris cover (for example Garhwal or Khumbu Himalaya, Bhambri et al., 2011b;Bolch et al., 2008).In addition, a negative mass budget and area loss favour the increase of debris cover on glaciers (Bolch et al., 2008;Stokes et al., 2007) but the investigated glaciers in the Shyok valley are more stable.
Complex image geometry of the Hexagon images is possibly the main source of uncertainty in the estimate for glacier change.However, Hexagon images from the 1970s are very helpful to identify surge events and provide insight into longterm behaviour of glaciers in remote and inaccessible areas in the absence of historic (since 1970s) aerial photographs and are useful to evaluate inventories based on topographic maps (Bhambri and Bolch, 2009;Bolch et al., 2010b).Further, uncertainty may arise due to the fact that the upper  glacier boundaries were considered constant in the present study and the adjustments were restricted to the lower part of the glaciers.This approach can be justified because the area changes in the upper glacier parts are usually minor and occur only when there is a significant thinning along with a clear negative mass budget (Pelto, 2010).However, in the case of balanced budgets as true for most of the glaciers in Karakoram (Gardelle et al., 2012(Gardelle et al., , 2013;;Kääb et al., 2012) insignificant or only very little area changes occur in the accumulation area (Schwitter and Raymond, 1993;Pelto, 2010).In addition, the visible change in the higher elevation area is mainly due to variation in the seasonal snow cover, which makes the identification of real change in the highelevation regions impossible.This is a common problem with multi-temporal glacier inventories.Hence, keeping the upper glacier boundaries constant will result in less uncertainty than including variation that are likely caused by varying snow conditions on the available satellite scenes.Therefore the uncertainty was treated for each period separately using a buffer around the entire glacier in the present investigation, and was not restricted to the glacier parts where the visible changes occur.
The results confirm that there are evidences of stability or slight glacier area gain during the 2000s also in eastern Karakoram as found for other parts of the mountain range (Scherler et al., 2011;Minora et al., 2013).The gain in glacier area is statistically insignificant based on our uncertainty assessments.However, a growth is clearly visible for many glaciers on the satellite images (Fig. 4).Our study indicates that the area gain has not only occurred during the first decade of the 21st century but also the last decade of the 20th century.Unfortunately, the change analysis is restricted to a subsample of the glaciers due to unfavourable images.It is therefore not certain whether a possible area gain would also hold true for the entire region.However, it is unlikely that the inclusion of all glaciers will lead to a clear area loss.About 70 % of the glacier area increase between 2002 and 2011 can be attributed to surge-type glaciers,  while the remaining area increase is associated with slow advancement of non-surge glaciers.This slow non-surge glacier area gain may be a result of cumulative positive mass balance over several decades due to high precipitation in the Karakoram region (Fowler and Archer, 2006).Recently, Gardelle et al. (2012) and Kääb et al. (2012) reported balanced budgets or a possible slight mass gain in the central Karakoram for the first decade of the 21st century.Positive mass balance or balanced budgets were also measured in the Pamir north of our study region (Yao et al., 2012;Gardelle et al., 2013).Thus there is a need to carry out further mass balance studies to understand impact of climate change in the upper Shyok valley as well as the greater Karakoram region.Previous studies in Garhwal Himalaya (Dobhal et al., 2008;Bhambri et al., 2011b;Mehta et al., 2011) and Himachal Himalaya (Kulkarni et al., 2007) have reported an increase in the number of glaciers during the last few decades due to disintegration, whereas this study found that some glaciers coalesced after surging processes in the eastern Karakoram region.These findings confirm that the greater Karakoram

Surging glaciers in Shyok Valley
The present investigation results show heterogeneous patterns of glacier advances in upper Shyok valley (Tables 2  and 3).Eighteen (18) surge-type glaciers were identified, which cover ∼ 34 % of the entire glacierized area of upper Shyok valley.Hence, the advances of surge-type glaciers influence the glacier area statistics.For instance, the possible slight total glacier area increase between 1989 and 2002 can be attributed to the surge-type glaciers, whereas on average non-surge-type glaciers may have shrunk slightly during this period (Table 3).Out of all surge-type glaciers, only four glaciers (Rimo, Chong Kumdan, Kichik Kumdan and Aktash) in the eastern Karakoram have been investigated by previous studies using a limited number of satellite images (from 1990(from to 2009(from , e.g. Tangri et al., 2013)).These four glaciers cover 65 % of the area of the total 18 surge-type glaciers.The present study confirms the fluctuations of these glaciers.Surge-type glaciers flow at a substantially higher rate during a short period (varying from 1 to 15 yr; active phase) and with a long period of stagnation (10-100 yr; quiescent phase) in between (Jiskoot, 2011).Kichik Kumdan and Aktash actively surged for two and three years, respectively.These glaciers thus have rather short active phase duration similar to previously studied western Karakoram surge-type glaciers (Copland et al., 2011).For instance, Chiring Glacier, tributary of the main Panmah Glacier, Karakoram actively surged for three years (1994)(1995)(1996) (Jiskoot, 2011).However, the duration of the active phase is significantly longer (3 to 10 yr) for Svalbard glaciers than for surge-type glaciers in Karakoram (Dowdeswell et al., 1991).Almost all the surge-type glaciers showed advancements during the study period.Only one glacier has been identified as surge-type which continuously retreated during the investigated period.This glacier showed typical folded moraines and heavily crevassed surfaces at the glacier front in the 1973 Hexagon image so we suspect that the quiescent phase of this glacier is probably longer than the present period of investigation.Barrand and Murray (2006) found that surge-type glaciers are more likely to be longer and larger in the Karakorum region.This is in line with our study and only two surge-type glaciers have been identified which are smaller than 5 km 2 .
Investigations of long-term glacier changes in the Karakoram helps us to understand the complete picture of glacier fluctuations.The finding of a retreat of Rimo Glacier between 1930 and 1989 confirms previous observations by Hewitt (2005): most of the glaciers of the central and western Karakoram strongly diminished from the 1920s to the early 1990s except for some short-term advances and surges in the 1970s.Grant and Mason (1940) predicted that Chong Kumdan Glacier was likely to obstruct Shyok River and form a large lake in 1970s, and could again advance by 2013.Fortunately, this glacier did not block the Shyok River during the period of the present investigation.Five surge tributaries have also been identified out of the 18 surge-type glaciers; those detached from or joined to the main glaciers (e.g.North Shukpa Kunchang) during the study period.Similarly, Hewitt (2007) also noted 12 tributaries which have been separated from or joined to main glaciers in the central Karakoram since the 1860s that are not known to surge.
Several noteworthy studies have reported periodic surging during the last two centuries in central and western Karako-ram while only three glaciers (Kumdan glaciers) of the eastern Karakoram were investigated in detail till date (Mason, 1930;Hewitt, 1982Hewitt, , 2005Hewitt, , 2011;;Raina and Srivastava, 2008).In the present investigation, the existing published data of Kumdan glaciers were confirmed and extended and, in addition, five surge-type glaciers east of main stream of Shyok River were identified (Fig. 1).Surge events of Kumdan glaciers have affected the downstream areas during last two centuries (Mason, 1930;Raina and Srivastava, 2008;Hewitt and Liu, 2010).The lakes formed due to the blockage of Shyok River, periodically burst and discharged huge amounts of debris-laden water to downstream areas (Mason, 1929;Mason et al., 1930).Over the last 200 yr, 62 glaciers partially or completely blocked Karakoram rivers (Hewitt, 2011).Hewitt and Liu (2010) reported historic events of about 90 ice-dammed lakes associated to dozens of glaciers and outburst floods in the entire Karakorum, but no flood event was reported in upper Shyok valley after the 1930s.However, the studied Aktash Glacier touched the Shyok River during 2002/2003 and crossed the river in 2009 (Fig. 4).Fortunately the river tunneled the glacier front from its bottom and continued to flow without any blockage.Copland et al. (2011) reported 90 surge events in the entire Karakoram from 1960s onwards.During the present investigation, only 3 glaciers surge events during 1973-1989 and 15 surge events between 1989 and 2011 have been observed.This appears to be in cohesion with the previous observation that the number of new surges almost doubled after 1990 (Copland et al., 2011), which could be due to hydrological (Quincey et al., 2011) and/or thermal controls (Copland et al., 2011).

Conclusions
Our inventory of 2123 glaciers of the upper Shyok valley in eastern Karakoram (covering an area of 2977.9 ± 95.3 km 2 ) and detailed change analysis for the period 1973-2011 using Hexagon KH-9 and Landsat imagery will not only fill gaps in the GLIMS database but also support further detailed studies related to volumetric changes, surging mechanism and lake outburst modelling of ice-dammed lakes.Fresh snow and cloud cover hampered the generation of multi-temporal glacier inventories for the entire region and only a subset of 136 glaciers covering an area of about 1609.7 ± 51.5 (2002) could be studied in detail.The analysis indicates that on average, a slight area loss between 1973 and 1989 was followed by a slight gain.Overall, a possible slight net area gain was found (area 1973: 1613.6 ± 43. in upper Shyok valley is similar to other western Karakoram glaciers.However, the complex glacier surge mechanism and possible area and mass gain need to be further investigated by field investigations, geodetic mass budget estimations and glacier flow studies, in order to better understand the response of the Karakoram glaciers to climate change.

Fig. 1 .
Fig. 1.Location of upper Shyok valley, northeast Karakoram with four glaciers studied in detail.Glacier outlines are based on 2002 August Landsat ETM+ imagery.The colour-coded elevation is based on SRTM3 data.

Fig. 2 .
Fig. 2. Rimo Glacier, northeast Karakoram on different time spans; (A) 1930 extent based on Dainelli (1932) map and (B) 2002 situation based on Landsat ETM+ (2 August 2002).Map published with the permission of the Royal Geographical Society (with IBG).See Fig. 1 for location in study area.

Fig. 3 .
Fig. 3. Diagram showing the glacier size class vs. number of glaciers in upper Shyok valley.
while more recently (2002-2011) the glacier slightly reduced its area (−1.3 ± 3.7 % or −0.14 ± 0.4 % a −1 , Table 4).The central and south branch of Rimo Glacier were connected in the 1930 map but were separated in the 1974 image.The central branch receded −1300 ± 600 m during 1930-1974 and further by −800 ± 56 m until 1989.In contrast, this glacier exhibited an advance of +800 ± 52 m thereafter (1989-1999) while the glacier retreated again until 2011 by −400 ± 52 m (Fig. 2).No satellite image was available between 1989 and 1998.Therefore it is difficult to determine if the advancement of Central Rimo Glacier was gradual or sudden between 1989 to 1998.Kichik Kumdan Glacier retreated −1049 ± 56 m during 1974 to 1998, while it illustrated a sudden advance of ∼ 975 ± 52 m in the one-year period from 1998 to 1999 (Figs. 4 and 5).The Aktash Glacier had almost the same extent during 1974-1989 and terminated ∼ 200 m away from the Shyok River.This glacier touched the river during 2002/2003 and crossed it in 2009 (Fig. 4).An analysis of satellite data of 1973, 1989 and 1998 to 2011 (annual satellite data, 14 scenes) suggest that the left limb of Chong Kumdan Glacier began surging in 2002 and intersected with the right limb of the glacier between 2005 and 2006.

Fig. 5 .Fig. 6 .
Fig. 5. Cumulative length records of Central Rimo and Kichik Kumdan and Aktash glaciers during 1973-2011.Dotted line for Central Rimo Glacier indicates that there were no observations between 1930 and 1973 in our study.

Fig. 6 .
Fig. 6.Example of glacier changes in part of upper Shyok valley from 1973-2011.Background based on Landsat ETM+ (2 August 2002).See Fig. 1 for location in study area.
6 km 2 ; area 2011: 1615.8 ± 35.5 km 2 ) for the glaciers in the upper Shyok valley since 1973.Individual glaciers in the upper Shyok valley exhibited significant area changes and clearly visible fluctuations at their termini.Detailed investigations for length changes of a few surge-type glaciers revealed irregular behaviour as of 1973 (and for Rimo Glacier as of 1930).The heterogeneous response of glaciers to climate change www.the-cryosphere.net/7etal.: Heterogeneity in glacier response in the upper Shyok valley

Table 1 .
Details of the USGS satellite data used in the present study.

Table 2 .
Area (km 2 ) and number changes of surging and non-surging glaciers in upperShyok valley during 1973Shyok valley during  -2011.   .

Table 3 .
Changes in glacier area in the study areabetween 1973, 1989, 2002 and 2011based on Hexagon and Landsat images.

Table 4 .
Length and area fluctuation of Kichik Kumdan, Aktash and Rimo glaciers during study period.The central and south branch of Rimo Glacier were combined in 1930.