Spatial distribution of debris thickness and melting from remote-sensing and meteorological data, at debris-covered Baltoro glacier, Karakoram, Pakistan

Abstract

A distributed surface energy-balance study was performed to determine sub-debris ablation across a large part of Baltoro glacier, a wide debris-covered glacier in the Karakoram range, Pakistan. The study area is ~124km². The study aimed primarily at analyzing the influence of debris thickness on the melt distribution. The spatial distribution of the physical and thermal characteristics of the debris was calculated from remote-sensing (ASTER image) and field data. Meteorological data from an automatic weather station at Urdukas (4022ma.s.l.), located adjacent to Baltoro glacier on a lateral moraine, were used to calculate the spatial distribution of energy available for melting during the period 1–15 July 2004. The model performance was evaluated by comparisons with field measurements for the same period. The model is reliable in predicting ablation over wide debris-covered areas. It underestimates melt rates over highly crevassed areas and water ponds with a high variability of the debris thickness distribution in the vicinity, and over areas with very low debris thickness (<0.03 m). We also examined the spatial distribution of the energy-balance components (global radiation and surface temperature) over the study area. The results allow us to quantify, for the study period, a meltwater production of 0.058 km³.

1. Introduction

Calculating the distribution of ice melt over large areas of debris-covered glaciers represents a challenge in assessing water resources from glacier melt in remote areas such as Karakoram, Pakistan, and the Himalaya. Remote-sensing data are now available and can be used to analyze the spatial distribution of surface energy fluxes, along with traditional field studies that can be used to calibrate the new techniques.

While quite a few energy- and mass-balance studies have been performed on debris-free glaciers, studies on debris-covered ice are not numerous, especially for distributed energy- and mass-balance studies. Recently, Nicholson 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). Han 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 decade, a few papers have focused on debris-covered glaciers in the Himalaya and Karakoram (e.g. Hewitt and others, 1989; Mattson and others, 1989, 1993; Young and Hewitt, 1993; Nakawo and Rana, 1999; Kayastha and others, 2000; Nakawo and others, 2000; Takeuchi and others, 2000). Some studies have utilized remote-sensing data to analyze the spatial distribution of surface temperature to calculate the energy available for melting (Nakawo and others, 1993; Rana and others, 1997; Nakawo and Rana, 1999). Unfortunately these studies only calculate melt over small areas and short time-spans. One of the crucial input parameters for a successful application of melt models over large areas is a reliable debris-cover distribution. This can be obtained from remote-sensing data but still represents a challenge (Mihalcea and others, in press).

Several studies have dealt with the debris-covered area of Baltoro glacier and its tributaries in the past (Desio, 1954; Desio and others, 1961; Pecci and Smiraglia, 2000; Diolaiuti and others, 2003). More recently Mayer and others (2006) and Mihalcea and others (2006) analyzed the glaciological and meteorological characteristics of Baltoro glacier (glacier velocities, ablation rates and meteorological data from field expeditions), providing the observational and experimental basis for the present study.

Here we present the results of the application of a distributed energy-balance model calculation of sub-debris melt on a large Karakoram debris-covered glacier. The model predicts the magnitude of buried-ice melt using a debris thickness distribution derived from remote-sensing data and field meteorological data from a local automatic weather station (AWS), Urdukas (4022ma.s.l.), at Baltoro glacier.

The study time frame and area were determined by the 2004 Italian scientific–alpine expedition on Baltoro glacier ‘K2 2004 – 50 Years Later’. The study period is quite short (1–15 July 2004), but different data sources were available during this time span, permitting intercomparisons (i.e. satellite data and field investigations, including an AWS on the glacier surface) and model validation. In addition, our approach for calculating debris thickness and melt distribution had not been applied before on Karakoram glaciers, avoiding duplication of existing literature.

2. Study Site: The Baltoro Glacier System

Baltoro glacier, one of the world’s largest debris-covered glaciers, drains the south flank of the Karakoram range (35^˚35^'–35^˚56^' N, 76^˚04^'–76^˚46^' E) as a tributary to the Indus basin (Desio and others, 1961). The glacierized basin extends from 3370ma.s.l. at the terminus to 8611 ma.s.l. (K2 summit), comprising a drainage area of 1500 km². The actual glacier area (i.e. ice-covered, neglecting rock walls and steep firn slopes) is ~524km², and the ablation area in 2004 was ~372km². The longest flowline of Baltoro glacier reaches 62 km (Mayer and others, 2006). The main glacier tongue, which forms below Concordia at 4600ma.s.l, is debris-covered and oriented in an east–west direction. The thickness of the supraglacial debris ranges from a few mm to more than 1 m of continuous debris cover below 3900 ma.s.l. and towards the terminus.

