Energy balance to compute debris thickness (model F12)

The F12 approach is described briefly below, as it was used and further developed in view of its application to Himalayan glaciers. Debris thickness d is mapped for every pixel of the spaceborne debris surface temperature T _s by solving the energy balance for the debris-covered glacier surface with d as a residual. The energy-balance equation describes that the change in heat stored within a defined volume of debris is equal to the sum of all fluxes towards the debris volume, given by

where ΔD is the rate of change in heat stored within debris, S _net is net shortwave radiation, L _net is net longwave radiation, H is the sensible turbulent heat flux and G is the conductive heat flux. Energy contributions provided by rain and latent heat are neglected, assuming that the debris layer is dry at the time of satellite image acquisition (Reference Brock, Rivera, Casassa, Bown and AcuñaBrock and others, 2007, Reference Brock, Mihalcea, Kirkbride, Diolauiuti, Cutler and Smiraglia2010). Fluxes are considered positive when directed towards the debris layer.

Energy transfer within debris is reduced here to conductive heat flux G, neglecting other transfer energy sources and sinks (e.g. convection, advection, phase changes or ventilation). The one-dimensional conductive heat flux G is estimated from Fourier’s law as a function of the effective conductivity K of the debris material and the temperature gradient within the layer dT/dz. Assuming a linear temperature gradient within the debris surface, and temperature of the ice–debris interface to be at melting point,

where T _s is debris surface temperature and d is debris thickness. Conductivity K is strongly influenced by the moisture content of the debris surface. Here we assume that the uppermost level of the debris layer is dry and use a spatially constant value for the conductivity of dry material.

The rate of change in heat stored, ΔD, at the time of the image acquisition results from the morning warming of the debris layer. This heat flux depends on the density and specific heat of the debris layer, the rate of temperature change and debris thickness. Stored heat cannot be determined precisely in this approach as satellite thermal band imagery has a daily resolution and does not represent a sequence of images. Reference Foster, Brock, Cutler and DiotriFoster and others (2012) therefore assumed that the rate of change in heat stored is a fraction F (hereafter referred to as ‘stored heat factor F’) of the conductive heat flux in the form

Using these assumptions, debris thickness is given by

Net shortwave radiation S _net is calculated from S _in(1 − albedo), using a spatially constant value for albedo, according to Reference Rounce and McKinneyRounce and McKinney (2014). Reference Foster, Brock, Cutler and DiotriFoster and others (2012) proposed taking a spatially constant value for albedo for Miage glacier due to minor measured spatial variation (Reference Brock, Mihalcea, Kirkbride, Diolauiuti, Cutler and SmiragliaBrock and others, 2010). Incoming shortwave radiation S _in is computed using the equations of Reference Strasser, Corripio, Pellicciotti, Burlando, Brock and FunkStrasser and others (2004) for a clear-sky situation. Incoming longwave radiation can be estimated through parameterizations as described by Reference BruntBrunt (1932) or Reference BrutsaertBrutsaert (1975). Outgoing longwave radiation L _out is computed based on the Stefan–Boltzmann law, using surface temperature T _s from the satellite images and spatially constant emissivity. The sensible heat flux H is calculated using the bulk aerodynamic approach

where u is wind speed, T _a and T _s are air and debris surface temperatures, respectively, and D is a bulk exchange coefficient. D can be calculated using a stability correction based on the bulk Richardson number (Reference BrutsaertBrutsaert, 1975; Reference OkeOke, 1987).

Air temperature T _a distribution over glaciers is commonly extrapolated from station data by applying a spatially constant lapse rate, which may not hold for debris-covered glaciers. Reference Foster, Brock, Cutler and DiotriFoster and others (2012) proposed that air temperature is strongly related to surface temperature T _s in the morning, with strong insolation and low wind speed. At this time of the day, under strong insolation, above debris-covered surfaces, sensible as well as longwave heat fluxes are directed upwards transferring heat to the near-surface air (Reference Reid and BrockReid and Brock, 2010). Based on data from Miage glacier, Reference Foster, Brock, Cutler and DiotriFoster and others (2012) computed an empirical relationship in the form

