Debris thickness of glaciers in the Everest area (Nepal Himalaya) derived from satellite imagery using a nonlinear energy balance model

Debris thickness is an important characteristic of debris-covered glaciers in the Everest region of the Himalayas. The debris thickness controls the melt rates of the glaciers, which has large implications for hydrologic models, the glaciers’ response to climate change, and the development of glacial lakes. Despite its importance, there is little knowledge of how the debris thickness varies over these glaciers. This paper uses an energy balance model in conjunction with Landsat7 Enhanced Thematic Mapper Plus (ETM+) satellite imagery to derive thermal resistances, which are the debris thickness divided by the thermal conductivity. Model results are reported in terms of debris thickness using an effective thermal conductivity derived from field data. The developed model accounts for the nonlinear temperature gradient in the debris cover to derive reasonable debris thicknesses. Fieldwork performed on Imja– Lhotse Shar Glacier in September 2013 was used to compare to the modeled debris thicknesses. Results indicate that accounting for the nonlinear temperature gradient is crucial. Furthermore, correcting the incoming shortwave radiation term for the effects of topography and resampling to the resolution of the thermal band’s pixel is imperative to deriving reasonable debris thicknesses. Since the topographic correction is important, the model will improve with the quality of the digital elevation model (DEM). The main limitation of this work is the poor resolution (60 m) of the satellite’s thermal band. The derived debris thicknesses are reasonable at this resolution, but trends related to slope and aspect are unable to be modeled on a finer scale. Nonetheless, the study finds this model derives reasonable debris thicknesses on this scale and was applied to other debris-covered glaciers in the Everest region.


Introduction
Debris-covered glaciers are common in the Everest area of the Himalayas.The debris cover has a large impact on the sub-debris ablation rate and hence the evolution of the glacier.A thin debris layer may enhance ablation by reducing the albedo, causing the surface to absorb more radiation compared to clean ice, while a thicker debris layer will insulate the glacier, causing the ablation rate to decrease.The critical thickness at which the debris cover reduces ablation is around 2 cm (Ostrem, 1959;Mattson et al., 1993;Kayastha et al., 2000).Field studies have supported these results, showing that beyond this critical thickness, the melt rate greatly decreases (Nakawo and Young, 1981;Conway and Rasmussen, 2000;Nicholson and Benn, 2006;Reid and Brock, 2010;Reid et al., 2012).The role that debris cover has on the evolution of glaciers in the Everest area is summarized well by Benn et al. (2012).In short, the debris cover increases towards the tongue of the glacier, where the slopes are gentler.The spatial variation of debris cover causes the ablation to be predominately concentrated in areas of thinner debris behind the tongue of the glacier.The differential melting causes the tongue of the glacier to become stagnant and promotes the development of supraglacial lakes.
The sub-debris ablation rate is controlled by the debris thickness, the thermal properties of the debris, and meteorological conditions.The debris thickness may be measured by surveying exposed ice faces (Nicholson and Benn, 2012) or via manual excavation (Reid et al., 2012).Surveying exposed ice faces greatly reduces the amount of labor involved in measuring the debris thickness, but it may not be representative of the entire glacier and is limited to regions with significant differential melting.Due to the labor-intensive nature of this work, few other surveys of debris thickness have been performed in the Everest area (Nakawo et al., 1986).The thermal property associated with describing the debris cover is the effective thermal conductivity.Studies have found the thermal conductivity of debris cover in the Khumbu to range from 0.85 to 1.29 W m −1 K −1 (Conway and Rasmussen, 2000;Nicholson and Benn, 2012).The water content and lithology of the debris cover may partly explain the variation in thermal conductivity as the water content will change the effective thermal conductivity of the debris (Nicholson and Benn, 2006) and the lithology will influence the bulk volumetric heat capacity, which is used to derive the thermal conductivity (Nicholson and Benn, 2012).
In addition to the properties of the debris cover, the meteorological conditions will affect the sub-debris ablation rate.The net solar radiation has been found to be the main source of energy responsible for ablation on debris-covered glaciers (Inoue and Yoshida, 1980;Kayastha et al., 2000;Takeuchi et al., 2000); however, the turbulent heat fluxes are still significant (Brock et al., 2010).Many studies have modeled the energy balance on debris-covered glaciers with varying levels of success (Nakawo and Young, 1982;Nakawo et al., 1999;Han et al., 2006;Nicholson and Benn, 2006;Mihalcea et al., 2008b;Reid and Brock, 2010;Reid et al., 2012;LeJeune et al., 2013).These models integrate meteorological data from automatic weather stations with knowledge of the debris cover to solve for the surface temperature of the debris, which may then be used to calculate the sub-debris ablation rates.These models are limited by their knowledge of how the debris cover varies over the glacier, or they require a great deal of site-specific information.
This has led other studies to use satellite imagery to derive the properties of the debris cover.These studies use surface temperature data from satellite imagery in conjunction with an energy balance model to solve for the thermal resistance, which is the debris thickness divided by the thermal conductivity (Nakawo and Rana, 1999;Nakawo et al., 1999;Suzuki et al., 2007;Zhang et al., 2011).If the thermal conductivity of the debris is known, the model can solve directly for debris thickness (Foster et al., 2012).Mihalcea et al. (2008a) used a different approach by deriving debris thickness from linear relationships between surface temperature and debris thickness for different elevation bands.
One problem associated with the studies that solved for the thermal resistance is that, while the spatial distribution of thermal resistances typically agreed well, the actual values of thermal resistances were significantly lower than those derived from field studies.Suzuki et al. (2007) attributed their low thermal resistances to the mixed-pixel effect, which refers to the pixels in the satellite imagery comprising supraglacial ponds, ice cliffs, and bare ice areas.Nakawo and Rana (1999) also commented on areas with exposed ice cliffs reducing the surface temperature of the pixel, thereby lowering the calculated thermal resistances.Zhang et al. (2011) did not address the low values of thermal resistances, but did attribute the small disagreement between modeled and observed melt rates to the unknown variations in meteorological conditions caused by altitude, aspect, and shading in different areas, as well as the unknown nature of water content in the debris.The mixed-pixel effect and the spatial variation in meteorological conditions may reduce the thermal resistances, but it is unlikely to cause the satellitederived thermal resistances to be 1 or 2 orders of a magnitude lower than those found in the field.Foster et al. (2012) is the first study, to the authors' knowledge, that accurately derives debris thickness from satellite imagery.The model uses a digital elevation model (DEM) generated from an airborne lidar survey and compares the results of a sloped model, which accounts for variations in topography, and a flat model.The sloped model resulted in thicker debris areas when compared to the flat model, but also identified some pixels as having unrealistically high or negative debris thicknesses.These errors occurred in pixels with steep slopes and high surface temperatures and were replaced with the values from the flat model.Unfortunately, the model is difficult to transfer to other glaciers because a great deal of site-specific data were used.Their modifications to their energy balance include the addition of a heat storage term that is a fraction of the ground heat flux and an empirical relationship between the surface temperature and air temperature.The relationship between the surface and air temperature was used to reduce the values of sensible heat flux, since initial results using an instability correction were found to be unrealistically high.
We report a method for deriving the debris thickness of debris-covered glaciers using an energy balance model with Landsat7 Enhanced Thematic Mapper Plus (ETM+) satellite imagery and apply the method in the Everest region of Nepal.The performance of various models is assessed via comparison with field data.First, the use of an approximation factor that accounts for the nonlinear temperature gradient of the debris cover is investigated.This simple nonlinear energy balance model is then used to compare a flat model with a sloped model, which accounts for the variations in topography.The affect of the quality of the DEM is then explored by comparing DEMs of different resolutions.Lastly, the model is applied to other glaciers in the Everest region.

