1. Introduction

Since the 1960s, many efforts have been made to investigate the effect of supraglacial debris on ice ablation, due to its importance in evaluation of water resources and the morphological evolution of debris-covered glaciers. Observations of ice ablation under a debris layer suggest that a very thin debris cover may accelerate melting, while an increasing thickness of debris, when in excess of a threshold thickness of approximately 0.02 m, inhibits melting (Reference ØstremØstrem, 1959; Reference LoomisLoomis, 1970; Reference FujiiFujii, 1977; Reference Mattson, Gardner and YoungMattson and others, 1993). For numerical estimation of ice ablation under a supraglacial debris cover, Reference Nakawo and YoungNakawo and Young (1981, Reference Nakawo and Young1982) proposed an energy-balance model (herein referred to as Nakawo’s model), which takes the meteorological variables and physical properties of the debris material into consideration, to calculate the ablation under a thin debris layer. By computing the variables in an energy-balance equation and employing an empirical parameter of effective thermal resistance of the debris layer, Nakawo’s model can predict the melting rates beneath a thin debris cover. To date, this model and its variations have been successfully employed to estimate the rates of ice ablation in debris-covered areas on Peyto Glacier, Rocky Mountains, Canada (Reference Nakawo and YoungNakawo and Young, 1981, Reference Nakawo and Young1982), Langtang Valley, Nepal, Himalaya (Reference Rana, Nakawo, Fukushima and AgetaRana and others, 1997) and Khumbu Glacier, Nepal, Himalaya (Reference Kayastha, Takeuchi, Nakawo and AgetaKayastha and others, 2000; Reference Takeuchi, Kayastha and NakawoTakeuchi and others, 2000). Despite its good results in estimating ablation, the energy-balance model encounters some difficulties in evaluating the ice ablation over large debris-covered areas. One of the difficulties is that it might be unable to give a satisfactory prediction of ice ablation under a thick debris layer, when debris depth is greater than 0.5 m. In Nakawo’s model, rates of ablation are estimated with the assumption that the variation of heat stored in the debris layer can be neglected. Thus, ground heat flux is directly employed to calculate the rate of ice ablation under the debris cover without weakening during downward conduction. For many debris-covered glaciers (e.g. continental-type glaciers in the Tien Shan, China), areas with >0.5m thick debris cover account for more than 70% of the total debris-covered area. In such areas, calculation of ice ablation without taking the heat storage of the debris layer into account will inevitably give an overestimate. To date, few studies have considered the ice ablation under a thick debris cover. Reference KonovalovKonovalov (2000) proposed a numerical model to evaluate ablation under various thicknesses of debris cover, but we are not aware of any reports where his model has been tested against observed data. Another difficulty of Nakawo’s model is linked to the collection of meteorological data. External variables (e.g. wind velocity, air and surface temperature, relative humidity and solar radiation) are needed to calculate the net radiation flux, sensible heat flux and latent heat flux, and thus conductive heat flux from the energy-balance equation. However, the collection of these variables is difficult and time-consuming.

For the reasons discussed above, this paper aims to establish a simple and robust model for estimating ice ablation under various thicknesses of debris cover, especially under thick debris layers. The methodology of the model is based on the theory of heat conduction. Given surface temperature, depth and thermal properties of the debris layer, the proposed model can estimate mean temperature at interfaces of different debris layers, and then the heat flux for ice ablation. From this, rates of ice ablation under the debris layer can be calculated.

2. Numerical Solution To 1-D Debris Heat Transfer

One-dimensional (1-D) heat transfer in vertically homogeneous debris can be represented by (e.g. Reference SellersSellers, 1965)

(1)

where Tis the debris temperature (K), K ₀ is a constant thermal diffusivity (m² s^–1), T₀ is the initial debris temperature at the beginning of the calculation (K), t is the time (s) and z is the depth of debris (m). Note that Equation (1) is valid only in the absence of other heat sources or sinks, such as melting or freezing. Here the theory of fractional derivative/integral (Reference Miller and RossMiller and Ross, 1993) is employed to solve Equation (1) with the given boundary conditions. Following Reference Wang and BrasWang and Bras (1999), the vertical gradient of the debris temperature may be conveniently expressed in terms of the half-order time-derivative of T. Thus the (kinematic) debris heat flux Q(z,t) (Wm^–2), defined as positive downward, can be written as (Reference Wang and BrasWang and Bras, 1999):

