Changes in glacier surface cover on Baltoro glacier, Karakoram, north Pakistan, 2001–2012

ABSTRACT

The presence of supraglacial debris on glaciers in the Himalaya-Karakoram affects the ablation rate of these glaciers and their response to climatic change. To understand how supraglacial debris distribution and associated surface features vary spatially and temporally, geomorphological mapping was undertaken on Baltoro Glacier, Karakoram, for three time-separated images between 2001–2012. Debris is supplied to the glacier system through frequent but small landslides at the glacier margin that form lateral and medial moraines and less frequent but higher volume rockfall events which are more lobate and often discontinuous in form. Debris on the glacier surface is identified as a series of distinct lithological units which merge downglacier of the convergence area between the Godwin-Austen and Baltoro South tributary glaciers. Debris distribution varies as a result of complex interaction between tributary glaciers and the main glacier tongue, complicated further by surge events on some tributary glaciers. Glacier flow dynamics mainly controls the evolution of a supraglacial debris layer. Identifying such spatial variability in debris rock type and temporal variability in debris distribution has implications for glacier ablation rate, affecting glacier surface energy balance. Accordingly, spatial and temporal variation in supraglacial debris should be considered when determining mass balance for these glaciers through time.

KEYWORDS:

• Supraglacial debris
• Baltoro glacier
• Karakoram
• debris-covered glacier
• surface energy balance
• debris distribution
• glacial geomorphology

1. Introduction

Debris-covered glaciers, where the majority of the ablation area is covered in rock debris, are ubiquitous in mountainous regions globally (Benn & Lehmkuhl, 2000), including the Southern Alps of New Zealand, European Alps, the Andes and in High Asia. The Himalaya and Karakoram ranges contain the largest area of glacial ice outside of the Polar Regions, with one-fifth of the world’s population depending on the run-off from these glaciers as a water resource (Bolch et al., 2012). Approximately 23% of glacier area in the Himalaya is covered by debris, increasing to 33% in the Everest region (Scherler, Bookhagen, & Strecker, 2011; Thakuri et al., 2014). The input of debris to glaciers in the Karakoram is driven by high rates of rock uplift in the mountain range and correspondingly high rates of erosion, ranging between 0.06 and 2.5 mm a⁻¹ (Seong et al., 2007). Debris distribution is controlled by the delivery of rock debris and debris-transport by ice flow, which slows as the glacier loses mass, resulting in a constantly changing distribution of supraglacial debris (Bolch, Buchroithner, Pieczonka, & Kunert, 2008; Rowan, Egholm, Quincey, & Glasser, 2015; Thakuri et al., 2014). Understanding the nature of these changes in debris distribution is currently limited, due to the lack of data describing how supraglacial debris layers evolve through time. However, our knowledge of these glaciers needs to progress to determine their response to current climatic change, and the impact these changes will have on the people who rely on them.

Supraglacial debris distribution is important when calculating glacier mass balance, as such debris layers attenuate the ablation of the underlying ice depending on the debris thickness (Evatt et al., 2015; Østrem, 1959). A thin debris layer (up to around 0.02 m thick), when compared to clean ice, enhances melt by increasing the albedo of the glacier surface, causing it to absorb more solar radiation and melt more quickly than debris-free ice. At a critical debris thickness, the melt rates of the debris-covered ice and the clean ice are equal. Above the critical thickness, which differs between glaciers, the ablation rate of underlying ice decreases exponentially with increasing debris thickness (Kayastha, Takeuchi, Nakawo, & Ageta, 2000; Nicholson & Benn, 2006; Østrem, 1959). In addition to meteorological parameters, the resultant ablation rate of ice under a debris layer is affected by debris characteristics, including lithology and porosity of the debris layer. The thermal capacity and conductivity of the debris clasts is partially dependent on the lithology of the debris and so will affect how much thermal energy is held in and transferred through the debris layer (Conway & Rasmussen, 2000; Nicholson & Benn, 2013). Incorporating the evolution of debris distribution, thickness and lithology into calculations of the mass balance of debris-covered glaciers will result in better constraint of the role of debris in modifying the response of these glaciers to climatic change.

