Slope glaciers on Kilimanjaro (ca. 5000–6000 m MSL) reached their most recent maximum extent in the late nineteenth century (L19) and have receded since then. This study quantifies the climate signal behind the recession of Kersten Glacier, which generates information on climate change in the tropical midtroposphere between L19 and present. Multiyear meteorological measurements at 5873 m MSL serve to force and verify a spatially distributed model of the glacier’s mass balance (the most direct link between glacier behavior and atmospheric forcing). At present the glacier is losing mass (522 ± 105 kg m−2 yr−1), terminates at 5100 m, and the interannual variability of mass and energy budgets largely reflects variability in atmospheric moisture. Backward modeling of the L19 steady-state glacier extent (down to 4500 m) reveals higher precipitation (+160 to +240 mm yr−1), higher air humidity, and increased fractional cloud cover in L19 but no significant changes in local air temperature, air pressure, and wind speed. The atmosphere in the simulated L19 climate transfers more energy to the glacier surface through atmospheric longwave radiation and turbulent heat—but this is almost entirely balanced by the decrease in absorbed solar radiation (due to both increased cloudiness and higher surface albedo). Thus, the energy-driven mass loss per unit area (sublimation plus meltwater runoff) was not appreciably different from today. Higher L19 precipitation rates therefore dominated the mass budget and produced a larger glacier extent in the past.
Unlike many climate proxy data that are linked to biotic processes, the behavior of glaciers follows physical laws exclusively. The immediate impact of weather and climate governs the glacier volume, through mass and energy exchanges between glacier surface and atmosphere, while the change in glacier extent is the delayed response to these exchanges through ice flow dynamics. A correct understanding of the glacier–climate interaction therefore allows one to infer climate change from recorded glacier changes, which date back to preinstrumental periods in many parts of the earth.
Glaciers in the tropics have a particular relevance in light of this possibility. At present they exist in the tropical sector of the American Andes, equatorial East Africa, and in Indonesia (Kaser and Osmaston 2002). They occur exclusively at high elevation where precipitation falls as snow, which is the precondition for glacier formation. This high-altitude location (mostly >5000 m MSL) relates to the vitally important question of climate change in the midtroposphere; a difficult issue since free-atmosphere observations have basic problems for trend analyses (e.g., Karl et al. 2006) and observations and climate model results may disagree on whether surface trends are amplified in the free troposphere (see Pepin and Lundquist 2008, for a discussion). Intensified field research on tropical glaciers during the past decade in South America (e.g., Kaser et al. 2003; Favier et al. 2004) and Africa (see below) can therefore contribute to understanding climate change at very high altitude.
Compared to other tropical glaciers the ice fields on Kilimanjaro (Tanzania, East Africa) are related to several peculiar circumstances, which require particular care for the interpretation of glacier–climate interactions (see section 2 and Kaser et al. 2004). Nonetheless, situated on top of a massive freestanding mountain (5895 m MSL), these glaciers have the advantage of being exposed freely to the midtropospheric flow, which minimizes local climate forcings that may dominate glacier behavior in complex terrain. We have been studying these glaciers intensively during the last few years in terms of both microscale processes that control the glacier fluctuations (Mölg et al. 2003b, 2008, 2009b; Kaser et al. 2004; Mölg and Hardy 2004; Cullen et al. 2006, 2007a) and meso- and large-scale climate dynamics that have set the boundary conditions for glacier shrinkage on Kilimanjaro during the past ∼120 yr (Mölg et al. 2006, 2009a, hereafter MCGC).
This paper builds on all of these studies to infer quantitatively the change in climate between the late nineteenth century (L19) and the present from respective glacier extents (Fig. 1). This provides unique information on climate change for a place and time that are not covered by instrumental records. Section 2 gives the details about the glaciers on Kilimanjaro and their suitability for climate reconstruction, while section 3 presents the methods (measurements and modeling). Section 4 extends insight into local climate forcing of glacier volume at present, and section 5 explores the required change in climate forcing to enable the L19 glacier extent. A brief discussion of associated changes in large-scale climate drivers (section 6) precedes the conclusions (section 7).
2. Recent changes to glaciers on Kilimanjaro
The recent behavior of glaciers on Kilimanjaro, which are confined to the central peak Kibo, can be understood only when distinguishing between two glacier systems (Geilinger 1936; Kaser et al. 2004). First, the tabular-shaped glaciers on the rather flat summit plateau show pronounced, almost vertical ice cliffs at their margins (Fig. 1), where snowfall cannot accumulate. Thus, these cliffs lose mass rather constantly because of absorbed energy from radiation (Mölg et al. 2003b), which leads to a near-constant decrease in the plateau glaciers’ surface area (the variable depicted in Fig. 1; Cullen et al. 2006). Extracting climate change details from areal changes of plateau glaciers is therefore tricky, particularly since the mechanical sensitivity of ice cliffs to climate fluctuations is not yet understood (Kaser et al. 2004). Second, slope glaciers on Kibo below the summit plateau have shown a variable decrease in area during the past century (Cullen et al. 2006) because they convert changes of mass gain and loss—the sum of which is the mass balance of a glacier—into changes of glacier extent in the typical fashion (see the introduction). Thus, we focus on the slope glaciers for the reconstruction of past climate. Geothermal heat due to Kilimanjaro’s volcanic nature may impact ice loss of plateau glaciers to a very limited extent (e.g., Fig. 6 in Kaser et al. 2004), but this is not at all evident for the slope glaciers. Hence, the recent historical glacier changes on Kibo (Fig. 1) are actually forced by climate (Downie and Wilkinson 1972; Kaser et al. 2004).
Kersten Glacier is the largest remaining slope glacier (Fig. 1) and is located on the southern flank of the mountain. Several fresh terminal moraines were recorded for the most recent glaciation of Kibo (also below Kersten Glacier; Fig. 1), which are large enough to indicate a period of standstill with equilibrium between mass loss and gain on the glaciers (Osmaston 1989). The implications of this situation for the modeling approach are detailed in section 3e. When the first explorers reached the glacierized area of Kibo between 1887 and 1889, the glaciers were still in contact with these moraines or had retreated just a short distance (Meyer 1900). Similar evidence comes from the other two glacierized mountains of East Africa, Rwenzori and Mount Kenya, where glaciers commenced to shrink from a maximum extent during the 1880s and 1890s (Hastenrath 1984; Kaser 1999). Prominent rises of lake levels in East Africa, which started between 1820 and 1850, also came to an end between 1880 and 1900 (e.g., Hastenrath and Kutzbach 1983; Verschuren et al. 2000; Nicholson and Yin 2001), so the most recent maximum extent of glaciers on Kibo and the underlying climate conditions are associated with the late nineteenth century.
Climate in the Kilimanjaro region follows a bimodal distribution of rainfall (Coutts 1969), which reflects the wet seasons from March to May (MAM) and October to December (OND) in equatorial East Africa (e.g., Camberlin and Okoola 2003). A long dry season prevails from June to September (JJAS), and a moderately dry season in January and February (JF). So far, the application of our high-altitude measurements to physically based mass balance models for point locations on the northern ice field (station 1 in Fig. 1; Mölg and Hardy 2004) and the Kersten Glacier (station 3 in Fig. 1; Mölg et al. 2008) demonstrated that mass balance fluctuations primarily depend on precipitation variability, consistent with our early observations (e.g., Hardy 2003). Here these preliminary studies are extended to an entire glacier and different extents, but the point results suggest that the climate of the L19 glacier extent most probably featured more precipitation, which agrees with the history of lake levels (above) that are also controlled by precipitation (Hastenrath and Kutzbach 1983; Nicholson and Yin 2001).
a. Mass balance model
The mass balance of a glacier is the most direct link between atmospheric forcing and glacier response. The specific mass balance b (kg m−2) is the sum of accumulation c (mass gain processes) and ablation a (mass loss processes) over a specified unit of time. The employed mass balance model is described completely in Mölg et al. (2008) who verified the model at a point location (at AWS3, see Fig. 1). We utilize, however, the spatially distributed version of the model, which is expressed in the descriptions below as well as in section 3b. Similar models [known as two-dimensional (2D) models] were applied to extratropical mountain glaciers (e.g., Klok and Oerlemans 2002; Reijmer and Hock 2008) but not to an African glacier before.
Our model computes b at an hourly step as the sum of solid precipitation (csp), melt (msfc), sublimation (ssfc) and frost or dew deposition (cdep) at the surface, and englacial accumulation (cen) by refreezing of meltwater in snow or at the snow–ice interface. The first governing equation therefore reads
QL is the turbulent latent heat flux of water vapor, and QM signifies the latent energy flux of melting; they are connected to b through the latent heats of sublimation/deposition (LS) and melt (LM). Mass fluxes that imply a mass loss for the glacier column have a negative sign. The last three right-hand terms can be understood from the surface energy balance (SEB) of the glacier (in W m−2), which leads to the second governing equation:
where S↓ is incoming shortwave radiation, α is the spectrally integrated surface albedo, L↓ and L↑ are incoming and outgoing longwave radiation, QS is the turbulent sensible heat flux, QC is the conductive heat flux in the subsurface, and QPS is the energy flux from shortwave radiation penetrating to the subsurface. An energy flux has a positive (negative) sign when it induces a heat gain (sink) at the surface. The sum of these fluxes yields a resulting flux F, which represents QM only if surface temperature (Tsfc) is 273.15 K. The model’s SEB module interacts through Tsfc with the subsurface module, which prognoses the englacial temperature (Ten) on a numerical grid (e.g., Fig. 4 in Mölg et al. 2008).
The radiation scheme of the model calculates S↓ and L↓ from an effective cloud cover fraction (neff, see section 3b), air temperature (T), humidity, and pressure. Mölg et al. (2009b) optimized the scheme for Kibo and give a full description. The S↓ is corrected for slope and aspect of the model grid cell (e.g., Mölg et al. 2003a), so S↓ equals the global radiation (G) merely for a horizontal surface. Shaded grid cells only receive the diffuse component of S↓ (Mölg et al. 2009b). The albedo scheme determines α as a function of snow depth and time difference to the last snowfall event (Oerlemans and Knap 1998). Since snow depth results directly from the first three right-hand terms in Eq. (1) (and indirectly from cen through snow density), it is influenced by all SEB terms in Eq. (2) and, therefore, by all atmospheric variables that force the SEB (see section 3b). The control parameters in the albedo scheme derive from an optimization through measured shortwave radiation and snow depth (Mölg et al. 2008), which spans 824 days (section 4c). The fourth radiation term L↑ is obtained from the Stefan–Boltzmann law, assuming a surface emissivity (ɛsfc) equal to unity.
The turbulent heat flux computation follows the Monin–Obukhov similarity theory, using analytical expressions of the “bulk” method (Mölg and Hardy 2004). Surface-air gradients of humidity (for QL) and temperature (for QS), as well as air pressure (p), horizontal wind speed (V), surface roughness lengths for momentum (z0m) and heat (z0h; Cullen et al. 2007a), and stability of the surface boundary layer control the magnitude and direction of turbulent fluxes.
In the subsurface, Ten is solved from the thermodynamic energy equation on 14 layers in the vertical, including heat release by QPS and the englacial refreezing process (Mölg et al. 2008). Subsurface layers are centered at 0.09-, 0.18-, 0.3-, 0.4-, 0.5-, 0.6-, 0.8-, 1.0-, 1.4-, 1.8-, 2.2-, 2.5-, and 3-m depth. The resulting Ten profile serves to calculate QC. From the refreezing meltwater, a fraction fsui (0.3) refreezes at the snow–ice interface rather than within the snow (e.g., Cheng et al. 2006) and forms superimposed ice. The factor Cref (in units per time step) determines if the refreezing process is allowed to reach completion (Illangasekare et al. 1990). Here, Cref varies linearly with terrain slope (from 0 to 1 h−1 between 90° and 0° inclination), which is a simple representation of runoff time scales (e.g., Reijmer and Hock 2008). Meltwater that does not refreeze is assumed to run off.
Surface temperature governs a number of surface processes (L↑, QS, QL, QC) and also forces the thermodynamic energy equation at the upper boundary. An iterative Tsfc scheme was employed earlier for point locations (Mölg and Hardy 2004; Mölg et al. 2008) but tends to overestimate melting (Mölg et al. 2008). We therefore follow the scheme of Klok and Oerlemans (2002) for the 2D modeling. It converts the net flux F into a temperature change over a defined surface layer thickness (dsfl), which approximates the depth where the diurnal temperature amplitude diminishes to 5% of the amplitude in Tsfc (Klok and Oerlemans 2002). In the model runs of Mölg et al. (2008) this occurs at roughly half a meter below the surface, so dsfl = 0.5 m. The performance of the Tsfc scheme is validated explicitly in section 4a. An extensive sensitivity analysis investigates the impact of all key parameters involved in the computation of b (sections 4c and 4d).
b. Model input requirements
The mass balance model requires two types of inputs: (i) a digital terrain model (DTM) that provides altitude, aspect, slope, and sky-view factor for each grid cell (lower boundary condition); and (ii) meteorological data that drive the model over a certain time period (atmospheric forcing).
The DTM is constructed from data of the Shuttle Radar Topography Mission (e.g., Rabus et al. 2003). Figure 2 depicts the 2003 extent of Kersten Glacier in the DTM at a horizontal resolution of 50 m, which meets the real situation very well (Fig. 1). The glacier is orientated rather uniformly toward south-southeast, flat in its upper part and steep in the tongue area (Fig. 2). The reference model run (REF) uses a coarser resolution of 100 m, which, for the simple shape of Kersten Glacier, does not impact results to any large degree (see sensitivity analysis in section 4c).
The atmospheric forcing comprises six variables at hourly resolution: T, relative humidity (RH), neff, p, V, and precipitation (P). The mass balance model is run over a period of 1079 days, for which the forcing data are available from automatic weather stations (AWSs) on Kibo (Fig. 1), as detailed in the next two sections. The variable neff can be calculated from measurements of L↓ (all day) or G (for daytime), a method developed in Mölg et al. (2009b). The neff series for REF rely on G for the daytime window (1000–1600 LT) and L↓ for the remaining day. To extrapolate the point measurements to the entire glacier, a T gradient of −6.5 K km−1 (Osmaston 1989) at constant RH is prescribed, which produces T and water vapor pressure (e) for each grid cell. A vertical gradient in V of 1.3 m s−1 km−1 is assumed, as evident in the data analyzed by MCGC. A negative gradient in P on Kilimanjaro was measured by Røhr and Killingtveit (2003) between ∼2200 and ∼4600 m MSL but is unknown for the glacierized regions. Thus, for REF we assume the same P at all altitudes of Kersten Glacier, while a threshold of 2.5°C (Tpsl) separates solid from liquid P (Mölg et al. 2008). Model sensitivities to changing gradients are explored in section 4d. Table 1 summarizes the links between forcing data and mass balance (and related) terms.
c. Measurements and observations
Three AWSs currently collect meteorological data in the summit region of Kibo (Fig. 1): AWS1 since February 2000 (operated by Massachusetts) and AWS2 and AWS3 since February 2005 (operated by Innsbruck). AWS3 on Kersten Glacier at 5873 m MSL provides the data to drive and validate the physically based mass balance model, as demonstrated by Mölg et al. (2008) for the year 2005. Relevant measurements comprise incoming and outgoing radiative fluxes (relative to a horizontal surface) in the shortwave spectrum 305–2800 nm and longwave spectrum 5000–50 000 nm with a Kipp and Zonen CNR1 net radiometer (accuracy ±10%), T and RH with a Rotronic MP100A (protected with radiation shields) at initially 1.75 m above the ground (±0.2°C, ±2% units), V and wind direction with a RM Young 05103–5 at initially 1.85-m height (±0.3 m s−1, ±5°), p with a Setra 278 (±0.2 hPa), and distance to the surface with a Campbell SR50 sonic ranging sensor (±0.4%), which gives continuous records of (i) surface accumulation (decreasing distance) and ablation (increasing distance), (ii) the above instruments’ heights, and (iii) snow depth estimates. All measurements are sampled in one-minute intervals and then stored as half-hourly means on a Campbell CR-23X datalogger. For more details about the AWS measurements on Kibo we refer to our preceding studies cited in section 1 (e.g., Cullen et al. 2007a). Similar setups of AWSs have proven to provide reliable data on glaciers around the world (e.g., Oerlemans and Knap 1998; Van den Broeke et al. 2005; Cullen et al. 2007b), and for this paper we can rely on a high-quality dataset from AWS3 (see next section)—which is exceptional in view of the harsh weather conditions and working environment at almost 6000 m MSL. The net b at AWS3 from 9 February 2005 to 23 January 2008 (1079 days), determined by the SR50 sensor and manual measurements in the field, amounts to −654 ± 21 kg m−2 (which equals mm water equivalent as used in glaciology). It consists of 0.84 m (actual height) ice loss with an assumed ice density of 850–900 kg m−3 (which causes the above uncertainty) and a net gain of 0.22 m snow with a measured snow density (ρs) of 367 kg m−3.
Initially, measurements of b were planned along the centerline of Kersten Glacier between its tongue and AWS3 by traditional ablation stakes, which ameliorates validation of the mass balance model for locations other than the station (e.g., Reijmer and Hock 2008). The bittersweet reality is that we failed to establish this network so far because of the difficult accessibility of Kersten Glacier’s lower parts (accessible only from above) and the unavoidable physical exhaustion of the participants in field expeditions. A few measurements of snow depth and surface height change at 5690 and 5760 m MSL are nonetheless available for model validation (see section 4a). Further, extensive terrestrial photography as well as aerial photography in 2005 (Cullen et al. 2006) document that Kersten Glacier has thinned in the uppermost parts (emergence of rock), which strongly suggests that this area is an ablation zone (b < 0 over multiyear period) and agrees with the point-measured b above. The existence of an ablation zone at the upper end of slope glaciers on Kilimanjaro was already noted by early researchers and attributed to more effective solar insolation and/or lower precipitation than on midslopes (see Osmaston 1989 for a concise review). A decrease or stagnation in annual b toward the upper end was also measured during the field program on Lewis Glacier of nearby Mount Kenya (Fig. 3). Thus, the characteristic shape of the vertical b profile and its interannual variability, as observed by Hastenrath (1984) for a similar African glacier (Fig. 3), provide one additional chance to assess the spatial performance of the mass balance model.
d. Quality control of measurements
After installation, AWS3 was revisited in July 2005, February and July 2006, January and October 2007, and January 2008 for routine work (data retrieval, platform checks to ensure that instruments are level, and manual inspection of instrument heights to support measurements obtained from the SR50). Gap-free records exist from all sensors, except for the SR50, which failed from days 825–883 and 919–1079. The P as model input for these days is taken from the SR50 at AWS1, where P hardly differs from that at AWS3 (Mölg et al. 2008). Postprocessing of all data aims to detect measurement errors (see the appendix). Associated corrections only apply to a minority of data, but we deem the descriptions in the appendix nonetheless appropriate to demonstrate data were examined thoroughly.
Figure 4 presents daily means of five selected variables, which are used to force the mass balance model (Table 1). The T shows very little fluctuation, and even day-to-day variability remains almost entirely within the amplitude of the mean diurnal cycle. Monthly means vary much less than the diurnal cycle (“thermal homogeneity”: Kaser and Osmaston 2002). The P, e, and cloudiness exhibit the characteristic annual cycle of wet and dry seasons discussed earlier, while V tends to minima in the JJAS dry season. Here T, e, and V are corrected to screen level (2 m) in Fig. 4 from their actual measurement height by log-linear functions (Oerlemans and Klok 2002). This has a negligibly small impact, but is done for clarity and for the employment of the radiation scheme (Mölg et al. 2009b). Most noticeable for the period depicted in Fig. 4 are unusually dry conditions in OND 2005 and particularly strong precipitation in OND 2006 (MCGC). This suggests that measurements capture an important aspect of interannual climate variability at our site, since wet seasons in OND are generally more variable than in MAM (e.g., Nicholson 2000). Further measurements from AWS3 are presented in the model validation (section 4a).
e. Backward-in-time modeling
For a glacier with a rather constant area–altitude distribution (like Kersten) it can be assumed that the mass balance along the glacier centerline at steady state, between glacier top and terminal moraine, is zero (e.g., Kruss 1983; Yamaguchi et al. 2008). Hence, the necessity for a glacier flow model (see introduction) cancels for such a condition. This means the atmospheric forcing can be changed until the glacier reaches the L19 steady-state extent as indicated by the terminal moraine (Fig. 1). For clarity these changes in forcing are described after the present-day results (section 4), which provides a context. Figure 1 portrays the model grid points of the centerline, over which a steady-state profile must be valid for L19 conditions. Mean b over the chosen grid points mimics mean b area integrated over Kersten Glacier very precisely, which is demonstrated in conjunction with the results (section 4d).
The approach taken still evokes the question of how representative our time slice (3 yr) is of the present climate over Kilimanjaro. Figure 5 compares climate variables between our measurement period and the most recent 30-yr period from the National Centers for Environmental Prediction–National Center for Atmospheric Research (NCEP–NCAR) reanalysis (REA) data (Kalnay et al. 1996). The only systematic deviation arises in the T plot. This should not be interpreted as recent warming but as multiyear variability, which appears throughout the past 60 yr in the REA for our location, with 2005–07 being a warm episode of this variability (cf. Fig. 4 in Cullen et al. 2006). In addition, such small changes in T at subfreezing conditions do not affect b on Kibo’s glaciers importantly (Mölg et al. 2008). Hence, our 3-yr measurement period does not represent anomalous conditions within the twentieth-century climate.
4. Results and discussion: Present
a. Model validation (reference run)
Figure 6 compares measured and modeled incoming radiation at the AWS3 site, which yields very high correlation coefficients (r). Moreover, the root-mean-square difference (RMSD) is clearly smaller than the nominal accuracy of the CNR1 (see section 3c) for both components. The radiation scheme was optimized by Mölg et al. (2009b) for the first 700 days of data from AWS3, so only the last 379 days of our time slice are considered for the validation in Fig. 6.
Figure 7 turns to the mass fluxes over the entire 1079 days and demonstrates that the model reproduces surface height change and the associated net b well. An overestimation by the model is apparent after mid-July 2006 (roughly one month), which results from an overestimation of snow depth (Fig. 7, middle). By contrast, an underestimation of snow depth occurs in February and March 2007. Snow depth derived from the SR50, however, does not allow explicit detection of superimposed ice (SI), so this estimate is the sum of snow and SI depth. The model indicates SI formation in the same time span (which agrees with the observation in October 2007), which helps account for the discrepancy between the curves. Generally, the model differs most notably from measurements in the three main phases of decaying snow depth. Snow density is a pivotal variable for the decay speed (light snow ablates faster), but the few ρs measurements so far make validation of the simulated ρs difficult. Still, the simulated maximum (497 kg m−3) agrees with measurements of old snow in the summit region (300–500 kg m−3). The bottom panel of Fig. 7 shows the validation of Tsfc. “Measured” Tsfc is determined from longwave radiation data (e.g., Van den Broeke et al. 2005). A typical measurement uncertainty of 10 W m−2 leads to a RMSD of 2.3 K from the Tsfc measurement in Fig. 7. Hence, the 1.7 K between model and measurement are acceptable, especially since the model captures the variability well. The prominent minima of Tsfc in June and July were analyzed by Mölg et al. (2008) for 2005 and result from a period dominated by negative net radiation. In summary, Fig. 7 suggests a very good model performance, even if model values reflect a gridcell average and not a true point measurement like AWS3.
As mentioned earlier, there is a lack of validation data for the lower regions on Kersten Glacier. Figure 8 depicts the ones that exist and does not indicate a basic model deficiency. Surface height change and two different snow conditions of downslope increase (July 2006) and constant depth (January 2008) are reasonably captured. Interannual variability in modeled b profiles also seems realistic (shown below). To counteract the lack of validation options for lower regions, the mass balance of Kersten Glacier is eventually discussed from an ensemble of model runs (section 4d).
b. Causal factors of mass balance variability
The modeled vertical profiles of annual b are depicted in Fig. 9 (top left). They do show the maximum slightly below the upper end, and exhibit a small (big) vertical variation under dry (wet) conditions as in the year 2005 (2006). This is tied to a larger mass turnover in wet conditions that is mainly due to higher precipitation over all altitudes and more melting at the lower elevations (Fig. 9, left). Variability in b profiles therefore follows a similar pattern like measured on Mount Kenya (Fig. 3). The vertical variation of b, however, is expectedly larger for the lower-lying Lewis Glacier, where the atmosphere is moister and warmer (Hastenrath 1984). Sublimation occurs prominently at all altitudes on Kersten Glacier, while deposition as well as englacial accumulation (apart from the highest regions) are small (Fig. 9, left).
The S↓ and L↓ are the only ground-independent energy sources for a glacier surface. All other SEB components are only partly independent from the glacier surface temperature and/or the snow conditions (e.g., QS, QL) or entirely passive responses (e.g., L↑), as visible in Table 1. The shape of the b profiles in the upper glacier region is linked to the pattern of S↓ (Fig. 9, right), which indicates the highest efficiency of solar insolation (relative to the G curve) at the top where terrain slope is small (Fig. 2) and the sky-view factor is high (little shading by terrain). Thus, b maximizes below the top where the first minimum in energy input by S↓ (at ∼5700 m MSL) occurs. This favors a number of self-amplifying feedback processes, which involve a relatively long snow cover duration (292 days per year at 5750 m, 235 at 5859 m), a high α and small net shortwave radiation, Snet = S↓ (1 − α), and minimum Tsfc (expressed by L↑) (Fig. 9, right). In the lower glacier region, increased L↓ and QS, as well as a reduction of the heat sink by QL, overcompensate the effects from the second S↓ minimum, which develops toward the tongue. This leads to more residual energy and, thus, more frequent melting (892 h yr−1 at 5028 m, 340 at 5750 m), and drives the large mass losses at the lower elevations. Table 2 relates the total (area integrated) mean net b of Kersten Glacier (btot) to atmospheric forcings, and to total mass and energy balance components. It is evident that only a few watts per square meter make the difference between rather different ablation conditions (cf. Fig. 9), so the interaction between slope glaciers and climate on Kilimanjaro is a very sensitive system. The total ground heat flux (Table 2) approximates zero also at all elevations, and is thus not included in Fig. 9. In a dry year (2005) sublimation can be greater than melting even glacierwide (Table 2), not just in the upper parts (Fig. 9).
Of further interest is that a number of profiles during 2006 in Fig. 9 are shifted along the x axes relative to the year 2005, but 2007 occasionally breaks this pattern in the upper glacier region (e.g., msfc, ssfc, α). This can be interpreted as a “memory” effect from the wet year 2006, which saw snowfall mostly in OND (Fig. 4). The year 2007 therefore began with a well-pronounced snow cover that lasted late into the year (Fig. 7), leading to a high mean albedo and low mean Tsfc but also greater meltwater retention in the snowpack (Fig. 9). Lower Tsfc, however, may increase heating by QS, and the QS increase of 2007 was higher in the hours from 1400 to 1600 LT (the most likely daytime with Tsfc = 0°C) than the L↓ increase of 2006, which mostly affected nocturnal hours. Thus, surface melt in 2007 was slightly greater in the upper regions than in 2006 (Fig. 9, bottom left). This underlines the need to look into diurnal cycles of SEB components on tropical glaciers (Mölg and Hardy 2004), where dominant phases of ablation may occur year-round on short subdaily periods (e.g., Mölg et al. 2009b), not over entire seasons like on extratropical mountains (e.g., Klok and Oerlemans 2002).
The importance of snowfall in OND 2006 and subsequent snow cover duration support the main message from Fig. 9 and Table 2: Interannual variability in btot and the b profiles, as well as in the driving energy and mass fluxes behind, displays moisture variability over Kilimanjaro’s summit. Dry (wet) years typically feature strong (weak) sublimation, low (high) accumulation, and low (high) meltwater fluxes. The only adverse effect on b during wet conditions is the increase in surface melt (due to less negative QL and reduced sublimation), as found on other tropical glaciers (Wagnon et al. 1999). Its importance on Kersten Glacier, however, is outpaced by mass gains from more precipitation and refreezing and limited by less absorption of solar energy. This implies that the interannual variability in Snet (and, thus, albedo) is most important to the mass balance from a quantitative view of energy sources (Table 2), which once more explains the strong precipitation dependency of Kibo’s glaciers (Hardy 2003; Mölg and Hardy 2004; Mölg et al. 2008).
c. Sensitivity to well-known model parameters
The reference run can rely on a number of parameters that were carefully optimized for Kersten Glacier (Mölg et al. 2008, 2009b). Table 3 still investigates the sensitivity of results to such parameters to learn more about the model behavior. Parameters that control Snet (first seven lines) represent a relatively sensitive model component. This is expected since Snet is the most variable energy source (see discussion above), and it also appears in other 2D models (e.g., Reijmer and Hock 2008). Parameters are varied in only one direction if the preceding optimization yielded almost identical results for a certain range. For instance, setting the depth scale (d*) to 30–36 cm or the time scale (t*) to 5.4–8.0 days hardly affected the performance of the albedo scheme in optimization mode. These scales, and albedos of fresh and old (aged) snow, differ from Mölg et al. (2008) who did the optimization solely for the dry year 2005. The updated values (αfrs = 0.89, αols = 0.51, d* = 36 cm, t* = 5.4 days) are optimized over dry and wet years, so are certainly more robust for longer-term studies. The last four lines in Table 3 target the grid structure and subsurface module. Doing the more expensive runs for 50-m resolution (Fig. 2) or for 30 subsurface layers (Mölg et al. 2008) affects modeled btot marginally. This also applies to the constant bottom temperature (Tbot) at 3-m depth and the parameter that controls QPS in ice (ζi), both of which were optimized in Mölg et al. (2008). To sum up, individual sensitivities in Table 3 remain within 10% of btot, which is smaller than the estimated uncertainty (next section).
d. Sensitivity to less known model parameters
Model parameters from the literature or estimated ones (e.g., dsfl) are varied in Table 4, which leads to 17 runs from which the uncertainty in modeled btot is assessed. From the parameters that govern the model physics (1–12), the density of solid P (ρp) is the most sensitive component. Density of fresh snow is already high on tropical glaciers, 150–300 kg m−3 (Sicart et al. 2002), and on Kibo significant P events contain graupel, which can increase density even further (MCGC). Thus, 285 kg m−3 seems like a good compromise between measured density of fresh snow in January 2006 (228 kg m−3) and of rather fresh P with graupel-like elements in January 2008 (367 kg m−3). Nonetheless, ρp is a more important parameter in 2D models than has been stated so far and deserves examination (e.g., Sicart et al. 2002). Parameters that control englacial refreezing (8–9) are not well known at all but fortunately prove very insensitive, as also reported in another recent study (Reijmer and Hock 2008). For vertical gradients, the one tied to P exhibits the expected large sensitivity (run 14). The vertical T gradient change (run 15) entails a slower increase of e downslope (because of the constant RH requirement), which strengthens sublimation and weakens melt (see section 4b), and therefore induces an appreciable positive btot response. Although representative roughness lengths (z0m = z0h = 1.7 mm) were obtained for two days on the northern ice field by eddy correlation measurements (Cullen et al. 2007a), which was pioneering for that altitude, they are not constant year-round (Brock et al. 2006). Run 17 therefore employs time-varying roughness: z0m = z0h = 1 mm for fresh snow (Brock et al. 2006), which increases linearly to 4 mm (after t*) for aged snow (Wagnon et al. 1999), and z0m = 20 mm and z0h = 10 mm (Winkler et al. 2009) for likely “penitentes” conditions (daily dewpoint T < −25°C: Mölg et al. 2008). This experiment shows a moderate sensitivity in Table 4. The ensemble mean for Kersten Glacier yields btot = −522 ± 105 kg m−2 yr−1 (one standard deviation uncertainty). This is not drastically out of equilibrium, since btot of receding extratropical mountain glaciers can drop below −1000 kg m−2 yr−1 (e.g., Klok and Oerlemans 2004). Figure 10 delivers the vertical profile of the ensemble b and suggests the likelihood that (i) a small accumulation zone exists between 5700 and 5800 m, but (ii) the present glacier top will vanish in near future if current climate persists. Further, the reference run and its centerline points (btot = −542 and −545 kg m−2 yr−1, respectively) represent the ensemble mean well (Fig. 10 and estimate above).
e. Sensitivity to idealized climate change
Standard tests of the mass balance sensitivity to climate perturbations alter one atmospheric forcing variable (typically T or P) and so do not exhibit real climate change effects except the basic sensitivity of a glacier (e.g., Klok and Oerlemans 2002; Mölg and Hardy 2004). On extratropical mountain glaciers, an idealized 35% change in P amount has roughly the same effect on b as a 1-K change in T (Klok and Oerlemans 2004). Given the results above (section 4), it is not surprising that such a P effect on the total Kersten Glacier is ∼25% stronger than the T effect (Table 5) and appears most strongly (>100%) in the upper glacier region (Mölg et al. 2008). It can even be strengthened easily if the precipitation frequency (Pfr) also increases (Table 5), here by simply shifting each one-third of the daily maximum P in a month to the first two days with P = 0 in the same month. The little impact of T on accumulation in determining the phase of P, as Kersten Glacier lies above the mean freezing level (Fig. 4, run 13 in Table 4), and the largely T-independent Snet as the most variable energy flux (Table 2) explain this behavior. Table 5 additionally targets the role of the coexistence of sublimation and melt. It demonstrates that the existence of sublimation helps the glacier reduce mass loss tremendously, because sublimation is energetically 8.5 times more expensive than melting (LS ≈ 8.5LM). Without sublimation (and associated negative QL), more energy can be expended by the cheap process of melt [Eq. (2)]. This was noted in SEB studies of point locations before (Wagnon et al. 1999) but deserves quantification for an entire tropical glacier (Table 5).
5. Results and discussion: L19
a. Scenarios and resultant glacier behavior
The scenarios of the L19 climate assume that certain conditions in the past climate (e.g., wet and dry spells) were not anomalous but just occurred at a higher or lower frequency than in the present climate. This was demonstrated for moisture changes in tropical Africa in general (Nicholson 2000) and for East Africa specifically (Hastenrath 2001; Mölg et al. 2006; MCGC). Changes in tropical air temperature also show these frequency changes (see below). We consider both a moister or colder climate in the L19.
This means that the conditions (T, RH, P, neff, V, p) on dry (warm) days are replaced by those of wet (cold) days in the model forcing data for moisture (temperature) changes in each concerned month separately. Alternatively or additionally, a certain period (e.g., season) may be overwritten by a different period of the same length. The approach leads to changes in all atmospheric forcing variables, which relies on the measured intervariable dependencies. The order of “wet days” in a month is based on the daily P sum and that of “dry days” on the daily mean RH (see Mölg et al. 2009b; MCGC), while for T changes the order is given by daily mean T. Table 6 summarizes the scenarios. For the temperature scenario E, we attempt to change the pattern of the frequency distribution in T similar to the observed pattern between longer-term cold and warm periods in the REA data (Fig. 11; the choice of which 30% are used does not influence the model results critically). All model runs now use T and P gradients favorable for glacier health (options 14 and 15 in Table 4), so the estimate of the change between L19 and present climate is probably conservative rather than extreme.
Table 7 presents the results of the backward modeling. Three times more frequent anomalous OND precipitation over East Africa in the L19 as in the twentieth century (scenario A), as indicated by conditions in the Indian Ocean (Mölg et al. 2006), induces a strong response of Kersten Glacier but further moisture changes are needed (scenario B) to eventually allow the glacier to grow to its L19 moraine. Hence, higher moisture in the L19 was probably not concentrated to one climatological season, as also supported by scenarios C and D. Colder conditions (scenarios E and F) cannot account for the L19 glacier extent, since they coincide with reduced P on Kibo summit (as measured at AWS3; Table 7). Thus the P effects overwhelm the T effects (cf. section 4) in colder conditions. We tried many further scenarios for the backward modeling, but atmospheric conditions that allow Kersten Glacier to reach the L19 extent always fall within the range shown in Table 7.
b. Causal factors of L19 mass balance
Figure 12 illustrates the driving processes of the L19 steady-state glacier. The equilibrium-line altitude (ELA) at ∼4900 m marks the transition to mass losses that culminate to about 2500 kg m−2 yr−1 at the tongue. A similar magnitude was measured at the tongue of Lewis Glacier on Mount Kenya at the same altitude (Fig. 3). In the moister L19, surface melting increases, but this change is compensated for by reduced sublimation and more refreezing of percolating meltwater (due to more persistent snow cover). Hence, the actual mean mass loss per unit area changes little compared to the present, so the stronger precipitation dominates the L19 mass budget. The higher surface ablation originates from increases in L↓, QS and the ground heat flux, and from a reduction of heat sinks by L↑ and mainly QL (Fig. 12, right). Decreasing Snet, due to both a cloudier sky (ΔS↓ = −14 W m−2) and a higher albedo (0.65 instead 0.54 at present), almost succeeds to balance the higher energy input alone, but an increase in the available energy for melt (QM) of 3.5 W m−2 remains. Hence, as for so many components of the climate system, the probable processes on Kersten Glacier in the past (Fig. 12) resemble those that occur today on shorter time scales, like interannual variability (Table 2, Fig. 9).
c. Evaluation of L19 climate and mass balance
Theoretical considerations of tropical glaciers as well as previous reconstructions of L19 climate in East Africa allow us to evaluate the results of the backward modeling. Regarding the former, the reconstructed L19 profile of b (Fig. 12) meets the theoretical calculations of Kaser (2001) nicely, with a characteristic large (small) vertical gradient in b in the ablation (accumulation) zone below (above) the ELA. Kaser and Osmaston (2002) extended to determining the theoretical accumulation area ratio (AAR; accumulation zone to total glacier surface area) for steady-state glaciers in tropical conditions. They concluded that AAR in the tropics must be greater than in midlatitudes (∼0.67). The b profile in Fig. 12 yields AAR ≈ 0.74 (assuming a fairly constant area–altitude distribution, see section 3), which provides mutual support for theory and measurement-based modeling in this study.
To our knowledge, two quantitative reconstructions of L19 precipitation climate exist. Nicholson and Yin (2001) employed a physically based water balance model to infer rainfall over Lake Victoria from lake level changes. For the period 1858–78 they obtained P = 1941 mm yr−1, which is 116% of the twentieth-century mean determined in Yin and Nicholson (2002). Kruss (1983) modeled the L19 extent of Lewis Glacier, which was a pioneering study but based on much fewer measurements than our computations here. He obtained ΔP ≈ +150 mm yr−1 in the L19, or 119% of the estimated mean annual P of 800 mm on Lewis Glacier during the 1960s to 1980s (Hastenrath 1984). L19 P on Kersten Glacier in Table 7 is 155–239 mm yr−1 higher than at present (128%–144%). These comparisons suggest that the relative P change is smaller in the lowlands (Lake Victoria) than in the summit regions of the mountains. MCGC found the same pattern for seasonal P variability in the Kilimanjaro region and provided an explanation through the mesoscale atmospheric dynamics. Thus, our results (Fig. 12, Table 7) do not contradict other evidence and provide a reference for climate change at 500 hPa over equatorial East Africa.
6. Possible large-scale climate forcings
The initial (rather rapid) L19 moisture drop in equatorial East Africa (EEA) coincides with changes of sea surface temperature (SST) and zonal wind patterns in the Indian Ocean, which appear both in historical observations (Hastenrath 2001) and a paleoclimate simulation with an atmosphere–ocean coupled general circulation model (AOGCM) (Mölg et al. 2006) and lead to the strongest observed retreat rates of Kibo’s glaciers in the early twentieth century (Cullen et al. 2006). A recent study of observed SSTs found an anomalously high frequency of “negative dipole” events in the Indian Ocean during 1880–1919 (Ihara et al. 2008), which fosters precipitation deficit in EEA and agrees with the AOGCM results (Mölg et al. 2006). Rapid P changes in equatorial Africa from SST-forced circulation changes were also observed in the 1960s for the western continent (Nicholson and Webster 2007), so they appear not to be unusual.
However, negative dipoles in the Indian Ocean have become less frequent in the most recent decades (Ihara et al. 2008) but glaciers on Kibo are still shrinking (Cullen et al. 2006). Vegetation changes on Kilimanjaro (Hemp 2005) as well as problems for nomadic life forms (Mwangi and Desanker 2007) and agricultural capacity (Funk et al. 2008) also point to present drying in EEA. There is increasing evidence that global warming (like in most recent decades) affects regional P changes in the tropics, where convecting zones generally get wetter and nonconvecting zones get drier (e.g., Held and Soden 2006; Chou et al. 2009). The boundary zone (the convective margin) is expected to shift inland for a land region with moisture inflow from a neighboring ocean (Lintner and Neelin 2007), a precondition largely met for Kilimanjaro (Chan et al. 2008).
Figure 13 compares precipitation of the most recent 10 years with 1979–88 along 3.75°S in EEA. MAM precipitation appears to have decreased at all longitudes, which conforms to later onsets of MAM rainfall detected in station data (Camberlin and Okoola, 2003). Funk et al. (2008) identified an anticyclonic circulation over EEA in this season as the main cause (driven by Indian Ocean warming), which disrupts moisture advection from the east. OND rainfall, on the other hand, does exhibit a structure in the vein of the convective margin theory (cf. Fig. 4 in Lintner and Neelin 2007), with drying confined to the land >28°E (Fig. 13). While this very preliminary analysis neglects other involved complexities (e.g., meridional shifts, changes of moisture source regions), it indicates at least that a convective margin shift could contribute to the present drying in EEA. Mentioned changes in large-scale moisture advection are likely to feed back with local evaporation in EEA, but the latter’s contribution to precipitation (“recycling”) is relatively small in general (cf. maps in Trenberth 1999) and more dependent on changes in local surface characteristics (e.g., Hemp 2005).
Kersten Glacier shows an annual net mass loss in the present climate, which explains its observed shrinkage in recent decades (Cullen et al. 2006). The year-to-year variability in the total mean specific mass budget is strongly influenced by variability in atmospheric moisture over Kilimanjaro’s summit (precipitation, air humidity, cloud cover fraction). Dry (wet) years typically coincide with little (much) mass accumulation, high (low) sublimation, and low (high) surface melt, which produces a small (large) mass turnover. The most variable energy flux between such contrasting conditions is net shortwave radiation. These processes clearly demonstrate the dependence of Kersten Glacier on moisture changes, in particular precipitation.
Several scenario calculations suggest that higher precipitation (160–240 mm yr−1), higher relative and absolute humidity (2%–6% units and 0.1–0.3 hPa, respectively), and increased fractional cloud cover (2%–4% units), which are all mutually linked, enable Kersten Glacier’s steady-state extent of the late nineteenth century. Changes in local air temperature, air pressure, and wind speed appear to be less important. In terms of causal factors, the changes between the late nineteenth century and present climate resemble the variability seen during the 3-yr study period. However, the refreezing of meltwater in the more persistent snow cover of the past climate limits the effect of increased surface melt more strongly than during wet years of the present, which makes solid precipitation even more dominant in the mass budget of the past.
Our results complement other evidence of drying in East Africa during the twentieth century (e.g., Hastenrath 2001; Hemp 2005; Funk et al. 2008; MCGC) and point to an amplified drying at very high altitude like the summit of Kilimanjaro. While probable causes of the initial moisture drop are known (Hastenrath 2001; Mölg et al. 2006), the most recent drying merits further investigation in particular for the October–December wet season. A convective margin shift under present global warming (e.g., Lintner and Neelin 2007) represents one hypothesis. However, if atmosphere–ocean dynamics over the Indian Ocean become increasingly favorable to moist conditions in East Africa in the twenty-first century, as simulated by GCMs in response to future warming (Vecchi and Soden 2007), monitoring glacier behavior on Kilimanjaro could provide a further key to understanding implications of large-scale climate change for regional climate change in the tropics at high altitude.
This study is funded by the Austrian Science Foundation (FWF; Grants P17415-N10, P20089-N10). DRH is supported by the National Science Foundation and NOAA Office of Global Programs, Climate Change Data and Detection Program (Grant No. 0402557), and by the U.S. Global Climate Observing System. Local support for Kilimanjaro is provided by the Tanzania Meteorological Agency (TMA; with special thanks to Tharsis Hyera), the Commission of Science and Technology (COSTECH), Tanzania and Kilimanjaro National Park Authorities (TANAPA and KINAPA), and the Tanzania Wildlife Research Institute (TAWIRI). Gridded data were provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, from their Web site at http://www.cdc.noaa.gov/.
Postprocessing of Data
Potential radiative heating of the unventilated temperature sensor is corrected by the method described in Mölg et al. (2008). Only 0.79% of all data are affected. Because of the subfreezing T (Fig. 4), RH series are rescaled to account for saturation with respect to ice rather than liquid water (e.g., Cullen et al. 2007b). Regarding radiation data, the upward-looking sensors are particularly prone to measurement errors (Van den Broeke et al. 2004). For shortwave data, day-to-day changes in daily accumulated albedo αacc (Oerlemans and Knap 1998) are initially compared to snow depth changes from the SR50 record (for the first 824 days). Whenever αacc rises (drops) by >0.2 but snow depth decreases (increases) by >0.01 m (which contradicts theory, see section 3a), G data were reexamined. All suspicious days indeed show RH close to 100%, which indicates a situation favorable for icing of the upward-looking radiation dome (Van den Broeke et al. 2004). The diurnal cycle in G of affected days was then multiplied by a correction factor to produce a daily mean of G that equals the mean over days with an average RH > 95% and within ±2 days of the suspicious day. In that way, the recalculated αacc does change as one would expect from theory. However, only a few G data are affected (0.65%). We screen longwave data for window heating offset and riming of the sensor (Van den Broeke et al. 2004), but they don’t show obvious errors. In the wind speed data, nine periods of zero wind longer than one day occur (2.1% of total). These coincide with RH ∼100% and largely uniform (spurious) wind direction, again indicative of icing and therefore most likely freezing of the anemometer (as observed during the January 2008 visit). Affected periods were filled with hourly data obtained from the mean diurnal cycle over the two adjacent days. Air pressure data do not exhibit any apparent errors.
Corresponding author address: Dr. Thomas Mölg, Department of Geography, University of Innsbruck, Innrain 52 6020, Innsbruck, Austria. Email: email@example.com