Area and mass changes of Siachen Glacier (East Karakoram)

Abstract

Here, we present a comprehensive assessment of Siachen Glacier (East Karakoram), in terms of its area and elevation change, velocity variations and mass budget, utilizing different satellite datasets including Landsat, Hexagon, Cartosat-I, Shuttle Radar Topography Mission, Envisat Advanced Synthetic Aperture Radar and Japanese Advanced Land Observing Satellite Phased Array-type L-band SAR. The total areal extent of Siachen Glacier did not change significantly between 1980 and 2014; however the exposed-ice area decreased during that period. The terminus of the glacier has experienced substantial downwasting (on average 30 m) over the period of 1999–2007, followed by a retreat of the transition between exposed and debris-covered ice by a distance of 1.3 km during the short span 2007–14. The spatial patterns of the elevation difference and velocity are heterogeneous over the large areal extent of Siachen Glacier. The average velocity of the entire glacier, as computed between 11 December 2008 and 26 January 2009, was 12.3 ± 0.4 cm d⁻¹, while those estimated separately for the accumulation and ablation regions were 9.7 ± 0.4 cm d⁻¹ and 20.4 ± 0.4 cm d⁻¹, respectively. The mass budget of Siachen Glacier is estimated to be –0.03 ± 0.21 m w.e. a⁻¹ for the period of 1999–2007.

1. INTRODUCTION

Glaciers in the Karakoram mountain range are an important source of freshwater for one of the most densely populated river basins (Indus Basin) in the world (Immerzeel and others, 2013). The volume of discharge generated by snow and glacial melt far exceeds the volume of the discharge generated in the downstream areas of the basin (Immerzeel and others, 2010). The status of glaciers in the Karakoram has a direct impact on some of the key components of the global water cycle, global mean sea-level rise (Gardelle and others, 2012), freshwater availability and the risk of Glacial Lake Outburst Floods and droughts. Besides, continuous monitoring of glaciers is crucial for identifying the regional impacts of changing climate (Bolch and others, 2012), particularly where long-term climatic measurements are rare.

Components of glacier mass budget have been monitored throughout the world for more than six decades (WGMS, 2008; Zemp and others, 2009). As a result, various measurement methods have evolved with their own sets of merits and demerits (Hoinkes, 1970). For example, field-based glaciological methods (Mayo and others, 1962; Kaser and others, 2002) are widely used but are extremely manpower intensive. Inevitably, field-based methods in the Himalaya are biased towards short-term assessment of certain glaciers, such as Chhota Shigri (Wagnon and others, 2007; Azam and others, 2014), which are easy to access and smaller in size. Similarly, the accumulation area method suggested by Kulkarni and others (2011) uses empirical, field-based relationships and can only be extrapolated to glaciers that lie within a certain climatic region (Kulkarni and others, 2004). Geodetic methods, comparing DEMs from at least two different times, have been extensively implemented for various regions around the globe, including Karakoram glaciers (Gardelle and others, 2012), Lahaul Spiti glaciers (Berthier and others, 2007), glaciers in Khumbu Himalaya (Bolch and others, 2011; Nuimura and others, 2012), the Tien Shan (Bolch, 2015; Pieczonka and Bolch, 2015) and in general for Pamir-Karakoram-Himalayan glaciers (Gardelle and others, 2013). In contrast, the hydrological mass budget method involves subtracting hydrological outputs (i.e. evaporation and runoff) from hydrological inputs (i.e. mass gain from snowfall) for a glacier in order to determine its mass budget. This method has been rarely used in the Himalaya since Bhutiyani (1999).

The Karakoram Range, comprising ~18 800 km² of glacier area (Bolch and others, 2012), accounts for nearly 3% of the total area of ice outside the ice sheets in Greenland and Antarctica (Cogley, 2012). Even though Karakoram glaciers have shown irregular behavior over longer time periods, they have retreated on average over most of the 20th century (from 1920 to 1990) (Hewitt, 2005). Further, earlier reported expansion of high-relief glaciers in the Karakoram (Hewitt, 2005) was later negated by Bhambri and others (2013) and Minora and others (2016), who reported heterogeneous glacier behavior but on average no significant change in the Upper Shyok Valley, NE Karakoram since 1970 and in Central Karakoram since 2000. Owing to the political sensitivity and logistic difficulties of the Karakoram region, studies of glacier dynamics (Copland and others, 2009; Scherler and others, 2011) and mass budget (Gardelle and others, 2012; Kääb and others, 2012) have been primarily based on satellite observations. Earlier studies of glacier dynamics in the Eastern Karakoram ranges have reported increases in surface velocity (Heid and Kääb, 2012) with average stable termini (Scherler and others, 2011; Bhambri and others, 2013). Also, a large number of surge-type glaciers have been reported in the Karakoram since the 1860s (Barrand and Murray, 2006; Copland and others, 2011). And perhaps most interestingly, surge activity has increased in recent years (e.g. Copland and others, 2011; Quincey and others, 2011; Rankl and others, 2014). Gardelle and others (2012), based on remotely sensed geodetic measurements over 5615 km² of glacierized area in the central Karakoram, computed a region-wide mass budget of 0.11 ± 0.22 m w.e. a⁻¹ for the period of 1999–2010. More specifically, the mass budget of surging and non-surging glaciers was reported to be 0.11 ± 0.31 m w.e. a⁻¹ and 0.10 ± 0.19 m w.e. a⁻¹, respectively. Thus, the Karakoram glaciers were either in equilibrium or gaining mass between 1999 and 2010 (Gardelle and others, 2012, 2013). Like many other Himalayan and Karakoram glaciers, Siachen Glacier has not been monitored adequately owing to remoteness and years of military conflict in the area. Bhutiyani (1999) made an independent study of hydrologic mass budget for Siachen Glacier during 1986–1991 and obtained a large negative mass budget of −0.51 m w.e. a⁻¹. However, Zaman and Liu (2015) made corrections to the catchment area of Siachen Glacier and found mass budget values with upper and lower bounds of +0.22 m w.e. a⁻¹ and −0.23 m w.e. a⁻¹, respectively. However, since the hydrological method involves differencing of hydrologic variables with large magnitudes such as runoff, evapotranspiration and precipitation over the entire catchment, it is prone to large relative errors making the method less reliable (Zaman and Liu, 2015). The study by Gardelle and others (2013) also covered Siachen Glacier. The authors reported a slightly positive mass budget of +0.14 ± 0.14 m w.e. a⁻¹ during the period 2000–10. This estimate is however uncertain because of incomplete coverage of the glacier.