Geomorphological mapping of debris distribution using satellite imagery is a useful tool to identify temporal changes in debris distribution across an entire glacier surface, reducing the need for expensive and often challenging fieldwork. By identifying the sources of debris, and the interaction of different debris units and glacier flow, controls on spatio-temporal variations in debris distribution can be determined. A greater understanding of the interaction between hillslope and glacial processes can also be developed with regard to the input of supraglacial debris into the glacial system. Baltoro Glacier is highly appropriate for supraglacial debris mapping due to its location in a region of high erosion rates and concurrently extensive debris cover, and the large number of tributary glaciers which form the glacier system and lead to a complex flow regime and spatially varied debris distribution.

2. Study site

Baltoro Glacier is located in the Karakoram mountain range in northern Pakistan. The glacier and its tributaries cover an area of ∼660 km², with a centre-line length of ∼62 km (Figure 1). The second highest mountain in the world, K2, is located ∼12 km north of the main glacier tongue. The main tongue of Baltoro Glacier flows from Concordia (4600 m) in an east to west direction, to its terminus at ∼3410 m, and has 12 tributary glaciers along its lengths (Tributary Glaciers (TGs) 3–14, Main Map). Above Concordia Godwin-Austen Glacier (TG 1) and Baltoro South Glacier (TG 2) converge to form the main tongue of Baltoro Glacier. Baltoro Glacier has an extensive debris cover, with ∼38% of its area covered by supraglacial debris (Mayer, Lambrecht, Belò, Smiraglia, & Diolaiuti, 2006). Debris thickness has a general pattern of increasing thickness downglacier, reaching around 3 m in thickness at the glacier terminus (Mihalcea et al., 2008).

Figure 1. The regional extent of the Baltoro Glacier system and relevant locations, including Concordia and named tributary glaciers.

3. Previous work

Baltoro Glacier is relatively well studied despite its location in a remote and politically unstable region. Much of the existing research literature is almost exclusively based on remote sensing datasets (e.g. Collier et al., 2014; Quincey et al., 2009a), although some studies have been complemented by fieldwork (e.g. Mayer et al., 2006; Mihalcea et al., 2008; Searle et al., 2010). Supraglacial debris thickness for the entire glacier was calculated by Mihalcea et al. (2008) and Minora et al. (2015) for 2004 and 2011, respectively, but these maps only consider debris thickness. These maps therefore disregard lithology and the potential impact of spatial differences in thermal energy stored in the debris due to different lithologies across the glacier surface. Lack of segregation between debris lithology may cause errors in these debris thickness output maps, as thermal capacity of rock debris influences debris surface temperature, the parameter used to calculate debris thickness in these two studies (Steiner & Pellicciotti, 2015). Collier et al. (2013) incorporated a debris layer with a simulated thickness into a distributed energy balance model for the glacier, but did not consider spatial variations in debris lithology either. Omitting lithological variation is a potential error in the resultant ablation rates calculated by the study, as differing thermal capacity and rates of heat transfer for different debris lithologies were not incorporated into calculations. Other studies focused on different aspects of the Baltoro Glacier system; Quincey et al. (2009a) investigated the changing velocity of the glacier between 1993 and 2008 and linked their observations to climatic variability, while both Diolaiuti, Pecci, and Smiraglia (2003) and Belò, Mayer, Smiraglia, and Tamburini (2008) identified the surging nature of the confluent Liligo Glacier (TG 7) (Figure 1) and reconstructed its recent dynamic history.

4. Methods

