2.1. Site description

Chhota Shigri glacier (32.2° N, 77.5° E) is a valley-type glacier located in the Chandra-Bhaga river basin of Lahaul and Spiti valley, Pir Panjal range, western Himalaya. This glacier extends from 6263 to ~4050 m a.s.l., is ~-9km long and covers an area of 15.7 km². Its snout is easy to locate from one year to the next because it is well defined, lying in a narrow valley and giving birth to a single proglacial stream. The main orientation of the glacier is north, but its tributaries have a variety of orientations (Fig. 1). The lower ablation area (<4400ma.s.l.) is partly covered by debris representing ~3.4% of the total surface area. The glacier is located in the monsoon-arid transition zone and is influenced by two atmospheric circulation systems: the Indian monsoon during summer (July-September) and the Northern Hemisphere mid-latitude westerlies during winter (January-April) (Reference Singh, Jain and KumarSingh and others, 1997;Reference Bookhagen and BurbankBookhagen and Burbank, 2006;Reference Gardelle, Arnaud and BerthierGardelle and others, 2011).

Fig. 1. Map of Chhota Shigri glacier with the measured transverse cross sections (lines 1-5), the ablation stakes (dots) and the accumulation sites (squares). Also shown are longitudinal sections (lines A-E) used to calculate thickness and ice velocity variations (see Section 4.2). The map (contour lines, glacier delineation) was constructed using two stereoscopic pairs of SPOT5 (Systeme Probatoire pour l’Observation de la Terre) images acquired on12and13 November 2004 and 20 and 21 September 2005 (Reference WagnonWagnon and others, 2007). The map coordinates are in the UTM43 (north) World Geodetic System 1984 (WGS84) reference system.

2.2. Mass balance

The first series of mass-balance measurements on Chhota Shigri glacier was performed between 1987 and 1989 (Reference Nijampurkar and RaoNijampurkar and Rao, 1992;Reference Dobhal, Kumar and MundepiDobhal and others, 1995; Reference KumarKumar, 1999). The bedrock topography and surface ice velocity were also surveyed over the same period by gravimetric and stake displacement methods respectively (Reference Dobhal, Kumar and MundepiDobhal and others, 1995;Reference KumarKumar, 1999). We reinitiated the mass-balance observations in 2002. Since that year, annual surface mass-balance measurements have been carried out continuously on Chhota Shigri glacier at the end of September or the beginning of October using the direct glaciological method (Reference PatersonPaterson, 1994). Ablation was measured through a network of ~22 stakes distributed between 4300 and 5000 m a.s.l. (Fig. 1), whereas in the accumulation area the net annual accumulation was obtained at six sites (by drilling cores or pits) between 5100 and 5550ma.s.l. (Reference WagnonWagnon and others, 2007). In the accumulation area, the number of sampled sites is limited due to difficulty of access and the high elevation. The glacier-wide mass balance B_a is calculated according to

(1)

where b_i is the mass balance of the altitudinal range i (m w.e. a^-1), of map area s_i , and S is the total glacier area. For each altitudinal range, b_: is obtained from the corresponding stake readings or net accumulation measurements.

2.3. Surface velocity

Annual surface ice velocities were measured at the end of each ablation season (September-October) by determining the annual stake displacements (~22 stakes) using a differential GPS (DGPS). These geodetic measurements were performed in kinematic mode relative to two fixed reference points outside of the glacier on firm rocks. The accuracy of x (easting), y (northing) and z (elevation) at each stake position is estimated at ±0.2 m depending mainly on the size of the hole in which the stake was set up. Thus the surface ice velocities measured from stake displacements have an accuracy of ±0.3 m a^-1.

2.4. Ice thickness

GPR measurements were conducted in October 2009 to determine ice thickness on five transverse cross sections (Fig. 1) between 4400 and 4900 m a.s.l. A pulse radar system (Icefield Instruments, Canada) based on the Narod transmitter (Reference Narod and ClarkeNarod and Clarke, 1994) with separate transmitter and receiver, was used in this study with a frequency centred near 4.2 MHz and an antenna length of 10 m. Transmitter and receiver were towed in snow sledges along the transverse profile, separated by a fixed distance of 20m, and used to record measurements every 10 m. The positions of the receiver and the transmitter are known through DGPS measurements, within an accuracy of ±0.1 m. The speed of electromagnetic wave propagation in ice has been assumed to be 167 m ms^-1 (Reference Hubbard and GlasserHubbard and Glasser, 2005). The field measurements were performed in such a way as to obtain reflections from the glacier bed located more or less in the vertical plane with the measurement points at the glacier surface, allowing the glacier bed to be determined in two dimensions. The bedrock surface was constructed as an envelope of all ellipse functions, which give all the possible reflection positions between sending and receiving antennas. Ice thickness was measured along four transverse profiles (profiles 1-4) on the main glacier trunk and one (profile 5) on a western tributary (Fig. 1).