In total, about 38% of the glacier area is debris-covered (Mayer and others, 2006). Below 5000 m, meltwater ponds and superficial streams exist, but only account for a few per cent of the glacier area. Therefore these areas of surface water are not treated explicitly in this study. Accumulation zones are generally found above 5200–5800m elevation (Young and Hewitt, 1993). The analyzed area extends from 3650 to 5400 ma.s.l., covering an area of ~124km² (equal to approximately 72% of the ablation area of the main Baltoro glacier excluding tributaries). This is the area where results are of most significance due to the abundance of debris cover and very high melt rates.

3. Data

For this study, a multidisciplinary approach was applied, combining meteorological data, field measurements and remote-sensing information. Meteorological data were available for the 2004 summer period from Urdukas AWS (Mihalcea and others, 2006). Urdukas AWS (SHARE network operated by Ev–K²–CNR Committee) was installed on a moraine ridge close to the left glacier margin and has run continuously since 18 June 2004. At Urdukas the meteorological parameters (air temperature, wind speed and direction, air pressure, humidity, precipitation, global radiation) are collected as hourly mean values.

Ablation data were collected during the 2004 summer season using 3m long stakes drilled in the ice across the Baltoro ablation area. The influence of varying elevation and debris thickness was taken into account in positioning the stakes. In particular, a specially designed experiment was performed on the glacier at site SF (Stake Farm) (Fig. 1), located in the central part of the glacier at about 4190 ma.s.l. Here ablation rates were measured for varying debris thickness and surface conditions (orientation and slope). For further details see Mihalcea and others (2006).

Fig. 1. Location of the study area, and schematic map of the Baltoro basin outlining the two major glacier branches where investigations were conducted.

The period of available field data is quite short, but measurements were carried out during the peak ablation season. This allows us to calibrate the model from low melt rates, during cloudy and cool days, to maximum melt rates. Therefore this period should be representative, as long as the continuous measurements of climate data are available as driving force. The still ongoing monitoring of climate parameters at Urdukas will allow an improved evaluation of changing climatic conditions within the next few years.

The debris temperature at different thicknesses close to the ablation stakes was measured by 20 single-channel 8-bit data loggers (Gemini Tinytag Plus). These 10 kΩ negative thermal coefficient thermistors operate over a –30 to +50˚C range. The surface temperature (T _s) was determined from debris temperature measurements by assuming a linear gradient through the debris layer at six locations spread over the analyzed area.

An Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) level 1B granule (14 August 2004, 05:46GMT, 1046 h local time) was processed for level 2 surface kinetic temperature (NASA processing-on-demand service (Global Land Ice Measurements from Space (GLIMS) project)). The spatial resolution of this product is 90 m, the temperature resolution is 0.5 K and the relative and absolute accuracies are 0.3 K and 4 K, respectively. A glacier mask was then used to restrict the subsequent calculations on the respective glacier area, where the debris thickness distribution and thermal resistance maps of the ablation area of Baltoro glacier were derived (Mayer and others, 2006).

All datasets were spatially referenced to the local Universal Transverse Mercator coordinate system (UTM zone 43N, WGS84) and were incorporated as grids into a Geographical Information System (GIS) database for analysis.

4. Methods and Results

This study aims at quantifying buried ice melt by evaluating the energy flux at the debris–ice interface through the analysis of remote-sensing data and field meteorological information. For our study we follow an approach mainly based on the use of debris surface temperature (T _s) and debris effective thermal resistance (R) to determine the energy available for ice melting (Nakawo and Young, 1981; Nakawo and Takahashi, 1982). The study of supraglacial debris conditions (T _s and R) was performed using remote-sensing (ASTER) and field (meteorological and reconstructed T _s) data, extended to the whole debris-covered ablation area of the main Baltoro glacier (without tributaries: 124 km²; see Fig. 1). The time-span we analyzed to calculate buried ice melt is 15 days (1–15 July 2004), since ablation rate data from field surveys were available for this period. This permits a comparison of modelled and measured ablation data to evaluate the reliability of our computations.