To our knowledge, the present study presents the most comprehensive assessment of Siachen Glacier in terms of its mass budget, velocity and area. The specific aims of this study include: (1) analysis of changes in the area of Siachen Glacier (2) improving understanding of how surface mass budget affects the surface velocity pattern and thickness variations of the glacier; of particular interest, in this regard, are differences between the elevation changes of exposed and debris-covered ice (3) quantitative characterization of mass budget for Siachen Glacier utilizing high-precision, remote sensing-based geodetic techniques.

2. STUDY AREA

Siachen Glacier is the largest glacier in the entire Karakoram Range and extends between the latitudes of 35°10′N-35°42′N and longitudes 76°46′E-77°25′E (Fig. 1). It spans a length of ~74 km and its width varies between 1 and 8 km, covering an area of ~936 km² with an estimated mean thickness of ~300 m (Frey and others, 2014). This makes Siachen the largest glacier not only of the Karakoram but also of all of High Mountain Asia. This glacier is located near the political boundary between India and Pakistan, in the eastern part of the Karakoram Range.

Fig. 1. Location of Siachen Glacier. Inset shows complete area of Siachen Glacier. The red line shows the ELA of 5250 m.

Siachen Glacier lies within the altitude range of 3400–7300 m a.s.l. with an average elevation of 5500 m a.s.l. Westerlies are an important source of precipitation for the Karakoram including Siachen Glacier. Nearly two-thirds of the annual snowfall in these areas occurs under the influence of extra tropical cyclones known as western disturbances, primarily during the winter season (Bhutiyani, 1999).

3. DATA

3.1. Cartosat-I

The Cartosat-I (IRS-P5) sensor is ideally suited for DEM generation in glaciological terrains, offering a base to height ratio of 0.62, a radiometric resolution of 10 bits and a spatial resolution of 2.5 m (Tiwari and others, 2007; Bolch and others, 2011; Pieczonka and others, 2011). The time difference between acquisitions of the members of stereo image pairs is 52 s, ensuring minimal temporal aliasing, which is particularly important for glaciological applications. In this study, three Cartosat-I images (Table 1) with an overlapping area of 98 km² were used to get an almost complete coverage of Siachen Glacier. The Cartosat-I ortho-pair was radiometrically corrected and coarsely geolocated using the path–row referencing system, which gives a horizontal accuracy of 170–300 m (Gianinetto and Fassi, 2008; Pieczonka and others, 2011). Rational polynomial coefficients (RPCs), also provided for each image of the stereo pair, were used to define the mathematical relationship between pixel and ground coordinates more accurately (Titarov, 2008). Since no field study could be conducted and no topographical maps are available, we have processed the Cartosat- I DEM without ground control points (GCPs).

Table 1. Details of the optical image datasets used

Also included is the SRTM.

* Used for glacier mapping in places of non-coverage of Cartosat-I scenes.

3.2. Landsat

In this study we utilized Landsat L1T products, which provide improved radiometric, geometric (by incorporating GCPs) and topographic accuracy (by employing DEMs) and a gap-filled Landsat scene obtained on 12 August 2009 from www.landsat.org. For glacier mapping, nearly cloud free scenes (cloud cover <30%) acquired at the end of the ablation period (starting from July to August) with minimal seasonal snow were used (Table 1).

3.3. KH-9 hexagon mapping camera (MC)

The KH-9 program (code-named Hexagon) of the US military, operational from 1973 to 1980, was aimed at reconnaissance and map making. In 2011, KH-9 images of the frame-MC with a resolution of 20–30 ft (~6–9 m), a frame of 23 × 46 cm² and a focal length of 30.5 cm (Surazakov and Aizen, 2010; Padmanabha and others, 2014) were declassified. The KH-9 scene used in this study (Table 1) covers an area of ~161 × 241 km², is mostly cloud free and has good image contrast for glacier mapping.

