Southwest-facing slopes control the formation of debris-covered glaciers in the Bhutan Himalaya

To understand the formation conditions of debriscovered glaciers, we examined the dimension and shape of debris-covered areas and potential debris-supply (PDS) slopes of 213 glaciers in the Bhutan Himalaya. This was undertaken using satellite images with 2.5 m spatial resolution for manual delineation of debris-covered areas and PDS slopes. The most significant correlation exists between surface area of southwest-facing PDS slopes and debriscovered area. This result suggests that the southwest-facing PDS slopes supply the largest quantity of debris mantle. The shape of debris-covered areas is also an important variable, quantitatively defined using a geometric index. Elongate or stripe-like debris-covered areas on north-flowing glaciers are common throughout the Bhutan Himalaya. In contrast, south-flowing glaciers have large ablation zones, entirely covered by debris. Our findings suggest that this difference is caused by effective diurnal freeze–thaw cycles rather than seasonal freeze–thaw cycles, permafrost degradation, or snow avalanches. In terms of geographic setting, local topography also contributes to glacier debris supply and the proportion of debris cover on the studied glaciers is suppressed by the arid Tibetan climate, whereas the north-tosouth asymmetric topography of the Bhutan Himalaya has less influence on the proportion of debris cover.


Introduction
Glaciers in the Himalayas and other mountain regions are valuable as seasonally sustainable water reservoirs for downstream agricultural regions (e.g.Immerzeel et al., 2010;Kaser et al., 2010), and they are recognized as sensitive indicators of climate change (e.g.Kaser et al., 2006;Bolch et al., 2012).Many of the larger Himalayan glaciers have ablation zones covered by a mantle of unconsolidated sediments, ranging in size from dust to boulders.These are derived from rockfalls and avalanches on mountain slopes, or from subglacial erosion at the ice-bedrock interface (e.g.Moribayashi and Higuchi, 1977;Fushimi et al., 1980;Nakawo et al., 1986;Hambrey et al., 2008).Experimental field studies have revealed that thinner debris layers accelerate ice melting, whereas thicker layers potentially reduce melting through insulation (e.g.Østrem, 1959;Nakawo andYoung, 1981, 1982;Mattson et al., 1993).For instance, thick debris coverage induced by a landslide can drastically change glacier mass balance and flow velocity (Shugar et al., 2012).However, topographically rugged surfaces are also known to cause considerable ice melting in debris-covered areas, where numerous ice cliffs and supra-glacial ponds effectively absorb and transport the solar radiation for ice melting (Sakai et al., 2000(Sakai et al., , 2002)).Although stable fronts of heavily debris-covered glaciers have been identified (Scherler et al., 2011a), regionally averaged surface-lowering rates of Himalayan glaciers are similar for debris-covered versus debris-free areas (Kääb et al., 2012).Significant lowering rates on debris-covered surfaces have been observed during other remote sensing studies (Berthier et al., 2007;Bolch et al., 2011;Nuimura et al., 2012) and in situ measurements (Nuimura et al., 2011).To understand regionally heterogeneous shrinkage of glaciers throughout the Himalayas (Fujita and Nuimura, 2011), it is therefore necessary to recognize and understand the coexistence of these aspects on debris cover.
Spatially heterogeneous debris cover also has potentially important implications for risk assessment of glacial lake outburst floods.Glacial lakes are larger and more numerous in the eastern Himalaya, such as Nepal and Bhutan, compared with the western Himalaya (Gardelle et al., 2011).In the Nepal Himalaya, glaciers with thinner debris cover typically have lakes at their termini (Suzuki et al., 2007).Inclination and surface lowering of the debris-covered areas are probably associated with lake formation (Sakai and Fujita, 2010;Salerno et al., 2012).In Bhutan Himalaya, glacial lakes on the southern side expanded faster than those on the northern side (Komori, 2008), which is probably associated with topography (Kääb, 2005).
To understand the formation conditions of debris-covered glaciers, therefore, it is necessary to examine their topographic setting.Scherler et al. (2011b) suggested that the ratio of debris-covered area to total glacier area is related to the steepness of the ice-free zone above the snow line.This is based on the hypothesis that steeper ice-free areas tend to cause avalanches and therefore supply more debris to the glacier.On the other hand, it has been pointed out that asymmetric shape of terminal moraines were strongly correlated with aspects of source walls (Benn, 1989).In the Kangchenjunga Valley of the Nepal Himalaya, frequency of rockfall events is potentially influenced by slope orientation and diurnal fluctuations in surface temperatures; fewer rockfalls occur on north-facing slopes because small diurnal fluctuations in surface temperatures minimize the freeze-thaw cycle (Regmi and Watanabe, 2009).Previous studies of debriscovered glaciers have paid little attention to environmental controls on the extent and shape of debris-covered areas; however, recent findings by Regmi and Watanabe (2009) suggest that debris supply rate is controlled by surface conditions of the headwall where the debris mantle originates.
In this study, we examine the formation conditions of debris-covered glaciers by focusing on the relationship between the size of debris-covered areas and topographic variables in the areas surrounding glaciers in the Bhutan Himalaya.Previous research has undertaken direct measurements of rockfall events or temperature variation in small glacial basins (e.g.Rapp, 1960;Hewitt, 1968;Thorn, 1979;Shiraiwa, 1992;Regmi and Watanabe, 2009).In this study, satellite-derived remote sensing data are used to statistically analyse geomorphic information across the entire mountain range.Specifically, we analyse topographic and thermal features of slopes surrounding debris-covered glaciers to understand the supply of debris to the glacier and its influence on the extent of the debris-covered area.