The pattern and distribution of supraglacial debris over the Baltoro glacier ablation area was derived from the surface kinetic temperature (ASTER, 14 August 2004, 1046 h local time) according to the approach introduced by Mihalcea and others (in press). Thus, a debris thickness map with 90 m pixel resolution was obtained.

To calculate the required effective thermal resistance of the debris layer, an empirical approach was applied: the relation between the 2004 field debris thickness (DT) and the effective thermal resistance (R), calculated from measured debris thickness, ice ablation and surface temperature (T _s) data, was evaluated (Mihalcea and others, 2006). A linear relation between these two parameters was found and was applied to the debris thickness map (derived from ASTER 2004) to calculate R over the ablation area on a 90 m grid.

The mean daily glacier surface temperature (T _s) was modelled starting from the available meteorological information, in particular global radiation (G) data recorded at Urdukas AWS. The daily mean of the local value of G was spatially distributed according to the approach introduced by Oerlemans (2001), to reconstruct the pattern and temporal evolution of this parameter over the whole glacier surface. An empirical relation was derived between G and the average daily T _s at the locations of the field measurements. The resulting equation was applied to the entire G dataset to calculate a daily mean T _s over the whole debris-covered ablation area (124km₂).

In order to employ a simplified heat flux model with the assumption of a stable and linear temperature gradient within the debris cover, a 24 hour interval was applied (Nakawo and Takahashi, 1982; Nicholson and Benn, 2006).

Starting from the maps of mean daily T _s data and R over the whole debris-covered ablation area, the energy available for melting (Q _m) at the debris–ice interface was calculated for each of the 15 days analyzed, and then the theoretical ablation value. In addition, for the glacier sectors without debris cover, a daily ablation rate measured in the field and adjusted according to the pixel elevation was applied. Within the GIS system a map could then be produced showing the accumulated ablation, over 15 days, of the whole Baltoro glacier debris-covered area on a 90m raster. The ablation data were compared with those collected in the field during the same period to verify the reliability of our method.

4.1. Debris thickness distribution derived from remote-sensing data

The debris thickness distribution is one of the key variables needed to calculate ablation over wide glacier areas. Field measurements are time-consuming and often very difficult in remote areas, so new remote-sensing techniques applied to derive the debris pattern are very useful. Studies have confirmed a good correlation between DT and T _s for a relatively thin debris layer (0.01–0.4 m) (Mihalcea and others, 2006, in press). When the supraglacial cover is thicker than 0.4m of continuous rock debris, the influence of the glacier ice on T _s ceases, which makes it difficult to calculate thick (>0.4 m) DT from remote-sensing data (Taschner and Ranzi, 2002; Ranzi and others, 2004).

For this study, the Baltoro T _s distribution was applied from the surface kinetic temperature data (ASTER, 14 August 2004, 1046 h local time) on 200m elevation bands. Our analysis showed that elevation dependence is not an important factor in the spatial distribution of T _s on Baltoro glacier, at least at the time of image acquisition (1046 h, 14 August 2004, clear-sky conditions).

Two datasets of debris field measurements on Baltoro glacier were available, which were used independently. One set was analyzed with respect to ASTER T _s, to derive an equation for calculating a debris thickness map. The second was applied to check and validate the calculated debris thickness map.

To derive the DT distribution over the ablation area of Baltoro glacier, an exponential relation between field measurements of DT and T _s data from ASTER was found for individual pixels and subsequently applied on the entire ablation area. The equation was obtained by using (i) the minimum and the maximum DT (measured in the field on extended areas at pixel scale) and (ii) several (40) DT field measurements performed in areas with homogeneous debris cover. These data were analyzed with respect to the corresponding T _s from the ASTER image. To derive the DT map, a threshold was applied for debris-free areas with T _s≤273.15 K. The approach we followed was the same as applied by Mihalcea and others (in press) on Miage glacier, Mont Blanc, Italian Alps, except that on Baltoro glacier it turned out that the separate treatment of elevation bands was not necessary. On Miage glacier, several linear functions (one for each 100 m elevation band) were found to best fit the T _s–DT relation, whereas in the case of Baltoro glacier the data analysis suggested an exponential function as a best fit (r = 0.9):

(1)

Equation (1) was applied to the ASTER T _s data for the entire study surface (124 km²) to estimate the debris cover thickness distribution over the area.