Meteorological data
The energy balance model uses meteorological data from an automatic weather station, Pyramid Station (27.959 • N, 86.813 • E, 5035 m a.s.l), which is located next to the Khumbu Glacier (Fig. 1  speed, incoming shortwave radiation, and incoming longwave radiation from October 2002 to December 2009.All the meteorological data are assumed to be constant over the Khumbu region, with two exceptions.The air temperature was adjusted based on the elevation of each pixel using a lapse rate of 6.5 K km −1 .Furthermore, in the sloped model, the incoming shortwave radiation term was corrected for the effects of topography, altitude, and shading, similar to the methods of Hock and Noetzli (1997).The hillshade tool in ArcGIS was used to determine if any pixels were shaded from the surrounding terrain based on the position of the Sun at the time the satellite images were taken.The flat model also corrects for the effects of shading due to surrounding terrain, but assumes each pixel has a slope and aspect of zero.

Remotely sensed data
Landsat7 ETM+ (hereon referred to as Landsat7) satellite imagery over the same period as the meteorological data were used to derive the thermal resistance of the debris.All clear-sky images from the same period of time that meteorological data are available in the melt season were used.The melt season was defined as 15 May to 15 October, which is the time period where the daily mean temperature in the debris was above freezing (Nicholson, 2005).Twelve Land-sat7 images met this criterion (Table 1).Figure 2 shows an example of the Landsat imagery and the corresponding meteorological data from Pyramid Station from 4 October 2002.
All scenes were downloaded from the NASA Land Processes Distributed Active Archive Center (NASA LP DAAC, 2011).The processing level of the Landsat7 images were all L1T, indicating the images were all geometrically rectified using ground control points (GCPs) from the 2005 Global Land Survey in conjunction with the 90 m global DEM generated by the Shuttle Radar Topographic Mission (SRTM).Land-sat7 satellite imagery comprises eight different bandwidths with various resolutions.The two bands of interest here are the thermal band (Band 6) and the panchromatic band (Band 8).The thermal band has a resolution of 60 m, but is automatically resampled to 30 m and was used to derive surface temperature (Fig. 2) according to NASA (2011).It was atmospherically corrected using the methods described by Coll et al. (2010).The required meteorological data for the MODTRAN 4 model used by Coll et al. (2010) was taken from Pyramid Station.The image-to-image coregistration for Landsat7 is 7.3 m, and the uncertainty of the derived surface temperature data is estimated to be ± 1.0 K (Barsi et al., 2003;Coll et al, 2010Coll et al, , 2012)).The panchromatic band has a horizontal resolution of 15 m and was used to coregister the images.
The high-resolution DEM used in this study was generated by Lamsal et al. (2011) from Advanced Land Observing Satellite (ALOS) Pranchromatic Remote-sensing Instrument for Stereo Mapping (PRISM) images.The generated DEM has a horizontal resolution of 5 m and relative error of ± 4 m.In order to coregister the DEM with the panchromatic band from the Landsat7 imagery, a shaded version of the DEM was generated using the Hillshade tool in ArcGIS 10.3.The swipe visualization tool in PCI Geomatica 2013 showed that the images were properly coregistered without any further processing.The coarser-resolution global DEM used in this study was the ASTER GDEM, which is composed of automatically generated DEMs from the Advanced Spaceborne Emission and Reflection radiometer (ASTER) stereo scenes acquired from 2000 to present (METI/NASA/USGS, 2009).Nuth and Kaab (2011) found the accuracy of the ASTER GDEM to be similar to the validation summary (METI/NASA/USGS, 2009) when applied to debris-covered glaciers in New Zealand.They found the ASTER GDEM to have biases up to 10 m and root mean square error (RMSE) of 5-50 m.The horizontal resolution of the ASTER GDEM has been found to be better than 50 m (Fujisada et al., 2005).The swipe visualization tool in PCI Geomatica 2013 was used with a shaded version of the ASTER GDEM to confirm that the images were properly coregistered.While residual anomalies and artifacts may exist in this experimental/research grade product, it has been used in this study to develop an understanding of how the quality of DEM will affect the thermal resistances.