3.1. Glacier-wide mass balance and mass-balance profile

The annual glacier-wide mass balance and cumulative mass balance of Chhota Shigri glacier between 2002 and 2010 are plotted in Figure 2. The glacier-wide mass balance was negative except for three years (2004/05, 2008/09 and 2009/10). It varies from a minimum value of −1.42 m w.e. in 2002/03 to a maximum of +0.33 m w.e. in 2009/10. The cumulative mass balance of Chhota Shigri is -5.37 m w.e. between 2002 and 2010, while the glacier-wide mass balance averaged over the same period is −0.67 mw.e. a^-1.

Fig. 2. Cumulative (line) and annual glacier-wide mass balances (histograms) of Chhota Shigri glacier during 2002–10. The value of annual glacier-wide mass balance (mw.e.) is also given with the corresponding histogram

The quantitative uncertainty associated with the glacio- logical mass balance requires a distinction between the accumulation zone and the ablation zone. In the accumulation zone, the surface mass-balance measurements were obtained from shallow boreholes (auger). Therefore, they are based on core length and density determination. In the ablation zone, the measurements have been carried out from ablation stakes. The overall errors (standard deviation) on point measurements are estimated at 0.30 and 0.15 m w.e. in the accumulation and ablation zones, respectively. The overall error comes from a variance analysis (Reference Thibert, Blanc, Vincent and EckertThibert and others, 2008) applied to all types of errors (ice/snow density, core length, stake height determination, liquid-water content of the snow, snow height). Although conducted on a glacier in the European Alps, the analysis of Reference Thibert, Blanc, Vincent and EckertThibert and others (2008) can be generalized to other glaciers because it is based on measurement errors that are similar on every glacier when using the glaciological method. However, only six sites are sampled in the accumulation zone (11.6 km²), and 22 sites in the ablation zone (4.1 km²). The uncaptured spatial variability of surface mass balance may cause systematic errors in the glacier-wide mass balance. In the accumulation zone, the spatial variability remains unknown and is probably very high as observed for other glaciers (e.g. Reference Machguth, Eisen, Paul and HoelzleMachguth and others, 2006). In the ablation zone, stakes set up at the same altitude show similar values except on the terminal tongue which is debris-covered (0.54 km²). Consequently, the overall uncertainties in mass-balance profile have been assessed at 0.5 mw.e. in the accumulation zone, 0.25 m w.e. in the white ablation zone and 0.5 m w.e. in the debris-covered area of the glacier. The surface area estimation also causes systematic error. The estimated uncertainty in the surface area calculated for each altitudinal range is 5%. Combining these errors at different altitudinal ranges using Eqn (1), the uncertainty in the annual glacierwide mass balance is 0.4mw.e.a^-1. As revealed by other studies (e.g. Reference VincentVincent, 2002;Reference Thibert, Blanc, Vincent and EckertThibert and others, 2008;Reference Huss, Bauder and FunkHuss and others, 2009), this estimation confirms that the glacio- logical method needs to be calibrated by a volumetric method over a long period of monitoring (>5 years) in order to limit the systematic errors and improve the accuracy of absolute values of mass balance. Note that the uncertainty of relative changes in mass balance from year to year is smaller than those inherent in annual mass balances, as the influence of systematic errors can be reduced.

We also calculated the mass-balance profile between 2002 and 2010 (Fig. 3). For each altitudinal range, we computed the average of all available measurements. Figure 3 shows that melting in the lowest part of the ablation area (<4400ma.s.l.) is reduced by ~1 mw.e.a^-1 irrespective of its altitude. This is due to the debris cover (approximately 5-10 cm thick debris mixed with isolated rocks) which reduces the melting in this region (Reference Mattson, Gardner and YoungMattson and others, 1993;Reference WagnonWagnon and others, 2007). Moreover the lower part of Chhota Shigri glacier flows in a north-south- oriented deep and narrow valley (Fig. 1), causing the glacier tongue to receive less solar radiation due to the shading effect of the steep valley slopes.