The DT distribution map obtained in this manner (Fig. 2) provides information about 90 m×90m averaging of debris thickness at pixel scale. Maximum values of 1–3m of DT were found close to the terminus, while low values of a few centimetres of debris occur on the upper part of the glacier. On Godwin Austen Glacier, Baltoro south and where tributaries join the main tongue, there are larger debris-free areas. On the DT map, a continuous layer of debris covers the glacier from Concordia downstream, with debris of 0.01 m up to a few metres thick. Medial moraines which form below Concordia, on Godwin Austen Glacier and Baltoro south, are also evidenced in the DT map (Fig. 2).

Fig. 2. Debris thickness distribution derived from ASTER surface temperatures.

Calculated debris thickness is generally greater than measured thickness. The average difference between calculated and measured values for 56 points is ~0.026 m. This difference seems acceptable considering the wide area covered by our analysis and the limited effect due to this small DT overestimation. These differences between measured and calculated DT can probably be explained by the DT spatial variability with respect to the pixel size. The DT variability at pixel scale was evaluated from SF stake-farm ablation measurements and DT data (Mihalcea and others, 2006). The SF location corresponds to a few individual and identified pixels. Survey points located within the same pixel present different debris thickness (Table 1).

Table 1. SF data for buried ice ablation (column 5) and DT (column 2). The measured debris data (point data in column 2 and average for each pixel in column 3) were compared with DT derived from ASTER T _s (column 4)

The average of field measurements within the same pixel yields a better result than single points, when compared to calculated debris (Table 1). Also there is spatial variability of the debris thickness within a pixel. At site SF11 the calculated debris matches the measured debris value.

4.2. Global radiation distribution

The global radiation (G) distribution was calculated from daily mean global radiation data collected at Urdukas AWS. The radiation distribution was calculated as the daily mean value for the period 1–15 July 2004 and also for the exact time of ASTER acquisition (05:46GMT+ 5 h, 14 August 2004). In order to distribute the G measured at Urdukas AWS over the studied area, a relation for elevation dependence as proposed by Oerlemans (2001) was applied:

(2)

where we assume that the daily mean atmospheric conditions are constant across Baltoro glacier. G is the global radiation calculated at each pixel, G _U is the daily global radiation at Urdukas, h is the elevation difference (pixel elevation – Urdukas elevation, 4022 m). The second term on the righthand side of Equation (2), 1, was used instead of Orlemans’ 0.79 because the value 0.79 is referred to sea level and to higher latitudes. In our study, we also use measured global radiation from Urdukas at lower latitude (~35.5^˚ N), i.e. at higher solar elevation.

To analyze topographic influences on the radiation budget, shading maps were produced for every 2 hours during daylight over the Baltoro glacier area to quantify shading and aspect influence on the G distribution. It was found that <5% of the study area is influenced by these two factors: at 0600 h, 6% of the analyzed area is shaded, at 1800 h only 3%. Thus, for the 24 hour mean values of the G distribution, we only took into account the elevation factor. At smaller time-steps (hourly), these factors must be considered, as a previous study (Mihalcea and others, 2006) demonstrated by comparing shortwave incoming (SW_in) radiation at two AWS located in the Baltoro glacier area (Base Camp AWS (5033ma.s.l.) and Urdukas AWS (4022ma.s.l.)). In the morning and afternoon, after sunrise and before sunset, SW_in at Base Camp is less than at Urdukas due to shading. However, during the hours of full exposition (1100–1500 h) the elevation factor prevails. Therefore when calculating the G distribution at 1100 h on 14 August 2004, only the elevation factor was considered.

Solar radiation values are in agreement between the two AWS (Mihalcea and others, 2006) during clear-sky hours, but become more complicated, with an irregular, patchy cloud cover; for 24 hour mean values, however, this effect can be neglected. Most small-scale shading effects are not resolved in the spatial resolution of our study (90m×90 m) and are therefore not included.

4.3. Surface temperature

Surface temperature (T _s) is a key element in calculating the energy available for melting of debris-covered ice (Nakawo and Takahashi, 1982; Nicholson and Benn, 2006; Mihalcea and others, in press). Previous studies have shown that T _s correlates with global radiation (G) on debris-covered areas (Mihalcea and others, 2006). We regressed T _s on G (local data from Urdukas AWS), at six sites spread over the debris-covered area and for the entire period of measurements (140–150 points), to provide linear relationships for each location. All the correlations (Pearson r values) are high: 0.7–0.93 (Fig. 3). This encouraged us to use distributed G values for estimating the T _s distribution.

