1. Introduction

The largest glacierized region outside the Arctic and the Antarctic is High Mountain Asia (HMA), which covers an area of 118 200 km² (Reference GardnerGardner and others, 2013). Changes in glacier extent and volume in this region are spatially heterogeneous and poorly known (Reference BolchBolch and others, 2012). Indeed, recent studies revealed that most of the northwestern Himalaya have experienced less glacier shrinkage than the eastern parts of the same mountain range (Reference Bhambri and BolchBhambri and Bolch, 2009; Reference BolchBolch and others, 2012; Reference Kääb, Berthier, Nuth, Gardelle and ArnaudKääb and others, 2012). In the western and central Karakoram region, glaciers showed long-term irregular behavior with frequent advances, and possible slight mass gain in the last decade (Reference CoplandCopland and others, 2011; Reference HewittHewitt, 2011; Reference BolchBolch and others, 2012; Reference Gardelle, Berthier and ArnaudGardelle and others, 2012, Reference Gardelle, Berthier, Arnaud and Kääb2013; Reference Kääb, Berthier, Nuth, Gardelle and ArnaudKääb and others, 2012; Reference MinoraMinora and others, 2013; Reference SonciniSoncini and others, 2015). Reference Gardelle, Berthier and ArnaudGardelle and others’ (2012, Reference Gardelle, Berthier, Arnaud and Kääb2013) recent studies demonstrate how, in contrast to widespread global glacier retreat, glaciers in the Karakoram region as a whole have exhibited a general mass-balance stability (the so called ‘Karakoram anomaly’; Reference HewittHewitt, 2005, Reference Hewitt2011). Advances of individual glaciers have also been reported in the Shyok valley (eastern Karakoram) during the last decade (Reference Raina and SrivastavaRaina and Srivastva, 2008). These individual advances and mass gain episodes could be attributed to surging (Reference Barrand and MurrayBarrand and Murray, 2006; Reference HewittHewitt, 2007; Reference CoplandCopland and others, 2011; Reference Quincey, Braun, Bishop, Hewitt and LuckmanQuincey and others, 2011), temperature drops (Reference Shekhar, Chand, Kumar, Srinivasan and GanjuShekhar and others, 2010) and increased solid precipitation in the accumulation areas (Reference Fowler and ArcherFowler and Archer, 2006; Reference Bocchiola and DiolaiutiBocchiola and Diolaiuti, 2013). The Karakoram glaciers are a strategic resource for Pakistan, because they provide fresh water for civil use, hydropower production and farming. The glacierized Karakoram is therefore a key area for studying the effects of ongoing climate change on present and future meltwater discharge.

This study focuses on the glacier ablation areas within the Central Karakoram National Park (CKNP), with the aim of assessing the magnitude and rate of ice ablation and evaluating the derived meltwater amount. For this purpose, we applied a distributed model able to describe ablation in debris-covered and debris-free conditions (Reference Pellicciotti, Brock, Strasser, Burlando, Funk and CorripioPellicciotti and others, 2005; Reference MihalceaMihalcea and others, 2008a). Indeed, a significant portion of the glaciers in the CKNP is covered by a supraglacial debris layer, modulating the magnitude and rate of ice ablation (Reference Nakawo and YoungNakawo and Young, 1981; Reference Nakawo and TakahashiNakawo and Takahashi, 1982; Reference Nicholson and BennNicholson and Benn, 2006; Reference MihalceaMihalcea and others, 2008a,Reference Mihalceab; Reference Reid and BrockReid and Brock, 2010). This debris layer must therefore be accurately considered in distributed modeling of ice melt.

While quite a few energy- and mass-balance studies have been performed on debris-free glaciers, studies including debris-covered ice are not numerous. In the recent past, some authors have focused their attention on debris-covered ice only, and at single-point sites. For example, Reference Nicholson and BennNicholson and Benn (2006) presented a modified surface energy-balance model to calculate melt beneath a debris layer from daily mean meteorological data on two European debris-covered glaciers (Ghiacciaio del Belvedere, Italy, and Larsbreen, Norway). Reference Han, Ding and LiuHan and others (2006) proposed a simple model to estimate ice ablation under a thick debris layer by using surface temperature and debris thermal properties on Koxkar glacier, Tien Shan, China. During the last two decades, a few papers have focused on debris-covered glaciers in the Himalaya and Karakoram (e.g. Reference Hewitt, Wake, Young and DavidHewitt and others, 1989; Reference Mattson and GardnerMattson and Gardner, 1989; Reference Mattson, Gardner and YoungMattson and others, 1993; Reference Young and HewittYoung and Hewitt, 1993; Reference Nakawo and RanaNakawo and Rana, 1999; Reference Kayastha, Takeuchi, Nakawo and AgetaKayastha and others, 2000; Reference Nakawo, Raymond and FountainNakawo and others, 2000; Reference Takeuchi, Kohshima, Yoshimura, Seko and FujitaTakeuchi and others, 2000; Reference Lejeune, Bertrand, Wagnon and MorinLejeune and others, 2013). Some studies have used remote-sensing data to analyze the spatial distribution of surface temperatures and calculate the energy available for melting (Reference Nakawo, Moroboshi and UeharaNakawo and others, 1993; Reference Rana, Nakawo, Fukushima and AgetaRana and others, 1997; Reference Nakawo and RanaNakawo and Rana, 1999). Unfortunately, these studies only provided melt data over small areas and short time spans. Reference MihalceaMihalcea and others (2008a) modeled debris-covered ice ablation over the whole Baltoro glacier ablation area by applying a distributed approach, based on computation of the conductive heat flux through the debris layer and requiring information on debris thickness distribution derived from ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) thermal data. This approach has also been used by Reference Zhang, Fujita, Liu, Liu and NuimuraZhang and others (2011) who applied it on Hailuogou glacier, southeastern Tibetan Plateau, and more recently by Reference Fujita and SakaiFujita and Sakai (2014) on the Tsho Rolpa glacial lake–Trambau glacier basin in the Nepal Himalaya. Reference FyffeFyffe and others (2014) developed a melt model, which calculates sub-debris melt rates using an existing debris energy-balance model (DEB-Model introduced by Reference Reid and BrockReid and Brock, 2010) and melt rates for clean ice, snow and partially debris-covered ice using standard energy-balance equations. This approach is more exhaustive (but also more complex) than that of Reference MihalceaMihalcea and others (2008a), though its application to a whole glacierized watershed or an entire glacier region is not simple, and requires input data featuring high spatial and temporal resolution, not always available in remote high-elevation glacier zones. Therefore, the results we present in this study were obtained for the entire CKNP debris-covered ice zone by applying the model developed by Reference MihalceaMihalcea and others (2008b). Furthermore, we assessed the contribution of debris-free ice melt to fresh water from the whole CKNP ablation area. Modeling of debris-free ice melt has been extensively analyzed in the recent past, and several attempts have been made to apply the degree-day approach or simplified energy budget computations to large and representative glacierized catchments worldwide. Nevertheless, on large and remote glacierized catchments, only a few attempts to model debris-free ice melt coupled with approaches estimating debris-covered ice ablation are available (e.g. Reference SonciniSoncini and others, 2015). In our study, we evaluate ice melt in debris-free conditions by applying an enhanced T-index model (following Reference Pellicciotti, Brock, Strasser, Burlando, Funk and CorripioPellicciotti and others, 2005), which also considers solar energy inputs in driving ice melt. This melt model was coupled with the debris-covered ice melt model, thus obtaining a suitable tool to be applied on large and remote glacierized areas featuring both debris-covered and debris-free ice conditions.

