Introduction

To assess the relation between climate and ice ablation (mass loss) on the mostly flat, upper surfaces of the remaining tabular-shaped glaciers on Kilimanjaro, Tanzania, it is necessary to have a thorough knowledge of the surface energy balance (SEB):

where M is energy available for melt (M = 0 if the surface temperature T _s < 273.15 K), S↓ and S↑ are incoming and reflected shortwave radiation, L↓ and L↑ are incoming and emitted longwave radiation, H and λE are the turbulent fluxes of sensible and latent heat (where λ is the latent heat of vaporization (2.501×10⁶ J kg^–1) when the surface is melting or the latent heat of sublimation (2.848×10⁶ J kg^–1) when the surface is frozen, and E is the evaporation (sublimation)/deposition rate), G is the subsurface conductive and radiative heat flux, S _N and L _N are net shortwave and longwave radiation and R _N is net all-wave radiation. The sign convention used is that all fluxes are positive when the surface is gaining energy.

The retreat of glaciers on the summit plateau of Kilimanjaro (≥5700 m) differs from that of other tropical glaciers because they have high, nearly vertical walls at their margins (Reference Hastenrath and GreischarHastenrath and Greischar, 1997; Reference Mölg, Hardy and KaserMölg and others, 2003). Retreat of these tabular-shaped ice bodies is governed mostly by S _N-induced ablation on their 10–30m high vertical walls (Reference Mölg, Hardy and KaserMölg and others, 2003). Glacier retreat along the margins of the vertical walls is irreversible, and no change in 20th-century climate has altered their continuous retreat (Reference Cullen, Mölg, Kaser, Hussein, Steffen and HardyCullen and others, 2006). Thus, to extract a climate signal from the tabular-shaped plateau glaciers on Kilimanjaro, and to fully understand their total mass balance, an investigation of energy and mass fluxes on their flat, upper surfaces is required.

SEB model results over a horizontal surface on the largest ice body on Kilimanjaro showed that R _N dominates energy exchanges at the glacier–atmosphere interface, governed by variations in S _N (Reference Mölg and HardyMölg and Hardy, 2004). Spatial and temporal variability in ablation is controlled mostly by changes in albedo (α), which is dependent on precipitation amount and frequency. Importantly, λE is almost always directed away from the surface as mass loss through sublimation, and is large enough to account for most of the observed ablation. If surface melting does occur on any given day, it is typically of short duration (a few hours). Mean sublimation over 19 months was estimated to be 0.92 kgm^–2 d^–1 (Reference Mölg and HardyMölg and Hardy, 2004), which is comparable to rates on other high mountain glaciers in the tropics (Reference Wagnon, Ribstein, Francou and PouyaudWagnon and others, 1999, Reference Wagnon, Sicart, Berthier and Chazarin2003; Reference Favier, Wagnon, Chazarin, Maisincho and CoudrainFavier and others, 2004b). Because sublimation is about eight times less effective than melt as an ablation process, its dominance on the horizontal surfaces of the plateau glaciers is important in limiting total mass loss on Kilimanjaro.

Determining the magnitude of sublimation is dependent on the accuracy of the calculation of λE, which is typically estimated using the bulk method in SEB models on high-altitude tropical glaciers (e.g. Wagnon and others, 2003; Reference Favier, Wagnon and RibsteinFavier and others, 2004a, Reference Favier, Wagnon, Chazarin, Maisincho and Coudrainb; Reference Mölg and HardyMölg and Hardy, 2004; Reference Sicart, Wagnon and RibsteinSicart and others, 2005). SEB models require automatic weather station (AWS) data as input, which for logistical reasons are often difficult to obtain accurately on tropical glaciers. There still remains much uncertainty about stability functions and surface roughness lengths that are used in the bulk method, particularly on tropical glaciers (e.g. Reference Sicart, Wagnon and RibsteinSicart and others, 2005).