Fig. 3. Surface temperature (T _s) vs global radiation (G) data and best-fit linear regression.

T _s is also dependent on debris thickness: higher DT increases T _s (due to low albedo and diminishing influence of the ice body underneath). Therefore a good estimation of the T _s distribution should incorporate both DT and G. T _s is also influenced by other factors such as air temperature, wind speed and material properties of the debris. A complete energy-balance study based on the surface energy equation may give better results of the melt distribution, but the spatial distribution of these parameters might be difficult to obtain. In this study we concentrated on factors that correlate well with T _s.

As mentioned above, the influence of DT on the magnitude of T _s needs to be considered. We analyzed the relation between the individual parameters of the linear equations, T _s vs G (for all six locations) and T _s vs local DT measurements. This resulted in a power law for the slope and a logarithmic law for the offset of the linear equation. The power law is not valid for the entire range of debris thicknesses; in fact, for areas with DT≥0.17 m, the calculated surface temperature was too high. Therefore, whenever DT≥0.17 m, a linear equation was applied for the slope, showing better agreement with the control values.

To test this method, the T _s distribution was calculated for the same moment as the ASTER image acquisition (05:46GMT, 14 August 2004) from distributed G values at 1100 h local time and the DT map. The resulting T _s map (Fig. 4) was compared with the ASTER surface kinetic temperature (T _ASTER) map.

Fig. 4. Calculated T _s distribution map (1100 h, 14 August 2004).

The differences (T _ASTER – calculated T _s) are in the range –6 to 3.76˚C, and the correlation between the two datasets is high (Pearson r value = 0.9). The average difference is 0.2˚C. For the most part, the differences are rather small, and the maximum deviation is found only in some areas: underestimated T _s values (2.5–3.7˚C) in zones with very thin debris (0.01–0.05 m), and overestimated T _s values (–3.0 to –1.3˚C) in debris-free areas and areas with thicker debris of 0.2–0.3 m. Differences T_ASTER – T _s2˚C account for 52% of the analyzed area. Both T _s maps show a similar spatial distribution pattern (Fig. 5).

Fig. 5. Map of differences between ASTER T _s map and calculated T _s map.

By following the method described above, daily mean T _s maps were calculated for each day during the 1–15 July period from mean daily G and from ASTER-derived DT distribution.

As the reflectivity of the supraglacial debris influences the amount of energy absorbed, it is important to state that the inferred surface temperatures and the model calibration are only valid for this specific glacier with its characteristic lithology. For application to other glaciers, field measurements are required which account for changes in surface reflectivity. This method for calculating T _s values yields good results on debris cover areas thicker than 0.02m (78% of the investigated area). For very thin debris layers at pixel scale, debris-free areas exist between the debris cover, and our method underestimates values of T _s mainly due to the cooling effect of the ice.

4.4. Effective thermal resistance distribution

To derive the effective thermal resistance (R) over the whole debris-covered glacier area, a linear relation was found by least-squares analysis between R and DT:

(3)

This relation is based on field data of ice melt, debris thickness and surface temperature sampled on Baltoro glacier during the 2004 summer (Mihalcea and others, 2006). Equation (3) was applied to the DT map in order to calculate an R distribution map (Fig. 6):

Fig. 6. Effective thermal resistance map of Baltoro glacier.

Maximum values up to 60 (10^–2×m^{2 ˚}CW^–1) are calculated at the terminus, where the debris is 1–3m thick. Lower values relate to thinner debris.

4.5. Estimation of energy available for melting (Q _M)

The main aim of this study is to calculate the spatial distribution of energy available for melting and accordingly the melt distribution on the debris-covered area of Baltoro glacier. This method has not been applied previously to such an extended area (124km²). The general energy-balance equation at the debris layer surface can be expressed as

(4)

where Q _s, Q _l, Q _h, Q _e and Q _c (Wm^–2) are net shortwave radiation flux, net longwave radiation flux, net sensible heat flux, net latent heat flux and conductive heat flux into the debris, respectively. All the terms are taken to be positive towards the debris surface (Nakawo and Young, 1981).

As a first approach to distribute the energy available for melting (Q _M which is given by the conductive heat flux Q _c in Equation (4)), a linear variation of temperature is assumed in the debris layer for daily mean conditions (Nakawo and Young, 1981):

