Study area

Our regional-scale investigation of debris-covered area changes is carried out in the Hispar and Shimshal sub-catchments of the Hunza River basin in the Karakoram Range, at about 36° N, 75° E in the Northern Territory of Pakistan (Fig. 1). The Hunza River is a tributary of the Gilgit River, which flows into the Indus River. The climate of the area is dominated by westerlies with maximum precipitation in winter, and glaciers are winter-accumulation type (Reference Gardelle, Berthier, Arnaud and KääbGardelle and others, 2013; Reference Ragettli, Pellicciotti, Bordoy and ImmerzeelRagettli and others, 2013), in contrast to the monsoon accumulation–ablation type to the east. The monsoon also influences the climate and is stronger in the south than in the north. High-elevation areas in the north of the Hunza River basin are significantly more arid than those in the south (Reference Ragettli, Pellicciotti, Bordoy and ImmerzeelRagettli and others, 2013).

Fig. 1. Hispar and Shimshal sub-regions of the Hunza River basin, Karokoram, northern Pakistan. Glacier area is shown in dark gray and debris geometry is black. Glaciers shown in light gray were not considered.

Glacier mass balances in this region have been shown to be stable or slightly positive (Reference HewittHewitt, 2005; Reference BolchBolch and others, 2012; Reference Gardelle, Berthier and ArnaudGardelle and others, 2012, Reference Gardelle, Berthier, Arnaud and Kääb2013; Reference Kääb, Berthier, Nuth, Gardelle and ArnaudKääb and others, 2012), suggesting an anomalous behavior in comparison to other High Mountain Asia regions; this has been termed the ‘Karakoram anomaly’ (e.g. Reference HewittHewitt, 2005; Reference Gardelle, Berthier, Arnaud and KääbGardelle and others, 2013). These studies, however, rely on observations from the past decade, and evidence is missing for previous periods. In an update of previous studies, Reference Gardelle, Berthier, Arnaud and KääbGardelle and others (2013) reported an average mass gain of 0.10 ± 0.16 m w.e. a⁻¹ for the region from 1999 to 2011. The balanced or slightly positive mass budget of Karokoram glaciers has been explained by the precipitation regime and topographic controls. Precipitation from winter westerlies can reach higher elevations than the summer monsoon (Reference Scherler, Bookhagen and StreckerScherler and others, 2011), and it has been suggested that precipitation is increasing. This, together with its maximum occurring between 5000 and 6000 m, has been proposed as an explanation for glacier expansion (Reference Archer and FowlerArcher and Fowler, 2004; Reference HewittHewitt, 2005, Reference Hewitt2011; Reference Fowler and ArcherFowler and Archer, 2006). As a result of these factors, glaciers in the Karokoram have high accumulation–area ratio (AAR) values (>60%), in contrast to values of 30–40% for the Nepal and Bhutan Himalaya (Reference Gardelle, Berthier, Arnaud and KääbGardelle and others, 2013).

A large number of glaciers in the Karakoram are surge-type (Reference CoplandCopland and others, 2011; Reference Gardelle, Berthier and ArnaudGardelle and others, 2012; Reference Rankl, Vijay, Kienholz and BraunRankl and others, 2013), with short active phases (of months to years) where glacier mass is transferred from a reservoir area to a receiving area and sometimes accompanied by a rapid advance of the terminus (Reference Meier and PostMeier and Post, 1969). An advance is followed by a quiescent phase (years to decades) with a stable front or retreat.

Many of the glaciers in the area are extensively debris-covered. The total glacierized area in the two sub-catchments was 1567 km² in 2014, and the elevation range covered by glaciers was 2339–7850 m.

Glacier delineation

