Brief communication: Glacier thickness reconstruction on Mt. Kilimanjaro

Glaciers on Kilimanjaro are unique indicators for climatic changes in the tropical mid-troposphere of Africa. The history of severe glacier area loss raises concerns about an imminent future disappearance. Yet, the remaining ice volume is not well known. We reconstruct thickness maps for 2000 and 2011 for the Northern Icefield (NIF) and Kersten Glacier (KG) that are informed by ground-truth thickness measurements and multi-temporal satellite information. For 2011, we find mean thickness values of 26.6 and 9.3 m, respectively. The existing consensus estimate for global glacier ice thickness shows 10 unrealistically thick values for KG in areas that are meanwhile ice-free.


Introduction
The importance of tropical glaciers at high elevations as unique climate indicator for the tropical mid-troposphere has been previously highlighted (e.g. Kaser 2001, Kaser et al. 2004, Mölg et al. 2009a. As one of few remaining tropical locations with existing glaciers, Mt. Kilimanjaro, a stratovolcano with an elevation of 5895 m.a.s.l., is located in East Africa close to the 15 Tanzania-Kenya border (3°04' S / 37°21' E) (Fig. 1,overview). In addition to the very high elevation, the free-standing nature of the mountain causes the glacier on top of the summit to be directly exposed to tropospheric flows, minimizing the forcing of local climate on the glacier and creating a unique opportunity to study the mid-troposphere climate.
The modern glacier recession on Kilimanjaro has been well documented and mapping approaches have shown that from an estimated ice extent of 11.4 km² in 1912, only 1.76 km² remained in 2011, constituting to a severe 85% reduction in glacier 20 area (Cullen et al. 2013). While glaciological research on Kilimanjaro has been focused on mapping glacier area and glacier retreat (Kaser et al. 2004, Cullen et al. 2013, as well as mass and energy balance studies (Mölg et al. 2003(Mölg et al. , 2008(Mölg et al. , 2009, the research on the ice thickness of different glaciers on Kilimanjaro has been comparably sparse (Bohleber et al. 2017). However, in light of severe glacier recession, an assessment on current glacier thickness is important to better determine future recession.
A recent effort to reconstruct the distributed ice thickness for all glaciers outside of Antarctica from a consensus of up to 5 25 models (Farinotti et al. 2019), provides the possibility to compare the global estimate, which generated ice thicknesses from an ensemble of 2 (Northern Icefield, NIF) and 3 (Kersten Glacier, KG) models ( Fig. 2 C), to results of this study. Uncertainties of the global consensus estimate, such as the separation of NIF into three separate glacier entities in the most recent version of the Randolph Glacier Inventory (RGI6.0), as well as the consensus estimate not being informed of local thickness observations, https://doi.org/10.5194/tc-2019-310 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License.
have already been recognized as deficiency of the global estimate (Farinotti et al. 2019) and the recently observed separation 30 of KG into two fragments casts further doubts on the high thickness values there.
Here, we present well constrained thickness maps for KG and NIF using a mass conserving reconstruction approach that readily assimilated thickness measurements for the first time (Fürst et al., 2017). In different experiments we test the influence of varying thickness input for the glacier state of 2000, where we rely on digital elevation data with global coverage (SRTM), pursuing a new calibration strategy that uses multi-temporal satellite information on geometric changes in absence of 35 observational data on KG. In a third experiment, we combine the very high resolution digital elevation model KILISoSDEM (Sirguey et al. 2014) with results from the previous experiments to improve the reconstruction of the ice cliffs at NIF.

Data
The two glaciers share the following input for the distributed surface mass balance (SMB) model and the ice thickness reconstruction approach: climate data measured by the automatic weather station (AWS) located on KG (

Mass balance modelling
The mean annual climatic surface mass balance fields were generated using version 2.4 of the distributed, physically-based 60 MB model by Mölg et al. (2008Mölg et al. ( , 2009a, using meteorological input from the aforementioned AWS (Suppl. Fig. 1). As the full MB model has only been tested on and verified for KG before, we slightly altered the model code so that it allowed refreezing of meltwater on a bare ice surface with a slope angle below 5 degrees. In this way, the model is capable of reproducing the observed surface height changes observed by a Sonic Ranger mounted to the AWS. Because of the low slope angles, meltwater cannot run off from the surface of NIF's planar top before refreezing sets in . 65

TDX processing
The detailed processing steps to create a DEM, as well as a surface height change layer, from the TanDEM-X image can be found in e.g. Braun et al. (2019). Surface elevation changes between 2000 and 2011 were inferred from DEM differencing with respect to the SRTM DEM (Suppl. Fig. 2). Positive values, which indicate a height gain in the TanDEM-X layer, were removed as a height gain outside the 2011 glacier extent implies an increase in glacier thickness from 2000 to 2011, which is 70 unlikely. In total we removed 92 of 602 grid cells with a mean height gain of 0.19 m/a for NIF and 14 of 254 grid cells with a mean height gain of 0.25 m/a for KG.

Margin thickness generation
For KG, no in-situ thickness measurements are available. Therefore, multi-temporal DEM and glacier outline information is used to infer past ice thickness. First, glacier retreat is delineated from outline information in 2000 and 2011 ( Fig. 1 hatched 75 area). In the nowadays ice-free area, contemporaneous elevation changes (2000-2011) then give information on past ice thickness.

Ice thickness reconstruction
A detailed description of the two-step ice-thickness reconstruction, of which we only used the first model step as surface velocities were not available, can be found in Fürst et al. (2017). The reconstruction approach is based on the principle of mass 80 conservation and computes a glacier-wide flux field from the difference between the surface mass balance (Section 3.2) and contemporaneous elevation changes. The flux solution is converted into thickness values using the Shallow Ice Approximation (SIA; Hutter, 1983). This conversion involves the ice-viscosity parameter B, which is a-priori unknown. This parameter stems from assuming a Glen-type flow law, linking deviatoric stresses to strain rate components ɛ̇ij via the effective viscosity the higher DEM quality in 2011, the resolution was increase to up to 2 m. The processing was conducted separately for NIF and KG. The coupling length parameter is defined as a multiple of the local ice thickness serves to smooth the surface slope 90 during the reconstruction. For KG, the parameter is set to 1, a typical value for valley glaciers (Kamb & Echelmeyer, 1986).
For NIF, it had to be reduced because otherwise the step in the elevation profile over ice cliffs was not imprinted in the thickness field. A compromise value of 0.3 was chosen to guarantee sufficient smoothing of the flux streamlines.

Experimental Setup
The general strategy is to reconstruct a thickness field for KG and NIF combining SMB, elevation changes and glacier 95 geometry with in-situ measurements of ice thickness for two points in time, namely 2000 and 2011. In Experiment 1 we reconstructed the glacier state for 2000 with the generated margin thickness data for both NIF and KG. At these thickness measurements, an appropriate viscosity value is inferred. As KG is rather small, we expect a homogeneous ice viscosity. In Experiment 2, we therefore average all viscosity values and impose a constant value for KG. In this way, generated thickness values are no longer reproduced but spurious spatial viscosity variations stemming from the generic data are suppressed. For 100 NIF in Experiment 2, we chose to use the thickness measurements from Bohleber et al. (2017) as input, to check how observational data influences the glacier wide ice thickness in comparison to only using margin thickness information. With Experiment 3 we reconstructed the glacier state in 2011, using again the mean ice viscosity of 2000 as input for KG and the observations from Bohleber et al. as input for NIF, as these thickness inputs produced the best results for the respective glaciers in the 2000 reconstruction. We decided to do a reconstruction of the 2011 glacier state, to test the influence of a higher 105 resolution DEM, as well as a higher model resolution on the reconstruction result. Table 1 summarizes the three different experimental setups (Tab. 1).

Results
In the following section we discuss the results of the three experimental setups (Section 3.5). Results show generally larger ice thickness for NIF than for KG in all three experiments. Experiment 1 (Fig. 2 A) shows thicker ice of up to 15 m at the flat top 110 parts of KG and thickness up to 7.5 m at the central area of the steep part, with patches of thinner ice towards the glacier margins and a mean ice thickness of 6.2 m. NIF is up to 40 m thick in its center, decreasing towards the glacier margins and towards CG and has a mean ice thickness of 13.7 m.
In comparison, results from Experiment 2 show a similar thickness pattern on KG (Fig. 2B). For NIF, the magnitude differs significantly (Fig. 2 B). Now one large part of NIF's flat area and two smaller parts of CG exceed 40 m. Moreover, the ice 115 thickness in the steeper western areas of NIF and CG has increased by a factor of 2. The mean ice thickness also increases to https://doi.org/10.5194/tc-2019-310 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. 23.4 m. KG shows a similar thickness, but a distribution is smoother as compared to Experiment 1 and the flat top part of KG is still thicker than the slope part. The mean ice thickness with 6.9 m is very similar to Experiment 1.
For Experiment 3 (Fig 1), KG is now split into two parts and shows an ice thickness of up to 10 m at the flat top part and most of the slope being between 5-7.5 m thick. KG's mean ice thickness is 9.3 m. NIF's thickness distribution is similar to 120 Experiment 2, with the thickest areas of over 40 m at its flat part on the plateau. For NIF, the decrease in thickness is less noticeable than the lateral retreat and decrease of glacier area. The mean ice thickness of NIF in Experiment 3 is 26.6 m.

Discussion
We will first discuss the reconstructions for the year 2000 ( Fig. 2): Generally, the consensus thickness map (Fig. 2C) shows less difference in thickness magnitude between KG and NIF as compared to our reconstruction. For KG, no ice thickness 125 measurements are available, and it is uncertain to what extent the generated thicknesses along the glacier margin (Section 3.3) are useful to inform the reconstruction. The margin values result in a spatially varying viscosity field, which is transmitted into the ice thickness field (Fig. 2 A). As no strong viscosity variations are expected for the small KG, a second run was conducted with constant viscosity (Fig. 2 B). Results indicate a thick central flow unit, as one might expect for a steep glacier. In this case, the result provides a smoother ice thickness distribution, with higher thickness in the center of the glacier and thinning 130 towards the margins. Whether this smoothed result, or the less smooth reconstruction generated from margin observations, is closer to the reality is unclear. Considering the input of thickness observations was only available on the glacier margin, calibrating the model, which is based on ice flow, with data where ice flow is generally small, puts the reconstruction to its limits. However, as the thickness of most glaciers on Earth is unsurveyed, the use of margin thickness information, generated from outline differences enabled a local glacier-specific viscosity tuning which might be preferential to an empirical 135 temperature relation (Huss and Farinotti, 2012). The consensus map shows a similar thickness pattern as Experiment 2. The NIF's peculiar geometry is difficult to reconstruct using generic thickness observations around the margin (Fig. 2A). The ice is much too thin in the interior (Fig. 2 A), which underestimates the ice core lengths from Thompson et al. by 48,52 and 71% for the core locations C1, C2 and C3 respectively (Fig 2. A). When the interior GPR measurements are assimilated (experiment 145 2), differences decrease to 10, 17 and 53%. The consensus estimate for NIF also underestimates the ice core lengths by 34, 38 and 72%, respectively. The consensus map is a weighted mean of two model results on NIF. The better model shows a difference of 10, 12 and 76% for C1, C2 and C3, respectively. C3 is located near a separation line between the three RGI https://doi.org/10.5194/tc-2019-310 Preprint. Discussion started: 3 February 2020 c Author(s) 2020. CC BY 4.0 License. glacier units for NIF. This might be a reason why the consensus estimate has problems to reproduce a larger value. The upper value of the error margin has to be assessed cautiously, as the location of the drill site of one of the cores coincides with the 150 margin of the three separate glacier entities on NIF, where model 1 and consequently the composite show a notably larger difference to the measurements than for the other two locations. While this location is the one where our reconstructions show the largest difference to the measurements as well, the proximity to an internal separation line should not influence our reconstruction. If we only consider the other two core locations (C2, C3), our reconstructions underestimate the measured thickness by 48 and 52% and 10 and 17% for Experiment 1 and 2, respectively. The consensus result (model 01; model 02) 155 underestimates the measurements by 34 and 38% (10-12%; 65-72%). This shows that one of the models is strongly influenced by the separation of NIF into three separate entities and in turn worsens the composite result. NIF's mean ice thickness of 21.5 m is similar to our results, especially to the mean of Experiment 1 and 2 (21.6 m).
Analyzing the three separate models of the consensus, none of them creates a large bias, which could influence the composite result, but the KG reconstruction is too thick in areas where the glacier has disappeared in 2011. Thus, we can confidently say 160 that the consensus estimate is too thick for KG, but we cannot confirm how well our reconstruction results are able to depict reality.
Experiment 3 repeats the reconstruction for the year 2011 at a very high resolution. The overall distribution of ice thickness is barely affected by the increase in resolution. This is a desired effect and stems from the coupling length that scales with the thickness. Further experiments with 10 and 5 m model resolution showed barely any difference in thickness distribution, 165 verifying this effect. The mean ice thickness for KG and NIF have increased in comparison to Experiment 2 to 9.3 m and 26.6 m, respectively. This increase likely stems from the difference in DEMs that were used for the reconstruction. While an increase for the period 2000-2011 in mean glacier-wide ice thickness is unlikely, results for NIF now match the thickness estimates from Bohleber et al. (2017), where the very high resolution KILISoSDEM is used as input as well. While Bohleber et al. (2017) have reconstructed parts of NIF with a linear extrapolation of the bedrock, the thickness distribution in the areas 170 located on top of the summit plateau matches the result of our reconstruction. However, the slope area is generally thinner, which is caused by the influence of several ice-free areas that could not be identified from the coarse 2011 Landsat image.
From this we can deduce that the increase in ice thickness could be caused by a difference in DEM and not increased model resolution.

Conclusion & Outlook 175
The aim of this study was first and foremost to accurately determine the volume and distribution of ice on Mt. Kilimanjaro.
For this purpose, thickness observations are generated from multi-temporal information on the glacier geometry. These generic values can be inferred in areas that became ice-free or were ice-covered in the past. We assess the utility of lateral thickness information in constraining glacier volume. For Kersten Glacier, we report significantly smaller thickness values as compared to the current consensus estimate. The latter reconstruction is shown to be inconsistent with the observed glacier separation.
For NIF geometry, the lateral retreat information seems less utile as central ice thickness is strongly underestimated. Reasons for this worse performance might be the complex topography and the dynamic inactivity of NIF. We therefore speculate that thickness information from retreat is most useful in dynamically more active areas.
Furthermore, the reconstructions reveal that if there are no thickness observations available, better results can be achieved with a mean viscosity value as input for ice thickness, instead of margin ice thickness generated from DEMs and glacier outline 185 difference. However, if the glacier topography is as complex and peculiar as for NIF, both mean viscosity values and generated margin ice thickness observations underestimate the overall glacier thickness. In the case of NIF, ice thickness based on observational data and ice thickness based on mean viscosity differed by a factor of 2. Based on results obtained for NIF and KG, the mean viscosity can be used in future studies to generate an estimate for the ice thickness of retreating glaciers without actual thickness observations. As ice thickness reconstructions on a global scale can only rely on thickness observations for a 190 small percentage of all glaciers, the approach of using generic margin thickness observations, presented in this study, can provide additional input for glaciers with visible retreat. The first validation of this generic-thickness approach on the glaciers on top of Mt. Kilimanjaro might have been influenced by the complex topography of NIF and the absence of ground-truth data on KG. However, ice thickness reconstructions were still comparable to results by Bohleber et al. (2017) and Farinotti et al. (2019), which supports the use of the approach with mean viscosity values as ice thickness input to adjust the viscosity locally. 195 We imagine that the inference of margin thickness values for retreating glaciers is readily transferrable to many glaciers worldwide and it would provide a mean for glacier-specific calibration of reconstruction approaches on regional or global scales.

Acknowledgement
CS received primary funding from the project BR2105 /14-1 within the DFG Priority Program "Regional Sea Level Change 200 & Society". JJF was funded by the German Research Foundation (DFG) under grant number FU1032/1-1. Results presented in this publication are based on numerical simulations conducted at the high-performance computing center of the Regionales Rechenzentrum Erlangen (RRZE) of the University of Erlangen-Nürnberg. We would also like to thank Nicolas Cullen and Pascal Sirguey for constructive discussions and for providing the KILISoSDEM.

Author contribution 205
CS led writing of the manuscript, in which she received support from all authors. The research aims and setup were developed in regular discussion with JJF, TM and MB.