Geomorphological features of Baltoro Glacier were mapped using three time-separated and orthorectified Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) images. The three images were acquired in August 2001 (29th), 2004 (14th) and 2012 (20th), at a resolution of 15 m, and were selected because of low levels of cloud cover (46%, 1% and 13%, respectively). Images from the same month in each year also allowed direct comparison of seasonally controlled aspects of the glacier geomorphology (e.g. supraglacial water bodies and mass movement scars). ASTER images were georeferenced and orthorectified in ENVI (v.5.0) prior to using them for mapping surface cover on Baltoro Glacier. Georeferencing was undertaken using the ASTER georefering tool in ENVI (v.5.0) and a Landsat 7 ETM+ image (acquired on 9 August 2002), which is georeferenced and orthorectified prior to downloading. Tie points between each ASTER and the Landsat 7 ETM+ images were created as part of the georeferencing process. Georeferenced ASTER images were then compared to the Landsat 7 ETM+ image and to each other to confirm that there was no pixel shift present between images. Orthorectification was then undertaken on each ASTER image using the ENVI (v.50) ASTER orthorectification tool, which corrects the effects of sensor tilt and terrain, and creates a planmetric image.

An additional Landsat 7 Enhanced Thematic Mapper (ETM+) image, acquired in August 2001 at 30 m resolution, was used to map upper regions of the glacier that were covered in cloud in the ASTER 2001 imagery. Additionally, Quickbird data, accessed through Google Earth (2016) were used to investigate specific features at a horizontal resolution of 2.4 m. Features were mapped from ASTER imagery using a false-colour composite (Bands 3N, 2, 1) at 15 m resolution and were manually digitised in ESRI ArcGIS 10.1. Geomorphological features were initially mapped at a glacier-wide scale, then the boundaries of features identified at a finer scale to increase mapping accuracy. Debris units were identified manually using colour differences between pixels on the glacier surface, aided by comparison of the red, green and blue spectral signal of these pixels with different colours (Lillesand, Keifer, & Chipman, 2014). The units were then traced up-glacier to their source area and lithology was identified by matching the source area with geological units identified by Searle et al. (2010) in the geology map of the Karakoram. Mass movement deposits were identified by the presence of two features: a scar, identified as a lighter, elongate feature on the valley side, and an associated debris fan deposit on or near the glacier surface (Figure 2(a,b)). Pixels containing supraglacial water were identified using a normal difference water index (NDWI) (Gao, 1996) in ENVI, (v. 5.0) calculated from ASTER Bands 3 and 4 using Equation (1).(1) $\mathrm{NDWI}=\left(\mathrm{NIR}\right(\mathrm{Band}3)-\mathrm{SWIR}\left(\mathrm{Band}4\right)\left)/\right(\mathrm{NIR}\left(\mathrm{Band}3\right)+\mathrm{SWIR}\left(\mathrm{Band}4\right)).$(1) Manual identification and digitisation of supraglacial features resulted in some misclassification of pixels, particularly where changes in debris distribution or supraglacial pond size were on the order of a few pixels. The extent of this error was reduced by repeated digitisation of each feature three times and the three outputs were then compared against the imagery to determine the most appropriate outline. Digitised margins of debris-covered glaciers in previous studies have been calculated to have errors between 2.5% and 4.5% (e.g. Bolch, Menounos, & Wheate, 2010; Paul, Frey, & Le Bris, 2011). Consequently, changes on the order of less than 4.5% of the feature area were disregarded in this study.

Figure 2. Examples of mass movement deposits (a,b) and supraglacial water bodies (c,d) in original ASTER 2012 images and mapped examples, respectively.

The map associated with this paper is the 2012 geomorphological map of Baltoro Glacier, with insets from 2001, 2004 and 2012 showing changes in surface cover through time. The map is in a Universal Transverse Mercator (UTM) WGS 1984 projected coordinate system with off-ice 250 m contours derived from the ASTER Global digital Elevation Model (GDEM) (2000, 30 m resolution) using the contour tool in ArcGIS. The glacier outline used is adapted from the Randolph Glacier Inventory (Arendt, Bolch, Cogley, Gardner, & Hagen, 2012; Bajracharya, Shrestha, Bajracharya, Maharjan, & Guo, 2014).