where T _a is air temperature, T _s is surface temperature and a and b are empirical parameters that can be derived based on measurements of air temperature and surface temperature on debris-covered glaciers. Table 2 shows values for the linear regression of air temperature against surface temperature for three sites in the European Alps and the Himalaya with different elevation and debris cover. Here we take values from Lirung and Miage glaciers to extrapolate air temperature across the glacier surface.

Extension of the F12 approach

In view of the different regional context and comparably limited number of existing field-based data, we modified the F12 approach to make it applicable to an extensively debris-covered glacier in the Indian Himalaya. The modifications, illustrated schematically in Figure 2, are mainly related to: (1) the vertical temperature gradient within the debris layer, which is possibly not linear within the debris layer in the morning; and (2) the energy balance, which has been closed in this application for the uppermost layer (uppermost 10%) rather than for the entire debris cover.

Fig. 2.

Schematic view of changes added to the F12 approach presented by Reference Foster, Brock, Cutler and DiotriFoster and others (2012) to map debris thickness with reanalysis data and thermal band surface temperature. Additions to the original approach are highlighted in red.

Observations have shown that temperature gradients within a debris layer are not linear at ∼10.00–12.00 local time, as illustrated in Figure 3 (Reference Conway and RasmussenConway and Rasmussen, 2000; Reference Han, Ding and LiuHan and others, 2006; Reference Nicholson and BennNicholson and Benn, 2006; Reference Reid and BrockReid and Brock, 2010; Reference Lejeune, Bertrand, Wagnon and MorinLejeune and others, 2013). Taking a linear gradient within the entire debris leads to a large underestimation of debris thickness. Recently published results by Reference Rounce and McKinneyRounce and McKinney (2014) also indicate that accounting for a nonlinear temperature gradient is crucial. Here instead of taking a linear gradient, we estimate the depth d ₀ at which the linear temperature gradient would reach 0°C. Based on field measurements, the 0°C depth is estimated to be half of the total debris depth (d ₀ = i _d · d). Conductive heat flux is then computed by

where T ₀ is 0°C and the factor i _d is 0.5.

Fig. 3.

Temperature profiles within the supraglacial debris at about noon local time as published in the literature.

We assume that F is a function of debris thickness and thus should be represented by a dependent stored heat factor F(d). For a thin debris cover, we assume that the surface is close to an equilibrium temperature and is not warming strongly in the morning. The energy from the surface is conducted through the debris layer and melts the ice beneath, leading to a linear temperature gradient and a small F (small fraction ΔD/G). Where debris is thick, the diurnal temperature maxima of the debris surface increases. Energy from the surface is used to heat the debris layer, resulting in a larger F (large fraction ΔD/G). Based on data found in the literature (Reference Conway and RasmussenConway and Rasmussen, 2000; Reference Nicholson and BennNicholson and Benn, 2006; Reference Reid and BrockReid and Brock, 2010; Reference Foster, Brock, Cutler and DiotriFoster and others, 2012; Reference Lejeune, Bertrand, Wagnon and MorinLejeune and others, 2013) this relationship can be expressed by a polynomial function of the order 2 with R ² = 0.97 (Fig. 4). However, we found that the iterated result converges to 0 for F < 1. Therefore, the assumption of a linear relationship with F = 1 where d = 0 is necessary. Debris thickness d is finally derived by rearranging the energy balance and Fourier’s law in order to get the result by iteration:

Fig. 4.

Stored heat factor F for the upper ∼10–50% of the debris layer at about noon as a function of debris thickness using vertical temperature profiles published in the literature. Note that F increases with increasing debris thickness.

