Satellite observations show no net change in the percentage of supraglacial debris-covered area in northern Pakistan from 1977 to 2014

Abstract

Spatial evolution of supraglacial debris cover on mountain glaciers is a largely unmonitored and poorly understood phenomenon that directly affects glacier melt. Supraglacial debris cover for 93 glaciers in the Karakoram, northern Pakistan, was mapped from Landsat imagery acquired in 1977, 1998, 2009 and 2014. Surge-type glaciers occupy 41% of the study area and were considered separately. The time series of debris-covered surface area change shows a mean value of zero or near-zero change for both surging and non-surging glaciers. An increase in debris-covered area is often associated with negative regional mass balances. We extend this logic to suggest that the stable regional mass balances in the Karakoram explain the zero or near-zero change in debris-covered area. This coupling of trends combined with our 37 year time series of data suggests the Karakoram anomaly extends further back in time than previously known.

Introduction

Debris-covered portions of a glacier alter surface energy fluxes relative to bare ice and can have a significant impact on total glacier melt (Østrem, 1959; Nakawo and Rana, 1999; Reid and Brock, 2010). During recent decades, mapping supraglacial debris cover and quantifying its effect on glacier melt has been identified as an important component of monitoring and modeling mass balances of mountain glaciers (e.g. Nakawo and others, 2000; Lejeune and others, 2013; Pellicciotti and others, 2014). Before the effect of debris cover on glacier melt can be included in a distributed glacier melt model, the location of the debris must first be known, and, in order to resolve glacier volume change through time, measured or simulated debris-covered area should also evolve with time (Jouvet and others, 2011).

While several methods exist for mapping debris cover (summarized by Shukla and others, 2010), the method most suitable for generating a regional-scale time series with sufficient longevity is given by Paul and others (2004) and has been used extensively for mapping debris cover at a single point in time (e.g. Scherler and others, 2011; Bolch and others, 2012; Frey and others, 2012; Kienholz and others, 2015). The method uses multispectral satellite images and a band ratio to discriminate debris-covered ice from clean ice and can be automated, allowing application to a large sample of glaciers.

The timescale needed to observe supraglacial debris-covered area change will vary with setting and data resolution, but is largely unknown. Anderson (2000) provides a theoretical description of medial moraine evolution, with a timescale of 150 years for a medial moraine to widen by a factor of three in a stable regional climate. Deline (2005) reconstructed debris evolution of three glaciers in the Mont Blanc range from 1770 to 1940 during which a ‘clean’ glacier effectively became a rock glacier. Jouvet and others (2011) presented a forward model of Grosser Aletschgletscher, Switzerland, where the current 4% debris cover expands to complete coverage of the terminus by 2080 using their highest debris propagation parameter. Landsat satellites have been acquiring data continuously since 1972 and it is plausible that this duration offers a sufficient timescale to map debris-cover evolution (Stokes and others, 2007), provided a method can produce debris-covered area maps of comparable quality over two or more scenes.

Several studies have measured debris-cover evolution (Kirkbride, 1993; Popovnin and Rozova, 2002; Deline, 2005; Stokes and others, 2007; Kellerer-Pirklbauer and others, 2008; Mazué and others, 2009; Quincey and Glasser, 2009; Shukla and others, 2009; Lambrecht and others, 2011) yet focus on only one or a few glaciers. Analysis at a regional scale provides spatial context and reduces the risk of extrapolating anomalous behavior, but studies of debris-cover evolution at this scale are very limited (Bhambri and others, 2011). All the observations and models of debris-cover evolution mentioned were conducted or constructed in a climate that promotes glacier shrinkage (with the exception of Anderson, 2000). Nearly every glacier that has been studied shows an increase in debris-covered area, which is expected within a negative regional mass-balance regime and reduced-flow, or ‘ablation-dominant’, conditions (Kirkbride, 2000; Kirkbride and Deline, 2013).