(2)

where k₀ is the thermal conductivity (Wm^–1 K^–1), Co is the volumetric heat capacity (J m^–3 K^–1; C₀ = k₀/k₀) and s is the integration variable. Equation (2) indicates that debris heat flux, Q(z,t), at any depth, z, is completely determined by the history of the debris temperature at the same depth, given the thermal properties of the debris. The integration starts from a time when the debris heat flux is zero, as implied by a uniform initial debris temperature profile (Reference Wang and BrasWang and Bras, 1999).

Debris heat flux through a relatively homogeneous thermal debris layer can also be estimated from the debris temperature gradient within the debris layer (e.g. Reference Liang, Lettenmaier, Wood and BurgesLiang and others, 1994):

(3)

where T(z ₀) is the debris temperature at depth z₀. Given a time series of T(z₀) and relevant debris properties, T(z) can be iteratively solved using Equations (2) and (3).

It should be mentioned that the iterative procedure is started from the second time-step, since Equation (2) needs at least two sequential debris temperatures as its input, and the first T(z) must be initialized carefully (see section 4.2). Our experience is that the iteration converges rapidly. Errors in initial conditions, i.e. the uniform initial temperature profile and the constant debris thermal properties, are found to decay rapidly and have only a small influence on the model performance.

In the case of a thick supraglacial debris cover, several debris layers with different thermal properties should be classified for practical prediction of the debris temperature. Given surface temperature and thickness of debris cover, the proposed model can compute the debris temperature at various depths by applying the iterative procedure from the upper debris layer down to the subsequent ones. In the deepest debris layer, which is bounded by glacier ice at the bottom (assuming that the ice surface temperature is 0°C), debris heat flux can be calculated directly from Equation (3) and taken as the heat used for ice ablation, Q_m, with the assumption that heat conduction into the ice beneath is neglected. Then the ablation rate, r (cmd^–1), is calculated from:

(4)

where L_f is the latent heat of the phase change of ice (334 × 10³ J kg^–1) and p_i is density of ice (900 kg m^–3).

Notably, the proposed model assumes that heat flow within the debris cover is dominated by conduction, and it is justified by the relatively linear variation of the average temperature measured in the temperature profiles (especially for the deeper debris layers) on debris-covered glaciers (Reference Conway and RasmussenConway and Rasmussen, 2000). Caution is needed during winter when temperature probably increases with depth and buoyancy forces could induce air mixing through the debris layer (Reference Harris and PedersenHarris and Pedersen, 1998).

3. Study Site and Data Collection

The data used to test the model were collected from 15 September to 6 October 2004 on the debris-covered area of Koxkar glacier, northwest China. The location and surface cover of the glacier are illustrated in Figure 1. Supraglacial debris covers an area of about 19.5 km², which accounts for 83% of the total ablation zone, with thicknesses ranging from <0.01m on the upper reach of the ablation zone and on ice cliffs to >2.0m near the terminus. Meteorological conditions during the measurement periods are presented in Figure 2, in which air temperature and debris surface temperature were measured at the experimental sites and net solar radiation, precipitation, relative humidity and wind speed were recorded by an automatic weather station (AWS) established on the debris cover about 700m from the sites (Fig. 1). Air temperature was measured by a temperature and humidity probe (MP101A, ROTRONIC Instrument Corp.) mounted 1m above the ground surface. A thermistor (made by the Key Laboratory of Ice Core and Cold Regions Environment, Chinese Academy of Sciences) with an accuracy of 0.02°C was positioned a few millimeters under the debris surface to record the surface temperature. Both temperature sensors were connected to a data logger (DT 600, Datataker Pty Ltd) to provide 30 min average data records during the measurement periods. The meteorological data from the AWS were collected using the instruments and procedures reported by Reference Han, Ding and LiuHan and others (2005). Figure 2 indicates that there was an apparent temperature drop during 16–19 September, with low solar radiation and high relative humidity during the same period. A dryer period occurred during 27 September to 1 October, when the solar radiation and temperature were higher. The total precipitation during the experiment was 10.8cm, of which 77.8% occurred during 16–21 September; all precipitation was rain. During the experimental period, there was extreme variation in the wind speed, ranging from 0.14 to 7.37 ms^–1.