2. Study Site

The CKNP is a protected area established in 2009 covering 12 162 km² in northeastern Pakistan, at the border with India and China (Fig. 1). The park protects major natural resources for the country, including >700 glaciers, with a total area of 4632 km², corresponding to ~30% of the overall glacier area of the Pakistani Karakoram (ICIMOD, 2013). Glaciers in the CKNP span a broad range of sizes, types (mountain glaciers, glacierets, hanging glaciers, compound-basin valley glaciers), and surface conditions (debris-free and debris-covered ice). In addition, the Karakoram is known to host several surge-type glaciers (Reference Diolaiuti, Pecci and SmiragliaDiolaiuti and others, 2003; Reference Barrand and MurrayBarrand and Murray, 2006; Reference HewittHewitt, 2007, Reference Hewitt2011; Reference Kotlyakov, Osipova and TsvetkovKotlyakov and others, 2008; Reference Gardelle, Berthier and ArnaudGardelle and others, 2012), displaying cyclically short active phases involving rapid mass transfer from high to low elevations, and long quiescent phases of low mass fluxes. Reference CoplandCopland and others (2011) reported 90 surge-type glaciers in the Karakoram mountains. In particular, in the 14 years after 1990, they found double the number of new surges (26) compared with those counted in the 14 years before (13 surges).

Fig. 1. Study area. The map shows the border of the Central Karakoram National Park (CKNP) in northern Pakistan (red line), the automatic weather stations at Askole, Urdukas and Concordia (yellow dots), the glacier boundaries (blue lines), and the glacier areas covered by supraglacial debris (orange). Names of the widest glaciers are shown, with the area of Baltoro glacier highlighted by a box with blue dashed line.

The ablation areas of many glaciers in the CKNP are heavily debris-covered, because of abundant rockfalls from steep walls, and intense avalanche activity (e.g. Reference BolchBolch and others, 2012). Supraglacial debris covers ~10% of the total glacier area in the Karakoram (Reference BolchBolch and others, 2012).

Baltoro glacier is one of the largest debris-covered glaciers in the Karakoram Range, with a maximum length of ~62 km and an area of ~524 km², including all connected tributaries (Reference Mayer, Lambrecht, Beló, Smiraglia and DiolaiutiMayer and others, 2006). The total drainage basin of Baltoro glacier is ~1500 km². It extends in the east–west direction on the south side of the Karakoram Range, lying in the region 35°35′–35°56′ N, 76°04′– 76°46′ E. The Baltoro glacier elevation ranges between 3370 m a.s.l. and the K2 summit (8611 m a.s.l.). Supraglacial debris of diverse lithology (Reference Desio, Marussi and CaputoDesio and others, 1961) occurs below 5000 m a.s.l., and covers ~38% of the glacier area. The grain size shows a large variability on the glacier, from sub-millimeter dust to boulders of a few meters diameter (Reference Mayer, Lambrecht, Beló, Smiraglia and DiolaiutiMayer and others, 2006), and this is typical of the whole Karakoram debris-covered area (Reference MihalceaMihalcea and others, 2008a; Reference Zhang, Hirabayashi, Fujita, Liu and LiuZhang and others, 2013). In the highest parts of Baltoro glacier, the debris occurs over medial moraines, then gradually spreads to a uniform cover across the entire surface, mantling most of the ablation zone. Near the terminus, the debris thickness may exceed 1 m (Reference Mayer, Lambrecht, Beló, Smiraglia and DiolaiutiMayer and others, 2006). Debris cover has been present on this glacier for at least the last century (historical data reported by Reference ConwayConway, 1894; Reference De FilippiDe Filippi, 1912).

Because of its size and the proportion of area covered by supraglacial debris (38% debris-covered and 62% debris-free), and its elevation (the equilibrium-line altitude (ELA) was reported to be at 5200 m a.s.l. by Reference Mayer, Lambrecht, Beló, Smiraglia and DiolaiutiMayer and others, 2006, and Reference BocchiolaBocchiola and others, 2011), Baltoro glacier can be considered paradigmatic of glacierized areas in the CKNP. It will therefore be our focus here for calibrating and validating the melt models.