Here the parameters m and n are taken from the linear regression in Figure 4. The debris thickness d is evaluated iteratively; it is first evaluated for a neutral F, this value is used to make an initial estimate of F, which is then used to recalculate d; the procedure is repeated until the value of d converges.

Satellite-derived surface temperatures

The F12 model is transferable in time (applicable to different satellite images), given that the pattern and distribution of computed d are consistent for different satellite images. Since d is calculated as a function of T _s, with increasing computed debris thickness for warmer surfaces, spatial patterns of T _s are also required to be consistent among different satellite images. Hence an important test for the model performance and transferability involves comparing the spatial distribution of remotely sensed surface temperatures between different satellite scenes.

Figures 5 and 6 show that the distribution and pattern of the standardized surface temperatures are highly consistent for different scenes and sensors (ASTER and Landsat 7), despite differences in mean surface temperature. Relatively high surface temperatures are observed for the medium moraine and the eastern edge of the glacier tongue, as well as across the lowest part of the tongue. Figure 7 shows that mean surface temperature is clearly correlated with modelled incoming shortwave radiation at the time of satellite image acquisition, with R ² = 0.75, and to air temperature from reanalysis data (for 32.5° N, 77.5° E at 600 hPa), with R ² = 0.74.

Fig. 5.

Histograms of surface temperature and standardized surface temperature for the debris-covered area of Bara Shigri Glacier, based on seven ASTER and six Landsat 7 satellite images taken in dry season conditions between 1999 and 2012 on clear-sky days. Date format is dd.mm.yyyy.

Fig. 6.

Maps of standardized surface temperature for Bara Shigri Glacier based on three ASTER thermal band satellite images. Relatively high surface temperatures are observed for the medium moraine and the eastern edge of the glacier tongue, as well as across the lowermost part of the tongue. Date format is dd.mm.yyyy.

Fig. 7.

Incoming shortwave radiation (using the approach described by Reference Strasser, Corripio, Pellicciotti, Burlando, Brock and FunkStrasser and others, 2004) and NCEP/NCAR air temperature versus mean ASTER (open symbols) and Landsat 7 (solid symbols) surface temperature of the debris-covered glacier area of Bara Shigri Glacier.

Reanalysis data versus station data

Reanalysis data, including wind speed, air temperature, relative humidity and precipitable water, are used to compute shortwave and longwave radiation and turbulent sensible heat flux. We compare reanalysis data for wind speed and air temperature from both ERA-Interim and NCEP/NCAR with the Keylong station data to check their representativeness within the Himalayan region.

700 hPa level air temperature from NCEP/NCAR and ERA-Interim has been used for the computation of downwelling longwave radiation and correlates well with measured air temperature from Keylong station (Fig. 8a). Reanalysis data tend to slightly underestimate local values for air temperature, probably because they represent free atmosphere and do not consider boundary layer effects.

Fig. 8.

Reanalysis (ERA-Interim and NCEP/NCAR) 700 hPa level data versus station measurements at Keylong for (a) air temperature and (b) wind speed.

Wind speed from reanalysis data, used for the computation of turbulent fluxes, is not correlated with measurements from meteorological stations (Fig. 8b). To overcome the low quality in input variables for wind speed, we conducted a statistical analysis of the Keylong and Pyramid station data. The analysis shows that wind speed for clear-sky days from September to November at 06.00 UTC is mostly 1–4 m s⁻¹ for Keylong station (Fig. 9a) and 3–7 m s⁻¹ for Pyramid station (Fig. 9b), whereas the distribution is much more scattered for the entire period including all time steps. Given the lack of correlation between reanalysis and measured wind speed, we used measured wind speed at the adjacent station and assumed that wind speed on Bara Shigri is close to 3 m s⁻¹ at the time of satellite overpass.

Fig. 9.

Histograms of wind speed measured at (a) Keylong and (b) Pyramid meteorological stations during all seasons and times (grey bars) as well as for September–November from 05.00 to 07.00 UTC for clear-sky conditions with incoming shortwave radiation >800 W m⁻² (black bars).