Study area
We examine debris-covered glaciers broadly distributed along the entire Bhutan Himalaya, at elevations of 4000 m to 7500 m above sea level (a.s.l.) (89 • 12 -92 • 00 E, 27 • 36 -28 • 30 N; Fig. 1).In this region, numerous (> 1000) debris-free glaciers are also distributed; however we do not deal with them in this study.Our study area includes the headwaters of the Brahmaputra River and several inland rivers in southern Tibet.Glaciers have been retreating throughout this region in recent decades (Karma et al., 2003).

Satellite data
Spatial outlines of debris-covered glaciers and related slopes were delineated using 2.5 m-resolution images produced by the Panchromatic Remote-sensing Instrument for Stereo Mapping (PRISM), housed on-board the Advanced Land Observing Satellite (ALOS).A total of 32 panchromatic images of PRISM (Level 1B1) were used in this study.These images were orthorectified with a PRISM-derived digital surface model using Ortho-image Generation Software for ALOS PRISM (Tadono et al., 2012), and they cover the entire study area, apart from several small glaciers on the Tibetan side (Fig. 2).To achieve the most accurate delineation, we compared PRISM images taken on different dates and selected the image with the least snow cover.When no PRISM images were available due to snow or cloud cover, we used images of the Advanced Visible and Near Infrared Radiometer type 2 (AVNIR-2) with 10 m spatial resolution on-board ALOS (Level 1B1), orthorectified by a digital elevation model derived from the Shuttle Radar Topography Mission (SRTM) (Fig. 2).All polygons of the delineated glaciers were overlaid on bird's-eye view images in Google Earth ™ , and delineation error was modified.
To obtain the topographic contour lines with which to delineate glacier outlines, we used the second version of the Global Digital Elevation Model generated from the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER GDEM-2).This model was released by the National Aeronautics and Space Administration of the United States of America, and the Ministry of Economy, Trade, and Industry of Japan (http://www.jspacesystems.or.jp/ersdac/GDEM/E/index.html) (Tachikawa et al., 2011).We chose to use this model instead of SRTM due to our focus on high and steep mountain localities where SRTM potentially has voids (e.g.Luedeling et al., 2007;Reuter et al., 2007).ASTER GDEM-2 data were created from multiple ASTER images between 2000 and 2010, and have a cell size of 1 arc second (∼ 30 m).
ASTER-derived spatial distributions of surface temperature (ASTER Level 2B03 product) were obtained for part of the study region.The Earth Remote Sensing Data Analysis Center has completed atmospheric correction and temperature/emissivity separation on these data.We used neighbouring scenes covering the Kulha Kangri Massif region (white rectangle in Fig. 1), with two continuous scenes taken in daytime (1 January 2008, 10:40 and 10:41 Bhutanese Standard Time (BST)) and one in nighttime (27 January 2005, 22:09 BST) (Fig. 2).To overcome topographically induced distortions in these data, two orthorectified thermal infrared (TIR) images of ASTER Level 3A01 product, which are produced from the same original ASTER data, were also utilized.Overlaid on the orthorectified images, topographically induced distortions on the daytime surface-temperature distribution were manually rectified by six tie-points on the mountain ridges with root mean square error (RMSE) of 147.7 m.It was not possible to apply the same procedure for the nighttime image because an orthorectified TIR image was unavailable.However, its distortion is considered to be on the same order of magnitude as the above-mentioned RMSE, less than 2 pixels of ASTER TIR sensors (180 m).

Delineation of glaciers, debris-covered areas, and potential debris-supply slopes
No promising algorithm has yet been developed for automated mapping of debris-covered glaciers (Racoviteanu et al., 2009;Paul et al., 2013), and although the Global Land Ice Measurements from Space (GLIMS) provides a global glacier database, its use of relatively low-resolution ASTER VNIR sensors (15 m) has resulted in identification of incorrect boundaries between glacier ice and ice-free surfaces covered by snow in the Bhutan Himalaya.As such, it was necessary to manually delineate glaciers on PRISM and AVNIR-2 images using the standard definition of glacier outlines, following Raup and Khalsa (2007) and Rastner et al. (2012).
Visual interpretation of higher resolution PRISM images enabled exclusion of snow-covered terrain, whilst gradational shading and crevasses were used to identify glacier ice.In the uppermost snowfield, there was less contrast between individual glacier flows, which made it difficult to identify glacier boundaries.Contour lines were thus derived from the ASTER GDEM2 and utilized to divide glaciers.Debris-covered areas were identified using PRISM images (2.5 m resolution) and Google Earth ™ high-resolution images.Rugged, small-relief supra-glacial ponds and ice cliffs are typical features for a relatively thick debris cover and are not found on lateral moraines or bedrock terrains around glaciers.The upper boundary of debris cover often varies in several PRISM images taken on different dates.In such cases, we adopted the maximum outline by considering temporal snow cover.
We define potential debris-supply (PDS) slopes that could supply the debris mantle.Such slopes were manually extracted as the continuous slope between the glacier margin and mountain ridges, downward to the glacier, using contour lines at a 20 m interval and a slope map derived from ASTER GDEM-2 (Fig. 3).More precise interpretation is given by the PRISM images.Rockfall-or avalanche-induced scars and traces are identified as rockfall possibility (Fig. 3).Slopes intercepted by hills or lateral moraines were excluded from the PDS slopes because rockfall would be prevented.Slopes with temporal snow cover were included to account for debris-supply induced by snow avalanches (e.g.Scherler et al., 2011b).Slopes of lateral or terminal moraines were also included if they inclined toward the glacier's surface.For further quality check, we imported delineated glacier polygons onto Google Earth ™ with the 0.8 m Quickbird-derived images taken in different dates.

Results
The study area contains 213 debris-covered glaciers with a total surface area of 950.2 km 2 (4.46 km 2 per glacier) on a horizontal projection (Fig. 1).The total debris-covered area is 208.9 km 2 (0.98 km 2 per glacier) and PDS slopes are present in 676.0 km 2 (3.17 km 2 per glacier).

Surface temperature and aspect of PDS slopes
Figure 4a shows the daytime surface temperature distribution on PDS slopes in the Kulha Kangri Massif region, recorded by ASTER TIR sensors on 1 January 2008.Surface temperatures clearly depend on the aspect of PDS slopes (Fig. 4b).Those exceeding 0 • C account for more than 60 % of eastfacing to southwest-facing PDS slopes and less than 20 % of north-facing slopes.The concentration of warmer surfaces on southeast-facing slopes may be due to the solar direction at the time of ASTER acquisition (10:40 and 10:41 BST).There is a clear contrast between nighttime and daytime surface temperatures in relation to altitude (Fig. 4c).There is a significant negative correlation between nighttime surface temperature and altitude (r = −0.72,p < 0.001), with a lapse rate of −6.0 • C km −1 and temperature range of −30.2 • C to −6.7 • C. In contrast, there is a weaker correlation between daytime surface temperature and altitude (r = −0.53,p < 0.001), with a wide temperature range of −23.9 • C to 27.8 • C (10.0 • C standard deviation).This temperature range is much larger than the nighttime range (4.1 • C).Moreover, the melting point is exceeded at numerous localities, even at high altitudes (4500-7500 m a.s.l.) in winter.

Debris-covered area and potential debris-supply slopes
To estimate the absolute area of PDS slopes, each one-pixel DEM area (900 m 2 ) was divided by the cosine of its slope gradient, yielding a significant positive correlation between the absolute area of PDS slopes and the debris-covered area (r = 0.83, p < 0.001; Fig. 5a).This implies that glaciers with larger PDS slope areas tend to have larger debris-covered areas.Two distinct groups of glaciers were identified based on the direction of glacier flow (Figure 5b): south-flowing (90-270 • clockwise from north; open circles in Fig. 5b) and northflowing (0-90 • , 270-360 • ; solid circles).The direction of glacier flow is defined as the orientation of a line connecting the centroids of the upper and lower halves of the glacier.Both groups show a linear relation between the absolute area of PDS slopes and the debris-covered area (Fig. 5a), indicating that north-flowing glaciers tend to have smaller debriscovered areas than south-flowing glaciers, for a given PDS slope area.This finding suggests that south-facing slopes supply more debris than north-facing slopes.We therefore focused on south-facing PDS slopes (90-270 • ) and compared these with debris-covered areas in the same manner (Fig. 5b).The correlation coefficient between south-facing PDS slope area and debris-covered area is slightly improved (r = 0.88, p < 0.001) compared with the case for the total PDS slope area (Fig. 5a), and there is little difference between the south-flowing and north-flowing groups.These results suggest a significant relationship between south-facing slopes and debris-covered areas.

Sensitivity of debris supply to slope aspect and gradient
To identify the PDS slope aspect and gradient that yield the greatest debris supply, we further analysed the sensitivity of correlation coefficients between debris-covered area and PDS slope area by adjusting the sampling aspect range (Fig. 6).In terms of slope aspect, high correlation coefficients between debris-covered area and PDS slope area (r > 0.8, p < 0.001) are only obtained for southwest aspects (210 • to 240 • ), even when the aspect range is narrowed; the coefficient becomes much lower for other slope aspects (Fig. 6a).The total area of PDS slopes is not dominated by southwestfacing slopes (dashed line in Fig. 6a); consequently, the highest correlation coefficient is not caused by a dominance of surface area on southwest-facing slopes.In addition to slope aspect, there is a strong positive correlation between debriscovered area and surface gradients of 60 • , for which the aspect is mainly to the southwest (225 • ) (Fig. 6b); this tendency is found for all ranges of aspect values considered here.For slopes steeper than 70 • , the correlation coefficient shows a rapid decrease with steepness, which may reflect a reduction in the slope area itself.PDS slopes consist of both ice-free and ice-covered surfaces.A relatively gentle slope is expected to be stable and covered by snow, while steeper slopes may supply more debris due to a higher frequency of rockfalls and/or avalanches.Although 26 % of PDS slopes have a surface gradient of around 40 • as the most frequent gradient (dashed line in Fig. 6b), no influence on the correlation coefficient is evident.This suggests that PDS slopes with gradients of around 60 • have the highest potential for debris supply, despite their relatively small total area compared with slopes with gradients of around 40 • .

Shapes of debris-covered areas
Figure 7 shows two debris-covered glaciers in Mount Kulha Kangri, flowing to the north and south, respectively.The north-flowing glacier has elongate or stripe-like debriscovered areas aligned in the direction of glacier flow, whereas the south-flowing glacier has a large ablation zone covered entirely by debris.To assess whether this feature (Fig. 7) is commonly observed in the Bhutan Himalaya, a geometric shape index of the debris-covered area (S d ) is proposed as follows:  where P d is the perimeter and A d is the surface area of a debris-covered area.This is a similar concept to an index proposed by Shugar and Clague (2011).This index increases for elongate or stripe-like debris-covered areas, and decreases for round areas.Because the shape of a debris-covered area may be influenced by the shape of the host glacier, the index is normalized by the shape of the glacier, as follows: where S n is the normalized shape index and S g is the shape index for the lower part of the glacier (i.e.below the maximum elevation of the debris-covered area).To examine the relationship between the normalized shape index and the glacial flow direction, the average of the normalized shape index is calculated for glacier groups that have been categorized according to their flow direction.
Figure 8 shows the normalized shape index for an aspect bin size of 40 • over the Bhutan Himalaya.Southflowing glaciers yield the lowest values (close to 1), suggesting that debris mantle generally covers the entire ablation zone.North-flowing glaciers yield higher values, suggesting debris-covered areas with a relatively elongated or stripe-like shape.The difference between north-flowing (0-40 • /320-360 • ) and south-flowing (180-220 • ) glaciers is statistically significant (Student's t test).Differences in the shape of debris-covered areas on the two glaciers shown in Fig. 7 are therefore representative of the Bhutan Himalaya overall.
Comparing the shapes of debris-covered area and the spatial distributions of south-and north-facing PDS slopes for these two glaciers (Fig. 7), dominant debris supply is expected from the south-facing PDS slopes rather than the north-facing ones.The south-flowing glacier has a southfacing PDS slope with a wide headwall immediately above the accumulation area (Fig. 7).Debris-covered areas are separated by debris-free surfaces on this slope, and the mantle in the accumulation zone is likely to be entrained and transported down the glacier, re-appearing in the ablation zone (e.g.Hambrey et al., 1999).On the other hand, the northflowing glacier has south-facing PDS slopes on either side of the glacier tributaries, with debris-covered areas starting at the foot of these slopes.The debris mantle appears to be primarily transported on the surface of the glacier, with less entrainment than that of the south-flowing glacier.The contribution of the south-facing PDS slope to debris supply suggested here is consistent with the result in Fig. 6a even though the approach is different.These differences suggest that the shape of debris-covered areas is also controlled by the spatial distribution of south-facing PDS slopes around a glacier, which is in turn coupled with the ice flow direction.

Discussion
Our analysis shows a significant positive correlation between the debris-covered area and area of southwest-facing PDS slopes.This suggests that the most active debris production occurs on southwest-facing slopes, and that there is limited debris production on north-facing slopes.North-south contrasts in surface temperature on PDS slopes (Fig. 4b) and in the shape of debris-covered areas (Figs. 7 and 8) are also recognized.Debris-production mechanisms on PDS slopes are thought to be diurnal/seasonal freeze-thaw cycles, permafrost degradation, and snow-and-ice avalanche (e.g.Rapp, 1960;Matsuoka et al., 1998;Noetzli et al., 2006;Scherler et al., 2011b).Here we discuss the contribution of these variables to the formation of debris-covered surfaces and to the geographic characteristics of the Bhutan Himalaya.

Diurnal freeze-thaw cycle
Frost shattering is caused by repeated temperature fluctuations across the melting point coupled with high moisture content in rock materials (e.g.Matsuoka, 1994;Sass, 2005;Andren, 2006).Fracturing begins immediately below 0 • C, with volumetric expansion and water migration occurring between 0 • C and −5 • C (Matsuoka, 1990a).An effective freeze-thaw cycle is defined as temperature fluctuations occurring in the range of −2 • C to +2 • C, where −2 • C is the critical value for fracture generation in porous or jointed rocks and +2 • C is a sufficient temperature for complete melting of frozen ice (Matsuoka, 1990b).In high mountains, freeze-thaw depths caused by diurnal freeze-thaw cycles have been recorded within a few centimetres to decimetres, not exceeding 50 cm depth (e.g.Fahey, 1973;Coutard and Francou, 1989;Matsuoka, 1994).In the French Alps and the Nepal Himalaya, the diurnal freeze-thaw cycle occurs up to 100 times per annum on southwest-facing slopes and less than 50 times per annum on north-facing slopes (Coutard and Francou, 1989;Regmi and Watanabe, 2009).Diurnal variations in surface temperature and rockfall frequency on slopes that effectively receive the solar radiation are also much greater than on north-facing slopes.These in situ observations indicate that debris supply is related to slope aspect and the intensity of the diurnal freeze-thaw cycle importantly, the latter occurs most frequently on south-facing slopes.Our results indicate the same north-south contrast in surface temperature in the Bhutan Himalaya (Fig. 4).Diurnal temperatures on south-facing slopes potentially span the effective freeze-thaw range (−2 • C and +2 • C).Although southeast-facing PDS slopes are warmest (Fig. 4b), this result reflects the angle of the sun at the time of ASTER data acquisition (10:40 and 10:41 BST); satellite images record discrete moments in time and their temporal representativeness is therefore uncertain.Nevertheless, we have identified a strong correlation between debris-covered areas and PDS on southwest-facing slopes (Fig. 6a).This relationship is probably attributable to solar radiation activating the freezethaw cycle during the afternoon, when air temperature peaks.These findings, together with those discussed above from other mountain regions, suggest that debris mantles are primarily produced by diurnal freeze-thaw cycles, which are most active on southwest-facing PDS slopes.

Seasonal freeze-thaw cycles
Seasonal freeze-thaw cycles also contribute to debris production.In the Nepal Himalaya, seasonal rockfall occurs mainly in the melting season, corresponding to the seasonal temperature change where seasonal freeze-thaw cycle requires sufficient and continuous freezing temperature in winter (Regmi and Watanabe, 2009).On north-facing slopes, daily maximum and minimum temperatures are lower than the melting point in winter, and exceed the melting point in summer (Fig. 5 in Regmi and Watanabe, 2009).On southfacing slopes, however, daily maximum temperatures exceed the melting point during most of the year.This suggests that www.the-cryosphere.net/7/1303/2013/The Cryosphere, 7, 1303-1314, 2013 north-facing PDS slopes are more susceptible to the seasonal freeze-thaw cycle than are south-facing PDS slopes.If this cycle is the main driver of debris production, northfacing PDS slopes should have higher correlation coefficients with debris-covered areas than south-facing PDS slopes, yet our results do not support this assumption.Instead, we have found that the size and shape of debris-covered areas on north-flowing glaciers are linked to the spatial distribution of south-facing PDS slopes, and that the diurnal freeze-thaw cycle is the main source of debris in the Bhutan Himalaya.

Permafrost degradation and slope failure
Permafrost in steep mountain bedrock has been investigated with respect to natural hazards and the recent period of climatic warming (e.g.Kääb, 2008;Harris, 2009;Haeberli et al., 2010).In high and steep mountain regions, permafrost degradation can lead to rockfalls (Noetzli et al., 2003;Haeberli et al., 2004;Fischer et al., 2006;Rabatel et al., 2008;Allen et al., 2009;Huggel, 2009).Furthermore, extreme warming is considered to be a trigger of large landslides in high mountains due to enhanced water production from melting snow and ice, and rapid thawing (Gruber et al., 2004;Huggel et al., 2010).Current climatic change may increase the frequency of large slope failures in steep glaciated and permafrost terrain; however, the relationship between climate change and rock instability through permafrost processes is poorly understood due to a lack of long-term investigations (Krautblatter et al., 2012).The permafrost base responds slowly to glacier retreat, at a rate of millennia (Wegmann et al., 1998).Discontinuous permafrost is thought to be distributed in the Bhutan Himalaya (Brown et al., 2001;Iwata et al., 2003).
In the course of our research, we recognized that several debris-covered areas had not reached the glacier termini (Fig. S1).Similar slope failures onto glacier surfaces have been reported in other regions (e.g.Blair, 1994;Reznichenko et al., 2011;Shugar et al., 2011Shugar et al., , 2012)).Although it is difficult to determine whether such failures are related to permafrost degradation, it is considered unlikely because this type of debris cover is rare in the Bhutan Himalaya.Most debris-covered areas instead appear to receive both regular and gradual supplies of debris mantle.Fischer et al. (2012) suggests that the location of active slope failures is consistent with the estimated permafrost boundary.However, slope failure events can occur without relation to slope aspect because the actual permafrost distribution in areas of steep topography is inhomogeneous and is strongly affected by local factors, such as topography and geology (Fischer et al., 2012).These studies suggest that the remarkable aspect-dependency of debris cover on southwest-facing PDS slopes identified here (Fig. 6a) is not caused by permafrost degradation.Thus, although permafrost degradation is expected to occur in the Bhutan Himalaya, its contribution to the formation of debris-covered areas is limited compared with diurnal freeze-thaw cycle.

Snow avalanches
Snow avalanches are generally expected to occur on steep mountain slopes with a gradient of 35-40 • (McClung and  Schaerer, 1993), contributing to the supply and transport of debris mantle onto glacier surfaces (Scherler et al., 2011b).Our findings demonstrate that although the maximum area of slope distribution in the Bhutan Himalaya is 35-40 • , the strongest correlation between debris-covered area and slope gradient occurs at 60 • (Fig. 6b).This suggests that snow avalanches are a limited source of debris in this region.In the Karakoram, on the other hand, avalanche-fed glaciers are well developed hosting the large area of debris cover, and their nourishment occurs in winter and summer simultaneously (Hewitt, 2011).Less possibility of freeze-thaw activity on the surrounding slopes is thus expected due to constant snow cover.In this region, therefore, snow avalanches are expected to have larger portions of debris supply in comparison to those in the Himalayas, where arid winter climate could cause the exposure of headwall.Comparison of glaciers in the different climates with the PDS-slope analysis will reveal regional differences of the debris-supplying mechanisms.

Influence of geographical setting on debris supply
It is also important to discuss the role of geographic environments in enhancing debris supply.Scherler et al. (2011b) found a positive relationship between steepness of accumulation zone and the ratio of debris-covered area to glacier area in the Himalayas.We compare the debris-covered ratio (i.e.debris-covered area divided by glacier area) and mean gradient of the PDS slope, which represent the local geomorphic features (Fig. 9a). Figure 9a shows that the relationship reported by Scherler et al. (2011b) is evident even in a relatively small region of the Bhutan Himalaya.A positive relationship (r = 0.40, p < 0.001) implies that the steeper the local relief around the glacier, the more developed the debriscovered area will be.No relationship has been found between the absolute area of debris cover and steepness of the PDS slope (Fig. S2a).
Topographic features of glaciated landscapes have been investigated by Kääb (2005) in the Lunana region in northern Bhutan.Steeper headwalls have developed on southern slopes in this region and these provide a more sustainable debris supply to glaciers than those on the gentler northern slopes.Elevation profiles of several glaciers clearly exhibit topographic differences between those on northern and southern slopes of the Himalayan main ridge (Fig. 9 in Kääb, 2005).In addition, lithological settings in the study region are categorized into metamorphic rocks in the southern area and sedimentary rocks in the northern area divided by the main Himalayan ridge (Long et al., 2011).Although their hardness is difficult to be evaluated quantitatively, the Himalayan uplift is likely to generate large areas of shattered bedrock in southern areas supported by the complex geological features shown in Long et al. (2011).To examine the influence of north-south asymmetric topography on debris supply, we plotted the mean gradient of the PDS slope against latitude (Fig. 9b) together with maximum and minimum elevations derived from ASTER GDEM-2 data (averaged for longitudes from 88.5 • to 92.5 • ).The maximum elevation profile peaks at 28.25 • E (i.e., the Himalayan main ridge), where the southern and northern slopes are steeper and gentler, respectively.Although this analysis produced the same asymmetric profile as Kääb (2005), there is no strong correlation between the mean gradient of PDS slope and latitude in the present study area.This means that the local relief of individual glaciers is not related to the asymmetric northsouth topography of the Bhutan Himalaya.
There are also contrasting climatic conditions along the Himalayan main ridge.Southern glaciers have rapidly retreated in recent decades (Karma et al., 2003).Glaciers affected by summer monsoon precipitation are sensitive to rising temperatures, which alter precipitation from snowfall to rainfall, consequently decreasing surface albedo and enhancing ice melting (Fujita, 2008).In the longer term, degrada-tion of glaciers since the Little Ice Age (LIA) should have released PDS slopes, and thus caused more preferable condition for debris supply in this region.To examine the influence of precipitation regime on debris production, we have compared their latitudinal relationship.Tropical Rainfall Measuring Mission (TRMM) data were used to determine the spatial distribution of monthly precipitation on a grid size of 0.25 • (http://mirador.gsfc.nasa.gov)(TRMM 3B43; Huffman et al., 2007).Figure 9c shows the debris-covered ratio and the highest and lowest annual mean precipitation (averaged for longitudes from 88.5 • to 92.5 • ) in relation to latitude.Annual mean precipitation decreases from south to north, as reported by Karma et al. (2003), whereas the debris-covered ratio shows a weak correlation with latitude (r = −0.16,p < 0.001).In terms of absolute area of debris cover, however, no relation with latitude is found (Fig. S2b).Frequent rainfall has caused slope instability in other mountain regions (e.g.Sass, 2005;Krautblatter and Moser, 2009), and permafrost degradation in steep bedrock could be strongly affected by water percolating through fractures (Gruber and Haeberli, 2007).On the Tibetan side of the present study area, where latitudes exceed 28.25 • , most glaciers have a debris-covered ratio of < 20 %, whereas many glaciers on the Bhutanese side have ratios of 20-90 %.The arid climate of Tibetan glaciers (latitudes of > 28.25 • ), as evident in Fig. 9b  and c, may suppress debris production and supply even with relatively steep PDS slopes around 40 • .On the Bhutanese side, however, sufficient precipitation might explain the lack of a latitudinal gradient in debris production.In this region, local relief of the PDS slopes, which has been exposed since the LIA would be the primary control of debris supply to the glaciers.

Conclusions
We performed a spatial analysis of potential debris-supply (PDS) slopes, assuming that these slopes are the main source of debris mantle on glaciers in the Himalayas.Previous studies on debris-covered glacier focused on debris supply via avalanches and did not discuss the contribution of rockfall activity which was affected by the slope aspect (e.g.Scherler et al., 2011b).Our research demonstrates that the dimension and shape of debris-covered areas are strongly controlled by the spatial distribution of southwest-facing PDS slopes.Surface temperatures on these slopes frequently fluctuate around the melting point, providing large amounts of debris under the influence of the diurnal freeze-thaw cycle.The relationship between PDS slope area and debriscovered area indicates that PDS slope aspect is an important variable controlling the rate of debris supply.The seasonal freeze-thaw cycle, permafrost degradation, and snow avalanches also occur in the study region.However, the wellcorrelated aspect and gradient tendencies of debris-covered areas do not support their remarkable contributions.In terms www.the-cryosphere.net/7/1303/2013/The Cryosphere, 7, 1303-1314, 2013 of geographic setting, local topography contributes to debris supply on glaciers regardless of latitude.The climate regime of the Bhutan Himalaya may influence on the smaller debris-covered ratio in the Tibetan side.The methodology adopted here is potentially transferrable to other regions such as Karakoram, the European Alps, and the Andes, and could therefore provide a global perspective on the formation conditions of debris-covered glaciers.The spatial distribution of PDS slopes may also be used to assess the formation conditions of growing glacial lakes, by providing a means of estimating the critical quantity of debris mantle required to trigger hazardous lake expansion.

Fig. 1 .
Fig. 1.Spatial distribution of 213 debris-covered glaciers in the Bhutan Himalaya.The background image is a mosaic of AVNIR-2 images.

Fig. 2 .
Fig. 2. Coverage areas of satellite data used in this study.

Fig. 3 .
Fig. 3. Visual interpretation of potential debris-supply slopes on a PRISM image (location indicated in Fig. 1).(1) Traces or scars show the possibility of debris supply.(2) Small hills identified by shadings and (3) depressions outside lateral moraine may prevent debris supply.

Fig. 4 .
Fig. 4. (a) Spatial, (b) aspect, and (c) altitudinal distributions of surface temperature on potential debris-supply slopes in the Kulha Kangri Massif region on 1 January 2008 (location indicated in Fig. 1).Background image is AVNIR-2.Aspect bin size in (b) is 10 • .Nighttime temperature distribution on 27 January 2005 is also shown and conpared in (c).

Fig. 6 .
Fig. 6.Correlation coefficients between the debris-covered area and potential debris-supply (PDS) slope area, plotted against (a) aspect (clockwise from north), and (b) gradient of the PDS slope.Correlation coefficients were calculated over four sampling aspect ranges (180 • , 90 • , 45 • , and 15 • ).Surface area distributions are calculated in (a) aspect bin of 30 • , and (b) gradient bin of 10 • .The central aspect is fixed to the southwest (225 • ) in (b), at which the highest correlation coefficient was found in (a).

Fig. 7 .
Fig. 7. Aspect distribution of potential debris-supply (PDS) slopes of north-flowing and south-flowing debris-covered glaciers in the Kulha Kangri Massif region.Background image is AVNIR-2.The location of this area is indicated in Fig. 4a.

Fig. 8 .
Fig. 8. Normalized shape index plotted against glacier flow direction.Averages and standard deviations for an aspect bin of 40 • are presented with the number of samples in each bin.

Fig. 9 .
Fig. 9. Relationships between (a) debris-covered ratio and mean gradient of potential debris-supply (PDS) slope, (b) mean gradient of PDS slope and latitude, and (c) debris-covered ratio and latitude in the Bhutan Himalaya.Averages with standard deviation of the debris-covered ratio are calculated for 5 • bins (red lines in (a)).In (b), the elevation range derived from ASTER GDEM-2 data shown by grey shading covers the longitudinal area of 88.5-92.5 • E. In (c), the mean annual precipitation is estimated from the TRMM 3B43 product (1998-2010; 88.5-92.5 • E), shown by grey shading.
Figure S1.Example of irregularly supplied debris cover (shading), probably supplied from the neighbouring PDS-slope (allow).Inset figure shows the location.

Figure S2 .
Figure S2.Relationships between a) debris-covered area and mean gradient of potential debris-supply slope, and b) debris-covered ratio and latitude in the Bhutan Himalaya.In (b), the mean annual precipitation estimated from the TRMM 3B43 product (1998-2010; 88.5º-92.5ºE) is shown by grey shading.