In addition, permanent automatic weather stations (AWSs) operate on the glacier (at Urdukas, on a lateral moraine at 3926 m a.s.l., and Concordia, on the glacier melting surface at 4700 m a.s.l.), as well as in Askole (3029 m a.s.l.), the first village down-valley from Baltoro glacier (Fig. 1). These AWSs are part of the SHARE (Stations at High Altitude for Research on the Environment) network (an international project launched by the Ev-K2-CNR Chartered Association). These stations, managed with the support and agreement of the Pakistan Meteorological Department, provided input data for developing the melt models.

According to the Köppen–Geiger classification, this area falls in a cold desert region with a dry climate, little precipitation, a mean annual temperature lower than 18°C, and a wide daily temperature range (this type of arid, desert and cold climate is identified with the symbol BWk; Reference Peel, Finlayson and McMahonPeel and others, 2007). The Nanga Parbat massif forms a barrier to the northward movement of monsoon storms, which intrudes little in Karakoram. The hydrometeorological regime is barely influenced by monsoons, while a major contribution results from snow and ice melt. Precipitation is gathered in two main periods, summer (July–September) and winter (January–March), i.e. the seasons of the monsoons and the westerlies, the latter providing the dominant nourishment for glaciers. Some studies indicate that the total annual rainfall is 200–500 mm, as generally derived from valley-based stations, less representative for the highest zones (Reference ArcherArcher, 2003). Estimates from snow pits above 4000 m a.s.l. range from 1000 to >3000 mm w.e. (Reference Winiger, Gumpert and YamoutWiniger and others, 2005). However, there is considerable uncertainty about the behavior of precipitation at high altitudes. Recent studies (e.g. Reference Bocchiola and DiolaiutiBocchiola and Diolaiuti, 2013; Reference MinoraMinora and others, 2013) have highlighted climate and snowcover trends within this area in the period 1980–2009. Even if no significant changes of precipitation amounts have been detected, the number of rainy days appears to have increased. Regarding air temperature, the maxima have generally increased, while the minima appear to have decreased only during summer.

3. Methods

Since snow depth data in the CKNP area are scant and spotty, our study focused on modeling ice melt only, neglecting snowmelt. Following previous studies (Reference Mayer, Lambrecht, Beló, Smiraglia and DiolaiutiMayer and others, 2006; Reference MihalceaMihalcea and others, 2008a; Reference BocchiolaBocchiola and others, 2011; Reference SonciniSoncini and others, 2015), we set the ELA at 5200 m a.s.l., and we then applied the model only to glacier areas below that level, i.e. in the ablation zone. According to this criterion, the ablation zones of glaciers cover 3138 km², or 67.75% of the total glacierized area, in the CKNP. In addition, our analyses covered the period 23 July–9 August 2011, corresponding to the season when the largest glacier melt occurs during the year. The choice of limiting the application of the models to areas <5200 m a.s.l. may lead to underestimation of the actual glacier melt, as melt can occur above this elevation threshold, however limited in this season. Indeed, Reference SonciniSoncini and others (2015) analyzed the hydrological regimes of the Shigar river, covering ~7000 km² in the upper Karakoram, and nesting ~2000 km² of glaciers (including Baltoro, Biafo, Chogo Lungma), and found that snowmelt contribution was limited during our time frame (i.e. <20% considering the entire basin and not only glacier areas).

We applied two distributed melt models, one for debris-covered and one for debris-free areas. Both models were calibrated using field data gathered during an expedition on Baltoro glacier performed during 2011 (see Table 2 further below).

Table 1. Source, acquisition date and code scene of each satellite image used for the assessment of debris thickness distribution. Site displayed by each image is also reported

Table 2. Dataset used to calibrate and validate melt models. Dataset indicates whether ablation recorded at that site was used to calibrate (C) or to validate (V) the models; the site was debris-covered (DC) or debris-free (DF); Elev: elevation (m a.s.l.); X and Y: projected coordinates (WGS84UTM zone 43N); DR: debris effective thermal resistance (m² °C W^–1); DT-res and M-res: debris-thickness and melt residuals (modeled minus observed values); err: melt residual (%). The period considered is from the end of July to mid-August

To model the ice melting amount in the whole CKNP glacier ablation area, we considered the following input data:

1. 1. The glacier boundaries: the CKNP glacier inventory was derived by Reference MinoraMinora and others (2013), who applied remote-sensing investigations. More precisely, Landsat Thematic Mapper (TM) and Enhanced TM Plus (ETM+) scenes of 2010 were processed and analyzed (Table 1).

2. 2. A digital elevation model (DEM) describing the CKNP area (derived from the Shuttle Radar Topography Mission, SRTM3).

3. 3. A supraglacial debris cover map: a map describing the occurrence of supraglacial debris was obtained by applying a supervised maximum likelihood (SML) classification to a 2011 Landsat false-color composite image (i.e. 543 bands) (Table 1). This map allowed the separation of the debris-free and debris-covered zones of each glacier.

4. 4. Meteorological input data: the daily mean air temperature (T _a-Askole; °C) and the daily mean incoming solar radiation (SW_in-Askole;W m^–2) were obtained from hourly data measured during summer 2011 by the permanent AWS installed at Askole (Fig. 1). We used T _a-Askole and SW_in-Askole to evaluate ice melt over debris-free areas, by applying an enhanced T-index approach (Reference Pellicciotti, Brock, Strasser, Burlando, Funk and CorripioPellicciotti and others, 2005). SW_{in - Askole} was also considered in the debris-covered ice melt model in order to estimate the surface debris temperature (T _S-point; °C) driving the energy available at the debris–ice interface for ice melt.

5. 5. Debris data: (a) a map describing the supraglacial debris thicknesses (DT_point; m) was derived from Landsat TM thermal band imagery from August 2011 (T _S-Landsat; K); (b) daily surface debris temperatures (T _S-point; °C) in each pixel of the DT_point map and for each day in our period were computed by considering both daily incoming solar radiation data (SW_in-point; W m^–2) and debris thickness values (DT_point; m); and (c) the debris effective thermal resistance (DR_point; m² °C W^–1) was evaluated from debris thickness values by applying an empirical relation developed by Reference MihalceaMihalcea and others (2008a). The data in (a), (b) and (c) are the main inputs in the debris-covered ice melt model because they allowed estimation of the conductive heat flux through the debris layer and, consequently, of the ablation rate.