Fig. 3. The 2002–10 average mass-balance profile and hypsometry of Chhota Shigri glacier. Altitudinal ranges are of 50m (e.g. 4400 stands for the range 4400–4450 m), except for 4250 and 5400 which stand for 4050–4300 and 5400–6250m respectively.

3.2. Ice thicknesses and cross-sectional areas

Thanks to clear reflections, the ice/bedrock interface was generally easy to determine in all profiles. Figure 4 provides an example of the radargram obtained at cross section 2. A radar wave velocity of 167 m ms^-1 was used to calculate ice thickness at all the profiles. The cross sections obtained from GPR measurements reveal a valley shape with maximum ice thickness greater than 250 m (Fig. 5). The centre-line ice thickness increases from 124 m at 4400ma.s.l. (cross section 1 in Fig. 1) to 270 m at 4900ma.s.l. (cross section 4). This confirms that the thicknesses obtained by gravimetric methods in 1989 (Reference Dobhal, Kumar and MundepiDobhal and others, 1995), twice as low as the present results, were underestimated as proposed by Reference WagnonWagnon and others (2007). The cross-sectional areas are given in Table 1. The accuracy of the calculated ice thickness is determined, in part, by the accuracy of the measurement of the time delays and the antenna spacing. Additional errors may arise because the smooth envelope of the reflection ellipses is only a minimal profile for a deep valley-shape bed topography, with the result that the ellipse equation will be governed by arrivals from reflectors located toward the side and thus not directly beneath the points of observation. Further errors may be introduced by assuming that all reflection points lie in the plane of the profile rather than on an ellipsoid. No errors associated with radar wave velocity variations between snow and ice have been accounted for because all cross sections were surveyed in the ablation zone or slightly above (with the firn/ice transition depth at the surface or <2 m deep). Hence, the radar wave velocity for ice (167 m ms^-1) was used to calculate all ice depths. The estimated overall uncertainty in ice thickness is ±15 m. Given that the uncertainty in ice surface coordinates is low (±0.1 m), the uncertainty in crosssectional areas mainly arises from the uncertainty in ice thickness. The uncertainties in cross-sectional areas are 16%, 9%, 10%, 10% and 15% for cross sections 1, 2, 3, 4 and 5 respectively.

Fig. 4. Radargram of cross section 2: radar signals plotted side by side from west to east in their true spatial relationship to each other (interval between each signal of 10 m). The x-axis gives the amplitude of each signal (50mV per graduation); the y-axis is the doubletime interval.

Fig. 5. Ice depth and surface topography of cross sections 1–5. The horizontal and vertical scales are the same for all cross sections. All cross sections are oriented from west to east except cross section 5 which is north–south oriented.

Table 1. Calculated ice flux, mean surface ice velocity and maximum ice depth at each cross section. The mean surface horizontal ice velocities are from DGPS measurements performed in 2003/04. The satellite-derived mean ice velocities are from the correlation of satellite images acquired on 13 November 2004 and 21 September 2005 (NA: not available)

3.3. Ice velocity

Annual surface ice velocities were also measured between 2002 and 2010. However, some data gaps exist due to discontinuous DGPS signal, or loss of stakes. The 2003/04 ice velocities were used in this study because they provided the most complete dataset (Fig. 6). The centre-line horizontal ice velocities at each cross section were calculated by linear interpolation along the centre line between the velocities measured immediately up- and downstream of the cross section (ablation stakes visible in Fig. 1). Mean cross-sectional velocities are required to compute the ice fluxes (see Section 3.4). A map of the surface ice velocity field was derived by correlating SPOT5 images acquired on 13 November 2004 and 21 September 2005 (Reference BerthierBerthier and others, 2005). Comparison of the satellite-derived velocities with 16 nearly simultaneous DGPS velocity measurements shows a mean difference of 0.2 m a^-1 and a standard deviation of 1.6 m a^-1. The ratio between the centre-line horizontal velocity and the mean surface velocity (all extracted from the satellite-derived 2004/05 velocity field) was found to be 0.80 and 0.78 for cross sections 2 and 3 respectively. Reliable velocity measurements could not be made with SPOT5 imagery for other cross sections. Using the mean value of 0.79, the mean horizontal velocity was calculated from the centre-line velocity for each gate cross section (Table 1).

Fig. 6. Measured ice velocities plotted as a function of the distance from the 2010 terminus position. Measurements were collected along the central flowline.