Fig. 1.

Location and map of study sites on Koxkar glacier, northwest China.

Fig. 2.

Meteorological conditions recorded by the AWS on days 259–280 (15 September to 6 October, 2004). (a) Air temperature (black line, observed at 1m above the ground surface) and debris surface temperature (gray line). (The gap in the temperature record from 1400 h on day 263 to 1800h on day 264 is due to measurement interruption when the experimental equipment was transferred from site 1 to site 2.) (b) Relative humidity (at 1m above the ground surface). (c) Net solar radiation (at 1.5m above the ground surface) and precipitation. (d) Wind speed (at 1m above the ground surface).

The debris temperature measurements were carried out at three experimental sites (Fig. 1) where the debris consists mainly of coarse sand (grain size ≈0.5 mm). The thicknesses of debris cover at the experimental sites were 1.2m (site 1), 2.0m (site 2) and 0.7m (site 3). Seven thermistors (the same as those used for surface temperature measurement) were used to measure the vertical temperature profile within the supraglacial debris, providing data for model input and model evaluation, as well as for thermal property determination. The traditional approach to observing ice ablation under a thin debris cover is to settle a number of stakes into the ice beneath and read the height variation of the debris cover periodically (the frequency ranging from a couple of days to several months). For ice ablation measurement over a period of a few hours, under a thick debris cover (e.g. >0.5 m), which we require to evaluate our model, this approach is inappropriate, as the ice ablation may be too small to be precisely detected by such a manual measurement and there will be disturbance from an uneven debris surface. We therefore placed a heat flux plate (CN-81, EKO; accuracy: 25.182 W m^–2 mV^–1) a few centimeters above the glacier ice to record the heat flux near the debris-ice interface, assuming that the overall heat flux at this depth was responsible for the ice ablation. The thermistors and heat flux plate were also connected to the data logger (DT 600, Datataker Pty Ltd). A 30 min average of debris temperature and debris heat flux with 10s samples was obtained during the measurements.

Measurements of density, porosity and water content of the debris at several depths for each site were conducted during the excavation for the installation of the sensors. Coupled with the temperature profile measurement within the debris cover, these data assisted in determination of the thermal diffusivity, the thermal conductivity and the volumetric heat capacity using the temperature profile procedure proposed by Reference Conway and RasmussenConway and Rasmussen (2000).

Owing to a shortage of instruments, the measurements were conducted over different periods of time for each site. The temporal coverage is 15–19 September for site 1, 24–28 September for site 2 and 2–6 October for site 3.

4. Application and Discussion

4.1. Layer classification within the debris cover

For the case study on Koxkar glacier, analysis of the temperature profile and variations of the thermal properties through the debris suggests a three-layer scheme for the model calculation. As indicated in Figure 3, layer 1 is characterized by a large diurnal variation in debris temperature due to strong impact from solar radiation, and the thermal conductivity of the layer is relatively small. The amplitude of the diurnal signal in layer 2 is decreased considerably, and the temperature varies linearly with a consistent positive temperature gradient. Layer 3, which is adjacent to the debris-ice interface, replicates the characteristics of layer 2 with lower temperature and more even temperature changes, except for a higher thermal conductivity (and volumetric heat capacity) due to high water content as a result of ice ablation. In this study, we set the thickness of layer 3 to be 20 cm uniformly for all the experimental sites, according to the observations of the temperature and water content as well as the related thermal properties. This is also justified by the good model predictions presented below. Figure 4 shows the schematic structure of the three layers, and Table 1 lists the depths (z₁, z₂ and z₃), the thermal conductivities (k₁, k₂ and k₃) and the volumetric heat capacities (C₁, C₂ and C₃) for the layers in each site that we used in the model calculation. Given a time series of ground surface temperature, T _s, and depth of debris cover, z₃, debris temperatures T₁ (at z₁) and T₂ (at z₂) can be solved using the procedure described above. Subsequently, heat used for ice ablation, Q_m, and ablation rate, r, could be calculated from the model. The time-step of the modeling is 30 min.