The total melting (M _TOT; m w.e.) in both debris-covered and debris-free ice zones was estimated as

(1)

where n and m are the total number of pixels (each pixel is 30 m × 30 m in size) of the digital image corresponding, respectively, to debris-covered and debris-free glacier areas, k is the length of the study period (days) and M _DC-point and M _DF-point are the melting rates over debris-covered and debris-free areas, respectively. M _DC-point and M _DF-point are fully described in Sections 3.3 and 3.4 respectively. The distribution of the meteorological parameters is reported in Section 3.1 and the evaluation of debris features is discussed in Section 3.2.

3.1. Distribution of the meteorological parameters

The meteorological data obtained from the AWS at Askole were corrected applying an altitudinal gradient to obtain estimated meteorological data on the whole glacier ablation areas of the CKNP. In particular, the daily mean air temperature (T _a-point; °C) was modeled by applying a constant lapse rate of –0.0075°C m^–1 (Reference Mihalcea, Mayer, Diolaiuti, Lambrecht, Smiraglia and TartariMihalcea and others, 2006). For the whole CKNP glacierized area, T _{a - point} was calculated as

(2)

where Δz is the difference in altitude (m) between the glacier pixel and Askole.

The daily mean incoming solar radiation (SW_{in - point}; Wm^–2) was estimated at each pixel based on the data gathered at Askole (SW_in-Askole) (Reference MihalceaMihalcea and others, 2008a):

(3)

3.2. Distribution of the debris thickness

There is no fully updated map of the glaciers in the CKNP. The first Park Glacier Inventory was only recently developed from Landsat imagery dating from 2010 (Reference MinoraMinora and others, 2013). The debris coverage within glacierized areas below 5200 m a.s.l. (i.e. the ablation zone) was assessed by the SML classification applied to Landsat TM images from 2011. This approach involved training the classification algorithm with a number of sites where the classification output (i.e. presence or absence of debris on the glacier surface) was known (Reference Brown, Lusch and DudaBrown and others, 1998). The SML algorithm assumes that values in each spectral band from Landsat TM are normally distributed and calculates the probability that a given image pixel is debris-covered or debris-free based on the values of all spectral bands. Each pixel is finally classified as debris-covered or debris-free according to the class that has the highest probability (Reference RichardsRichards, 1999). The details of the satellite images used are shown in Table 1. In particular, we used band combination 543 (as red, green, blue) of Landsat TM scenes to draw 20 regions of interest (ROIs) and trained the classifier. ROIs are sample areas that we know were covered by supraglacial debris in 2011. After training, the classifier was run on all the glacierized areas of the CKNP, assuming a probability threshold of 90% to separate debris-covered from debris-free pixels (i.e. a pixel was classified as ‘supraglacial debris-covered’ when the probability of a pixel belonging to this class was >0.9). The remaining pixels within glacierized areas and below the ELA were considered debris-free areas.

To map the thickness of supraglacial debris over the whole glacierized area of the CKNP, we used the method developed by Reference MihalceaMihalcea and others (2008b) for Miage glacier, Mont Blanc massif, Italy, and already applied to Baltoro glacier by Reference MihalceaMihalcea and others (2008a). This method is based on the relationship between surface temperature and supraglacial debris thickness (Reference Taschner and RanziTaschner and Ranzi, 2002). The input data are (1) debris thickness measured in the field on a wide and representative debris-covered glacier area and (2) satellite-derived surface temperatures. The empirical relationship between these data is a valuable tool for estimating debris thickness over unmeasured glacier zones (Reference MihalceaMihalcea and others, 2008a,Reference Mihalceab). This approach was initially developed on ASTER temperature data acquired on 14 August 2004 and applied to Baltoro glacier by Reference MihalceaMihalcea and others (2008a) (Table 1). Unfortunately, the ASTER images were not available for the whole CKNP area on the same date. We therefore modified the approach of Mihalcea and others (Reference Mihalcea2008a,Reference Mihalceab) to use Landsat TM images covering the entire CKNP area. To evaluate the suitability for debris assessment of Landsat TM images instead of ASTER ones, firstly we processed the Landsat image of the debris-covered portion of Baltoro glacier acquired on 14 August 2004, 05:18 GMT (10:18 local time), only 28 min before the acquisition of the ASTER image analyzed by Reference MihalceaMihalcea and others (2008a), and then we compared the results.

To assess surface temperature from Landsat images (T _S-Landsat; K), Landsat TM band 6 (i.e. thermal wavelength) digital numbers were first converted to radiance values (R _Landsat; Wm^–2 sr^–1 μm^–1) (Reference Coll, Galve, Sánchez and CasellesColl and others, 2010), and then T _S-Landsat was calculated applying the inverted Planck function:

(4)

where K ₁ and K ₂ are constant values (607.76 W m^–2 sr^–1 μm^–1 and 1260.56 K, respectively; Reference NASANASA, 2011), and " is the sky emissivity including atmospheric scatter (set to 0.95; Reference Barsi, Barker, Schott and SteinBarsi and others, 2003, Reference Barsi, Schott, Palluconi and Hook2005). The temperatures estimated using the two different images showed a good correlation (R ² = 0.91; mean, maximum and minimum temperature differences 2.1 K, 14.5 K and 0.0 K, respectively), thus supporting the use of Landsat data to describe supraglacial thermal conditions. Secondly, we used the same field data of debris thickness gathered in 2004 and used by Reference MihalceaMihalcea and others (2008a) to assess the best empirical function linking Landsat 2004 thermal data and debris thickness. The best-fitting function (R ² = 0.99) is

(5)