Parameterization of radiative fluxes

Longwave radiation fluxes are computed for the Khumbu region using the empirical relationships introduced by Reference BruntBrunt (1932) and Reference BrutsaertBrutsaert (1975). The parameterization by Reference BruntBrunt (1932) relates clear-sky emissivity to screen-level measurements of air temperature and humidity, using empirically determined coefficients. Reference BrutsaertBrutsaert (1975) developed a more physically based approach based on Schwarzschild’s transfer equations for simple atmospheric profiles of a certain air temperature and vapour pressure. Parameterizations of downwelling longwave radiation L _in are compared with meteorological data from Pyramid station under clear-sky conditions (S _in > 800 W m⁻²) from September to November at about noon. In addition, we compared ERA-Interim downwelling longwave radiation to the station measurements. ERA-Interim overestimates incoming longwave radiation by an average of ∼60 W m⁻² compared with station measurements (Fig. 10a). Both parameterizations are accurate for measured longwave radiation below ∼230 W m⁻², whereas an underestimation is found for the higher end of measured high longwave radiation. For the purpose of this study, the method of Reference BrutsaertBrutsaert (1975) seems to be more accurate, with a root-mean-square error (RMSE) of 5.9 W m⁻² (compared with RMSE = 11.2 W m⁻²) for situations with measured incoming longwave radiation below 230 W m⁻².

Fig. 10.

Modelled versus measured radiation fluxes for Pyramid station at 06.00 UTC between September and November 2002–08 where incoming shortwave radiation was above 800 W m⁻²: (a) downwelling longwave radiation and (b) incoming shortwave radiation. For longwave radiation, RMSE and R ² are given for measured incoming longwave radiation below 230 W m⁻².

Shortwave radiation was computed using the parameterization proposed by Reference Strasser, Corripio, Pellicciotti, Burlando, Brock and FunkStrasser and others (2004) accounting for the total transmittance of the atmosphere, given as the product of the individual transmittances. Precipitable water is taken from ERA-Interim reanalysis data, and a spatially constant value for incoming shortwave radiation is assumed. Direct solar radiation under clear-sky conditions is assumed to vary only slightly within the elevation range ∼3900–5000 m a.s.l. for Bara Shigri. In addition, the surface of debris-covered glacier tongues is assumed to be flat, as any improvement in model performance including the slope proved to be questionable, at least for the European Alps (Reference Foster, Brock, Cutler and DiotriFoster and others, 2012). As the local time of image acquisition is almost midday, spatial variations due to horizon obstruction and reflections are neglected. The parameterization of incoming shortwave radiation is compared with the Pyramid data for September–November for 06.00 UTC in order to validate the approach developed by Reference Strasser, Corripio, Pellicciotti, Burlando, Brock and FunkStrasser and others (2004) for clear-sky conditions at about noon in the Himalaya. As can be seen from Figure 10b, the parameterization of shortwave radiation in some cases tends to strongly over- or underestimate measured values. Nevertheless, the performance proves to be good (RMSE = 77.4, R ² =0.56) for the range of measured shortwave radiation of 800–1300 W m⁻² and corresponds to the approximate range of clear-sky situations.

Sensitivity analysis

In a next step we introduce a sensitivity analysis to see how the uncertainty of modelled debris thickness can be apportioned to energy fluxes, input variables and parameters. We therefore run the model 10³ times for each variable and parameter by varying one parameter while keeping the others constant. Input variables for wind speed were assumed to be normally distributed. Parameters were varied uniformly within a range of plausible values taken from the literature or normally distributed with estimated mean values and standard deviations (Table 3). Long- and shortwave radiation fluxes were computed with the parameterizations introduced by Reference BrutsaertBrutsaert (1975) and Reference Strasser, Corripio, Pellicciotti, Burlando, Brock and FunkStrasser and others (2004), respectively, and varied by the standard deviation, which was derived by comparing computations with station data. The sensitivity analysis was based on the satellite image taken on 23 September 2012, assuming that the uncertainty of the parameters and variables is transferable to other images. The uncertainty of the results is evaluated as a function of debris thickness.