5. Features related to debris distribution

5.1. Debris units

The schist, gneiss and granitic debris units are part of the Baltoro Batholith and associated metamorphosed rock units, while the metasediment debris unit originates from a small exposure of metamorphosed limestone near the headwall of TG2. The debris-covered glacier area has varied between 2001 and 2012 by ∼5%, but gneiss debris dominates the debris distribution on Baltoro Glacier throughout this period, covering ∼53% of the area (Table 1). The metasedimentary debris has the most variable area between 2001 and 2012 and does not reach the glacier terminus. It is likely that the increase in the area of the metasediment debris unit is due to it being a younger deposit than other debris units, and so it has not reached the terminus yet, or that it has reached the terminus but has been mixed or covered by other debris units before the terminus. All debris units form moraine structures, both lateral and medial, which are fed by debris input from the glacially incised valley walls and associated valley spurs at the confluences of glacier flow units through small but frequent landslide events, as seen in mountain valley environments (Rowan et al., 2015). These units are initially separated by debris-free ice, but converge below Concordia. The debris units have distinct margins in the upper and mid-sections of the main glacier tongue, converging to a massive debris unit in the lower 10 km of the glacier. The moraines increase in width with increasing distance from their source until they converge, due to increased debris input from englacial meltout and slumping of the moraines over time (Hambrey & Glasser, 2003). Little change in the main debris unit boundaries occurs in the mid- and lower sections of the main glacier tongue between 2001 and 2012, suggesting that once units converge their width is constrained by the other debris units.

Table 1. The total area of each debris unit type, based on lithology, for 2001, 2004 and 2012, and the percentage of each debris type as a proportion of the total debris cover for Baltoro Glacier and its tributary glaciers.

Year	2001	2001	2004	2004	2012	2012
Debris type	Area (km^2)	% of total debris	Area (km^2)	% of total debris	Area (km^2)	% of total debris
Gneiss	81.48	52.9	79.83	51.2	79.48	52.9
Metasediment	8.60	5.6	11.41	7.3	8.28	5.5
Schist	17.76	11.5	18.54	11.9	17.71	11.8
Granite	46.29	30.0	46.20	29.6	44.71	29.8
Total debris (m^2)	154.13		155.97		150.17

5.2. Supraglacial water bodies

The number of supraglacial water bodies, also referred to as ‘ponds’, increases from 234 to 570 over the study period, occupying 0.66 km² of the debris surface in 2001, and 2.04 km² in 2012, an increase in pond area of 209% (Table 2). Pond formation is associated with spatially varying ablation which results in the formation of hummocky topography where water accumulates in the depressions on the debris surface (Reynolds, 2000). Identifying areas of supraglacial water bodies is important for determining the state of a debris-covered glacier. Over decadal timescales, increasing supraglacial ponds can suggest build-up of debris cover and a lowering of the glacier surface gradient, resulting in a reduction of the driving stress and hence glacier stagnation (Quincey, Luckman, & Benn, 2009b; Wessels, Kargel, & Kieffer, 2002). In the case of this study, it is likely that increased pond number and area are responses to temporal variation in climatic conditions, specifically an increase in precipitation since 2000, which was also attributed to increasing surface velocity in 2005 (Quincey et al., 2009a). Such water bodies increase the local glacier ablation rate through increased absorption of incoming shortwave radiation (Sakai, Nakawo, & Fujita, 2002). Ice cliffs are often concomitant with supraglacial ponds, and are a focus for a large proportion of a debris-covered glacier’s net ablation (Juen, Mayer, Lambrecht, Han, & Liu, 2014; Reid & Brock, 2014; Sakai, Nakawo, & Fujita, 1998). Identifying the number of supraglacial ponds is therefore integral to understanding the controls on the ablation rate of a debris-covered glacier and the relative importance of these controls over time as the debris cover evolves. Although conditions may not be entirely comparable, a similar increase in pond number and area on Khumbu Glacier was thought to be a result of a period of recession following a period of heightened activity during the middle part of last decade (Quincey et al., 2009a), and similar conditions could also be true for Baltoro Glacier, but would need further investigation of variations in past surface velocity and mass balance.