where DT is debris thickness (m) and T _S-Landsat is the Landsat-derived surface temperature. This equation is similar to that found by Reference MihalceaMihalcea and others (2008a) and describes the nonlinear relation between debris thickness and surface temperature. Moreover, we compared DT values obtained applying the equation reported in Reference MihalceaMihalcea and others (2008a) to 2004 ASTER data against the values derived from Eqn (5) on 2004 Landsat data on the Baltoro glacier area. The results (Fig. 2) show a good correlation between the two datasets (R ² = 0.85) and suggest a similar performance of the two models.

Fig. 2. Scatter plot showing DT values derived from 2004 ASTER (using equation reported in Reference MihalceaMihalcea and others, 2008a), and the ones from 2004 Landsat (applying Eqn (5)); 8700 pixels were considered.

Hence, these preliminary tests support the suitability of Landsat-derived surface temperatures to describe supraglacial debris thickness. We therefore used the debris thickness dataset collected in the field on the surface of Baltoro glacier during an expedition in July–August 2011 (a total of 57 samples ranging from a few centimeters to 2 m at the tongue). Regarding the Landsat surface temperatures, a single image covering the whole CKNP was not available; therefore, we used two images acquired on 10 August 2011 05:18 GMT and on 17 August 2011 5:24 GMT (Table 1). The images selected were particularly useful for our analyses because they were taken during the same period as the field measurements, and they partly overlap; they both cover part of the Baltoro glacier tongue (where field DT data were sampled). These data allowed us to assess two empirical equations linking debris thickness measured in the field to surface temperatures derived from Landsat images. The best-fitting equation (R ² = 0.75) obtained from the image taken on 10 August 2011 (which covers the whole Baltoro glacier area) was

(6)

while the one (R ² = 0.91) from the image acquired on 17 August 2011 (covering part of the Baltoro glacier tongue) was

(7)

We then applied Eqn (6) to thermal data derived from the Landsat image acquired on 10 August 2011, and Eqn (7) to thermal data derived from the Landsat image acquired on 17 August 2011. For the area covered by both overlapping images, results from Eqn (6) applied to the 10 August image were preferred because Baltoro glacier was only partially covered by the 17 August image, while it was completely covered by the 10 August image. Thus, the use of results from the 10 August image provided consistent estimates of the supraglacial debris thicknesses over the whole ablation area of Baltoro glacier.

3.3. Melt over debris-covered areas

The amount of ice melt under a debris cover (M _DC-point; m w.e.) depends on the energy available at the debris–ice interface and can be estimated as

(8)

where G _point corresponds to the conductive heat flux (W m^–2), Δt is the time step, _ρi is the ice density (917 kg m^–3) and L _m is the latent heat of melting (3.34 ×10⁵ J kg^–1). According to Reference MihalceaMihalcea and others (2008a), G _point can be estimated assuming a linear temperature gradient from the top of the debris layer to the ice surface for mean daily conditions (Reference Nakawo and YoungNakawo and Young, 1981; Reference Nakawo and TakahashiNakawo and Takahashi, 1982; Reference MihalceaMihalcea and others, 2008a):

(9)

where T _i is the ice temperature (set to the melting point, 0°C; i.e. we neglected refreezing phenomena, which generally do not occur during the main ablation season; Reference Mihalcea, Mayer, Diolaiuti, Lambrecht, Smiraglia and TartariMihalcea and others, 2006, Reference Mihalcea2008a) and DR_point is the effective thermal resistance of the debris layer (m² °C W^–1).

To derive DR_point over the whole debris-covered glacier area, an empirical relationship was applied (Reference MihalceaMihalcea and others, 2008a):

(10)

DR_point can be assumed constant over an ablation season as it mainly depends on debris thickness, which is generally considered stable over short periods (1–2 months; Reference FyffeFyffe and others, 2014). To model the daily mean debris surface temperature at each pixel (T _S-point), we considered both daily incoming solar radiation (SW_in-point) and debris thickness (DT_point), because higher radiation and thicker debris lead to higher surface temperatures (Reference Mihalcea, Mayer, Diolaiuti, Lambrecht, Smiraglia and TartariMihalcea and others, 2006, Reference Mihalcea2008a,Reference Mihalceab; Reference MayerMayer and others, 2010). T _S-point was estimated according to the empirical function

(11)

with a root-mean-square error (RMSE) of 2.1°C. This relation was based on field data of debris thickness and surface temperature sampled on Baltoro glacier during summer 2011 and incoming solar radiation estimated in the same gridpoints. Finally, the daily ablation (M _DC-point; m w.e.) at each pixel of the CKNP debris-covered glacier area was modeled as

(12)

where Δt is the number of seconds in a day (8.64 × 10⁴).

3.4. Melt over debris-free areas

The daily ice melt at each pixel with debris-free ice (M _{DF- point}) was estimated by applying an enhanced T-index model (Reference Pellicciotti, Brock, Strasser, Burlando, Funk and CorripioPellicciotti and others, 2005):

(13)

where T _a-point is the daily mean air temperature (°C), αis the surface albedo, SW_in-point is the daily mean incoming solar radiation (W m^–2), and TMF (32.43 10^–4 md^–1 °C^–1) and RMF (0.79 × 10^–4 md^–1 W^–1 m²) are the temperature and radiative melting factors, respectively. These melting factors are assessed from ablation measured at some selected sites on Baltoro glacier (from 3939 to 5200 m a.s.l.) from 23 July to 7 August 2011 (Table 2). Melting factors estimated from field data are taken as constant in time and space (Reference HockHock, 1999). Albedo was estimated by analyzing incoming and outgoing solar radiation data recorded during 2012 by a net radiometer (CNR1, Kipp & Zonen) installed at the Concordia supraglacial AWS. The spectral range considered was 0.3– 3 μm. Data show a high debris-free ice reflectivity, with a mean value of 0.30. Since distributed outgoing solar radiation data are not available for Baltoro glacier, we used the mean albedo value observed at Concordia. A lower ice reflectivity may be due to the presence of water, dust, debris and organic matter increasing the absorbed solar radiation and therefore the melting rate. Consequently, our assumption of a constant albedo equal to 0.30 may have led to a slight underestimation of the amount of meltwater, i.e. if the glacier ice is not completely clean in some parts of its debris-free area.