The objective of this paper is to assess the uncertainty of calculating sublimation on Kilimanjaro using a SEB model. To achieve this we compare modelled (bulk method) turbulent heat fluxes to eddy covariance data collected over a brief period on the largest remaining ice body on Kilimanjaro. The eddy covariance method is the most direct procedure to measure turbulent heat fluxes and other key scaling parameters in the near-surface atmosphere (e.g. Reference MunroMunro, 1989; Reference Kaimal and FinniganKaimal and Finnigan, 1994). Scaling parameters obtained from the eddy covariance measurements are used to obtain site-specific surface roughness lengths, which in turn are used as input in the SEB model to improve estimates of sublimation. The uncertainty of modelled sublimation is further assessed by applying offsets to input data and model parameters over a series of SEB model runs, which enables us to describe with improved confidence the relative importance of sublimation in the ablation process on the remaining glaciers on Kilimanjaro.

Data Treatment and Model Description

Eddy covariance data used in this study were collected during the long dry season (June–September) in East Africa between 26 and 29 July 2005 on the tower of the existing Northern Ice Field AWS on Kilimanjaro (3˚04′ S, 37˚21′ E; 5794 ma.s.l.). The nearly horizontal surface at the site of the AWS has an average surface slope of about 2˚ and a fetch of several hundred metres in the prevailing wind direction. Single-level measurements of air temperature, relative humidity and wind speed from the AWS were used, as well as all four single components of R _N and changes in surface height and air pressure. Instrument type and accuracy are given in Table 1. Measurements of wind speed and all radiation components were sampled every 60 s, air temperature and humidity every 600 s, after which 60 min averages were calculated and stored on a Campbell Scientific Inc. (CSI) 23XB data logger. Surface height and air pressure were sampled and stored every 60 min. For additional information on AWS instrumentation and Northern Ice Field climate refer to Reference Mölg and HardyMölg and Hardy (2004). Mean atmospheric conditions over the measurement period are given in Table 2. Because data from our brief measurement period appear to reflect longer-term conditions well (Table 2), we feel confident that our temporally short eddy covariance dataset is representative not only of the long dry season but also of mean annual conditions.

A three-dimensional sonic anemometer (CSI, CSAT3) with a fine wire thermocouple (CSI FW05) and an ultraviolet hygrometer (KH20, CSI) were mounted at a height of 2 m on the AWS tower to obtain direct measurements of λE and H. Sensor separation between the CSAT3 and the KH20 was 0.1 m. Specifications of the eddy covariance instruments are given in Table 1. The instruments were sampled at 20 Hz using a CR5000 data logger (CSI). All high-frequency data were stored on a PCMCIA 512MB card inserted into the CR5000. Based on an estimate of the statistical uncertainty due to averaging time (Reference Lumley and PanofskyLumley and Panofsky, 1964), 30 min sampling periods were chosen to calculate all turbulence variables in this study. The total number of non-continuous (data gaps occurred because of rime and snowfall) 30 min intervals available using this sampling interval was 91 (total of 1.9 days of data). All averaging intervals were subject to a test of stationarity, using a run test as described by Reference Cullen, Steffen and BlankenCullen and others (2007), which showed that 85% of all cases fulfilled this requirement.

Before λE and H were calculated, a coordinate rotation of the velocity time series was performed in such a way that vertical and lateral wind vectors were equal to zero over each averaging interval (Reference Kaimal and FinniganKaimal and Finnigan, 1994; Reference Cullen, Steffen and BlankenCullen and others, 2007). Data were linearly detrended before turbulence statistics (including λE and H) were calculated in block averages over the chosen averaging interval. Flux loss resulting from limitations imposed by the physical size of the eddy covariance instruments was corrected using spectral transfer functions that account for path-length averaging and sensor separation (Reference MooreMoore, 1986).

Flux losses were sensitive to z/L (stability parameter), where z is height above the surface and L the Obukhov length, in stable atmospheric conditions (Fig. 1). Over almost all of the observed stability range (–0.1 < z/L < 0.25), flux losses were 2–4% and 7–14% for H and λE, respectively (Fig. 1). The correction for λE was larger than H because of the sensor separation between the KH20 and CSAT3 instruments.