Table 2. The total area and number of supraglacial waterbodies identified using an NDVI index for each time step between 2001 and 2012.

	2001	2004	2012
Number of water bodies	234	404	570
Area (km^2)	0.66	1.79	2.04

5.3. Areas of mass movement

Mass movement scars and deposits are visible along most of the tributary and main glacier margins, with an increasing number appearing between 2001 and 2012, particularly between 2001 and 2004. An increasing number of mass movement events may be related to a higher frequency of seismic activity or increased precipitation between 2001 and 2004 (Barnard, Owen, Sharma, & Finkel, 2001; Bookhagen, Thiede, & Strecker, 2005). The most evident mass movement deposits in the Baltoro Glacier system are found on TGs 4 and 5 (Inset B). Two mass movement deposits on TG5 and a series of smaller but more frequent protrusions on TG4 are interpreted as being sourced from intermittent mass movement events. The mounds of debris have then been transported with glacier flow to the north, at a rate of ∼20 m a⁻¹ (calculated from the distance the deposits have travelled between images, divided by the number of years between images). These periodic events, occurring on sub-decadal to decadal timescales, deliver an additional source of supraglacial debris to the glacier system, which are subsequently dissipated and result in thickening of the debris layer.

5.4. Tributary glacier convergence

In addition to Godwin-Austen and Baltoro South Glaciers, 12 tributary glaciers flow into the main tongue of Baltoro Glacier (Tributary Glacier inset, Main Map). These tributary glaciers, flowing perpendicular to the main tongue, have spatially and temporally varying covers of supraglacial debris, which affects the distribution of debris onto the main glacier tongue. TG4 in Inset B is predominantly debris free, with the exception of two longitudinal schistose debris bands. A series of supraglacial ponds have developed on the main glacier tongue to the east of the two debris bands on TG4 between 2001 and 2012. Pond development also occurs on TG3 and on the glacier surface below Concordia, where TGs 1 and 2 converge. Pond development in the convergence area below Concordia is likely due to the same process as found on clean ice glaciers where water ponding occurs in surface depressions (Banwell, Arnold, Willis, Tedesco, & Ahlstrøm, 2012). For TGs 4 and 3, the formation of the water bodies is likely due to a constriction of the main glacier tongue causing an undulant surface topography. These constrictions in glacier flow also manifested in debris thickening in the same areas, observed on the debris thickness map of Baltoro Glacier produced by Minora et al. (2015). TG8, flowing southeast into the terminus of Baltoro Glacier, is unique in the Baltoro Glacier system, as no obvious extensive debris cover is identifiable on its surface in satellite imagery. The lack of debris cover is attributed to its steep surface slope and high velocities likely to be associated with this, resulting in little time for debris to accumulate on the surface, as well as its relatively short length in comparison to all other tributary glaciers in the system, reducing the valley side area from which debris can be supplied (Young & Hewitt, 1993).

5.4.1. Surging tributary glaciers

The schistose debris unit on Liligo Glacier (TG7, Inset C) shows an apparent advance and subsequent retreat of the schistose debris unit extent between 2001 and 2004 and 2004 and 2012, respectively. Retreat of the debris unit is, however, unlikely, as movement of rock debris upslope would have had to occur. Therefore, the changes in debris unit boundary between 2001 and 2012 are attributed to dissipation and reorganisation of the debris, as the surface morphology changes following the active surge phase of Liligo Glacier between 1973 and 2001 (Belò et al., 2008; Diolaiuti et al., 2003). The distribution of the debris unit succeeds the ultimate advance of the glacier, likely to result from slumping of the debris into the depression in the proglacial area as the tributary glacier commences its quiescent phase and decouples from the main Baltoro Glacier tongue. The associated supraglacial pond identified by Belò et al. (2008) in 2001 persists until 2004, but had drained by 2012, likely due to downwasting of ice and redistribution of the debris in the area between Liligo Glacier and Baltoro Glacier. Changes in surface morphology associated with pond drainage may have also contributed to the shift in debris unit boundary.