For this study, we selected a subset of 93 glaciers within the Karakoram, Pakistan, where glaciers presently show a regional-average slight mass gain or stable mass balances (Gardelle and others, 2012, 2013; Kääb and others, 2012) and several glaciers exhibit surge-type behavior (Copland and others, 2011; Gardelle and others, 2012; Rankl and others, 2013). Clapperton (1975) identifies several factors that might cause a surge-type glacier to incorporate debris cover differently than non-surging glaciers, yet the interplay between surge behavior and debris cover is largely unknown.

The objectives of this study are to (1) propose a method that produces regional debris maps from different Landsat satellites, that are within an acceptable quality margin to resolve changes over time, (2) investigate how debris cover has evolved through time in a region with slightly positive or stable mass balances and (3) consider if surge-type glaciers exhibit unique behavior with respect to glacier-wide changes in debris cover.

Study Area and Methods

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 (Gardelle and others, 2013; Ragettli 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 (Ragettli 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 (Hewitt, 2005; Bolch and others, 2012; Gardelle and others, 2012, 2013; Kää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. Hewitt, 2005; Gardelle 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, Gardelle 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 (Scherler 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 (Archer and Fowler, 2004; Hewitt, 2005, 2011; Fowler 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 (Gardelle and others, 2013).

A large number of glaciers in the Karakoram are surge-type (Copland and others, 2011; Gardelle and others, 2012; Rankl 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 (Meier 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. Paul and others, 2004; Bolch and Kamp, 2006; Shukla 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. Nagai 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. Copland and others, 2011; Rankl and others, 2013; Quincey and Luckman, 2014) and elevation difference data (e.g. Gardelle and others, 2012, 2013). 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 Grant 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 Rundquist 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 (Paul 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.

Results

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.

Discussion

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 (Gardelle and others, 2012, 2013; Kää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 (Gardelle and others, 2012, 2013; Kää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? Shroder and others (2000) outline characteristics that may control the transition of a debris-covered glacier to a rock glacier. Similarly, Deline (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 Kirkbride (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 (Kirkbride, 1993; Popovnin and Rozova, 2002; Deline, 2005; Stokes and others, 2007; Kellerer-Pirklbauer and others, 2008; Mazué and others, 2009; Quincey and Glasser, 2009; Shukla and others, 2009; Lambrecht and others, 2011; Kirkbride 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 (Kirkbride and Deline, 2013). Both mechanisms will result in an increase in supraglacial debris cover.

Lambrecht 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. Kellerer-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. Mazué and others (2009) estimated a 1.25% a⁻¹ increase in debris cover for Estelette glacier, Italy. Shukla 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. Bhambri 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 Shukla and others (2009) and Bhambri 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 Shukla and others (2009) and Bhambri 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 (Jouvet and others, 2011; Pellicciotti 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 (Foster 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 (Raup 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.

Conclusions

We have presented a method that uses satellite imagery spanning three generations of Landsat sensors to investigate supraglacial debris-cover evolution on a regional scale. We applied this method to investigate debris-covered area change for 93 glaciers in the Karakoram. Our results show zero or near-zero change in debris-covered area, a finding that is unique to similar studies yet fits with the regional stable or slightly positive mass balances observed in the Karakoram. These findings suggest that, from a spatial perspective neglecting changes in debris thickness, both the input of rock material to, and its removal from, the glacier system have been in an equilibrium state from 1977 to 2014. The coupling of stable mass balances and stable debris-covered area trends suggests the Karakoram anomaly persisted further back in time than previously documented. Despite the zero change trend observed in the majority of the analyzed glaciers, a number of outliers showed abrupt changes due to particular events (e.g. surge events and landslides). Surge-type glaciers present a wider distribution of changes related to their unstable ice dynamic, but over the complete 37 year observation period we did not find significant differences between surging and non-surging glaciers in terms of total glacier area and debris-covered area trends. Further studies of this type in other areas in the world will reveal whether the stable debris-covered area changes measured here are unique or if regions of negative glacier mass balances can also exhibit similar values.

Author Contribution Statement

S.H. developed the method, performed all of the analysis and wrote the paper. F.P. initiated the study, provided guidance and feedback throughout the work and contributed to the writing. A.A. contributed to the writing and provided comments on both the method and the paper. C.K. helped write code and, together with A.C. and J.S., provided comments on the paper. A.S. helped secure funds for this project.