Fig. 1. Fraction of flux loss estimates for H and λE determined by eddy covariance measurements on Kilimanjaro. The vertical dashed lines show the typical stability range over the measurement period (92% of all cases).

To resolve all components in Equation (1), a SEB model was used, which did not account for the effects of snowdrift sublimation and heat transfers as a result of precipitation (both negligible on Kilimanjaro). As all SEB model equations are given by Reference Mölg and HardyMölg and Hardy (2004), they will not be reproduced here, but the following provides a brief description of the modelling approach used in this study. Importantly, the SEB model was run to best reproduce H and λE, which was achieved by using all available AWS data rather than parameterized inputs for S↑ and L↓. The latter approach is favoured when performing longer-term climate sensitivity studies.

As the SEB model was run in 30 min time-steps, all hourly AWS data were initially interpolated (cubic spline) to the smaller averaging interval. To calculate the SEB, measurements of S↓, S↑ and L↓ from the AWS were directly used. The bulk method, which depends on flux–gradient relationships between the surface and one atmospheric measurement level, was used to compute H and λE. Non-dimensional stability functions were expressed in terms of the bulk Richardson number as described by Reference Wagnon, Sicart, Berthier and ChazarinWagnon and others (2003). Surface roughness lengths for momentum (z _ov) and temperature (z _ot) (moisture (z _oq) equal to z _ot) were assumed to be constant with time. The approach to obtain site-specific roughness lengths using eddy covariance measurements is described in the following section. Heat conduction (G _co) and radiation transfers (G _s) into the snow/ice were numerically solved using the thermodynamic energy equation on grid levels spaced 0.01 m apart to a depth of 3 m. The bottom temperature was held constant at 268 K. Values for the thermal conductivity of snow and ice were taken from Reference PatersonPaterson (1994). G was calculated as the flux between the two uppermost model layers (surface and 0.01 m depth). The same subsurface module is described by Reference Bintanja and van den BroekeBintanja and Van den Broeke (1995).

The SEB model was initialized using an isothermal (sub)surface temperature profile (268 K). Essentially, the solving procedure of the SEB model then followed two steps every averaging interval: (1) calculation of L↑, H, λE and G, and (2) an update of the (sub)surface temperature field. The link between the two steps was modelled T _S, which was solved iteratively so that T _S converged toward a value to allow R _N to equal the sum of H, λE and G (Reference Mölg and HardyMölg and Hardy, 2004). If T _S was greater than 273.15 K, it was reset to 273.15 K and excess energy was used for melting. In each time-step, the SEB model converted λE and M to mass fluxes of sublimation and melting, respectively. Because the spin-up time of the (sub)surface temperature field was approximately 5–7 days (error from isothermal initialization removed), the SEB model was initiated on 16 July 2005 to allow model-derived H and λE to be compared to eddy covariance data (26–29 July 2005). For clarity, all H and λE values calculated using the bulk method are referred to as ‘modelled’, while those obtained from eddy covariance data are ‘measured’.

Results and Discussion

An important characteristic of climate on Kilimanjaro is that air temperatures at the altitude of the plateau glaciers remain below 0˚C (e.g. Reference Kaser, Hardy, Mölg, Bradley and HyeraKaser and others, 2004; Reference Mölg and HardyMölg and Hardy, 2004). The mean air temperature during the eddy covariance measurement period was –5.7˚C (Table 2), which is slightly higher than the longer-term air temperature mean over 19 months. Though the average surface temperature over our study period was lower than air temperature (Table 2), implying stable stratification of the atmosphere above the glacier surface, unstable conditions prevailed during the daytime, which was also observed by Reference Mölg and HardyMölg and Hardy (2004). This diurnal cycle of stability is not unique to Kilimanjaro but has been observed over other ice and dry-snow surfaces (e.g. Reference Bintanja and van den BroekeBintanja and Van den Broeke, 1995; Reference Wagnon, Sicart, Berthier and ChazarinWagnon and others, 2003; Reference Cullen, Steffen and BlankenCullen and others, 2007). As indicated in Figure 1, very stable conditions were infrequent over the measurement period. At times when the atmosphere was unstable (45% of all cases), temperature gradients were typically small (near-neutral conditions). This is important because it reduces one large uncertainty when calculating turbulent heat fluxes using the bulk method. The choice of stability functions for momentum and heat (moisture) is not critical because they tend to be similar (converge towards unity) in near-neutral conditions. This is discussed further at the end of this section.