Geomorphological features present on TG12 (Inset A) also suggest that a surge event has occurred recently and that the glacier is now in a quiescent phase. This tributary glacier has not been documented to surge previously. The sinuous nature of the debris unit margins at the confluence of TGs 12 and 13, and the development of a fold-type structure on the eastern schistose debris unit ∼1 km from the convergence of these glaciers (Figure 3), are both evidence for a surge event, where TG13 would have accelerated and its terminus advanced, folding the debris units downglacier from this confluence, analogous to folding of the debris bands on Susitna Glacier, Alaska (Meier & Post, 1969). Additional evidence for the occurrence of a surge event include heavy crevassing of the TG13, upglacier from the main fold in the schistose debris unit, while the area below the folded schist units are heavily ponded with a hummocky surface, suggesting downwasting of the ice below the surge front as the glacier recedes following a surge event (Evans & Rea, 1999).

Figure 3. The geomorphological features on Baltoro Glacier’s tributary glacier showing evidence of a surge event (Google Earth, 2016).

6. Conclusions

Supraglacial debris units were segregated into four lithologies, with each debris type, particularly metasediments, being composed of different minerals that would cause them to have different thermal capacities and albedo. Supraglacial cover varies both spatially and temporally on Baltoro Glacier. Temporal variations in supraglacial pond number and area are attributed to climatic conditions in this study, but should be monitored on Baltoro Glacier in the future to determine whether this trend continues and if it is linked to other aspects of glacier dynamics. Debris distribution is controlled by the input of material from frequent rockfalls in the upper reaches of the glaciers upglacier of Concordia, subsequently entrained into medial and lateral moraines. Larger and less frequent mass movement events in the ablation area add additional debris to the system and occur on sub-decadal to decadal timescales. Once debris has entered the glacial system, variability in its distribution is controlled by the interaction of flow units between tributary glaciers and the main Baltoro Glacier, complicated by the advance and recession of tributary glacier termini due to surge events. The changes in debris unit areas between 2001 and 2012 are attributed to error in identification of debris-covered pixels in the analysis of remotely sensed data, as the difference in glacier area is near to the assumed error value and so may due to error in manual digitisation. However, it is likely that variations in debris distribution do occur through time due to variations in glacier dynamics and input of debris to the glacier system. Consequently, analyses of variations in debris distribution over longer timescales are needed as spatial variation in debris unit lithology and layer thickness will affect the energy balance of the supraglacial debris layer, which would cause spatio-temporal variation in glacier ablation. It is therefore appropriate to undertake further investigations into rates of change in supraglacial debris distribution and debris thickness through time, as well as to determine the influence of debris lithology on heat transfer and subsequent ablation of underlying ice. Such spatial and temporal variations could then be incorporated into glacier surface energy balance modelling to refine subsequent calculations of ablation rates for debris-covered glaciers, to aid in determining the role of a supraglacial debris layer in the response of these glaciers to climatic change.

Software

ENVI (v. 5.0) was used to prepare the remotely sensed data for mapping. ESRI ArcGIS (v. 10.1) was used for mapping all features and for production of the Main Map and insets, alongside Google Earth (2016) to investigate and identify features in more detail. The final map was constructed in Adobe Illustrator (v. 15.1.0).

Acknowledgements

The authors would like to thank Dr Tom Holt and Stephen Jennings for consultation on glacial structures and associated processes. The reviewers are also thanked, in particular Dr Christoph Mayer, for providing useful comments on the manuscript, as these comments greatly improved the final map and associated paper.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

M. J. Gibson http://orcid.org/0000-0003-1707-5562