Table 3.

Standard deviation of estimated heat fluxes and ranges of input parameters and variables used for the sensitivity analysis. The right column explains for which energy-balance components the input values and parameters are used

The parameters introducing most uncertainty are related to the unknown temperature profile within the debris at the time of satellite image acquisition (Fig. 11a). The highest uncertainty is introduced by estimation of a linear temperature gradient within the debris (expressed by the 0°C depth factor i _d) and the relation of stored heat to conducted heat (expressed by the stored heat factor F). Moreover, thermal conductivity introduces considerable uncertainty, followed by surface temperature T _s. In this context, surface temperature estimation T _s from the satellite image leads to notable uncertainties for low debris thickness.

Fig. 11.

Coefficient of variation of modelled debris thickness as a function of surface temperature T _s (a) for different parameters (0°C depth factor i _d, stored heat factor F, conductivity K, surface temperature T _s) and (b) for different heat fluxes (incoming and outgoing shortwave (S _in, S _out) and longwave radiation (L _in, L _out) and turbulent sensible heat flux H).

Generally, the single energy fluxes present smaller uncertainty than the characteristics of the debris (Fig. 11). For low temperatures, L _out and H are important, whereas H and S _out (through albedo) are important for warm surfaces. For warm surfaces, the uncertainty of H becomes very high, mainly due to the unknown air temperature, followed by surface roughness z ₀ and wind speed u.

Debris thickness

The total uncertainty of model results that arises in estimating physical parameters and taking meteorological variables from reanalysis data was assessed by running the model several times and by varying randomly and simultaneously input variables and parameters. Figure 12 shows that uncertainty is small for cool surfaces but it increases significantly with increasing surface temperature, especially for temperature above ∼10°C for which debris thickness cannot be estimated with reasonable certainty. Even though i _d and F were held constant in Figure 12, the uncertainties are considerable. It seems that above ∼10°C, the method cannot be relied upon to give accurate results, although it should correctly identify the pixel as relatively thick debris. The spread in results is small for relatively cool surfaces, which means that thin debris can be identified with a lower uncertainty.

Fig. 12.

Distribution of the modelled debris thickness (varying the energy fluxes, input parameters and variables) for different surface temperatures as analysed for the ASTER scene taken on 23 September 2012 assuming constant 0°C depth (i _d = 0.5) and stored heat factor (F = 3).

A final debris thickness map was calculated for the entire debris-covered area, including each of the satellite imagery scenes (Fig. 13). The final debris thickness values of every pixel represent the mean of six satellite images, whereas the values of every image are from the median value of the 10³ runs. In their work conducted in the Alps and based on extensive field data, Reference Foster, Brock, Cutler and DiotriFoster and others (2012) concluded that the F12 model is able to properly identify areas with thick debris cover (>50 cm). By analogy and assuming comparable system behaviour, we can assume that debris cover passes this threshold at the medial moraine of the heavily debris-covered tongue of Bara Shigri, by visual inspection of SPOT satellite imagery and the analysis of photographs taken of the glacier tongue. However, the debris thickness values computed here remain below 30 cm in most cells and the model thus seems to underestimate debris thickness. The modelled result reveals relatively thick debris on medial moraines, the eastern edge and the lower tongue. The thick debris on the eastern edge is probably deposited from the steep slopes in the east, with likely greater erosion activity than on the opposite slopes. Thin debris of a few centimetres in depth is observed at the edges of the medial moraines, at the confluence and on the western edge of the tongue.

Fig. 13.

Mean debris thickness from three ASTER and three Landsat 7 scenes.