Study site
Field research was conducted in September 2013 on the debris-covered portion of Imja-Lhotse Shar Glacier (27.901 • N, 86.938 • E, ∼ 5050 m a.s.l.).Imja-Lhotse Shar Glacier refers to two debris-covered glaciers, Imja Glacier (the southeastern component; ∼ 4.5 km to the confluence) and Lhotse Shar Glacier (the northeastern component; ∼ 3.5 km to the confluence), that converge and terminate into Imja Lake (Fig. 1).The third glacier that is present south of Imja Lake is Amphu Glacier, which appears to no longer contribute to Imja-Lhotse Shar Glacier.Imja and Lhotse Shar glaciers are both avalanche-fed and extend from the calving front of Imja Lake (5010 m) up to elevations of 7168 and 8383 m for Imja and Lhotse Shar glaciers, respectively.The thickness of debris cover on Imja-Lhotse Shar Glacier increases towards the terminal moraine of the glacier and is primarily composed of sandy boulder gravel (Hambrey et al., 2008).The debris cover extends up to elevations of 5200 and 5400 m on Imja and Lhotse Shar glaciers, respectively.From the calving front of Imja Lake to the confluence of the glacier, the increase in elevation is less than 50 m, which is consistent with the findings of Quincey et al. (2007) that the tongue of the glacier is relatively stagnant with a surface gradient less than 2 • .The debris cover has a hummocky terrain with melt ponds and exposed ice faces scattered throughout.

Debris measurements
Debris thermistors (TR-52 ThermoRecorder, T&D Corporation) were installed at four locations (referred to as LT1, LT2, LT3, and LT4) at depths of 0, 5, 10, 15, 20, and 30 cm, and at the debris-ice interface.The debris thickness was 31, 47, 36, and 40 cm for LT1, LT2, LT3, and LT4, respectively.Holes were excavated to the debris-ice interface and, as the thermistors were installed, the debris was replaced in its original position as well as possible.The thermistors recorded temperature at hourly intervals from 13 September to 24 September.The first 48 h of data for each thermistor was discarded to allow the thermistors to equilibrate with the debris.One of the surface thermistors malfunctioned on 23 September, so the data from this thermistor beyond this date were discarded.
Debris thickness measurements were performed at 25 locations and were concentrated in one melt basin (27.901 • N, 86.938 • E, 5045-5055 m a.s.l.) that appeared to be formed by differential melting and backwasting (Fig. 1).The melt basin was selected as the focus area of this study because it appeared to be representative of the hummocky terrain on Imja-Lhotse Shar Glacier and was relatively easy to access.The melt basin was approximately 120 m long and 60 m wide with a topographic low in the center of the basin (5045 m a.s.l.) and a topographic high on the perimeter of the basin (5055 m a.s.l.).The elevation of the melt basin was only 10-20 m higher than the elevation of Pyramid Station.
To the best ability of the authors, the measurements were performed randomly throughout the melt basin.Measurements were conducted via manual excavation using a tape measure at 23 of the 25 sites.This process involved digging holes to the ice surface and measuring the perpendicular distance from the ice surface to the surface of the debris.The other two sites were the debris on top of an ice face, which was measured using a laser range finder (TruPulse 360B) because the ice face could not be accessed safely on foot.Twelve debris thickness measurements were also performed outside of the melt basin to understand if the melt basin was representative of the debris-covered glacier.More debris thickness measurements were unable to be made due to time and labor restraints.Furthermore, the maximum depth of excavation was 1 m because further excavation was too physically demanding.