Satellite images from the Landsat program were downloaded from http://www.earthexplorer.gov (Table 1). We preferentially selected images based on three criteria: (1) acquisition late in the melt season with the highest transient snowline elevation; (2) as few clouds in the scene as possible; and (3) a multi-year time interval so there is a greater likelihood that the debris-cover change signal is higher than the random error in the method. The glacier system was mapped manually and subdivided into individual glacier catchments based on visual identification of flow divides. Glacier outlines were manually digitized for 1977, 1998, 2009 and 2014 using the satellite scenes in Table 1. Errors are common for debris-covered glacier outlines, and while some solutions have been proposed for automated methods (e.g. Reference Paul, Huggel and KääbPaul and others, 2004; Reference Bolch and KampBolch and Kamp, 2006; Reference Shukla, Arora and GuptaShukla and others, 2010), manual delineation of debris-covered area is likely to be the most effective method to derive accurate and consistent glacier outlines (e.g. Reference Nagai, Fujita, Nuimura and SakaiNagai and others, 2013). Systematic errors in manual debris-cover mapping can be reduced if all of the digitization is done by the same person with a consistent definition of a glacier. Accumulation zones are challenging to map accurately due to the year-round presence of snow both on glacier ice as well as on the surrounding landscape. Because of this, we digitized glacier accumulation areas to the best of our ability using all of the scenes in Table 1. We then held the accumulation zone area constant for all four years, and only adjusted the glacier outlines in the ablation zone. True glacier area change within an accumulation zone is likely ambiguous and/or very slight.

Table 1. Landsat satellite scenes used for this paper, and the corresponding method and threshold value used to map debris cover. MSS: Multispectral Scanner; TM: Thematic Mapper; OLI: Operational Land Imager

Debris-covered area change was measured below an aggregate minimum transient snowline, so the full glacier area was used only when presenting debris-covered area as a percentage of total glacier area. At some locations it is difficult to identify a precise glacier boundary due to the resolution of Landsat sensors. For these cases, if there was no very obvious change in area between the years mapped we kept the same edge position with the intent of minimizing error and mapped only area change where we were confident there was a change in glacier geometry. Surge fronts also cause difficulties when generating glacier outlines because they can leave dead ice at their termini following a surge event. Glacier area change can therefore vary greatly, depending on the definition of a glacier used when generating outlines. Where the glacier front was ambiguous, we looked for evidence of glacier motion to define the terminus position. This is a more classical glacier definition but also will exhibit greater area changes that could be considered exaggerated in the situation of dead ice that is reactivated by a surge event.

Surge-type glacier identification

Several studies have identified surge-type glaciers in the Karokoram from optical data (e.g. Reference CoplandCopland and others, 2011; Reference Rankl, Vijay, Kienholz and BraunRankl and others, 2013; Reference Quincey and LuckmanQuincey and Luckman, 2014) and elevation difference data (e.g. Reference Gardelle, Berthier and ArnaudGardelle and others, 2012, Reference Gardelle, Berthier, Arnaud and Kääb2013). We used four Landsat scenes (Table 1) to inspect each glacier tongue and individual tributary branches for surge-type glacier characteristics and evidence of a surge event following surge-type behavior criteria from Reference Grant, Stokes and EvansGrant and others (2009). If a glacier (glacier tongue and tributaries) contained more than one surge classification, the classification with the greatest proportion of the total area was assigned to the entire glacier. This subjective method avoided the need to define glacier tributaries (surging and non-surging) for the entire region. Because surge behavior can influence changes in debris cover, we used the surge-type glacier inventory shown in Figure 2 to differentiate two populations for debris-cover change analysis. Glaciers identified as ‘possibly surge-type’ and in ‘quiescent phase’ were treated as not surge-type for this analysis.

Fig. 2. Distribution of surge-type glaciers with surge events constrained in time by the three observation periods used in this study. ‘Possibly surge-type’ glaciers may be advancing rather than surge-type.

Snowline and cloud map

Transient snowlines for all 93 glaciers were manually mapped from each of the four Landsat scenes listed in Table 1 and combined so that one aggregate snowline mask included snow-covered areas for all four scenes (Fig. 3). Similarly, cloud cover within the glacier domain was mapped manually for each of the four scenes and combined to create one aggregate mask to clip out of the debris-cover maps from each scene (Fig. 3). Both of these procedures produce a lower debris-cover area value than is physically present but allow an unbiased comparison between scenes. Figure 4 shows an idealized glacier with each of these components, and the resulting area that is applied to all glaciers considered in the study.

Fig. 3. Transient snowline from each image was manually mapped and combined to generate a single aggregate lowest snowline. The aggregate lowest snowline was used to define the upper-glacier edge of a spatial domain that enables a meaningful measure of debriscovered area change. Clouds, which would otherwise be erroneously automatically classified as debris cover, were manually mapped and the area of each cloud was removed from all scenes.