3.4. Ice fluxes from kinematic method

The ice flux Q (m³icea^-1) through each cross section was calculated using the cross-sectional area S _c (m²) and depth- averaged horizontal ice velocity U (ma^-1).

(2)

The depth-averaged horizontal ice velocity was derived from the mean surface ice velocity calculated in Section 3.3. Reference NyeNye (1965) gives ratios of depth-averaged horizontal ice velocity to mean surface ice velocity varying from 0.8 (no sliding) to 1 (maximum sliding). Here we assume a mean basal sliding, with a constant ratio of 0.9. The calculated ice fluxes and maximum depth at each cross section are given in Table 1. The flux through cross section 3 at 4750ma.s.l. is higher than the flux through cross section 4 at 4900 m a.s.l. This is due to the ice influx from the western part of the glacier (flux through cross section 5) which contributes to cross section 3 and not to cross section 4 (Fig. 1).

The largest uncertainty in the depth-averaged horizontal ice velocity results from the ratio between the depth velocity and the surface flow velocity. The estimated factor 0.9 and unknown variations in the basal sliding lead to an uncertainty of roughly ±10% in the calculated flux, which lies within the range of uncertainty of the other variables as discussed by Reference Huss, Sugiyama, Bauder and FunkHuss and others (2007). Consequently, we can infer that depth-averaged horizontal ice velocity at each cross section is known with an accuracy of 1.0-3.0ma^-1 depending on the cross sections. Combining these errors in the cross-sectional area and mean velocity, the uncertainties in the ice fluxes are 0.21, 0.69, 0.92, 0.89 and 0.38 × 10⁶ m³ a^-1 for cross sections 1, 2, 3, 4 and 5 respectively. We have considered the errors to be systematic, so these uncertainties are probably overestimated.

3.5. Ice fluxes obtained from surface mass balance

We also calculated ice fluxes using annual surface mass balance measured during 2002-10. Although dynamic changes are neglected here, this method allows us to estimate the ice fluxes for each section from mass-balance data according to

(3)

where Q is the ice flux (converted into m³ ice a^-1 using an ice density of 900 kg m^-3, hence the factor 1/0.9) at a given elevation, z, and b_i is the annual mass balance of the altitudinal range i of map area s,-. The altitudinal ranges taken into account in the calculation are located between z and the highest range of the glacier z _max (highest altitude of the glacier area contributing ice to the cross section). We assume that at each point of the glacier above z the surface elevation has remained unchanged from one year to the next.

The ice fluxes calculated from annual mass-balance data at the five cross sections each year are given in Table 2 (further below), while the average ice fluxes for the 8 years are given in Figure 7. The uncertainties in ice fluxes resulting from surface mass balance are directly derived from the mass-balance uncertainties (see Section 3.1) applied to areas contributing to each cross section.

Table 2. Ice fluxes (10⁶m³ ice a^–1), inferred at each cross section from annual mass-balance data

Fig. 7. Ice fluxes at every cross section derived from 2003/04 ice velocities and section areas (open squares), and from mass-balance method for a glacier-wide SMB = 0mw.e. (plain squares) or a glacier-wide SMB = –0.67mw.e. (triangles). The error range for mass-balance fluxes calculated from the mass-balance method assuming a steady state (±1 standard deviation: hyphens) is also given.

4.1. Null to slightly positive mass balance during the 1990s inferred from ice fluxes

The ice fluxes obtained by the kinematic method using ice thickness and 2003/04 ice velocities are much higher than the average fluxes derived from the 2002-10 SMBs, the latter often being negative (Table 2). Thus, to assess the mean state of the glacier corresponding to the ice fluxes obtained by the kinematic method, here we compare these measured ice fluxes to theoretical ice fluxes calculated from SMB assuming the glacier to be in steady state. The glacier-wide mass balance obtained by the glaciological method is -0.67 m w.e. a^-1 over the 2002-10 period. Consequently, the SMB needs to be increased by 0.67 m w.e. a^-1 for the glacier to be in steady state with the present surface area. For each year (2002-10), we calculated the theoretical ice flux from SMB at each cross section assuming the glacier was in steady state. For this purpose, every year, a theoretical SMB at each elevation has been calculated by subtracting the overall annual specific SMB of the same year. For instance, year 2002/03 was characterized by an annual glacier-wide SMB of-1.42 m w.e., so we calculated a new SMB profile by adding 1.42 m w.e. to the SMB at each elevation. In contrast, 2009/10 was characterized by an annual glacier-wide SMB of +0.33 m w.e., so we calculated a new SMB profile by subtracting 0.33 m w.e. from the SMB at each elevation. The resulting ice fluxes are reported in Table 3, together with the mean ice flux at each cross section over the 8 years and the corresponding standard deviations.