Energy balance model
The energy balance model used in this study is a steady-state surface energy balance for the debris cover similar to that developed by Nakawo and Young (1982).
where R n is the net radiation flux, H is the sensible heat flux, LE is the latent heat flux (assumed to be zero), and Q c is the ground heat flux (all in W m −2 ).The net radiation flux includes the shortwave radiation flux and the longwave radiation flux.
where S↓ is the incoming shortwave radiation (W m −2 ), α is the albedo (0.30), ε is the emissivity assumed to be 0.95 (Nicholson and Benn, 2006), L↓ is the incoming longwave radiation (W m −2 ), σ is the Stefan-Boltzmann constant (5.67 × 10 −8 W m −2 K −4 ), and T S is the surface temperature (K).For the sloped model, incoming shortwave radiation was corrected for the effects of topography, altitude, and shading similar to the methods of Hock and Noetzli (1997).The flat model assumes that each pixel has a slope and aspect of zero degrees.The incoming longwave radiation and surface albedo were assumed to be constant over the entire debris cover.The albedo used in this study (0.30) was the average albedo of the debris cover on Ngozumpa Glacier (Nicholson and Benn, 2012).The latent heat flux was assumed to be zero based on the assumption that the debris cover is dry.The sensible heat flux was tested using two approaches: (1) assuming a neutral atmosphere (Nicholson and Benn, 2006) and ( 2) correcting for the unstable atmosphere (Reid and Brock, 2010).It is likely that high debris temperatures and low air temperature would cause an unstable atmosphere above the surface of the debris; however, the incorporation of the stability correction caused the sensible heat fluxes to be unrealistically high, similar to the initial results of Foster et al. (2012).Therefore, this study calculates the sensible heat flux assuming a neutral atmosphere according to Nicholson and Benn (2006): where where ρ air is the density of air (1.29 kg m −3 ), P is the atmospheric pressure computed using the barometric pressure formula, P 0 is the atmospheric pressure at sea level (101 325 Pa), c is the specific heat capacity of air (1010 J kg −1 K −1 ), A is the dimensionless transfer coefficient, u is the wind speed at Pyramid Station, T air is the air temperature 2 m above the surface calculated using a lapse rate of 6.5 K km −1 , k vk is von Karman's constant (0.41), z is the height of meteorological measurements (2 m), and z 0 is the surface roughness length (assume z 0 = 0.016).
The ground heat flux is different for the linear and the nonlinear models.
Linear model : Nonlinear model : G ratio is the nonlinear approximation factor, k eff is the effective thermal conductivity (W m −1 K −1 ), and d is the debris thickness (m).In this study, the modeled results will derive debris thickness assuming an effective thermal conductivity from Sect.4.1.However, this model can also be applied to solve for the thermal resistance in the event that the effective thermal conductivity is unknown.The linear model assumes the temperature gradient in the debris is linear from the surface temperature to the debris-ice interface, which is assumed to be at 273.15 K.At the time that Landsat7 images are acquired (10:15 local time), this linear assumption is not accurate.G ratio is used to approximate the nonlinear temperature variation in the debris.A linear temperature gradient in the upper 10 cm of the debris is used to make this approximation.G ratio is therefore defined as the ratio of the linear temperature gradient in the upper 10 cm of the debris to the linear temperature gradient throughout the entire debris.
where T 0.1 m is the temperature 10 cm below the surface of the debris and T d is the temperature at the debris-ice interface (273.15K).As the Landsat7 images are acquired at 10:15 and the thermistors recorded hourly temperatures, the temperatures in the debris at 10:15 were computed by linearly interpolating between 10:00 and 11:00.These interpolated temperatures were used to compute G ratio .
4 Field results

Thermal conductivity
The effective thermal conductivity, k eff , of the debris cover was computed following the methods in Conway and Rasmussen (2000) assuming a density (ρ = 2700 kg m −3 ) and a specific heat capacity (c = 750 J kg −1 K −1 ) of rock.The average effective thermal conductivity was calculated to be 0.96 (± 0.33) W m −1 K −1 .This effective thermal conductivity agrees well with other thermal conductivities computed in this area, which range from 0.85 to 1.29 W m −1 K −1 (Conway and Rasmussen, 2000;Nicholson and Benn, 2012).
The thermal conductivity was greatly influenced by depth as the average values above and below 10 cm were 0.60 and 1.20 W m −1 K −1 , respectively.The drastic difference in thermal conductivity above and below 10 cm is likely due to the amount of water content in the debris.Nicholson and Benn (2006) found the thermal conductivity of wet debris (assuming the pores were saturated with water) to be 2-3 times larger than dry debris.These results indicate that the top 10 cm of the debris is dry, while 15 cm and lower is wet.This is consistent with observations in the field, where approximately the top 10 cm of the debris was dry and below this depth the debris was wet.The observed moisture in the debris was mainly the wetted surface of the grains.