Fig. 4. For this study, debris cover is defined as the area mapped as debris cover that lies outside a composite of every cloud mapped within the satellite images from 1977, 1998, 2009 and 2014, and the lowest aggregate transient snowline for all four years.

Debris-cover map

Supraglacial debris cover was mapped using three different sensors from the Landsat program, each with slightly different spectral ranges (Table 2). Because the Landsat 2 Multispectral Scanner (MSS) did not cover wavelengths greater than 1.1 μm, high-contrast band ratios that can eliminate shadows and discriminate bare ice could not be constructed for the 1977 image. We experimented with many MSS band combinations and found |Band7–Band5| provides the greatest contrast for discriminating bare glacier ice and debris-covered areas, similar to Reference Rundquist, Collins, Barnes, Bussom, Samson and PeakeRundquist and others (1980). The Landsat 5 Thematic Mapper (TM) band ratio Band4/Band5 does eliminate the effect of shadowing and is routinely used to discriminate debris cover (Reference Paul, Huggel and KääbPaul and others, 2004). The Landsat 8 Operational Land Imager (OLI) has slightly different spectral ranges than the TM sensor, but the ratio of OLI Band5/Band6 is very similar to TM Band4/ Band5. When a threshold is applied to these band combinations, debris cover is automatically discriminated; however, an appropriate threshold value must be selected for each image.

Table 2. Landsat bands used in this study, and their corresponding wavelengths

Virjerab Glacier (Fig. 1) was selected to derive debris map threshold values for all four scenes. It was selected because of its large size relative to the surrounding glaciers, the presence of complete debris cover at its terminus, and a network of medial moraines that vary in width. As a moraine narrows, it will demonstrate the limitations of the method due to the spatial resolution of the satellite data. Debris cover on Virjerab Glacier was mapped manually from each of the four Landsat scenes and was, for the application of this method, considered the true debris-covered area. Automated debris-covered area was computed for a wide range of threshold values, and a function was fit to the point results. Applying the manually generated area as the independent variable in each of the four functions, an optimized threshold value was found (Fig. 5). These individual-scene threshold values were then applied to the entire study region and inspected for quality.

Fig. 5. Automatically derived debris cover for different threshold values (black dots), and manually derived debris cover (black line) for Virjerab Glacier. By fitting a function to the automatically derived debris maps, we found an optimized threshold value that produced an equal percentage of debris cover to the manual debris-cover map. For 1977 (a) we used the optimized threshold value derived here for some glaciers but not all due to the limited spectral range of the sensor (Fig. 6). A second-order polynomial was used to fit the 1998 and 2009 data ((b) R ² = 0.99 and (c) R ² = 0.99, respectively), and a third-order polynomial fit was used for the 2014 data ((d) R ² = 0.99). The threshold values derived from these plots are given in Table 1.

We considered the optimized threshold value found at Virjerab Glacier for Landsat 5 TM and Landsat 8 OLI data to produce successful results when applied to the surrounding glaciers, but while the optimized Landsat 2 MSS threshold was adequate for Virjerab Glacier, visual examination of the value applied to other glaciers showed many errors. Short of manually digitizing the debris cover for every glacier, we took a less rigorous approach than we applied for Virjerab Glacier and selected an acceptable threshold value visually for each glacier. Figure 6 shows the threshold values we applied, with a few examples illustrating these choices.

Fig. 6. The use of one threshold for all of the glaciers in the 1977 image produced poor results for some glaciers. To resolve this, we manually selected the best value of a range of threshold values for each glacier. (b) shows Virjerab Glacier where both the threshold derived in Figure 5 (55) and a threshold of 45 do reasonably well. However, (a) and (c) show examples where debris cover is better delineated using a threshold of 45. The threshold distribution used for the remainder of this paper is shown.

Automated debris-mapping techniques can produce a high number of mapped debris patches that are only one to a few pixels in size. In principle, if the pixel(s) value falls within the defined debris-cover threshold, it should be considered debris, yet these small and isolated patches are often from error within the glacier outline (e.g. a bedrock pixel included within the glacier outline). Additionally, enclave patches of bare ice within the automated debris map can be misclassified by the presence of supraglacial lakes. For our large-scale area change study we elected to fill these holes and simply consider supraglacial lakes as debris-covered area. Debris patches and holes with an area of 2700 m² (three 30 m pixels) or less were removed or filled automatically, and holes larger than 2700 m², yet clearly sourced from supraglacial lakes, were filled manually.