Table 3. Ice fluxes (10⁶ m³ ice a ¹), obtained at every cross section, using steady-state mass-balance assumption for every surveyed year

These ice fluxes are close to the 2003/04 ice fluxes obtained by the kinematic method (Fig. 7), indicating that the dynamic behaviour of the glacier in 2003/04 is representative for steady-state conditions. This suggests that during the one to two decades preceding 2003/04, the glacier-wide mass balance of this glacier has probably been close to zero and that, in 2003/04, the ice fluxes had not adjusted to the last 3-4 years’ negative SMB.

This result is also supported by other observations. First, the ice velocities measured in 1987/88 (Reference Dobhal, Kumar and MundepiDobhal and others, 1995) are very close to the 2003/04 values (Fig. 6), suggesting that the dynamic behaviour of the glacier did not change much between 1988 and 2004. Second, the terminus fluctuation measured between 1988 and 2010 shows a moderate retreat of 155 m, equivalent to only 7ma^-1, in agreement with conditions not far from steady state. Given that Reference Berthier, Arnaud, Kumar, Ahmad, Wagnon and ChevallierBerthier and others (2007) observed a glacier-wide SMB of Chhota Shigri glacier of approximately -1 mw.e. a^-1 during the period 1999-2004, the glacier is likely to have experienced a null to slightly positive mass balance between 1988 and the end of the 20th century.

4.2. Glacier dynamics starting to adjust to 21st- century negative SMB

In theory, the response of ice fluxes to surface mass balance is immediate (Reference Cuffey and PatersonCuffey and Paterson, 2010, p.468), but observations show a 1-5 year delay (Reference Vincent, Vallon, Reynaud and Le MeurVincent and others, 2000, Reference Vincent, Soruco, Six and Le Meur2009;Reference Span and KuhnSpan and Kuhn, 2003). For instance, Reference Span and KuhnSpan and Kuhn (2003) found synchronous decrease in ice velocity between eight glaciers in the European Alps, which are driven by the same mass-balance changes (Reference Vincent, Le Meur, Six and FunkVincent and others, 2005). Consequently, the recent dynamic behaviour of Chhota Shigri glacier should be affected by the negative mass balance since 1999. However, the stake network on Chhota Shigri glacier, originally designed for SMB measurements, is not best suited to accurately compare either the ice velocities or the thickness variations because the measurements have not been performed at exactly the same location every year and are mainly restricted to the ablation area.

In spite of the above limitation, an attempt has been made to compare ice velocities and elevations from the available stake network. For this purpose, stakes measured at the beginning and end of the series have been selected on five short longitudinal cross sections (A-E in Fig. 1) along the centre line of the glacier where the network is most dense. The elevations in 2003 and 2010 and the ice velocities in 2003/04 and 2009/10 have been reported on these longitudinal cross sections to deduce thickness and velocity changes in the ablation area (Fig. 8;Table 4). Although the accuracy of the results is affected by the distance between the point measurements, we conclude that the part of the glacier below 4750 m a.s.l. is in strong recession. First, the thickness has decreased by 0.7-1.1 ma^-1 over the last 7 years. Second, the ice velocities have decreased by ~7m a^-1 between 2003 and 2010, resulting in a 24-37% decrease in the ice fluxes since 2003. Despite an improvable monitoring network, it may be surmised that the ice fluxes have been affected by the negative glacier-wide mass balance during (at least) the last 8 years, and the glacier dynamics are progressively adjusting to the negative SMB. Consequently, we expect terminus retreat to accelerate in forthcoming years. If the SMB remained equal to its 2002-10 average value in the future, the terminus would retreat by 5.6 km to reach 4870 m a.s.l. (the altitude where the ice flux is equal to zero) (Table 2).

Fig. 8. Elevation (dots) and surface ice velocity (triangles) between 2003 (solid lines) and 2010 (dashed lines) along the longitudinal sections A-E shown in Figure 1.

Table 4. Thickness and surface velocity changes between 2003 and 2010 on five longitudinal cross sections (NA: not available)