(5)

where T _s is the debris surface temperature relative to melting and R is the effective thermal resistance of the debris layer (m^{2 ˚}CW^–1)

The energy used for ice ablation Q_M is calculated from

(6)

where L _f is the latent heat of phase change of ice (334×10³ J kg^–1), ρ _I is density of ice (~900 kgm^–3) and a is ablation rate in ice thickness (m s^–1). The Q _M distribution was calculated as daily values for the period 1–15 July 2004, from T _s and R of the debris layer. In the Q_M calculation, the variable that drives the daily result is the T _s distribution (and hence G), while the other variables remain fixed (DT and R).

The ablation amount (A, value in m) over the whole investigated debris-covered area (n pixels) in 1 day (day number k) was evaluated as the sum of the ablation amount calculated for each of the n pixels with debris cover and the ablation value of each of the m pixels with bare ice condition:

(7)

T _sk is the surface temperature evaluated for day k on a specific pixel n, R_n is the pixel thermal resistance and c is a constant equal to L _f ρ _I. The number of seconds in a day is 8.64×104, a _d is daily ablation rate measured on bare ice and F_h is elevation factor (found by analyzing field ablation data with respect to elevation, which permits us to adjust the ablation rate with respect to the pixel elevation). This model was applied for each day of k in the investigated period, and the sum of the 15 values obtained permits us to obtain the total ablation amount (Fig. 7; values in m).

Fig. 7. Calculated total ablation map (1–15 July; values in m).

For the debris-covered area, the maximum value of total ablation during the 15 days was 0.74 m in areas with relatively thin covers (10–12 cm). Minimum values are found on the upper part of the investigated area (due to lower temperatures and thus less melt energy) and in areas with bare ice or with very thin debris cover (<1 cm). The approach applied on debris-covered ice underestimates melt rates for very low debris thicknesses. The thinner the debris layer, the lower is the resulting T _s and consequently the energy available for melting (Fig. 8). In the case of a thick debris layer, our method yields better results. For areas with a debris cover >1m (at the glacier terminus), the ablation amount was still 0.15–0.17 m, which shows that also for very thick debris there is energy available for melting in a 15 day period.

Fig. 8. Calculated total ablation compared to measured total ablation at 56 different sites.

In the areas without debris cover (i.e. debris layer <0.01 m), according to the approach applied for bare ice areas the ablation amount during the 15 days varies from a minimum of 0 to a maximum of 1.37m at 5200 and 3900m elevation respectively.

Summarizing, the results obtained from Equation (7) allowed us to calculate for the study period an ablation value of 0.058 km³w.e. over an area of 124 km², resulting in a mean thickness change of ~–0.47 m. This value is obtained by summing the 0.055 km³ of water melt under the debris layer >0.03m (73% over the investigated area) and the 0.003 km³ of water melt on bare ice (i.e. debris thinner than 0.01 m over 15% of the study area). Calculated total ablation results were compared with measured ablation in 2004 at different sites (56 points). Total ablation calculated at each site at different periods (not all the ablation stakes were measured at the same time) was analyzed and evaluated with respect to the measured values. The Pearson correlation between the two series of data was 0.7 (Fig. 8).

We compared our calculated ablation rates (m d^–1) with stake-farm ablation values (Table 1). The comparison resulted in an underestimation (mean value –0.016 m) of buried ice ablation under a thin debris layer (≤0.03 m). In addition, in areas with very thin debris layer a high variability of debris-cover pattern may occur, strongly influencing the ablation rates (Table 1).

5. Discussion and Conclusions

The spatial resolution obtained for the melt calculations is limited by the pixel size (90m×90 m) and is acceptable given the extent of the analyzed debris-covered area (124km²). The supraglacial debris cover estimated from the analysis of ASTER surface temperatures shows a pattern of increasing thickness towards lower elevations. This pattern corresponds well with debris thickness distribution derived from field measurements (Mihalcea and others, 2006).

The main limitation is due to the fact that the supraglacial debris thicknesses derived from ASTER T _s represent an average value at the pixel scale. The approach does not consider meltwater ponds, supraglacial lakes and sectors with crevasses and ice seals covering glacier areas less than a pixel. In addition, a small overestimation of DT occurs, where the calculated DT in fact is slightly thicker than the measured values, and the average difference between measured and calculated differences at 56 points was ~–0.026 m. Nevertheless this difference is acceptable considering the wide area covered by our analysis.