Debris thickness
The debris thickness in this melt basin ranged from bare ice (0 cm) to depths greater than 1 m (Fig. 3).The average debris thickness, assuming a maximum thickness of 1 m, was 0.42 (± 0.29) m.These debris thicknesses are consistent with the debris thickness of other debris-covered glaciers in the Everest region (Nakawo et al., 1986;Nicholson and Benn, 2012).Debris thicknesses greater than 1 m were found in the bottom of these melt basins where the debris had likely accumulated over time due to differential melting and backwasting of the debris cover.Areas of thin debris cover were located on the slopes of the melt basin.These trends are identical to those found by Nicholson and Benn (2012) and were also observed at the 12 other sites where debris thickness was measured outside of the melt basin.There did not appear to be any trends in debris thickness with respect to aspect.Ideally, debris thickness would be sampled over the entire debriscovered glacier to derive a debris thickness map that could be used to validate the modeled results.As this was not feasible due to restraints on time and labor, the modeled results from within this melt basin and the melt basin's adjacent cells will constitute the focus area of the satellite imagery that will be compared to the measured debris thicknesses to qualitatively assess the reasonableness of the modeled results.

Nonlinear approximation factor, G ratio
G ratio was computed from all the temperature profiles based on the interpolated temperatures at 10:15. Figure 4 shows the temperature profiles at site LT3 and a schematic of the temperature gradients used to compute G ratio .Figure 4 shows the temperature at the debris-ice interface was 0 • C for each day at site LT3 at 10:15.The temperature remained at 0 • C for all the sites throughout the entire study period.The average value of G ratio for the melt basin was 2.7 (± 0.4).The depth used to approximate the nonlinear temperature gradient was found to greatly influence G ratio , where G ratio decreased as the depth increased (Fig. 5).This relationship was expected because as the depth used to calculate G ratio approaches the thickness of the debris, G ratio will approach a value of 1. Ideally, the nonlinear temperature gradient would be approximated by a linear temperature gradient in the upper 1 cm of the debris or less.However, these measurements could not practically be performed in the field.Figure 5 shows the values of G ratio derived using depths of 5 and 10 cm (2.9 ± 0.8 and 2.7 ± 0.4, respectively) are similar; however, the standard deviation of those derived from 5 cm is much larger than those derived from 10 cm.As the two values were similar, G ratio derived from a depth of 10 cm was used in this study.
The other trend that may be expected is for G ratio to increase as the surface temperature increases.However, Fig. 5 shows that there is little correlation (R 2 = 0.29) between surface temperature and G ratio based on the data from all four sites.When the G ratio values are removed from site LT4, there appears to be a stronger relationship (R 2 = 0.74), but there is no physical justification for removing these data.Therefore, this study uses an average value of G ratio of 2.7.Conway and Rasmussen (2000) is the only other study, to the authors knowledge, that has measured temperature profiles in the Everest area with a small enough spacing (maximum 10 cm) between thermistors to compute G ratio .The values derived from their temperature profiles at Everest Base Camp on 21-23 May 1999 were 3.0, 2.6, and 2.6, respectively.This good agreement lends confidence to the use of G ratio in the Everest area.Foster et al. (2012) noted that the variation in surface temperature is very sensitive to small changes in debris thickness less than 0.5 m, while its change is very gradual for thicker debris.This caused them to have a high level of uncertainty in mapping debris thickness greater than 0.5 m.The debris thickness of the four sites used to derive G ratio ranged from 0.31 to 0.47 m, which approaches the upper limit of 0.5 m.This lends confidence to the use of G ratio for mapping debris thickness greater than 0.31 m.However, for regions with much smaller debris thicknesses, e.g., Miage Glacier (Reid et al., 2012), where a significant amount of the debris is less than 10 cm thick, it is possible that G ratio will be different.Future work should determine the value of G ratio for thin debris layers in addition to determining how G ratio varies over the melt season and how G ratio is influenced by temperature.