A greater source of uncertainty in the SEB model calculation of turbulent heat fluxes is related to surface roughness lengths. Methods of obtaining z _ov and z _ot from atmospheric measurements are described by Reference SunSun (1999), while a theoretically based model proposed by Reference AndreasAndreas (1987) infers that z _ov/z _ot is a function of the roughness Reynolds number Re _^* = u ^* z _ov v ^–1, where u ^* is friction velocity and v the kinematic viscosity of air (1.46× 10^–5m² s^–1). To calculate site-specific roughness lengths on Kilimanjaro we used u ^* and T ^* (temperature scale) from our eddy covariance measurements to solve the following equations (e.g. Reference CalancaCalanca, 2001):

where U is the mean wind speed, ΔT is the difference between air and surface temperature and Ψ_M and Ψ_H are empirical functions of the stability parameter z/L. To avoid the effects of stratification on the curvature of the mean wind speed, cases with high stable and unstable values were not used to calculate roughness lengths (–0.1 < z/L < 0.1 were excluded, or 24% of all eddy covariance data). The stability functions used were the Businger–Dyer formulations as described by Reference StullStull (1988). As snowdrift was negligible, which is known to increase z _ov (Reference Van den Broeke, van As, Reijmer and van de WalVan den Broeke and others, 2005), no threshold for u ^* was applied.

The median value for z _ov during our eddy covariance measurements was 1.7×10^–3m (Table 2), which is comparable to other values over snow and ice (Reference MorrisMorris, 1989). The value is one to two orders of magnitude larger than those reported by Reference Van den Broeke, van As, Reijmer and van de WalVan den Broeke and others (2005) for snow-covered sites in Antarctica. However, this is not too surprising, as the surface during our observations had very small penitente-type features, which we suspect were responsible for the higher z _ov values. Penitente-type structures are often observed on the horizontal surfaces of all plateau glaciers on Kilimanjaro. Estimates of z _ot yielded a median value of 1.1×10^–4 m, but the spread of z _ot was much larger than z _ov (Table 2). This larger spread in estimates of z _ot reflects our uncertainty of surface temperature, which was only measured with an accuracy of 0.5˚C (Table 1).

Though we feel confident that our brief eddy covariance measurement period is representative of longer-term conditions, we cannot rule out that a temporal variation in roughness lengths may exist. Importantly, no diurnal variation was observed, and no step change in roughness length estimates was observed after a small snowfall event, which resulted in 2.1 kgm^–2 accumulation during the late evening and early morning of 27–28 July 2005. However, there did appear to be a dependence of z _ov/z _ot on Re^* (Fig. 2), which is similar to that modelled by Reference AndreasAndreas (1987). The model proposed by Reference AndreasAndreas (1987) predicts a decrease in z _ot (and z _oq) with increasing Re^*. This is observed in the Kilimanjaro data, but estimates of z _ot at higher Re^* values appear to be larger than modelled data. This results in the majority of binned data of z _ov/z _ot being above the solid line (Fig. 2) at times when Re^* is greater than 1. Because of measurement uncertainties, in particular the estimate of surface temperature, we cannot go as far as to suggest our results validate the Reference AndreasAndreas (1987) model, but we do believe that it may be appropriate for use over snow/ice surfaces such as those found on Kilimanjaro. The alternative, which is perhaps equally valid given our uncertainties (Fig. 2), is to assume that the magnitude of z _ov is equal to z _ot and z _oq.