In any case, it is worth noting that other authors applying enhanced T-index approaches have also used constant albedo values (e.g. Reference Pellicciotti, Brock, Strasser, Burlando, Funk and CorripioPellicciotti and others, 2005).

4. Results

Modeled meteorological variables (Eqns (2) and (3)) agree well with those measured at Urdukas in 2011 and at Concordia in 2012. RMSEs regarding air temperature datasets are found equal to 1.2°C (for Urdukas) and 1.3°C (for Concordia), indicating that the local gradient by Reference Mihalcea, Mayer, Diolaiuti, Lambrecht, Smiraglia and TartariMihalcea and others (2006) can be considered accurate (Fig. 3). Modeled incoming solar radiation values resulted in a good match with the measured ones (Fig. 4), with RMSE values of 39 and 125 W m^–2 for Urdukas and Concordia, respectively.

Fig. 3. Daily mean temperatures recorded by the AWS installed at Urdukas during 2011 (x-axis) vs modeled daily mean temperatures (y-axis) obtained by applying a constant local lapse rate of –0.0075°C m^–1 to Askole temperatures (open box). The same analysis was performed for the Concordia dataset during 2012 (solid diamond).

Fig. 4. Daily mean incoming solar radiation recorded by the AWSs installed at Urdukas during 2011 and at Concordia during 2012 (x-axis) vs the modeled values (y-axis) derived from Askole data.

Supraglacial debris covers 518.47 km² (16.5%, 576 072 pixels) of the ablation zone of all the CKNP glaciers, while the extent of the debris-free area was 2619.61 km² (2 910 672 pixels). An example of debris occurrence is shown in Figure 5, where Panmah glacier (located in the central part of the CKNP, northwest of Baltoro glacier) is displayed. As regards the supraglacial debris thickness (Fig. 6), a mean value of 0.23 m was found, with maxima of ~3 m.

During the 2011 ablation season, we collected 29 measurements on Baltoro glacier (both debris-covered and debris-free conditions). We divided this dataset into two subgroups: one for calibrating our melt models and the other for validating them. Table 2 reports the two sub-datasets used to calibrate and validate the models.

Fig. 5. (a) Landsat TM image with RGB = 543 of portion of Panmah glacier. (b) The same image with pixels classified as covered by supraglacial debris highlighted in white.

Fig. 6. Map showing the supraglacial debris thicknesses over CKNP area, and a zoom on Baltoro glacier (lower left box).

The validation indexes display the performance of our models for estimating debris-free and debris-covered ice melt. In particular, we found a mean error of +0.01 m w.e. (corresponding to 2%) and a RMSE equal to 0.09 m w.e. (17%). In addition, we assessed any error due to the methodology applied for distributing the meteorological variables. For this purpose, we calculated the melt amount at selected debris-free (C-DF1, C-DF2, C-DF3) and debris-covered (C-DC1, C-DC2, C-DC3, C-DC4) ice field points varying the meteorological model inputs (T _a, T _S and SW_in) by their maximum RMSE (i.e. 1.3°C,± 2.1°C and ±125 W m^–2, respectively). Changing T _a and SW_in, the debris-free ice melt variations range from ±10% to ±25% (at higher altitudes); debris-covered ice melt instead shows differences around ±30% when changing SW_in, while variations in T _S drive a lower alteration around ±15%, not particularly influenced by elevation. Thus, the debris-covered ice melt model is more sensitive to the errors in the meteorological input data. However, debris-covered ice melt accounts for only 11% of the total melt. Moreover, these error tests were made considering the worst cases (maximum RMSE).

The debris-covered and debris-free ice melt models were therefore applied to the whole glacierized area of the CKNP below the ELA. The results are shown in Figure 7.

Fig. 7. Ablation map of CKNP glaciers below the ELA in the period 23 July–9 August 2011, and a zoom on Baltoro glacier (lower left box).

Given that the solar radiation was used to estimate debris surface temperatures, affecting in turn conductive heat fluxes, melt in debris-covered areas (M _DC) was largely linked to incoming solar radiation (SW_in). Indeed, the minimum and maximum daily melt (0.005 and 0.089 m w.e., respectively) occurred during days with the lowest and highest incoming solar radiation (respectively, 112 and 371 W m^–2, in Askole; Fig. 8a). Conversely, melting in debris-free areas showed extreme daily values (0.009 and 0.110 m w.e.) in days with extreme air temperatures (respectively +14.1°C and +22.7°C recorded at Askole; Fig. 8b). Overall, the greatest ablation occurred on 5 August, when incoming solar radiation was high, but not the highest, while the minima occurred on days (28–31 July) with minimum radiative input.

Fig. 8. Daily meltwater production from 23 July to 9 August 2011 from all the CKNP glaciers over the debris-free (DF) and debris-covered (DC) areas and the total (DC + DF). Same data are presented with (a) daily incoming solar radiation (SW) and (b) daily mean air temperature (Ta) recorded at Askole. Date format is mm/dd.

These findings indicate that (1) melt from the debris-covered parts of the glaciers (M _DC) is mostly influenced by the incoming solar radiation, since it depends on the conductive heat flux, and (2) melt of debris-free parts of the glaciers (M _DF) is more sensitive to air temperature.

On debris-covered areas of the whole CKNP, the daily average ablation was 0.024 m w.e. d^–1, while on debris-free areas it was 0.037 m w.e. d^–1. Considering both debris-free and debris-covered areas in the whole CKNP and the entire analyzed period, we estimated a total melt of 0.63 m w.e., corresponding to an average ablation of 0.035 m w.e. d^–1. Hence, over the period we considered, melting of the debris-covered parts of all the glaciers in the CKNP produced 0.223 km³ of meltwater (total M _DC), with a daily average of 0.012 km³ w.e. d^–1. The total meltwater from the debris-free parts (total M _DF) was 1.740 km³, with an average of 0.097 km³ d^–1. The total ice melt from the CKNP was thus equal to 1.963 km³ w.e., with a daily average of 0.109 km³ w.e. d^–1. This water volume equals ~14% of the reservoir capacity of the Tarbela Dam, a very large dam on the Indus River that plays a key role for irrigation, flood control and the generation of hydroelectric power for Pakistan (Reference ThompsonThompson, 1974). Table 3 shows a summary of the model results.