3.4. Envisat Advanced Synthetic Aperture Radar (ASAR), Japanese Advanced Land Observing Satellite Phased Array-type L-band (ALOS PALSAR) and Shuttle Radar Topography Mission (SRTM)-DEM

The ASAR is onboard the ENVISAT satellite launched in 2003. The ASAR operates in the C-band frequency range with a wavelength of 5.66 cm. At this wavelength, cloud cover has no impact on the reflected radar signals (See Table 2 for ENVISAT ASAR scenes used). We have used Envisat ASAR images to delineate debris-covered tongue using coherence mapping.

Table 2. Details of Envisat ASAR images used as pairs for the computation of coherence images

* IM refers to Image Mode.

For generating velocity maps, we utilized wintertime ALOS PALSAR data, acquired on 11 December 2008 and 26 January 2009 in fine beam single (FBS) mode and summertime data acquired on 8 June 2007 and 24 July 2007 in fine beam dual (FBD) mode. All acquisitions were along orbit 523 (Table 3).

Table 3. List of ALOS PALSAR datasets used for generating velocity maps

SRTM data have been used as a reference DEM in a number of glaciological studies to generate elevation difference maps by comparing them with recent (Berthier and others, 2007; Paul and Haeberli, 2008) and historical DEMs (Surazakov and Aizen, 2006; Schiefer and others, 2007; Pieczonka and others, 2013) and with repeat-track ICESAT data (Kääb and others, 2012; Neckel and others, 2014). Here we used the SRTM3 V2 without void filling, which was provided by NASA in 2002 (http://dds.cr.usgs.gov/-srtm/version2_1/SRTM3/) (Table 1).

4. METHODOLOGY

4.1. Estimating the equilibrium-line altitude (ELA)

For estimating the ELA of Siachen Glacier we have assumed that the late-summer snow line altitude (SLA), i.e. the elevation of the snowline at the end of the hydrological year, can be considered as representative of the ELA (Rabatel and others, 2012). Thereafter, we identified the SLA visually on a number of Landsat images near the end of the ablation period between the years 2000 and 2014 (Table 4). To facilitate the identification of the SLA, we created false-color composites from the short wave infrared (SWIR), near infra-red (NIR) and one visible band (GREEN) and mapped the boundary between bright snow and darker ice (Rabatel and others, 2012). Since Siachen Glacier has a number of tributary glaciers with different SLAs, we computed the mean SLA from different tributaries. Finally, we averaged the SLAs from the different images to obtain the best estimate of the ELA of 5250 m a.s.l. The accumulation-area ratio was 0.60, which is similar to the value 0.58 mentioned by Hewitt (2011).

Table 4. List of Landsat Images used for ELA determination

4.2. Glacier mapping

For creating exposed-ice boundaries using Landsat data products available for 1990, 2000 and 2014, the ratio of RED to SWIR digital numbers is computed, followed by a selection of an appropriate threshold value of 2.1 (Paul and Kääb, 2005; Bolch and others, 2010). To overcome noise and isolated pixels in an automated manner, a 3 × 3 median filter was applied. In addition, a few misclassified pixels were edited and deleted manually in order to refine the boundary. The final binary image was converted to a vector shapefile. Cartosat-I stereo scenes with a spatial resolution of 2.5 m, orthorectified on the basis of the Cartosat-I DEM generated in this study, were used to delineate the exposed-ice boundary manually for the year 2007. As the glacier was not entirely covered by the Cartosat-I scenes, we used the gap-filled Landsat ETM+ scene of the year 2009 for delineating the remaining parts.

For mapping Siachen Glacier using Hexagon KH-9 imagery, ~100 tie points (TP) were collected to co-register the image to the Cartosat-I/Landsat ETM+ scene. TPs were obtained both automatically, from the orthorectified Landsat ETM+ metadata file, and manually, from an orthorectified and pansharpened SWIR-NIR-RED (6-4-3) false-color composite from the OLI scene of 2014. A spline transformation, which optimizes for local accuracy, was used to reference the images based on the TPs. Utilization of this technique with a large number of well-distributed TPs yielded a good accuracy (in the order of 3–4 m) for the entire study area. Thereafter, this georeferenced Hexagon image was orthorectified using the SRTM DEM of 30 m resolution with the Create Ortho-Corrected Dataset module in ArcGIS in order to correct for geometric distortions.

Debris-covered ice was mapped manually using Landsat color composites (SWIR, NIR, RED) for the years 1990, 2000 and 2014. High-resolution orthorectified Cartosat-I (along with the gap-filled Landsat image) and Hexagon imagery was used to manually delineate debris-covered ice boundaries for the years 2007 and 1980, respectively.

In order to support the identification of the debris-covered tongue and to study changes in snout position, we have used coherence mapping, which relies on high-precision interferometric techniques utilizing SAR data. This method has been effectively used to delineate glaciers automatically (Atwood and others, 2010), study changes in snout position (Saraswat and others, 2013) and identify debris-covered portions of glaciers (Frey and others, 2012).

Coherence (γ) is a measure of spatial correlation between two SAR acquisitions. Better coherence results in the generation of high quality interferograms, while loss of coherence results in decorrelation. Changes in vegetation growth, permafrost freezing and thawing, soil moisture content and glacier motion between two SAR acquisitions are some of the pertinent causes of decorrelation.

To assess the changes in the position of the debris-covered snout, coherence maps were generated for the period 2004–08 by processing Envisat ASAR images with JPL/Caltech ROI_pac software using the methodology detailed in Simons and Rosen (2007). Saraswat and others (2013) estimated that there is a bias of ~20% in the coherence estimates especially in glacierized regions, which have low coherence values. Since the members of the interferometric pairs used in this study were acquired one year apart (Table 2), the glacierized region is affected by motion-related decorrelation whereas the ice-free area is free from this effect. Thus the area near the glacier snout, where the ground is exposed, shows better coherence in contrast to the pixels of the surrounding glacier (Atwood and others, 2010; Frey and others, 2012). The stable areas ahead of the debris-covered snout have coherences in the range 0.5–0.9 while the glacierized areas have coherence in the range 0.0–0.5. We used a coherence between 0.50 and 0.55 (depending on the scene) to delineate the debris-covered tongue. Comparison of the boundary generated using coherence mapping with that created manually reveals good fit between the two (Figs 2a–d). Hence, this method can be particularly important to quantify changes in terminus position for heavily debris-covered glaciers in a semi-automated manner, since manual mapping of debris-covered glacier tongues requires high-resolution datasets, which can be difficult to obtain.

Fig. 2. Coherence images from Envisat ASAR data used to validate the manual mapping and to delineate debris-covered tongue for different years. Shown by the black line is the debris-covered tongue as delineated manually. (a–d) Represent the interferometric pairs (Table 2).

4.3. Mapping uncertainty

The precision of automatic delineation of exposed ice is usually within half a pixel (Bolch and others, 2010; Paul and others, 2013). Since no higher resolution datasets are available to manually correct the automatically derived boundaries, mapping uncertainty for the Landsat scenes of 1990, 2000 and 2014 was calculated by using a buffer distance of 15 m around the automatically derived boundaries. For manually delineated boundaries using the Hexagon scene of 1980, we assumed a mapping uncertainty of ~8 m (i.e. 1 pixel) (Pieczonka and others, 2013). In the case of 2007 glacier boundaries, delineated using Cartosat-I images, a mapping uncertainty of 2.5 m (1 pixel) is assumed. For the area not covered by the Cartosat-I image, an error estimate of 15 m was used, similar to that of Landsat images. Overall, this resulted in an uncertainty of 1.82% of the total exposed-ice area of Siachen Glacier mapped for the year 2007 (Table 5).

Table 5. Changes in area and length of Siachen Glacier, 1980–2014

4.4. Generation of velocity map

The offset field between pairs of ALOS PALSAR satellite data acquired with a 46-day interval was used for the estimation of glacier velocities. Offset tracking of L-band SAR images is a robust and direct technique for estimating glacier motion (Rignot, 2008; Strozzi and others, 2008; Pohjola and others, 2011). This is particularly useful when Differential Interferometric Synthetic Aperture Radar (DInSAR) is limited by the loss of coherence due to long time intervals.

Offset tracking measurements made along the slant-range (incidence angle is ~35° for ALOS PALSAR observations) and azimuth directions are combined to retrieve horizontal displacements on the ground. Slant-range and azimuth offset estimation errors of ALOS PALSAR data are in the order of one tenth of a pixel. For a FBS with ground-range pixel spacing of 7.5 m and an azimuthal pixel spacing of ~3 m, two-dimensional (2-D) ice velocities can be computed with an expected error of 0.8 m in 46 days or ~10 m a⁻¹. For a FBD with ground-range pixel spacing of ~15 m, 2-D ice velocities can be computed with an expected error of 20 m a⁻¹. The displacement maps were geocoded to geographical coordinates with the WGS84 horizontal datum and the EGM96 vertical datum at 90 m resolution, using the SRTM Digital Elevation Database v4.1 (Jarvis and others, 2008).

4.5. DEM generation and co-registration

All Cartosat-I scenes of the year 2007 were pre-processed using PCI Geomatica OrthoEngine 10.2. Data pre-processing involved radiometric enhancements using a locally-adaptive Wallis filter, which adjusts the brightness of greyscale images in local areas as opposed to other global filters (Pieczonka and others, 2013). This adjustment was very critical for shadow regions and accumulation areas of glaciers since they lack sufficient contrast, causing errors in parallax determination.

The next step in DEM generation was to compute the satellite stereo model, in order to determine the ground position of each point in the stereo scene. Cartosat-I scenes are provided with RPCs (Titarov, 2008). Even though these RPCs can be improved further using GCPs, sub-pixel horizontal accuracy can only be achieved when the GCPs are at least two or three times more accurate than the spatial resolution of the satellite image (Höhle and Höhle, 2009). For Cartosat-I data with 2.5 m resolution, the required horizontal accuracy of GCPs should be in the range 1–1.5 m, which for all practical purposes can only be achieved with in-situ measurements. Therefore we have not used any GCP during the preparation of the Cartosat-I DEM.

The Cartosat-I DEM was geocoded in the Universal Transverse Mercator (UTM) system; zone 43 N, with a grid size of 30 m. The SRTM DEM was projected into the UTM coordinate system and then resampled bilinearly from 90 m to 30 m resolution in order to make it conformable with the Cartosat-I DEM. About 7.5% of the glacier area in the accumulation zone is not covered by the Cartosat-I DEM (Fig. S1). To identify the outliers caused due to Cartosat-I DEM failure, hillshaded images of the Cartosat- I DEM and SRTM DEM were compared. The outliers on the glacier surface seen as black patches in the uppermost accumulation area in Figure S1 were removed before further processing.

4.6. Co-registration of DEMs

DEM differencing necessitates co-registration to ensure that corresponding pixels in the two DEMs represent the same ground location.

According to Nuth and Kääb (2011), elevation difference, slope and aspect of non-glacierized pixels are related by the following equations:

(1) $${\rm dh} = a \cdot \tan \alpha \cdot \cos \;(b - \psi ) + \; {\rm dh}{\rm ^{\prime}}$$

and dh, α and ψ represent elevation differences at individual pixels, terrain slope and aspect, respectively. The terms a, b, dh′ denote the magnitude of the horizontal shift, the direction of the shift vector and the overall elevation bias between the two DEMs.

The parameters a, b and dh′ are computed by least-squares optimization. These values were used to adjust the Cartosat-I DEM iteratively until a maximum of 10 iterations was reached or until the horizontal offset a was <1 m. Residual elevation differences over stable non-glacierized areas due to rotation and scale distortions were minimized by second-order spatial trend corrections (applied to the Cartosat-I DEM). In this process a second-order trend surface was calculated for all ice-free areas and added to the DEM (Pieczonka and others, 2011).

4.7. Accuracy assessments

4.7.1. Uncertainty of the mosaicked Cartosat –I DEM

We used the overlapping area (~98 km²) of the three Cartosat-I scenes to quantify the uncertainty in the final mosaicked DEM. Considering the differences between each individual scene to be independent and using the principle of error propagation, we computed a RMSE of ±1.2 m in the final mosaicked Cartosat-I DEM and this was considered as the uncertainty of the mosaicked DEM.

4.7.2. Relative vertical accuracy

In order to quantify the accuracy of the relative fit between the reference DEM (SRTM) and Cartosat-I DEM we computed the mean and standard deviation for ice-free stable areas (Table 6). The mean was taken as a quantitative estimate of the bias due to inaccuracies in the coregistration of the DEMs. Since DEM errors increase with increasing slope (Jacobsen, 2005) we restricted the accuracy assessment to regions with slopes <30° (Pieczonka and others, 2011). Elevation differences outside the range of mean ±3σ were considered as outliers (Gardelle and others, 2013). The statistics for the unadjusted DEM (Table 6) reveal that the removal of outliers significantly reduced both the bias and the standard deviation. The adjusted DEM shows significant improvement in the standard deviation (with and without outlier removal), which reflects a good fit of the two DEMs after co-registration. The bias after removal of outliers was subtracted from the elevation differences over glacierized terrain.

Table 6. DEM difference error statistics before and after co-registration

For the quantification of uncertainty we computed the normalized median absolute deviation and the 68.3% quantile for elevation differences over stable ice-free areas. Because of large differences in these two error estimates and the non-normal distribution of the elevation differences over stable areas we chose the 68.3% quantile as a conservative estimator for elevation difference errors (Höhle and Höhle, 2009). The 68.3% quantile was added to the error of the mosaicked Cartosat-I DEM to yield a final error estimate of ±2 m for the entire study area (Table 7).

Table 7. Estimates of elevation change and mass budget for Siachen Glacier for 1999–2007

In Scenario 1, a constant density of 850 ± 60 kg m⁻³ was considered while in scenario 2 densities for accumulation and ablation areas were 600 and 900 kg m⁻³, respectively.

4.7.3. Correction for C-band penetration and curvature correction

It is crucial to make adequate corrections for the penetration depth of the SRTM's C-band radar beams into snow and ice (Gardelle and others, 2012; Kääb and others, 2012). The penetration depth, however, depends mostly on the properties of the snow and ice. In the accumulation zone, C-band radar can easily penetrate through the snow, which accumulated on the surface of the previous ablation period and even into the different firn layers of the previous years, which can easily have a thickness of several tens of meters (Huss, 2013). On the contrary, C-band radar penetrates much less into the ice of the ablation region. Therefore it is necessary to make elevation-dependent corrections to the SRTM DEM in order to obtain the most precise estimates of glacier surface elevation. Moreover, since C-band radar penetrates well through snow, we have considered the SRTM DEM to be representative of the ablation surface at the end of the previous melt season (i.e., October 1999, Berthier and others, 2006; Paul and Haeberli, 2008). Even though the penetration of X-band radar (with a wavelength of 3.1 cm) into snow may not be negligible (Seehaus and others, 2015), it is significantly less in comparison to C-band radar (with a wavelength of 5.6 cm).

Thus, a comparison between SRTM-X and SRTM-C permits an estimate of the elevation dependency of the C-band penetration (Gardelle and others, 2012). Since SRTM-X did not cover Siachen Glacier, we computed elevation differences between SRTM-X and SRTM-C over nearby glaciers of the Eastern and Western Karakoram Range. A relationship between penetration correction and altitude was thus formulated (Fig. 3). The resulting C-band penetration depths, which increase from ~0.5 m to at least 6.5 m over the elevation range of the glacier, were subtracted on a pixel-by-pixel basis from the elevation differences (Cartosat-I DEM minus SRTM DEM) in order to get the “true” elevation differences. This approach is different from that of Kääb and others (2012), who used an average penetration correction of 2.4 ± 0.3 m for the Karakoram glaciers. Our average penetration correction of 3.2 m for Eastern and Western Karakoram glaciers is similar to the average penetration of 3.4 m computed over Karakoram glaciers by Gardelle and others (2013). We could not perform any curvature correction since Siachen Glacier is a gently sloping glacier (average slope of 13°) with most of the area having maximum curvature values between −1 and 0 (Gardelle and others, 2012).

Fig. 3. Smoothed difference between SRTM-X and SRTM-C (blue line), which was applied to correct the penetration of the C-band radar beam. Also shown is the ±1 standard deviation confidence region (red lines).

4.8. Calculation of elevation and mass changes

The elevation differences obtained after DEM and penetration correction were averaged within 100 m altitude intervals (from 3400 to 7300 m), and the average elevation change for the entire Siachen Glacier was computed as an area-weighted average of elevation change in each altitude range. The same procedure was followed to calculate the average elevation changes for the ablation and accumulation regions.

The average elevation change obtained can then be multiplied by suitable densities to obtain mass change estimates. In this study, we have used two density scenarios: (i) density of 850 ± 60 kg m⁻³ (Huss, 2013) for the entire glacier; (ii) densities of 600 and 900 kg m⁻³ (Hagg and others, 2004; Schiefer and others, 2007; Kääb and others, 2012) for the accumulation and ablation areas, respectively.

5. RESULTS

5.1. Glacier area changes

The total area of the Siachen Glacier system was 936.9 ± 21.1 km² in the year 2000. The area experiences only insignificant variations throughout the study period (Table 5). The exposed-ice area has decreased by 1.3% from 1980 to 2014 while the debris-covered area has increased by 10.5% during the same period (Table 5). Although it is evident that the debris-covered tongue has remained stable over the study period, fluctuations of the boundary between exposed and debris-covered ice are much more prominent. The distal part of this boundary retreated by nearly 1.3 km during the short-time span of 2007 to 2014, even though it remained practically unchanged prior to that (Fig. 4).

5.2. Velocity patterns

Velocity information for almost the entire glacier area could only be obtained for the 11 December 2008–26 January 2009 period, since large parts of the accumulation area suffer from decorrelation during the summer period (8 June 2007–24 July 2007). This temporal decorrelation is most likely due to snowfall between June and July 2007. In general, glacier velocity during the summer acquisition (8 June 2007–24 July 2007) period is greater than the velocity during the winter acquisition (11 December 2008–26 January 2009) period (Fig. 5). However, the general pattern reveals decreasing velocities towards the tongue and very low velocities on the frontal part of the tongue. The greater velocity of Siachen Glacier during summer can be attributed to basal lubrication caused by sub-glacial water flow during the ablation period (Copland and others, 2009). The average velocity for the entire glacier during 11 December 2008–26 January 2009 is 12.3 ± 0.4 cm d⁻¹. The average velocities for the accumulation and ablation regions are 9.7 ± 0.4 cm d⁻¹ and 20.4 ± 0.4 cm d⁻¹, respectively. We studied the profiles of glacier velocity for the winter season (11 December 2008–26 January 2009) over the main trunk glacier and the western tributary glaciers, which have shown anomalous behavior (see Table 8). On the main trunk profile AA′ (Fig. 6), velocity in the accumulation region and upper parts of the ablation area (section A) exceeds 40 cm d⁻¹ (Profile AA′, Fig. 7a) and drops steadily in the ablation area (section B) from a maximum at the end of Section A. However, the lower ablation area and areas near the snout (section C) are characterized by a drastic decrease in velocity. All the western tributary glaciers (except the one with profile FF′) have velocities lower than average. The tributary with profile BB′ shows near-zero velocity in most of the accumulation area and lower portions of the ablation area but higher velocity in the middle. While the tributary with profile DD′ shows decreasing velocities throughout its length, those with GG′ and EE′ show relatively low velocities persistent through the entire length of the accumulation and ablation areas.

Fig. 4. Extent of exposed-ice on the lower tongue for 1980–2014.

Fig. 5. Velocity pattern over Siachen Glacier for two different seasons: (a) 8 Jun 2007–24 Jul 2007 (b) 11 Dec 2008–26 Jan 2009.

Fig. 6. Map displaying the traces of profile lines used to construct elevation change and velocity patterns for different glaciers. The letters (A–H) in the figure refer to the trace of the particular profile line; for example, letter A refers to the trace of the profile line AA′. Profile lines AA′ to GG′ run from accumulation to ablation area while profile lines HH′ to JJ′ are cross-profiles running from left to right (looking down-glacier).

Table 8. Velocity and elevation change for profiles along the main trunk of Siachen Glacier and selected tributary glaciers

The velocity estimates are for the period 26 January 2009–11 December 2008. Elevation changes, averaged along each profile after correction for C-band penetration are for the period 25 November 2007–11 to 20 October 1999. Uncertainties are ±2.0 m for elevation changes and 0.004 m d^–1 for velocities.

5.3. Elevation and mass changes

Elevation changes observed here are spatially heterogeneous with accumulation regions at higher altitudes showing thickening, while ablation zones at lower altitudes show characteristic thinning (Fig. 8). The three tributary glaciers with an area >100 km² show on average thickening because their accumulation areas are large in proportion to their ablation areas. Over smaller tributaries with an area >50 km² the average elevation change varies in dependence upon their range of altitude. Among these tributary glaciers those that terminate at altitudes > 5000 m have positive elevation change, while those terminating at lower elevations (<5000 m) show negative or zero elevation changes. The western tributary glaciers, particularly those joining the main trunk at lower elevations, show anomalous behavior in terms of their average velocity, which is lower than that of the main trunk. To get a clear picture of elevation changes across the whole glacier, profiles of elevation differences (after penetration correction) were constructed along the main trunk and tributaries of Siachen Glacier. For section A of the main trunk (Fig. 9a), extending from an elevation of 5550 m (point A) to 4850 m, insignificant elevation changes (−5 to +5 m) are observed during the study period. Section B covers most of the ablation region and lies in an altitude range 4850–3750 m. In the upper parts of the ablation region (between 4850 and 4000 m), slight thinning (up to −10 m) is observed. The thinning becomes more prominent in the lower reaches of section B (4000–3750 m) where thinning is in the range 10–25 m. Section C lies in the lower reaches of the ablation zone, near the debris-covered tongue. Here the thinning is more prominent, reaching a maximum of 35–40 m. Thus, the debris-covered part down-glacier from the exposed-ice regions shows considerable thinning, which suggests rapid downwasting in these regions of low flow velocity. The mean thinning for this terminus zone is 30 m for the period 1999–2007. It is interesting to note that the variations in elevation change for the smaller western tributary glaciers do not follow a fixed pattern (Figs 9b–g). While profiles CC′ and GG′ show insignificant elevation changes, profile DD′ shows thinning in the accumulation area and major portions of the ablation area. Profile EE′ shows minor thickening in the accumulation region, but major thinning (−20–70 m) in the ablation area, which is probably due to the blockage of its flow by other tributary glaciers. On the contrary, profile FF′ registers thinning throughout its length. Similarly, profile BB′ has shown thinning near the point B′, where it meets the glacier trunk.

Fig. 7. Velocity patterns along profiles AA′ to GG′ for the period 11 Dec 2008–26 Jan 2009 (Fig. 6). The curves have been smoothed using a moving average filter with a span of five observations.

Fig. 8. Elevation difference map for Siachen Glacier for the period 1999–2007.

The mass budget of Siachen Glacier is computed to be −0.03 ± 0.21 m w.e. a⁻¹ for 1999–2007 (Table 7) considering the density scenario 1 (uniform density) and −0.05 ± 0.18 m w.e. a⁻¹ for the density scenario 2 (different densities in accumulation and ablation regions). Hence, the density scenarios hardly affect the result and indicate a more or less balanced mass budget.

5.4. Effect of debris cover

Debris constitutes a key component of a glacier system and affects the glacier in different ways. Usually thin debris enhances melt, while thick debris cover reduces melt (Nicholson and Benn, 2006). Gardelle and others (2012) reported that, on average, the elevation loss of exposed and debris-covered areas for the glaciers of the Eastern Karakoram are similar. To test this for Siachen Glacier, elevation differences over exposed and debris-covered ice were computed separately for lower elevation regions where debris-covered ice occurs (Fig. S1). However, it is important to note that exposed and debris-covered ice are dynamically coupled; hence, to observe possible subtle differences between the two, elevation changes in exposed and debris-covered regions have been classified into similar altitude intervals (Fig. 10). Here, we have randomly selected pixels over exposed and debris-covered ice such that the sampled pixels have similar frequency distribution with respect to their altitude. Thinning over debris-covered ice is slightly less than or comparable with that over exposed-ice areas for the same altitude range. One probable reason for this interesting observation is that debris-covered parts lying near the glacier margins have lower velocities compared with the exposed-ice parts located in the middle and having higher velocities (Fig. S2). In the lower ablation region (3600–4000 m), both exposed and debris-covered ice show strong thinning. In the elevation range 4000–4500 m, exposed ice shows equal or higher rates of thinning compared with debris-covered ice, although both rates are generally less than those below altitude of 4000 m. The thinning in the altitude band 4500–5000 m is smaller still, and debris-covered and clean ice show similar elevation changes. Among the debris-covered areas, the region between exposed ice and the debris-covered snout (with an elevation range of 3600–3700 m) demonstrates the highest thinning with a mean value of −30 m. One probable cause of the downwasting might be the presence of glacial lakes, which can be seen along the main trunk and along the lower western tributary glaciers. These lakes could, along with ice cliffs, enhance the ice melt at the debris-covered portions of the glacier (Ragettli and others, 2016).

Fig. 9. Elevation changes along profiles AA′ to GG′ (Fig. 6). The curves have been smoothed using a moving average filter with a span of five observations.

Fig. 10. Elevation changes over exposed and debris-covered ice as a function of altitude. Each elevation change has an uncertainty of ±2.0 m.

6. DISCUSSION

6.1. Elevation changes and velocity patterns

The average velocity for Siachen Glacier, 12.3 ± 0.4 cm d⁻¹ during winter, is in close agreement with the average velocity of 13.7 cm d⁻¹ reported by Copland and others (2009) for the nearby Baltoro Glacier for the period 26 July 2006–27 June 2007. Annual velocity estimates for other debris-covered glaciers such as Gangotri Glacier (Scherler and others, 2007; Bhattacharya and others, 2016) and Khumbu Glacier (Quincey and others, 2009; Bolch and others, 2008) lie in the range 5.5–11 cm d⁻¹, which is lower than the winter velocity of Siachen Glacier reported here. However, spatial heterogeneity needs to be taken into account when assessing how mass budget might affect velocity patterns and vice versa. Heid and Kääb (2012) have shown that Siachen Glacier has negligible velocity changes in most of the lower ablation region, minor acceleration in its middle parts and large acceleration in its upper parts during the years 2000–08. Our study suggests that, in the upper ablation areas of the main trunk, climatic and topographic conditions resulted in high ice flux and insignificant elevation changes. However, low ice flux into the lower reaches of the ablation area is probably one reason for thinning. The region between snout and exposed-ice records thinning of 30 m for the period 1999–2007 and exhibit low velocity conditions. On the western tributary glaciers, profiles BB′, CC′, DD′, EE′ and GG′ (Fig. 6) showed significantly lower velocities than the average velocity of Siachen Glacier. These western tributary glaciers are characterized by an overall thinning and low velocities (Table 8). Profile FF′ (Figs 7f, 9f), located on the only western tributary glacier, which has velocities comparable with the average velocity of Siachen Glacier, shows elevation losses throughout its length. Profile DD′ shows elevation loss even in the accumulation regions and elevation gain in parts of the lower ablation zone. This is a hint for a recent surge of this tributary; tributary surges are quite common in the Karakoram (Belò and others, 2008; Paul, 2015). The velocity information from 2007 to 2008 for this tributary shows no abnormal behavior and hence indicates that the active phase of the possible surge was over by then.

6.2. Comparison of mass budget and area change estimates

The area of Siachen Glacier has not changed significantly since 1980. Similar results are obtained for the change in elevation of the glacier, which is equal to −0.23 ± 2.00 m for the 8 years between 1999 and 2007 or −0.03 ± 0.25 m a⁻¹. Moderate thickening at high altitude is evidently balanced by significant thinning at lower altitude (Fig. 8). Moreover, for estimation of mass budget, we had considered two density scenarios, the results from which were indistinguishable. There exists only one other recent mass budget estimate for Siachen Glacier. Our estimate (−0.03 ± 0.21 m w.e. a⁻¹) is roughly comparable with the +0.14 ± 0.14 m w.e. a⁻¹ for the period 2000–08 reported by Gardelle and others (2013). However, their estimate did not include contributions from the lower portions of the ablation region, where we found significant elevation loss. On the other hand, some portions of the accumulation region, which are potentially gaining mass, are unaccounted for in the present study. It seems that the balanced mass budget of Siachen Glacier is not a recent phenomenon since Zaman and Liu (2015) reported a mass budget with upper and lower bounds of +0.22 and −0.23 m w.e. a⁻¹ during the period 1986–1991. The mass budget results computed by Zaman and Liu (2015) using hydrological method with probable overestimation of mass loss due to evaporation seem to be consistent with results of this study. Our study shows that while the snout has been stable during the period 1980–2014, the transition between exposed and debris-covered ice advanced down-glacier between 1980 and 2007 and retreated rapidly up-glacier (by ~1.3 km) thereafter (2007–14). However the changes we describe above is in reference to the main trunk of Siachen Glacier only. Overall, exposed-ice area has decreased steadily, though slightly, over the period 1980–2014 (Table 5).

7. CONCLUSION

This study presents area change from 1980 to 2014 and mass budget from 1999 to 2007 of Siachen Glacier in the eastern Karakoram. The area of Siachen Glacier has shown no significant variations during the study period. However, the area of the exposed ice decreased steadily between 1980 and 2014. Also evident is the sudden drastic retreat (~1.3 km) of the distal part of the transition between exposed and debris-covered ice. This is likely to be caused by downwasting of ice in the lower ablation reaches near the terminus of the glacier. These observed changes are despite the fact that a near-zero mass budget for Siachen Glacier (−0.03 ± 0.21 m w.e. a⁻¹) has been obtained for 1999–2007.

The average velocity of Siachen Glacier for December 2008–January 2009 is 12.3 ± 0.4 cm d⁻¹. The velocity distribution for the main trunk of Siachen Glacier corresponds to that of a typical debris-covered glacier in the Himalayas, with low flow velocity at the head, increasing velocities in the upper accumulation area, almost constant velocities in the lower accumulation area and rapid decline towards velocities near zero near the terminus. The western tributary glaciers and lower ablation area of Siachen Glacier are sites of low ice flux. Debris-covered ice shows elevation losses comparable with those of exposed ice occurring at similar altitudinal ranges. The debris-covered ice located mainly near the margins of the glacier suffers from higher elevation losses. Thus, the spatial heterogeneity is an important component of the Siachen Glacier system and it should be analyzed in greater detail for future studies.