Nonlinear sloped model -5 m pixels
The nonlinear sloped model accounts for the nonlinear temperature gradient in the debris cover using the G ratio approximation factor and corrects the incoming solar radiation at each pixel for the effects of topography and shading.Initially, the model was applied and solved for every individual pixel at a resolution of 5 m.The average debris thickness in the focus area was 0.27 (± 0.19) m, with 17 of the 1080 pixels being undefined.All of the pixels that were undefined were pixels that were shaded by the surrounding terrain, causing their net radiation to be very low.These debris thicknesses agree fairly well with the debris thickness observed on Imja-Lhotse Shar Glacier (Fig. 3), although there are fewer pixels that have very high debris thickness (> 0.5 m).This is likely due to the lack of variation in surface temperature when the debris is greater than 0.5 m as Foster et al. (2012) discussed.
One problem with the modeled results was there were clear trends with respect to slope and aspect (Table 2).Table 2 shows the average debris thickness increases with increasing slope, which is opposite of the trend observed in the field.Furthermore, there was a clear trend with respect to aspect where east-facing pixels had the smallest debris thickness, followed by north-and south-facing pixels, while westfacing pixels had the highest values of debris thickness.This trend was not observed in the field.These trends occur from the corrections that are applied to the sloped model to account for the topography and surrounding terrain.Pixels that have steeper slopes are angled such that the incoming solar radiation is reduced.This reduces the net radiation, which lowers the net energy flux (net radiation and turbulent heat fluxes) used to derive the debris thickness.In some cases, the net energy flux is negative, which causes a pixel to be undefined.Otherwise, the net energy flux is positive, but small, which results in large values of debris thickness.A minimum threshold for the net energy flux of 10 W m −2 was set, such that unrealistically high debris thicknesses would be classified as undefined.The same changes in incoming solar radiation occur with respect to aspect, where east-facing pixels are oriented towards the Sun, thereby increasing the amount of incoming solar radiation compared to west-facing pixels where the opposite occurs.Therefore, despite the modeled debris thickness yielding reasonable results, the model is incapable of capturing the fine local variations.  of the thermal band.The thermal band has a resolution of 60 m (automatically resampled to 30 m), which causes the surface temperatures over the 60 m pixel to be combined.This is referred to as the mixed-pixel effect.Conventionally, the mixed-pixel effect has been used to explain how bare ice faces reduce the surface temperature of the pixel, causing the derived debris thickness or thermal resistance to be low.While this may be true, the mixed-pixel effect also explains how local variations in surface temperature are not properly accounted for.Table 2 reveals that the average surface temperature in the focus area is almost constant and does not vary with respect to slope or aspect.A higher resolution thermal band would show higher surface temperatures on south-and east-facing slopes, since their orientation allows them to receive more incoming shortwave radiation throughout the morning.North-and west-facing slopes that do not receive radiation would have lower surface temperatures, which would reduce the derived debris thickness.The mixed-pixel effect explains why the modeled results agree well with the average measured values, but do not capture the local variations.One way to overcome this problem is to apply the corrections to the incoming solar radiation and then resample these values of incoming solar radiation to be consistent with the surface temperature pixels (30 m).Ideally, the values of incoming solar radiation would be resampled to a 60 m resolution, since this is the resolution of the thermal band.However, this cannot be performed since the Land-sat7 L1T product automatically resamples the thermal band to 30 m.The rest of the results shown in this study compute the incoming solar radiation, air temperature, and pressure for each pixel and then average these values to a resolution of 30 m.The debris thickness is then derived at this 30 m resolution.

Nonlinear sloped model -30 m pixels
The debris thickness derived using the nonlinear sloped model with 30 m pixels on Imja-Lhotse Shar Glacier are shown in Fig. 6a.Only two pixels on the entire glacier are classified as undefined.These pixels were almost completely shaded by surrounding terrain and had predominately northand west-facing aspects.This is a substantial improvement to the model's performance and shows the importance of resampling the corrected incoming solar radiation.The modeled results also reveal that the debris is thick on the terminal moraine, with most pixels having values greater than 0.40 m.
Behind the calving front of Imja Lake the debris is also thick, especially in the center of Imja-Lhotse Shar Glacier.This region of thick debris in the middle of Imja-Lhotse Shar Glacier is likely a result of the slope being gentler, thereby  allowing debris to accumulate.Furthermore, this portion of the glacier receives debris from both Imja and Lhotse Shar glaciers since it is downstream for their confluence.The latter likely explains why the debris is thinner towards the lateral moraines and thicker in the center.Further up glacier on Imja and Lhotse Shar glaciers, the debris cover thins, which is expected because the slopes are steeper and they are approach-ing areas of clean ice on each glacier.These trends, in which the debris is thicker on the moraine and thins upglacier, are consistent with debris thickness surveys performed on the Khumbu Glacier (Nakawo et al., 1986) and Ngozumpa Glacier (Nicholson and Benn, 2012).
The modeled debris thickness in the focus area had an average debris thickness of 0.29 (± 0.13) m.These debris thicknesses agree fairly well with those measured on the glacier (Fig. 3) with the exception of the high debris thicknesses (> 0.5 m) being underrepresented as previously discussed.The lack of ability in modeling debris thickness greater than 0.5 m is problematic for estimates of total volume or mass of debris on the glacier.However, if the debris thickness is used to estimate ablation rates, then this limit of 0.5 m is not a problem because there is little change in the ablation rates for debris greater than 0.5 m thick.If a limit of 0.5 m is applied to the field measurements such that any debris thickness greater than 0.5 m is given a value of 0.5 m, then the average debris thickness would be 0.34 (± 0.18) m.If the same limit is applied to the modeled debris thickness, the average is 0.28 (± 0.11) m, which agrees quite well.

Linear sloped model -30 m
The debris thicknesses derived using the linear sloped model (Fig. 6b) are significantly smaller than those using the nonlinear model (Fig. 6a).The same trends indicating thicker debris on the moraine and behind the calving front with thinner debris upglacier are apparent.The only difference is the nonlinear approximation factor is not included in the linear model, which causes all the debris thicknesses to be 2.7 times smaller than the nonlinear model.The average debris thickness in the focus area is 0.11 (± 0.05) m.These results indicate that the linear model severely underestimates the debris thickness and show that it is critical to account for the nonlinear temperature profile.