A good estimation of the T _s distribution should incorporate both DTand G. T _s is also influenced by other factors such as air temperature, wind speed and material properties of the debris. As a control, the calculated T _s were compared to the ASTER temperatures for the day and the time of ASTER acquisition. The two datasets are highly correlated (Pearson r value = 0.9), and both maps (calculated T _s and ASTER surface temperatures) show a similar distribution pattern (Fig. 5).

Regarding the ablation, an underestimation occurs for areas covered by very thin mean debris thickness (i.e. 1 < DT< 3 cm), especially at lower elevation (e.g. ice sails, ice cliffs and crevassed areas). This is mainly because under thinner debris our model calculates lower T _s values which produce less ablation in the model. For these areas (equal to 12% of the whole study area) a different approach should be applied and further investigations are required.

In order to evaluate the sensitivity of our model to global radiation variability, a test was performed comparing ablation rates calculated during different days (i.e. with clear-sky and cloudy-weather conditions) and different debris thicknesses (varying from 3 to 50 cm). The radiation values varied between 96Wm^–2 (the minimum daily value in the time-span we considered) and 406Wm^–2 (the maximum daily value in our dataset). The mean daily ablation obtained by our method, considering all the debris thickness classes, was 0.8–4.7 cm (for 96 and 402Wm^–2 respectively). Thus a variation in mean daily global radiation of ~24% may affect the buried ice ablation in the same order of magnitude. Furthermore, during cloudy days the variability of ice ablation with debris thickness was smaller than during sunny days. These results agree with the principal approach of our model which permits ablation to be evaluated on the basis of the mean daily energy input and with a linear dependence of T _s on G.

In addition, our model does not consider the effect of liquid precipitation on buried ice ablation, which on cloudy days could increase ice melting, thus giving ice losses larger than those predicted by our model.

For areas with scarce or absent debris cover (DT < 0.01 m), ablation was evaluated by only considering field data (i.e. daily ablation rates on bare ice) collected during the study period (Mihalcea and others, 2006) and adjusting them with respect to their elevation. While this second approach permits these areas (where only 5% of the entire ablation occurred) to be taken into consideration, it does not consider the energy input causing ablation and/or the meteorological factors driving it. Usually a distributed degree-day method works well for these debris-free areas. Our model proved valuable and trustworthy for ~88% of the investigated area (i.e. the sector debris-covered with DT> 0.3, and the bare ice zones, 73% and 15% of the study area respectively).

Furthermore the comparison of calculated ablation vs the measured values yields a correlation value of 0.7. This demonstrates that the method we used for calculating melt distribution from remote-sensing and meteorological data, despite numerous simplifications, gives reliable results. A comparison of single points of measured ablation vs the mean calculated ablation at pixel scale is difficult, due to scale difference (1–3m² in the field against 8100 m² in the model). In a previous study (Mihalcea and others, 2006), the high variability of the ablation rate over small areas was outlined; in fact, at the SF location, with the same DT over areas of ~20–30m², different ablation rates were measured (Table 1, column 5).

In this study, the buried ice melt was quantified by evaluating the energy flux at the debris–ice interface through analysis of remote-sensing data and field meteorological information. Over the study period, an ablation of 0.058 km³w.e. over an area of 124 km² was calculated. This amount is equal to a mean thickness change of ~–0.47m and was obtained by summing the 0.055 km³ of melt under the debris layer >0.03m (73% over the investigated area) and the 0.003 km³ of melt on bare ice (i.e. debris thinner than 0.01 m equal to 15% of the analyzed area).

The approach developed in this paper could be extended to several other debris-covered glaciers in the high Asian mountains (Karakoram, Himalaya and Pamir) when satellite remote-sensing data and only limited field information (mainly meteorological data) are available.

A complete energy-balance study based on the surface energy equation may give better results for the melt distribution, but the spatial distribution of the involved parameters might be difficult to obtain. The method developed represents a further step in our research and will be used to analyze 1 year of the Urdukas AWS data.

Another important goal is to evaluate the accumulation amount on Baltoro glacier upper sectors; this information together with the energy balance will permit us to evaluate the glacier mass balance, giving important information on the yearly changes occurring on a large debris-covered glacier and helping us to understand its response to ongoing climate change.