Surface temperature

The aim of the energy-balance method is to account for all effects determining surface temperature so as to consider exclusively the effect of debris thickness. Figure 7 shows that mean surface temperature is clearly correlated with modelled incoming shortwave radiation and air temperature from reanalysis data. This finding is in line with previous work in the Alps where T _s was shown to correlate well with global radiation on debris-covered areas (Reference Mihalcea, Mayer, Diolaiuti, Lambrecht, Smiraglia and TartariMihalcea and others, 2006). For instance the effect of air temperatures and incoming shortwave radiation on surface temperature should be considered and subtracted by the energy-balance approach. It is controlled by the denominator of Eqn (8) via S _in, L _in and H. For situations with low air temperature and low solar radiation, the denominator is lower, thus resulting in higher debris thickness for a pixel of a certain surface temperature T _s than for a situation with high temperatures and high radiation. The high correlation between T _s and shortwave radiation, as well as T _s and air temperature, thus supports the use of an energy-balance approach to estimate debris thickness based on remote sensing.

The spaceborne surface temperature is provided, with a spatial resolution of 90 and 30 m for ASTER and Landsat 7 imagery, respectively, and every pixel represents a mean temperature over an assumable inhomogeneous surface. The surface temperature at a certain point depends on various factors that are neglected in the model (e.g. shading, aspect, elevation, humidity of debris material, snow cover on debris, cooling during the night and influence of supraglacial ponds and ice cliffs). Where debris is patchy or characterized by a mix of fine material and large boulders, pixels will yield a mean temperature, which may not be representative of adequate heterogeneous surface conditions. However, for large glaciers with continuous debris cover such as Bara Shigri, spaceborne satellite data may well represent surface conditions.

Where the surface temperature of debris is below 0°C, the F12 approach cannot be applied properly, as pixels with negative surface temperatures will be classified as debris-free. The presence of snow or morning shadows as well as the strong cooling during the night may lead to relatively low surface temperatures at the time of overflight and therefore to an underestimation of debris thickness. It is therefore important to consider only satellite images with clearly positive surface temperatures across the debris-covered area.

Parameterizations of radiative fluxes

The computation of downwelling longwave radiation is based on an input value for air temperature. The parameterization of incoming longwave radiation as defined by Reference BrutsaertBrutsaert (1975) performs well and shows higher correlation than the values derived from the approach of Reference BruntBrunt (1932) or the ERA-Interim data. The underestimation of longwave radiation by these two parameterizations for situations with relatively high downwelling longwave radiation is probably due to variations in cloudiness. In the Bara Shigri region and elsewhere in the Indian Himalaya, convective clouds typically form around noon, thereby leading to increased downwelling longwave radiation. Without measurements of cloud cover, it is not possible to model the increase in radiation due to cloudiness. Nevertheless, at the time of satellite image acquisition and the scenes selected, the sky is cloud-free and we can expect longwave and shortwave radiation to be modelled reasonably. The uncertainty in outgoing longwave radiation was higher than for downwelling longwave radiation, which stems from the error in surface temperature as derived from satellite imagery and estimated to be in the order of 1–1.5 K.

The uncertainty introduced through the parameterizations of incoming shortwave radiation is larger than for downwelling longwave radiation. Outgoing shortwave radiation is a function of surface albedo, which adds additional uncertainty to this heat flux. We conclude that for a large glacier like Bara Shigri the spatial variable value for albedo adds uncertainty to the model result and should be assessed in more detail using for example multispectral Landsat Thematic Mapper (TM) satellite data or geological maps.

Turbulent sensible heat flux