Error estimation

The quality of both the glacier outlines and debris-cover area is critical for meaningful results, especially when data from different sensors with different spectral ranges and resolutions are required to be compatible in accuracy. The method we used to solve for a debris/bare-ice threshold for Virjerab Glacier optimizes on total debris-covered area alone, not its spatial distribution. However, for determining the accuracy of the automated debris geometry, the spatial distribution is important. By design, the method will contain two error quantities that will be effectively equal in magnitude: debris-covered area that the automated routine missed and bare-ice area that the automated routine classified as debris-covered. These quantities will be equal because the routine overcompensates for area missed by including additional area with only slight variations due to the inclusion or exclusion of bifurcated pixels. By summing the two terms and dividing by the spatial domain used during automated classification (with area above the transient snowline and clouded area removed) a percent error term can be applied to the surface-type classification method. Dividing the summed error terms by the total glacier area gives a percent error term for the glacier-wide debris map. For the measure of debris-covered area change between two images, the glacier-wide error estimates were summed in quadrature to find debris-covered area change error. We extrapolated these error values with the assumption that the errors observed at Virjerab Glacier are representative of the other glaciers included in this study.

Debris-covered area change

Glaciers in the Shimshal and Hispar valleys (Fig. 1) were 20.9 ± 4.1% debris-covered in 1977, and 21.1 ± 1.8% debris-covered in 2014. Total glacier area remained nearly constant, changing from 1573 to 1567 km² from 1977 to 2014.

Figure 7 shows the resulting debris maps as a change in percent debris cover for each glacier (%_{deb time2} − %_{deb time1}, as defined in Fig. 4) and as the actual debris-covered area geometry from both of the years that are compared. Where the debris is shown as light-gray or black, there has been a removal or addition, respectively, of debris at that location. There are several mechanisms that could produce a change from debris-free to debris-covered, or from debris-covered to debris-free for a given location on a glacier: (1) a translated feature that is not parallel to the glacier flowline (e.g. down-glacier progression of a landslide that was deposited onto the ice or a surge loop); (2) new debris accumulating at the surface either through englacial melt-out or extraglacially sourced mass transportation; (3) glacier dynamics where longitudinal strain can cause extension or compression of a debris cover or the reactivation of debris-covered dead ice (e.g. during a surge event); or (4) debris that is removed from the surface either through evacuation at the terminus, tumbling off the side or englacial transportation (e.g. by falling into a crevasse or removal by a hydrological process down a conduit). While translated features that are not parallel to the glacier flowline (e.g. a wavy moraine produced from glacier flow instabilities) reveal dynamic information, our analysis relies on summed debris-covered area change over a single glacier, eliminating any debris-covered area change from features that are translated without a change in area. Debris-covered area change from extensional or compressional strain associated with feature translation is accounted for.

Fig. 7. Supraglacial debris-cover area change over 37 years, derived from four Landsat satellite images. Glaciers that are white have no data (see Fig. 10). Boxes a–d in the 1977–98 map are detailed in Figure 11.

Figure 8 presents glacier-wide debris-covered area and total glacier area change. Non-surging glaciers (59% of total glacierized area) and surge-type glaciers (41% of total glacierized area) are considered separately for each observation period. Our results had little variance between glaciers of different sizes, so we did not further subdivide the presentation of our data based on glacier surface area. The debris-cover change results are presented in three ways: (1) as a change in the percentage of debris-covered area for each glacier; (2) as a rate of change in the percentage of debris-covered area for each glacier; and (3) as a rate of change of debris cover (km²) for each glacier. Figure 9 shows how glacier-wide error estimates were obtained for Virjerab Glacier which were assumed to be applicable to the other glaciers in this study and are included in Figure 8a. Error in the 1977 debris map nearly doubles the error estimated for the other satellite scenes, but this is likely explained by the 60 m, rather than 30 m, satellite image pixel resolution.