Fig. 2. The ratio z _ot/z _ov as a function of the roughness Reynolds number (Re_*) using geometrically bin-averaged values (grey squares). The median values in each calculated bin are also shown (black squares). Each bin represents about 10% of available data after being sorted by magnitude (Re_*). The error bars are one standard deviation in each bin category and are only shown for median values (same for geometric mean). The solid line is the relationship predicted using the Reference AndreasAndreas (1987) model.

To assess the accuracy of the bulk method used in our SEB model, the surface roughness lengths given in Table 2 were used to calculate turbulent heat fluxes, and then compared to eddy covariance measurements (Fig. 3). For both H and λE, modelled (bulk method) values agree well with the (eddy covariance) measurements (r = 0.82 for H and 0.77 for λE), with root-mean-square differences (RMSD) of 3.1 and 4.8Wm^–2 for H and λE, respectively. Importantly, the sign of H was well reproduced by the SEB model, with the diurnal variation in stability captured by the bulk method. Differences in the sign of H always occurred during transition periods from one stability regime to the other (13% of all cases), but no bias towards any particular stratification occurred (not shown). Figure 3 shows that the magnitude of modelled (bulk method) H in both stable and unstable conditions was not as large as measured (eddy covariance). Because of a cancelling effect imposed by changes in stability, averages of both modelled and measured H were small.

Fig. 3. Comparison of modelled (bulk method) and measured (eddy covariance) H and λE. Stable (unstable) cases as defined by the eddy covariance data are shown as black (grey) squares. Correlation coefficients (r) and RMSDs are also displayed.

Importantly, values of λE obtained from eddy covariance measurements confirm that mass loss due to sublimation was continuous in both stable and unstable conditions (Fig. 3). Though the scatter between modelled (bulk method) and measured (eddy covariance) λE was greater than for H, the diurnal amplitude of the former is reproduced more accurately. Total measured (eddy covariance) sublimation was 2.74 kgm^–2 over the brief observation period (1.9 days), which was 14% greater than modelled (bulk method) sublimation. Because of time gaps in the eddy covariance data, it was difficult to use surface height data to estimate ablation over the same interval, which would have provided a measure of total mass loss and allowed us to assess the relative contribution of sublimation to that loss. Though changes in surface height were available at the same temporal resolution as measured (eddy covariance) and modelled (bulk method) λE, we are of the belief that surface height data can only be used to estimate daily totals of ablation. The accuracy of the sonic ranging sensor (CSI SR50) is 0.01 m (Table 1), but the resolution is 0.1 mm, which is about 10% of daily mass loss due to sublimation.

To overcome the temporal limitation of our eddy covariance dataset, we extended our analysis to the entire second week of our SEB model run (25–31 July 2005) to briefly assess the contribution of modelled sublimation (bulk method) to total mass loss. This also gave us the opportunity to more accurately assess the sensitivity of modelled sublimation to changes in data and model parameter inputs. A useful first step in assessing the performance of the SEB model over a longer time interval was to compare modelled with measured surface temperatures (Fig. 4). Though modelled surface temperatures reproduced the measured data well (r = 0.91), the RMSD between the two temperatures was larger than expected (RMSD = 1.8 K). The major cause for this difference is that modelled surface temperatures reached melting point during the daytime, on average for 3.2 hours d^–1, while measured temperatures always remained below 273.15 K. This discrepancy is important, as modelled surface temperature determines all surface fluxes except S _N, G _s and L↓, and affects the proportion of energy available for melt. Data obtained from the instrument measuring surface temperature (Table 1) are sensitive to how measurements are linearized, such as the choice of offset. This issue must be partly responsible for the difference between modelled and measured surface temperatures, as observations in the field indicated that the surface reached melting point during the daytime.

Fig. 4. Comparison of modelled (black line) and measured (grey line) surface temperatures during the week eddy covariance measurements were made (25–31 July 2005).