Reference Brock, Mihalcea, Kirkbride, Diolauiuti, Cutler and SmiragliaBrock and others (2010) showed that turbulent sensible heat flux is a very large and crucial component of the surface energy balance of a debris-covered glacier and should not be neglected. We suggest that sensible heat flux may be the most uncertain component of the energy balance and constitutes an important source of model error. Uncertainty is mainly due to the uncertainties in surface temperature for surfaces close to 0°C and uncertainties in estimated 2 m air temperature, followed by surface roughness and wind speed. Whereas Reference Foster, Brock, Cutler and DiotriFoster and others (2012) found that model results are sensitive to surface roughness, we found that the uncertainty introduced as a result of unknown surface roughness is relatively small.

For the estimation of turbulent sensible fluxes, a spatial extrapolation of air temperature is needed. Reference Foster, Brock, Cutler and DiotriFoster and others (2012) showed that the relationship between surface temperature and 2 m air temperature is well correlated. This correlation was confirmed by measurements on Haut Glacier d’Arolla and Lirung Glacier. The value of b in Eqn (6) is similar for all glaciers, whereas the value of a is much lower for Haut Glacier d’Arolla. Miage and Lirung glaciers are both extensively debris-covered, whereas Haut Glacier d’Arolla is only partially debris-covered. Based on the results of this study, we assume that the relation of air temperature to surface temperature may depend on the degree of debris cover of the glacier surface. Where a glacier surface is only partially debris-covered, air temperature is less governed by heating of the debris surface, which would be expressed by a lower parameter a in Eqn (6), as measured for Haut Glacier d’Arolla. On this glacier, air temperature might by influenced strongly by katabatic winds associated with low ice surface temperatures, which may explain the relatively low value for a. Since data only exist for one Himalayan glacier, clearly more measurements are needed to determine this relationship at other debris-covered glaciers and to understand how this relation is influenced by for example wind speed and direction, as well as the degree of debris cover.

Moreover, wind speed is an important input value for turbulent sensible heat flux. Reanalysis data are not correlated with local wind speed estimations, calling for improved near-surface wind speed estimations as input values to the model. Reference ZouZou and others (2008) observed that the wind over the debris-covered tongue of Rongbuk Glacier is dominated by thermally developed winds (i.e. valley and mountain winds) or katabatic winds. Local wind systems over glaciers in this region may be poorly influenced by synoptic winds as a result of topographic shielding effects. We therefore suggest that wind speed over Bara Shigri is dominated by local valley mountain winds in the morning on clear-sky days, rather than by synoptic-scale meteorology. As a consequence, the quality of input data for wind speed could be improved considerably by taking statistically derived wind speed instead of reanalysis data, and by selecting clear-sky situations in the morning from September to November. As we did not have in situ wind data from at least one nearby station, the application of the F12 approach can be seen as limited in this particular respect.

Vertical profile of temperature and rate of change in heat stored

Temperature gradients are assumed to be linear in the uppermost debris layers, with a gradient defined by the 0°C depth of 50% of the total debris thickness, which can be confirmed by data found in the literature and the results of Reference Rounce and McKinneyRounce and McKinney (2014). Also the assumption that the fraction F is dependent on d is supported by in situ data from other glaciers. Without these new assumptions, the estimated debris thickness would result in much lower values and consequently be strongly underestimated. However, the unknown local vertical gradient as well as the rate of change in heat stored (ΔD) lead to considerable uncertainties in model results, which cannot be improved unless field data of vertical temperature gradients become available for heavily debris-covered glaciers such as Bara Shigri.

Thermal conductivity

Conductivity depends on several factors such as debris composition, lithology, water saturation and density. Both the increase of water saturation or density of the material will result in increased thermal conductivity. As these characteristics are spatially heterogeneous across the glacier, conductivity may be highly variable in space, and a single value does not represent well the conditions of supraglacial debris. As the model does not account for the spatial variability of thermal conductivity K, we assessed the uncertainty of the model result introduced by the assumption of a spatially constant K. We found that uncertainty is important but is smaller than uncertainties caused by the unknown temperature profile of the debris. Again, field measurements are needed to assess the spatial distribution of thermal conductivity on a Himalayan glacier.