Fig. 8. Debris-cover and glacier area change for all glaciers. The data in each panel are presented as a pair of box plots for each time interval, the first for non-surge-type (n = 62; 59% of total glacierized area) and the second, colored blue, for surge-type glaciers (n = 10; 41% of total glacierized area). For each box plot, the red line is the median, the edges of the box are the 25th and 75th percentiles, the whiskers extend to 3 × IQR (interquartile range) and outliers are plotted individually. (a) Change in glacier percentage that is debris-covered (as defined in Fig. 4). The number in parentheses above each box plot containing surge-type glaciers is the number of active surges that took place during that time interval. The gray area shows the error range explained in Figure 9. (b) Change in percent debris cover per year. (c) Change in square kilometers per year. (d) Change in glacier area. Glaciers where <10% of the debris cover was mapped (glaciers that appear red in Fig. 10 (n = 21), which occupy ∼70 km² of the glacierized area) are excluded from this figure.

Fig. 9. Automated debris map error calculated for Virjerab Glacier for each of the four Landsat scenes used in this study. By the design of the method, the two error terms, (1) debris-covered area missed and (2) bare ice erroneously classified as debris-covered, will be roughly equal. The automated debris algorithm was applied below the transient snowline (gray line), and the percent error of the algorithm (ε _algorithm) is found by summing the two error terms and dividing by the area below the gray line. Because the area above the snowline can be positively classified as not debris-covered, a glacier-wide error term (ε _glacier-wide) can be found by dividing the summed error over the entire glacier area.

Our summed regional results (Table 3) suggest there was no significant change in glacier area or debris-covered area for surging and non-surging glaciers from 1977 to 2014.

Table 3. Region-wide changes in debris-covered surface area and glacier surface area between 1977, 1998, 2009 and 2014. These results present a summation of individual glacier values for all of the glaciers studied and pertain to 1502 km² of glacierized area (value from 1977) which is 95.5% of the total area initially considered. Glaciers where <10% of the debris cover was mapped (glaciers that appear red in Figure 10 (n = 21), which occupy about 70 km² of the glacierized area) were excluded. Our area change estimates are not confident beyond one decimal place but where a rate value would otherwise round to zero we show the result to 10⁻². Debris-covered area change errors were calculated in Figure 9 and assumed constant for all glaciers

Excluded area

Our analysis was limited to the glacier area below the aggregate minimum transient snowline and aggregate cloud mask shown in Figure 3 and explained in Figure 4. Figure 10 gives an estimate of how comprehensive our findings are for the entire region. The 1998 image had almost no clouds and the highest transient snowline for most glaciers, so a debris map generated from this image with no imposing domain restrictions from other scenes is likely close to incorporating the maximum debris cover exposed at the glacier surface below the true equilibrium-line altitude (ELA) for that year. The difference between the 1998 debris map and the restricted domain for the same image used in this study gives per-glacier estimates of the area either captured or missed for our debris analysis. Many glaciers had debris-covered areas that were only slightly reduced (blue in Fig. 10), suggesting that for these glaciers our methods provide meaningful glacier-wide estimates of debris-cover change. For some glaciers, transient snow covered the glacier completely during the acquisition of at least one of the images. Twenty-one glaciers, comprising ∼70 km² or 4.5% of the total glacierized area (red in Fig. 10), were limited to <10% of the debris-covered area after considering transient snow cover and were excluded from this study.

Fig. 10. The difference between total debris-covered area present in the late melt season, low cloud cover, 1998 Landsat 5 scene and the reduced spatial domain used for this study applied to the same scene. This illustrates a per-glacier estimate of the debris-covered area either captured or missed for our analysis.

Debris-cover evolution

Our results of little or zero change in debris-covered area from 1977 to 2014 mirror the stable or slightly positive mass balances observed in the Karakoram (Reference Gardelle, Berthier and ArnaudGardelle and others, 2012, Reference Gardelle, Berthier, Arnaud and Kääb2013; Reference Kääb, Berthier, Nuth, Gardelle and ArnaudKääb and others, 2012), suggesting that, in a stable mass-balance regime, debris-covered area may also remain stable. While our results provide an independent dataset that corroborates the Karakoram anomaly, we can further extend this coupling to suggest that the anomalous, stable or slightly positive mass balances within the region extend further back in time than is shown in previous work (Reference Gardelle, Berthier and ArnaudGardelle and others, 2012, Reference Gardelle, Berthier, Arnaud and Kääb2013; Reference Kääb, Berthier, Nuth, Gardelle and ArnaudKääb and others, 2012) which relies on observations from the most recent decade. It is interesting that the glaciers that surged, which are clearly not dynamically stable, also exhibited a long-term average of stable debris-covered area.