Nonlinear flat model -30 m
One method to fill in the undefined pixels from the sloped model is to use a flat model (Foster et al., 2012).Figure 6c shows the debris thickness map derived from the nonlinear flat model.The flat model captures the trends of greater debris thickness on the terminal moraine and behind the calving front with thinner debris upglacier.However, these trends are not as prominent as in the sloped model, and the debris thicknesses appear to be significantly smaller in comparison to the sloped model.The focus area reveals the flat model underestimates the debris thickness on the glacier as it has an average value of 0.19 (± 0.05) m.These debris thicknesses are greater than the linear model, but much lower than the nonlinear sloped model.Table 2 reveals that the large difference between the sloped model and the flat model is caused by differences in the net radiation.The flat model only corrects for pixels that are shaded and does not correct for the effects of topography, thereby causing the incoming solar radiation to be overestimated in most cases (the exception is for east-facing pixels, which are underestimated).The overestimation in incoming solar radiation causes the net radiation to increase, thereby decreasing the modeled debris thickness.This is important because, if the flat model values are used to fill in the undefined pixels in the sloped model, one must un-derstand that the debris thicknesses will be lower.A preferable alternative may be to use the average debris thickness from the sloped model similar to Table 2 based on the average slope and aspect of the 30 m pixel.However, resampling the incoming solar radiation to 30 m eliminates almost all of the undefined pixels.

Importance of DEM resolution
The requirement of a high-resolution DEM limits the ability to transfer these models to other regions where these data may not be available.The ASTER GDEM was used to assess the importance of the DEM resolution.Figure 6d shows the debris thickness map derived using the nonlinear sloped model with the ASTER GDEM.The average debris thickness in the focus area was 0.18 (± 0.02) m.These debris thicknesses also underestimate the measured debris thicknesses and are very similar to the values derived using the flat model.This underestimation likely occurs because the 30 m resolution of the DEM is too poor to capture local variations in the surface topography of the glacier.Hence, the variations in incoming solar radiation are not captured as well as those using a high-resolution DEM.However, the debris thickness map using the ASTER GDEM with the nonlinear sloped model appears to capture the trends associated with the debris thickness better than the flat model, despite the findings in the focus area.The thick debris region behind the calving front in the center of Imja-Lhotse Shar Glacier (Fig. 6d) more closely resembles the nonlinear sloped model (Fig. 6a) than the nonlinear flat model (Fig. 6c).Therefore, while a high-resolution DEM will yield the best results, we recommend the ASTER GDEM should be used instead of a flat model to get an estimate of the debris thickness in an unknown area.However, one must use caution with these estimates as the analysis here shows the debris thickness are likely to be significantly underestimated.

Debris thickness for glaciers in the Everest region
The ASTER GDEM was used to derive debris thickness maps for the debris-covered glaciers in the Everest region (Fig. 7).The Ngozumpa and Khumbu glaciers have been outlined using the Glacier Land Ice Measurements from Space (GLIMS) database to compare the derived thermal resistances with previous debris thickness measurements (Nakawo et al., 1986;Nicholson and Benn, 2012).The results on the Khumbu Glacier show relatively good agreement with the debris thickness map generated by Nakawo et al. (1986).The debris is thicker close to the terminal moraine and thins upglacier.Furthermore, the modeled debris thickness shows the central zones of thin debris, with thicker debris towards the lateral moraines, which Nakawo et al. (1986) observed.The debris thickness is slightly underestimated once again, likely resulting from the use of the ASTER GDEM as opposed to a high-resolution DEM as The Cryosphere, 8, 1317-1329, 2014 previously discussed.There are similar trends with respect to debris thickness on Ngozumpa Glacier as well, where the terminal moraine has thicker debris and the debris thins upglacier as measured by Nicholson and Benn (2012).Furthermore, Spillway Lake has been growing near the terminal moraine of Ngozumpa Glacier and can clearly be seen as the region of low debris thickness amongst the thicker debris in the terminal moraine.Once again, the debris thickness is underestimated compared to Nicholson and Benn (2012) due to the poor resolution of the DEM.Future work should seek to quantify how much the debris thickness is underestimated and develop a method that is able to accurately quantify the debris thickness from a poor-resolution DEM.This would require a detailed debris thickness survey such that the resulting debris thickness maps could be properly validated.