Table 3. Modeled melt rates over debris-covered (DC) and debris-free (DF) areas, and the total ablation in the period 23 July–9 August 2011

5. Discussion

Supraglacial debris thickness was derived from Landsat thermal data (60 m × 60 m pixel size), limiting the spatial resolution of debris-covered ice melt calculations. The results could, however, be acceptable given the extent of the analyzed debris-covered area (518 km²). The obtained DT values were cross-checked against a selection of field data, and a good fit was found (see Table 2). The main limitation comes from the fact that the supraglacial debris thicknesses derived from Landsat thermal data are average values at the pixel scale. The approach does not consider meltwater ponds, supraglacial lakes and sectors with crevasses and ice seals covering glacier areas smaller than the pixel size. Consequently, the model performs better in estimating debris layers thicker than 0.1 m (i.e. debris coverage is relatively continuous), while slight overestimation occurs for thin and sparse debris areas (<0.1 m; Table 2). The same limitation in DT modeling by means of remote sensing was found by Reference MihalceaMihalcea and others (2008a).

Mapping of debris thickness is fundamental for estimating debris resistivity, and therefore debris-covered ice melt. Other approaches have been proposed to produce debris thickness maps at higher resolution than ours (Reference Foster, Brock, Cutler and DiotriFoster and others, 2012), but they require meteorological data (including, among others, wind speed and direction and turbulent heat fluxes) on the glacier surface, as well as high-resolution DEMs (e.g. from lidar surveys), which were not available for glaciers in the CKNP area. Hence, our simple approach is suitable for investigating a wide and remote glacier area where high-resolution information is not available.

In the debris-covered ice melt model we assumed a linear vertical temperature profile within the debris layer (with 0°C at the melting ice surface), which is usually considered a good approximation for calculating daily melt rates (Reference Mattson, Gardner and YoungMattson and others, 1993; Reference MayerMayer and others, 2010). In the period under study, refreezing during the night can be considered negligible in the debris-covered areas, because hourly air temperatures at Urdukas were always positive. This is also in agreement with data collected on the tongue of Hinarche glacier, Bagrot valley, CKNP, at 2757 m a.s.l. from 26 July to 5 August 2008 (Reference MayerMayer and others, 2010), where air temperature (2 m above the debris surface) ranged between +14°C and +24°C, and never dropped below +9.9°C. When our model is applied to the onset and the end of the ablation season, it may overestimate meltwater discharge, and therefore may require further calibration.

We performed several sensitivity tests and evaluated model responses to varying input data at field survey sites (Tables 4 and 5) as well as over the whole CKNP ablation area (Table 6). First, we considered the debris-covered areas. We varied the daily incoming solar radiation by ±10% and ±20%. Then we studied the effect of varying the debris thickness upon melt results (±10%, 1 cm, ±5cm and ±10 cm with respect to the actual debris thickness values). The model response at field survey points (C-DC1 to C-DC4) is shown in Table 4.

Table 4. Sensitivity tests performed by applying different input data to the debris-covered ice melt model. We applied the model to four points where actual ablation data were collected in the field (the calibration points in Table 2) and calculated melt anomalies (ΔM) with respect to M _DC by modifying the incoming shortwave radiation and debris thickness. The reference modeled melt is given by M _DC mod

Table 5. Sensitivity tests performed by applying different input data to the debris-free ice melt model. We applied the model to three points where actual ablation data were collected in the field (the calibration points in Table 2) and calculated melt anomalies (ΔM) with respect to M _DF by varying the air temperature, the incoming shortwave radiation and the albedo. The reference modeled melt is given by M _DF mod

Table 6. Sensitivity test performed by applying different input data to both the debris-free and debris-covered ice melt models. The model results without input variation are shown in line 2 (M). We considered the whole CKNP ablation area the daily air temperature by ±0.1, ±1.0 and ±2.5°C with respect to the measured values. Finally, we investigated the effect of changing the albedo values by ±10%. Table 5 shows the model responses at field survey points (C-DF1 to C-DF3).

These tests suggest that changing the debris thickness or radiative input noticeably affects the debris-covered ice melt. In particular, this appears more evident in the presence of a thin debris thickness. Indeed, whenever shallow debris layers occur (see C-DC3 compared to C-DC1 in Table 4), even slight input variations entail evident changes in the underlying ice ablation, as the debris insulating effect is weaker.

Next, we considered the debris-free areas. We varied the daily incoming solar radiation by ±10%. Then we shifted the daily air temperature by ±0.1, ±1.0 and 2.5°C with respect to the measured values. Finally, we investigated the effect of changing the albedo values by ±10%. Table 5 shows the model responses at field survey points (C-DF1 to C-DF3).