Though we have shown that modelled sublimation (bulk method) was slightly underestimated compared to eddy covariance data, and modelled surface temperatures favoured melting, Figure 5a shows that λE was still the most important non-radiative energy-balance term during the SEB model testing period (25–31 July 2005). S _N clearly provided the energy available for sublimation and the other non-radiative processes (Fig. 5a). The importance of S _N and its variability as controlled by precipitation (changes in albedo) over 19 months is described in detail by Reference Mölg and HardyMölg and Hardy (2004) and will not be further discussed here. Of more relevance to this study is that sublimation accounted for 92% of observed ablation (Fig. 5b). Although energy available for melt (M) was large enough to result in ablation rates of 3.9 kgm^–2 δ^–1, this clearly did not occur, as observed mass loss was only 7.8 kgm^–2 week^–1 (Fig. 5b). The only plausible explanation for this is that meltwater refroze on the mostly flat surfaces of the plateau glaciers.

Fig. 5. Energy (a) and mass (b) fluxes calculated using the SEB model and in situ measurements for the period 25–31 July 2005. Accumulation and ablation in (b) were determined from surface height measurements, while sublimation (kgm^–2 week^–1) was calculated from modelled λE (bulk method). Error bars in (b) reflect surface height instrument resolution (accumulation and ablation), while uncertainty in mass loss by sublimation is defined in Table 3.

The error bars on total mass loss due to sublimation in Figure 5b were determined by performing a series of SEB model sensitivity runs (Table 3) following an approach described by Reference Greuell and SmeetsGreuell and Smeets (2001). A net uncertainty in modelled sublimation (bulk method) was determined by increasing or decreasing each individual input variable relative to the reference SEB model run according to its uncertainty (Table 3). The same approach was used for calculated roughness lengths, but the (maximum) uncertainty was assessed to be one order of magnitude. The effect of stability functions on the magnitude of modelled sublimation was assessed by entirely switching them off during one sensitivity run. For each input variable and model parameter an absolute average was obtained, and the total uncertainty (2.2 kgm^–2 week^–1) was calculated as the square root of the sum of the squares of the individual uncertainties (Table 3).

The total uncertainty in modelled sublimation was mostly influenced by the one-order-of-magnitude change in z _ov and z _ot (Table 3). By obtaining in situ values of both these parameters using eddy covariance measurements, we believe we have greatly reduced the uncertainty of obtaining reliable estimates of modelled sublimation (bulk method) on Kilimanjaro. This sensitivity analysis also confirms our previous claim that the choice of stability functions is not critical on Kilimanjaro. Uncertainties in input data also do not greatly affect total mass loss due to sublimation. Taken as a whole, we are confident that sublimation was well reproduced by our SEB model, and that it can clearly account for almost all of the observed surface mass loss on the horizontal surfaces of plateau glaciers on Kilimanjaro.

Conclusions

The uncertainties of obtaining reliable estimates of sublimation using a SEB model on Kilimanjaro have been assessed using eddy covariance data. Although the eddy covariance measurement period was small because of the logistical difficulties of spending time on the top of Kilimanjaro, this dataset provided us with a unique opportunity to evaluate the response of our SEB model to changes in input data and model parameters.

The magnitude of modelled (bulk method) turbulent heat fluxes was affected more by changes to z _ov and z _ot than any other data input or model parameter. By obtaining eddy-covariance-derived values of z _ov and z _ot, and using them in our SEB model to calculate turbulent heat fluxes, we have established more robust estimates of total modelled sublimation (bulk method). This study confirms previous assertions (Reference Mölg and HardyMölg and Hardy, 2004) that sublimation can account for the majority of mass loss on the horizontal surfaces of plateau glaciers on Kilimanjaro. Such a finding clearly illustrates the importance of having a thorough knowledge of the full SEB over a glacier surface. Without such knowledge it could be easy to have the false impression that other processes, such as air-temperature-induced melt, could be more important in contributing to total mass loss.