Sensitivity analysis
The model developed in this study relies heavily upon meteorological inputs, measured parameters, and assumed values associated with the debris cover.The meteorological inputs are subject to instrument error and may not be directly transferable from the site of the automatic weather station to the debris-covered glaciers.The particular meteorological parameters of interest are wind speed (u), air temperature (T air ), and incoming solar radiation (S in ).The parameters associated with the debris cover that may affect results are the surface roughness length (z 0 ), albedo (α), and effective thermal conductivity (K eff ).In addition, there is uncertainty associated with the nonlinear approximation factor (G ratio ) and the derived surface temperature (T s ) from the Landsat7 thermal band.A sensitivity analysis with respect to these parameters was performed on the focus area (Table 3) to identify how each affects the modeled debris thickness (d).
The sensitivity analysis reveals the model is most sensitive to the wind speed and incoming solar radiation.The assumption that the incoming solar radiation at Pyramid Station is the same as the incoming solar radiation over the study area is a reasonable assumption, since all the images used in this study had completely clear skies over both Pyramid Station and the study area.Therefore, it is unlikely that the incoming solar radiation would be ± 10% different.The model's sensitivity to wind speed is concerning because the automatic weather station is located 10 km away from the glacier.The model assumes the wind is the same on the glacier as it is at the automatic weather station, but no data on this exist.Future work should investigate the relationship between meteorological parameters at Pyramid Station and on the debriscovered glaciers.
The model is moderately sensitive to the surface temperature and air temperature on the glacier.Both these parameters affect the temperature gradient in the sensible heat flux term; hence an increase in the gradient will result in an increase in the sensible heat flux, thereby reducing the net energy flux and increasing the modeled debris thickness.Furthermore, the increased sensitivity to surface temperature is quite interesting when one considers the mixed-pixel effect.The mixed-pixel effect causes the actual surface temperature to be reduced due to surrounding bare ice faces.This reduction in surface temperature directly reduces the derived debris thickness, but also reduces the temperature gradient in the sensible heat flux term, thereby further reducing the debris thickness.The model's sensitivity to the surface temperature shows how much the mixed-pixel effect can alter the derived thermal resistances.
With respect to parameters associated with the debris cover, the model is very sensitive to the value of effective thermal conductivity.Therefore, it is important that the effective thermal conductivity is measured at multiple sites for various debris thicknesses throughout the debris-covered glacier such that the modeled debris thickness maps are accurate.In the event that there is large uncertainty in the effective thermal conductivity, the model can be solved for thermal resistance instead.The assumption of a constant albedo over the debris-covered glacier is another limitation of this model, especially since the model is sensitive to albedo.Methods exist to use other Landsat7 bands to estimate albedo (Liang, 2001); however, the authors had no way of validating these results.Therefore, the average value of albedo of 0.30 determined by a previous study on Ngozumpa Glacier (Nicholson and Benn, 2012) was used in this study.Future work should seek to derive the albedo from satellite imagery in conjunction with field measurements to validate these results.Lastly, the model was least sensitive to changes in G ratio , lending confidence to its use in this study.

Conclusion
The model described in this study allows the debris cover or thermal resistance on debris-covered glaciers to be derived from Landsat7 satellite imagery in conjunction with meteorological data from an automatic weather station nearby.The model was applied to glaciers in the Everest region and the resulting debris thicknesses were compared to field measurements.The model accounts for the nonlinear temperature gradient in the debris through the use of a nonlinear approximation factor, which yields reasonable results.Furthermore, the use of a high-resolution DEM greatly improves the results of the modeled debris thicknesses.In the event that a highresolution DEM is not available, the authors recommend using a lower-resolution global DEM to estimate debris thicknesses as opposed to using a flat model.It is important to use caution with these results if a lower-resolution DEM is used since the derived debris thicknesses are likely to be underestimated.A sensitivity analysis reveals that the model is most sensitive to wind speed and incoming shortwave radiation.With respect to debris cover parameters, the model is very sensitive to albedo and the effective thermal conductivity.In the event that the thermal conductivity is uncertain, the thermal resistances may be derived instead.Future work should explore how the meteorological conditions vary spatially and seek to derive the albedo from satellite imagery.The main limitation of this work is the poor resolution of the Landsat7 thermal band.The derived debris thickness must be resampled to the resolution of the thermal band before they are used in melt models or other applications.

Figure 1 .
Figure 1.Panchromatic band from Landsat7 on 4 October 2002 showing Imja Lake, the focus area on Imja-Lhotse Shar Glacier, and Pyramid Station amongst debris-covered glaciers in the Everest region.

Figure 3 .
Figure 3. Debris thickness measurements on Imja-Lhotse Shar Glacier inside and outside of the melt basin.

Figure 4 .
Figure 4. Temperature profiles at site LT3 at 10:15 with the temperature gradients used to compute G ratio identified.

Figure 5 .
Figure 5. G ratio as a function of depth (left) and corresponding G ratio values at 10 cm depth for each sensor (right).

Figure 6 .
Figure 6.Debris thickness (m) with 30 m pixels using the (A) nonlinear sloped model, (B) linear sloped model, (C) nonlinear flat model, and (D) using the ASTER GDEM with the nonlinear sloped model.

D
. R. Rounce and D. C. McKinney: Debris thickness of glaciers in the Everest area (Nepal Himalaya)

Figure 7 .
Figure 7. Debris thickness (m) using the ASTER GDEM and nonlinear sloped model for debris-covered glaciers in the Everest region with Ngozumpa, Khumbu, and Imja-Lhotse Shar glaciers, left to right, respectively, outlined using the GLIMS database.

Table 1 .
Overview of satellite imagery used in this study.

Table 3 .
Sensitivity analysis for select meteorological and model parameters.Bold numbers showing major changes to average debris thickness and dashes indicating the baseline parameter was used.