The debris-free ice model is very sensitive to variations in air temperature, and the ablation varied by ±45% with changes of ±2.5°C. Minor impacts derived from changing SW_in inputs, showing a maximum variation of only 6%. This is a consequence of applying an enhanced T-index model, which indeed gives a primary role to temperature in driving ice melt, and a complementary role to incoming solar radiation (see, e.g., Reference Pellicciotti, Brock, Strasser, Burlando, Funk and CorripioPellicciotti and others, 2005). Concerning ice albedo, our model assumes a constant value of 0.30 for the whole area, thus probably entailing an over- or underestimation of the actual ice melt. Common albedo values for snow and ice surfaces range from 0.20 to 0.85; the albedo therefore has a very large and important influence on the total shortwave radiation absorbed by the surface, SW_in (1 – σ), and hence on ablation. In the absence of direct measurements, albedo is often estimated from ‘typical’ published values for snow or ice (Reference Cutler and MunroCutler and Munro, 1996): a clean ice surface generally features an albedo of 0.30–0.46, while a debris-rich ice surface is characterized by an albedo of 0.06–0.30 (Reference Cuffey and PatersonCuffey and Paterson, 2010). Thus, the choice of albedo is a very critical issue in accurately estimating the ice melt. In this study, we adopted the mean value (i.e. ~0.30) obtained by incoming and outgoing solar radiation data gathered by the supraglacial AWS placed at Concordia (in a debris-free area of Baltoro glacier). In previous studies, some authors applied similar approaches using an albedo of 0.30 (e.g. Reference Pellicciotti, Brock, Strasser, Burlando, Funk and CorripioPellicciotti and others, 2005). Reference OerlemansOerlemans (2001) reported a mean albedo value for debris-free ice of ~0.30. So we followed these previous studies supporting the use of a constant albedo of 0.30. The sensitivity test at field survey sites showed that changing the albedo by ±10% may lead to melt change of up to ±9% on debris-free areas (Table 5).

In addition to these model sensitivity tests, we considered the whole CKNP area totally debris-free, obtaining a total melt of 2.22 km³ w.e., with an increase of 0.48 km³ w.e. (more than twice as much) with respect to that obtained on actual debris-free areas (Table 6). This suggests that the debris layer is thick enough (more than the local critical value; Reference Mattson, Gardner and YoungMattson and others, 1993) to constrain the ice melt rates on average. To assess the effects of albedo, we changed the albedo of debris-free areas by a factor of ±10%, finding only a moderate impact on total melt (± 5%). Similar results were obtained by changing SW_in by ±10%. Moreover, stronger impacts (±7%) are caused by changing air temperature by ±1.0°C. Finally, we investigated the impact of DT by changing its values by ±10%, ±50% and +100%. In spite of the small impact on the total melt amount (+3.9% with –50% of DT and –3.2% with +100% of DT), the applied changes largely affected debris-covered ice melt. As the overall mean DT we derived from Landsat imagery (0.23 m) is surely higher than the local critical value (~0.05 m on Baltoro glacier according to Reference Mihalcea, Mayer, Diolaiuti, Lambrecht, Smiraglia and TartariMihalcea and others, 2006), the model is more sensitive to reduction than to increases of the actual DT value. This agrees with the well-known nonlinear relation between debris-covered ice melt and DT (see also fig. 7 in Reference Mihalcea, Mayer, Diolaiuti, Lambrecht, Smiraglia and TartariMihalcea and others, 2006). Indeed, when DT was decreased by 50%, melt in debris-covered areas increased by up to +34%, while when it was doubled, melt decreased by –28% (see Table 6).

6. Conclusions

We applied a simple model to evaluate the meltwater from debris-free and debris-covered ice ablation in the CKNP area, a wide glacierized zone of Pakistan. Our model estimates melt below 5200 m a.s.l. by applying an enhanced T-index model over the debris-free areas, and computing the conductive heat flux through the debris layer on the debris-covered zones. We neglected snowmelt, since snow data in the study area are not systematically available. We then focused on the peak ablation season, from 23 July to 9 August 2011, when meltwater is largely derived from ice melt, with snow thaw playing a minor role (Reference SonciniSoncini and others, 2015). Glacier features (i.e. surface area, supraglacial debris occurrence and thickness) were estimated from remote-sensing analysis of recent satellite imagery (2010– 11). Meteorological input data were distributed starting from data acquired at Askole. The data distribution procedure was validated by comparing the results with data recorded by two AWSs within the SHARE network (Urdukas and Concordia). The modeled ablation data were in strong agreement with measurements collected in the field during 2011 on Baltoro glacier, which can be considered representative of CKNP glaciers.

Our model estimated 0.223 km³ (on average, 0.012 km³ d^–1; min–max 0.006–0.016 km³ d^–1) of meltwater from the debris-covered parts, and 1.740 km³ (on average, 0.097 km³ d^–1; min–max 0.041–0.139 km³ d^–1) from debris-free sectors of the CKNP glacier ablation zone from 23 July to 9 August 2011. The total fresh water from the ablation areas of CKNP glaciers during the same period was therefore 1.963 km³ (on average, 0.109 km³ d^–1), corresponding, for example, to 14% of the water contained in a large strategic dam along the Indus River, of which all CKNP glaciers are tributaries.

The present model requires only a small number of input data, such as air temperature and SW_in (recorded by most of the standard AWSs), a DEM, and debris thickness measurements collected in the field. The relatively simple models we developed should provide portability to other regions, even if adjustments of the parameters against field measurements are necessary. In particular, (1) the lapse rate to distribute the air temperature (see Eqn (2)) should be locally evaluated; (2) the use of a constant albedo of 0.30 might be invalid for areas with debris-free ice affected by dust and black carbon deposition (see Reference Azzoni, Senese, Zerboni, Maugeri, Smiraglia and DiolaiutiAzzoni and others, 2014), thus requiring dedicated analyses; (3) the debris effective thermal resistance (DR) estimation requires debris-covered ice ablation and debris surface temperature data collected in the field.

The sensitivity tests suggest that melting will increase largely if summer air temperature increases. Also, any increase in the extent of debris coverage (which will likely occur due to augmented macrogelivation processes and rockfall events) will affect melt depending on new debris thickness. Thus, it will be important to monitor debris cover variations in time to update these crucial input data. Finally, albedo variations have to be properly considered, because surface darkening is reported as a result of increasing amounts of fine debris (Reference Oerlemans, Giesen and Van den BroekeOerlemans and others, 2009; Reference Azzoni, Senese, Zerboni, Maugeri, Smiraglia and DiolaiutiAzzoni and others, 2014). A further improvement of our approach will be the spatial distribution of debris-free ice albedo by applying methods based on remote-sensing investigations (see Reference Klok, Greuell and OerlemansKlok and others, 2003).