If we hypothesize that during a stable mass-balance regime debris-covered area does not significantly increase, at what point did these glaciers become 21% debris-covered? Reference Shroder, Bishop, Copland and SloanShroder and others (2000) outline characteristics that may control the transition of a debris-covered glacier to a rock glacier. Similarly, Reference DelineDeline (2005) used historical data to reconstruct the transition from mostly debris-free glaciers during the Little Ice Age to what are essentially rock glaciers at the present time. Both corroborate the model proposed by Reference KirkbrideKirkbride (2000): an inefficient, ablation-dominated setting will promote an increase in debris cover. This suggests that during this study we observed either one or a combination of two things: (1) a plateau of debris-cover evolution where 21% debris cover is the steady sum of the transportation rate of rock material to the glacier surface and the level of debris-cover evacuation efficiency, or (2) 37 years is an insufficient window of time to detect supraglacial debris-cover changes at this location from Landsat imagery.

We hypothesize that the zero or near-zero change measured in the Hispar and Shimshal valleys is unique to an area that is in a stable glacier mass-balance regime, but this hypothesis can only be tested by applying this method to adjacent regions with negative mass balances and elsewhere in the world. Our analysis suggests the findings for one glacier might be specific and not representative of average glacier trends, so it is preferable to rely on regional analyses to provide a more representative picture.

Other studies have obtained an increase in debris area in regional settings of negative mass balances (Reference KirkbrideKirkbride, 1993; Reference Popovnin and RozovaPopovnin and Rozova, 2002; Reference DelineDeline, 2005; Reference Stokes, Popovnin, Aleynikov, Gurney and ShahgedanovaStokes and others, 2007; Reference Kellerer-Pirklbauer, Lieb, Avian and GspurningKellerer-Pirklbauer and others, 2008; Reference Mazué, Deline and KirkbrideMazué and others, 2009; Reference Quincey and GlasserQuincey and Glasser, 2009; Reference Shukla, Gupta and AroraShukla and others, 2009; Reference LambrechtLambrecht and others, 2011; Reference Kirkbride and DelineKirkbride and Deline, 2013). It is commonly assumed that as a glacier loses volume in a warming climate, the often steep valley walls that were once supported by the ice mass are more likely to become unstable and fail, increasing extraglacial rock deposition onto the glacier. Additionally, sustained melt not counterbalanced by accumulation will lead to an increased exposure of englacial debris (Reference Kirkbride and DelineKirkbride and Deline, 2013). Both mechanisms will result in an increase in supraglacial debris cover.

Reference LambrechtLambrecht and others (2011) obtained an increase in debris cover from 16% to 23% from 1991 to 2006 for glaciers in the Adyl-su valley in the Caucasus, but no changes between 1971 and 1991. Very small changes (from 6.2% to 8.1%) were obtained for the glaciers in an adjacent valley. Reference Kellerer-Pirklbauer, Lieb, Avian and GspurningKellerer-Pirklbauer and others (2008) found a small increase in debris-covered area for Pasterze glacier, Austria: from 5.4% in 1964 to 5.5% in 1981 and 7.3% in 2000. Reference Mazué, Deline and KirkbrideMazué and others (2009) estimated a 1.25% a⁻¹ increase in debris cover for Estelette glacier, Italy. Reference Shukla, Gupta and AroraShukla and others (2009) measured a 76.5% increase in debris cover from 2001 to 2004 on Samudratapu Glacier, Himachal Himalaya, but it is important to note that transient snow cover may not have been accounted for. Reference Bhambri, Bolch, Chaujar and KulshreshthaBhambri and others (2011) derived an increase in debris-cover area of 17.8 ± 3.1% and 11.8 ± 3.0% in two basins in the Garwhal Himalaya, India. These results from Reference Shukla, Gupta and AroraShukla and others (2009) and Reference Bhambri, Bolch, Chaujar and KulshreshthaBhambri and others (2011) cannot be directly compared to our results or those previously mentioned because they calculated debris-covered area change as a difference between debris-covered area at two times without considering the change in glacier geometry. In a warming climate, a decrease in glacier area is likely to coincide with an increase in debris-covered area. Thus, a more meaningful measure of debris-cover change would be the change in debris as a percentage of total glacier area. In this case, the values obtained by Reference Shukla, Gupta and AroraShukla and others (2009) and Reference Bhambri, Bolch, Chaujar and KulshreshthaBhambri and others (2011) would be smaller. In summary, evidence of the assumed increase in debris cover in a warming climate is still sparse, and limited to a few glaciers. Regional studies in areas of the world that differ in terms of climate and glacier characteristics will be important to establish sound evidence to support a hypothesis that could have a substantial impact on future glacier mass balance and runoff modeling (Reference Jouvet, Huss, Funk and BlatterJouvet and others, 2011; Reference Pellicciotti, Carenzo, Bordoy and StoffelPellicciotti and others, 2014).