Fig. 3.

Representative debris temperature profiles (1400h (line a), 0400 h (line b), 1000 h (line c) and 1800 h (line d)) (all local time) and estimated thermal conductivities (gray line with solid stars) at site 1, where variations determine the three-layer classification scheme adopted in the model.

Fig. 4.

Schematic representation of thermal debris layers in the case study of Koxkar glacier. The depths and thermal parameters for each experimental site are listed in Table 1.

Table 1.

Depths, z, average thermal conductivities, k, and average volumetric heat capacities, C, of each debris layer for the three sites

4.2. Temperature modeling

Figures 5 and 6 show the modeled (dashed line) and the observed (solid line) debris temperatures T₁ and T₂ for sites 1, 2 and 3. In general, the modeled debris temperatures agree well with the observations, except for an apparent phase shift. We believe that this error is caused by the unsynchronized variations between the debris temperature and the heat flux through the debris layer during the calculation. When the temperature at the top of the layer experiences rapid changes, the temperature difference through the layer changes correspondingly, while the heat flux through the layer changes less (or not at all) due to the buffering effect of the debris on heat conduction. As a result, a phase shift occurs between observed and modeled temperatures. The phase differences between variations of temperature and heat flux seem to be related to the magnitude of the time-step. We tested the model with time-steps of 30 min and 2, 6 and 12 hours, in order to examine the effects of time-step duration on the phase shift in the calculation. Figure 7 shows that the phase shift between observed and modeled temperatures can be partly removed by using a larger time-step in the modeling. The disadvantage of adopting a large time-step is that the detail in the temperature variation is reduced. For some applications (e.g. in some permafrost studies on diurnal variations of frozen soil) this can be important. In this study, with the objective of ice ablation modeling, we can reduce the influence of the phase shift in the temperature prediction on the ablation rate by removing it through a calibration procedure. The phase difference between modeled and observed debris temperature is determined in the calibration period and, for the subsequent prediction period, the phase lead can be applied directly to the model outputs to correct for phase.

Fig. 5.

Comparison of modeled (dashed line) and observed (solid line) debris temperature T₁ at Koxkar study sites for days 259–280 (15 September to 6 October 2004). (a) Site 1, at z₁ = 0.5 m; (b) site 2, at z₁ = 0.5 m; and (c) site 3, at z₁ = 0.3 m.

Fig. 6.

Comparison of modeled (dashed line) and observed (solid line) debris temperature T₂ at Koxkar study sites for days 259–280. (a) Site 1, at z₂ = 1.0 m; (b) site 2, at z₂ = 1.8 m; and (c) site 3, at z ₁ = 0.5 m.

Fig. 7.

Comparison of observed 71 with 71 modeled with different time-steps at site 1.

In addition to the phase shift indicated by Figures 5 and 6, the amplitude of the modeled temperature is larger than the observed temperature, particularly at sites 2 and 3. This may contribute to overestimation of ablation at these sites. Corresponding to Figures 5 and 6, Tables 2 and 3 summarize the statistics for observed and modeled temperatures at each site. It can be seen from the tables that the mean amplitude of the modeled temperature is greater than that of the observed temperature by 9–67%, especially at site 2, where these values are 42.9% and 66.7% for T ₁ and T ₂, respectively. Such a large deviation may be due primarily to systemic errors in the model. Firstly, the ability of the model to predict the temperature profile of the debris cover is based on the measurement of the time series of surface temperature. However, due to the large impact of solar radiation, the surface temperature often exhibits strong diurnal variation, which is passed to the modeled debris temperature through the iterative procedure, resulting in larger amplitudes of the modeled temperatures. Secondly, the observed T ₁ and T ₂ have less variation: for example, at site 2 T ₁ ranges from 2.14 to 3.45°C and T₂ from 0.55 to 0.81 °C. Thus any small fluctuation in the modeled temperature may lead to a large deviation between the amplitudes of modeled and observed temperatures. Despite this seemingly large deviation, it is found to have only a small influence on the estimation of ice ablation. Tables 2 and 3 also show that the differences in maximum, minimum and mean temperatures between modeled and observed temperatures fall into an acceptable range.

Table 2.

Summary of the statistics for modeled and observed debris temperature T ₁ at each site

Table 3.

Summary of the statistics for modeled and observed debris temperature 7 ₂ at each site

Other aspects that may also affect model accuracy include (i) debris heat flux estimation from time series of debris temperatures by the half-order derivative/integral relations and (ii) debris temperature initialization to start the iterative procedure.

Solution of the 1-D heat diffusion equation by the halforder derivative/integral method (Equation (2)) constructs a unique link between debris temperature and debris heat flux at a given location (Reference Wang and BrasWang and Bras, 1999). The validity of Equation (2) depends on two seemingly restrictive assumptions: firstly, that the thermal properties of the debris are uniform, independent of debris moisture and debris temperature, and, secondly, that there is a time when the debris temperature profile is uniform. In the first evaluation of the half-order derivative/integral method to estimate ground heat flux from surface debris temperature and vice versa (Reference Wang and BrasWang and Bras, 1999), the impact of these assumptions on the accuracy of the solutions was discussed in detail. It was concluded that the assumptions of constant thermal properties and uniform initial debris temperature profile would not adversely impact on the usefulness of Equation (2). The results of the Koxkar case study support this argument.

Initialization of debris temperature, T(z), in Equation (3) is found to affect model accuracy considerably. Care must be taken in choosing the initial value of T(z), which the case study shows should be taken as the observed debris temperature following the duration of a phase difference between the modeled and observed debris temperatures. For example, if the iteration starts at midnight and the determined phase shift during the calibration period is 9 hours, then the observed debris temperature at 0900 h is taken as its initial value.

4.3. Ablation rate modeling

Figure 8 shows the modeled ablation rate (dashed line) and that derived from the debris heat flux measurement (solid line). To present a better comparison, the phase shifts between observed and modeled data are removed through the procedure given above. Table 4 summarizes the statistics of Figure 8. It shows that, despite a slight difference between modeled and observed ablation rates, the diurnal variability of ablation rate is well captured by the proposed model. The apparent overestimate of ablation rates by the model, as seen in Figure 8b and c, may be associated with a slight overestimate of the debris temperature. According to the numerical procedure, Equation (3) is directly used to calculate the ablation heat, Q_m, and in turn the ablation rate, r, when the interface temperature T ₂ is determined. Therefore, the modeled ablation rate varies linearly with T₂, leading to the replication of the overestimate from T ₂ to r. A larger amplitude in the estimate of T₂ may similarly contribute to an overestimate. Another reason for the difference between modeled and observed ablation may be related to errors during the ablation heat measurement. Due to the impact of the large volumetric heat capacity from the high water content found at the bottom of the debris cover, the observed heat flux may be less than the incoming heat flow from the upper debris layer, since part of the incoming heat is absorbed as heat storage of the wet debris layer. This under-measurement of the ablation heat may result in a slight overestimate of the modeled ablation rate compared with the observed value.

Fig. 8.

Comparison of modeled (dashed line) and observed (solid line) ablation rates, r, at Koxkar study sites for days 259–280 (15 September to 6 October 2004). (a) Site 1; (b) site 2; and (c) site 3. The phase shifts have been removed by phase correction.

Table 4.

Summary of the statistics for modeled and observed ablation rate, r, at each site

Despite the dominant trend of overestimation in the modeled ablation rate at sites 2 and 3, Figure 8 and Table 4 also show an apparent underestimate at site 1, particularly during 16–18 September. It is inferred from the surface temperature variations of the debris cover (Fig. 2a) that this underestimate may be related to the large temperature drop in the same period, which is also reflected in the modeled ablation rate through the iterative procedure.

4.4. Comparison with Nakawo’s model