While this study has focused on measuring two-dimensional debris-covered area evolution, it is also critical to resolve the thickness of the debris cover and its evolution through time. The coupling of debris-covered area and thickness evolution both from measurements and further modeling is needed to enable meaningful simulations of mountain glacier geometry and melt. Further methodological advancements of satellite-based debris-cover thickness estimation to expand the spatial domain to include many glaciers and enable the observation of meaningful changes in thickness by successive measurements through time are a feasible next step (Reference Foster, Brock, Cutler and DiotriFoster and others, 2012).

Glaciers that exhibit rapid changes in debris cover

While the average change in debris-covered area for the whole region is zero or near zero, specific glaciers exhibited outlying, nonzero changes in debris cover. We selected four glaciers that host outlier behavior to serve as examples of some sources of rapid change in debris-covered area (Fig. 11). All of the glaciers shown in Figure 11 are surge-type but some of the phenomena are not strictly surge-related (e.g. high-volume, low-frequency rock avalanches or landslides onto a glacier surface). We hypothesize that during a surge event accumulated supraglacial debris cover is reduced either by rapid evacuation through the glacier terminus (e.g. Fig. 11a) or by englacial entrainment from extensive crevassing (e.g. Fig. 11d). After a surge, crevasses close, shallow englacial debris will quickly melt out to the surface and quiescent-phase glacier velocities will facilitate reaccumulation of debris cover from englacial and extraglacial sources and cause both these positive and negative changes to cancel out. Another critical element, illustrated in Figure 11c, is whether stagnant ice is included as part of a glacier. While stagnant ice is included in the Global Land Ice Measurements from Space (GLIMS) definition of a glacier (Reference Raup and KhalsaRaup and Khalsa, 2007), we found the classical definition of a glacier, requiring deformation from glacier flow, to be more desirable and less ambiguous to implement. However, problems arise with the classical definition when stagnant ice becomes reactivated and glacier mass is gained through accretion. Provided the glacier definition used is consistent and the affected glaciers are identified, interpretation of the data should be unobstructed.

Fig. 11. Four examples illustrate instances of rapid debris-cover change. (a) A classic surge event of an unnamed glacier where debris cover is almost completely removed from the surface, then debris cover begins to reaccumulate in 2014. (b) An unnamed glacier, showing the addition of supraglacial rock avalanche debris (identified with a white arrow). (c) Gharesa Glacier, showing a situation where stagnant ice not considered part of the glacier (identified with a white arrow) was reactivated by a surge event (black area = aggregate cloud coverage). (d) A tributary branch of Hispar Glacier, showing another surge-related phenomenon, where formerly debris-covered area becomes debrisfree as crevasses open and transports supraglacial rocks to an englacial environment. By 2014 the area is again debris-covered.

Debris-cover expansion into the accumulation zone

In order to maintain a measure of actual debris-cover expansion or reduction, we held the upper bounds of the spatial domain constant at the lowest minimum transient snowline. This method does not allow observations of debris-cover expansion through up-glacier migration of the ELA. This is certainly a relevant process in debris-cover evolution, but we did not have sufficient data to identify the ELA relative to the transient snowline and to map the debris cover up to the ELA for each time step. Conversely, if debris cover is mapped to the transient snowline on multiple images and compared, possibly very erroneous changes in debris cover will be found because it is unknown whether the changes are from a physical addition or reduction of rock material or whether the debris mapped from one image is also present in the other but undetectable from snow cover. By holding a constant lowest minimum transient snowline we might underestimate debris-covered area and we cannot account for debris expansion by means of accumulation zone contraction. However, Figure 10 suggests that, for many glaciers, our method considered 90% or more of the complete debris-covered area (including area above our transient-snowline restricted domain) and, due to the stable mass balances observed within the study region, ELA migration has likely been slight or negligible from the 1970s to the present.