In order to compare the performance of our model with that of Nakawo we tested Nakawo’s model against the observed ablation data. The detailed numerical solution of this model is found in the works of Reference Nakawo and YoungNakawo and Young (1981, Reference Nakawo and Young1982). Figure 9 shows the differences between the modeled subdebris ablation rates and the observed data for Nakawo’s model and for our model. It clearly shows that the difference between our model and the observed data is slight, fluctuating around 0 cmd^–1. In contrast, the curve generated by Nakawo’s model has large fluctuations, indicating a larger deviation between the modeled results and the observed data. This error is primarily due to poor knowledge of the heat storage within the debris cover, particularly when the debris cover is thick. In addition, in Nakawo’s model the thermal resistance, R, defined as the ratio of the surface temperature to the conduction heat flux (which corresponds to the ground heat flux in our model) through the debris layer, is assumed to be constant for a given depth of debris cover. Although the ground heat flux and the debris surface temperature generally exhibit similar trends in diurnal variations, the numerical relation between these two meteorological variables can rarely be described in a linear manner. One disadvantage of making this assumption is that the modeled ablation rate will inherit the strong diurnal variability of the surface temperature, as illustrated in Figure 9. Other factors that can contribute to the discrepancy between the model and the observations are fully discussed by Reference Nakawo and YoungNakawo and Young (1982).

Fig. 9.

Comparison of the differences between observed and modeled ablation rates, r, for our model (black line) and for the Nakawo model (gray line). (a) Site 1; (b) site 2; and (c) site 3.

The cumulative ablation curves, presented in Figure 10, may be of more value for practical applications, since workers are usually interested in ablation over the long term, often an entire season. Figure 10 indicates that the accumulated ablation predicted by Nakawo’s model is several times larger than observed data over the experimental period, particularly at sites 2 and 3. This implies that an even larger discrepancy will occur over a longer period. In comparison, the observed sub-debris ablation is well predicted by our model, and the accumulated error is less than 27.5% at the end of the calculation.

Fig. 10.

Comparison of accumulated ablation: observed (solid line), modeled by our model (dashed line) and modeled by Nakawo’s model (gray line). (a) Site 1; (b) site 2; and (c) site 3.

Nakawo’s model is a useful and effective tool for computing the ice ablation under a thin debris layer. Work by Reference Nakawo and YoungNakawo and Young (1981, Reference Nakawo and Young1982) and many others gives solid support for this statement. However, major improvements must be made for the model to perform well when modeling ice ablation under a thick debris cover.

5. Conclusion

The simple model proposed in this paper develops a new method of estimating ice ablation under a supraglacial debris cover. The key method employed is to establish a link between debris heat flux and debris temperature at a certain depth when the heat transfer in the debris layer is described by a diffusion process. This method has been successfully applied in the estimation of sensible heat flux (Reference Wang and BrasWang and Bras, 1998) and ground heat flux (Reference Wang and BrasWang and Bras, 1999). As demonstrated by the numerical solution, the model is essentially an effective debris temperature calculation procedure. Given a time series of surface temperature, debris thermal properties and relevant boundary conditions, the model can find the debris temperature at any depth using an iterative procedure. Subsequently, the ablation rate of a debris-covered glacier can be estimated from the temperature gradient to the bottom of the debris layer. The case study on Koxkar glacier shows that, in general, the proposed model gives good results for the prediction of debris temperatures and ice ablation rates. However, the periods of debris temperature measurement during the experiment are short and represent the minimum for adequately validating the model. Further efforts must be made through a series of studies to improve and justify the model’s ability to predict ablation under thick debris cover over a longer period, in particular, on the temporal scale from months to seasons, during which the accumulated errors from the modeling may be considerable and need to be handled properly.

The major advantage of the proposed model is that it requires only a few parameters to model debris temperatures and ablation rates for a debris-covered glacier. This could be particularly important for studies carried out in regions where there is limited availability of observed data.

The proposed model, with its ability to predict the temperature profile in supraglacial debris cover, may find further use in a wide range of studies, for example, research on the thermal stability of frozen soil. Although some problems (such as phase shift) have not yet been completely understood, further research will aim to overcome these problems and exploit the full potential of the model.