Debris-covered glacier energy balance model for Imja–Lhotse Shar Glacier in the Everest region of Nepal

Debris thickness plays an important role in regulating ablation rates on debris-covered glaciers as well as controlling the likely size and location of supraglacial lakes. Despite its importance, lack of knowledge about debris properties and associated energy fluxes prevents the robust inclusion of the effects of a debris layer into most glacier surface energy balance models. This study combines fieldwork with a debris-covered glacier energy balance model to estimate debris temperatures and ablation rates on Imja–Lhotse Shar Glacier located in the Everest region of Nepal. The debris properties that significantly influence the energy balance model are the thermal conductivity, albedo, and surface roughness. Fieldwork was conducted to measure thermal conductivity and a method was developed using Structure from Motion to estimate surface roughness. Debris temperatures measured during the 2014 melt season were used to calibrate and validate a debris-covered glacier energy balance model by optimizing the albedo, thermal conductivity, and surface roughness at 10 debris-covered sites. Furthermore, three methods for estimating the latent heat flux were investigated. Model calibration and validation found the three methods had similar performance; however, comparison of modeled and measured ablation rates revealed that assuming the latent heat flux is zero may overestimate ablation. Results also suggest that where debris moisture is unknown, measurements of the relative humidity or precipitation may be used to estimate wet debris periods, i.e., when the latent heat flux is non-zero. The effect of temporal resolution on the model was also assessed and results showed that both 6 h data and daily average data slightly underestimate debris temperatures and ablation rates; thus these should only be used to estimate rough ablation rates when no other data are available.


Introduction
Debris-covered glaciers are commonly found in the Everest region of Nepal and have important implications with regard to glacier melt and the development of glacial lakes.It is well understood that a thick layer of debris (i.e., > several centimeters) insulates the underlying ice, while a thin layer of debris (i.e., < several centimeters) may enhance ablation (Østrem, 1959;Nakawo and Young, 1981;Nicholson and Benn, 2006;Reid et al., 2012).Spatial variations in debris thickness, particularly where the debris layer thins upglacier, can also lead to reverse topographic and ablation gradients, glacier stagnation and, ultimately, the development of lakes (Benn et al., 2012).These glacial lakes and their surrounding bare ice faces also play a crucial role in glacier melt as they typically have ablation rates that are orders of magnitude greater than those observed beneath debris cover (Benn et al., 2012).The importance of debris thickness has led many studies to develop models in conjunction with knowledge of the surface temperature to derive debris thickness (Zhang et al., 2011;Foster et al., 2012;Fujita and Sakai, 2014;Rounce and McKinney, 2014).With knowledge of debris thickness, energy balance models may be used to model debris surface temperature, sub-debris ablation rate, and/or runoff downstream (Nicholson and Benn, 2006;Reid et al., 2012;Collier et al., 2014;Fujita and Sakai, 2014).The main factors affecting the performance of these models are the amount of knowledge of the debris properties, the spatial and temporal resolution of the meteorological data, and the assumptions/complexity of the model.
The properties of the debris typically required in debriscovered glacier energy balance models are the albedo, thermal conductivity, and surface roughness.The albedo of debris on glaciers in the Everest region has been found to range from 0.1 to 0.6 (Inoue and Yoshida, 1980;Kayastha et al., 2000;Nicholson and Benn, 2012;Lejeune et al., 2013).Specifically, Nicholson and Benn (2012) reported that 62 % of measured values ranged between 0.1 and 0.3.Similarly, Kayastha et al. (2000) showed that most values fall between 0.2 and 0.4.The thermal conductivity of debris in the Everest region has been found to range from 0.60 to 1.29 W m −1 K −1 (Conway and Rasmussen, 2000;Nicholson and Benn, 2012;Rounce and McKinney, 2014).The surface roughness, z 0 , is arguably the most difficult parameter to measure as it requires an eddy covariance instrument, horizontal wind speed measurements at multiple heights above the surface, or detailed microtopographic measurements (Brock et al., 2006).In the Everest region, Inoue and Yoshida (1980) estimated z 0 to be 0.0035 and 0.060 m for two sites, one consisting of small schist and bare ice and another comprising mainly large granite, respectively.Takeuchi et al. (2000) estimated a similar value of z 0 on the Khumbu Glacier of 0.0063 m.On Miage Glacier in the Italian Alps, Brock et al. (2010) measured z 0 to be 0.016 m on a debris-covered glacier.
In addition to the properties of the debris, the amount and source of meteorological data available may also greatly influence the model performance.In particular, knowledge re-lated to the latent heat flux on debris-covered glaciers is very limited.This has led previous studies to assume the surface is dry (Foster et al., 2012;Lejeune et al., 2013;Rounce and McKinney, 2014), assume it is dry unless the surface relative humidity was 100 % (Reid and Brock, 2010;Reid et al., 2012;Fyffe et al., 2014), assume a relationship between debris thickness and wetness (Fujita and Sakai, 2014), or use a reservoir approach to model the moisture in the debris (Collier et al., 2014).Collier et al. (2014) suggested that if the atmospheric surface layer is well mixed, then the water vapor partial pressure between the surface and the air may be assumed to be constant, thereby resulting in a latent heat flux based on the vapor pressure gradient.Fyffe et al. (2014) also commented that the lower portion of the debris near the ice interface was observed to be saturated indicating that there may be evaporation and condensation occurring within the debris, albeit small, even when the surface relative humidity is less than 100 %.The lack of knowledge of the moisture in the debris and at its surface makes it difficult to accurately model the latent heat flux term.These problems are further exacerbated in data scarce regions where automatic weather stations are not available.In these situations, reanalysis data sets must be used for all the required meteorological data (Fujita and Sakai, 2014).
This study develops a method to estimate z 0 using a microtopographic method in conjunction with Structure from Motion (SfM) photogrammetry techniques (Westoby et al., 2012).The z 0 values are used with measured values of thermal conductivity, and previously reported values of albedo to calibrate a debris-covered glacier energy balance model on Imja-Lhotse Shar Glacier.Temperature sensors installed at various depths at debris-covered sites were operated from May to November 2014 on Imja-Lhotse Shar Glacier and are used for calibration and validation of the model.Various methods for estimating the latent heat flux are investigated.Furthermore, sub-debris ablation rates are compared to ablation stake measurements to assess model performance; and the effects of temporal resolution are investigated.

Field data
Field research was conducted on the debris-covered portion of Imja-Lhotse Shar Glacier (27.901 • N, 86.938 • E; ∼ 5050 m a.s.l., Fig. 1) from May to November 2014.Imja-Lhotse Shar Glacier refers to both Imja Glacier and Lhotse Shar Glacier, which are avalanche-fed debris-covered glaciers that converge and terminate into Imja Lake.The debris primarily consists of sandy boulder gravel (Hambrey et al., 2008), with the debris thickness increasing towards the terminal moraine.A more detailed description of the glacier may be found in Rounce and McKinney (2014).The field expedition focused on 19 sites on the debris-covered portion of the glacier to determine how debris thickness and topogra- phy affect ablation rates.Four sites were used to analyze the surface roughness through the use of SfM and are referred to as Sites A-D (Fig. 2).These sites were selected to represent various grain sizes and mixes of debris that were observed on Imja-Lhotse Shar Glacier.Site A was relatively homogenous with the majority of debris being cobble and gravel ranging in size from 0.05 to 0.25 m.Site B comprised similar cobbles typically ranging in size from 0.15 to 0.25 m, with larger boulders lying on top of the cobble of up to 1 m.Site C had the finest debris, which primarily consisted of fines and gravel, with some cobbles on the surface up to 0.15 m in size.Lastly, Site D was the most heterogeneous site with boulders ranging up to 0.40 m overlying a surface of cobble of similar size to Site A mixed with the fine and gravel material found in Site C. Temperature sensors and ablation stakes were installed at 20 other sites; however, data could only be retrieved from 15 of the sites (Table 1) as sensors were lost due to large changes in the topography and some of the loggers failed.Sites 4-14 were located in a single area that appeared to have developed from differential backwasting over the years.This was the same focus area as described in Rounce and McKinney (2014) and was selected because it appeared to be representative of the debris-covered terrain on Imja-Lhotse Shar Glacier and was accessible.Sites 15-20 were located outside of the focus area in an adjacent melt basin to determine if the focus area was representative of other debris-covered areas.At each site, the debris thickness was determined following the methods described in Rounce and McKinney (2014)  with the exception of Site 4, where the debris thickness was greater than 1 m and therefore was estimated assuming a linear temperature profile from the mean temperatures over the study period similar to the extrapolation used in Nicholson and Benn (2012).The debris thickness of these sites ranged from 0.07 m to greater than 1 m.A debris thickness of 1 m was considered the maximum due to labor constraints.The slope was also approximated by measuring two points, one 0.5 m uphill from the site and the other 0.5 m downhill, using a total station (Sokkia SET520, ±2.6 mm 100 m −1 ).The slope at each site ranged from 17 to 37 • .The aspect of each site was measured using a compass (Table 1).
Temperature sensors (TR-42 ThermoRecorder, T&D Corporation) were installed and successfully retrieved at 10 sites.These sensors recorded data every half hour from 19 May to 9 November 2014.Each of the 10 sites had a sensor at its surface, which was considered to be installed 1 cm into the debris since debris was placed on top of the sensor.Sites 4, 11, 13, and 14 also had temperature sensors installed within the debris to capture the nonlinear temperature variations in the debris; and at three of the four sites the sensors were retrieved such that the thermal conductivity could be estimated (Conway and Rasmussen, 2000).
Ablation stakes were also installed at 14 sites.One site had a debris thickness greater than 1 m, so an ablation stake could not be installed.The ablation stakes were installed by excavating to the debris-ice interface, at which point the debris thickness was measured, and then a 2 in.diameter hole was drilled vertically approximately 1 m into the ice using a manual ice drill (Kovacs Enterprise).A 2 m piece of 1.5 in.PVC pipe was placed into the hole and the height from the top of the ice to the top of the pipe was measured to determine the exact length that the PVC pipe was inserted into the ice.A PVC end cap was then place on top of the pole to prevent anything from entering the hole through the pipe.The debris was then replaced in its approximate original position.

Meteorological data
The meteorological data used in the model calibration and validation were from an automatic weather station (AWS), Pyramid Station (27.959 • N, 86.813 • E; 5035 m a.s.l., SHARE Network operated by EV-K 2 -CNR), located off-glacier, next to the Khumbu Glacier, approximately 14 km northwest of Imja-Lhotse Shar Glacier.The meteorological data provided by Pyramid Station were unvalidated; i.e., minute measurements of air temperature, wind speed, relative humidity, global radiation, precipitation, and snow depth were used prior to their quality-control processing.The data were processed to be consistent with the halfhour debris temperature measurements on Imja-Lhotse Shar Glacier.The air temperature, wind speed, relative humidity, and global radiation data were reviewed and deemed plausible, so no adjustments were performed.The half-hour precipitation data were determined by summing the precipitation over each half-hour time step.A few of the minute measurements recorded negative precipitation, which were assumed to be zero as negative precipitation is not feasible.The halfhour snow depth data were processed to assume a snow depth of zero if snow was not recorded on the ground for the entire half-hour.The average snow depth over the half-hour was then computed and any average snow depth less than 0.001 m was considered to be zero.Wind speed data were collected at 5 m and adjusted to 2 m to be consistent with air temperature measurements for the turbulent heat fluxes, assuming a logarithmic dependence (Fujita and Sakai, 2014).The snow depth data were used to derive a snowfall rate, assuming a density of snow of 150 kg m −3 .The data were available from 31 May to 12 October 2014 with a few short gaps (11.9 % data missing).The first 2 days of meteorological data were used as start-up time for the model.
Long-wave radiation was not measured at Pyramid Station during this period; therefore, the downward long-wave radiation flux from NCEP/NCAR reanalysis data (Kalnay et al., 1996) was used with a minor modification.A comparison of the downward long-wave radiation flux from NCEP/NCAR and the incoming long-wave radiation flux at Pyramid Station from 2003 to 2010 (neglecting any data gaps) between the months of June and September revealed that NCEP/NCAR overestimated the incoming long-wave radiation by an average of 29 W m −2 (results not shown).Therefore, the NCEP/NCAR downward long-wave radiation flux was adjusted to account for this overestimation when being used in conjunction with the Pyramid Station data.This reanalysis data set provides 6 h meteorological data and was resampled using a linear interpolation such that the temporal resolution of the incoming long-wave radiation agreed with the half-hour debris temperature measurements.
Ablation rates were modeled over the same time period as the ablation stakes (18 May to 9 November).For days where no meteorological data were available, i.e., the data gaps, the ablation for that day was assumed to be equal to the daily ablation rate for that specific month.As the available meteorological data began on 31 May, the daily ablation rate for the month of May was assumed to be equal to the daily ablation rate of the first week of June.Temperature sensors revealed the debris was snow-covered from 26 May to 1 June, so the melting during these days was assumed to be zero.Temperature profiles also show the debris was snow-covered from 13 to 20 October and deeper thermistors revealed the temperature remained around freezing until the sensors were removed in November.Therefore, the melt rates after the 12 October were assumed to be zero.

Surface roughness (z 0 )
Structure from Motion (SfM) was used to derive fineresolution (i.e., centimetric) digital elevation models (DEMs) at four sites (Sites A-D) located on the debris-cover of Imja-Lhotse Shar Glacier (Fig. 2).In brief, SfM relies upon the acquisition of a series of overlapping images that capture the features of the terrain from a number of different vantage points.Computer vision techniques detect matching features between images using multiscale image brightness and color gradients, and a highly iterative bundle adjustment procedure is used to develop a three-dimensional structure of the sur-face (Snavely et al., 2008).Camera positions and orientations are solved simultaneously with surface geometry utilizing the high level of redundancy afforded by a large overlapping image set.Ground control points (GCPs), collected using a total station with an error less than 0.4 mm, are then used to transform the relative three-dimensional surface into an absolute coordinate system.The resulting point-cloud data are comparable in both density and accuracy to those generated by terrestrial laser scanning (Westoby et al., 2012) and can either be used as they are, or decimated (as in this study) to generate gridded elevation data.The use of SfM within geoscience is well reviewed by Westoby et al. (2012) and specific details of the mathematical operations involved can be found in Snavely (2008) and Szeliski (2011).Here, we therefore focus mostly on our field method and subsequent roughness analysis.
At each of our sites ∼ 40 photos were taken around a roughly 2 m × 2 m grid.Cones were placed in the four corners of the grid as GCPs and their location was measured using a total station with a local coordinate system.The GCPs and photos were processed using Agisoft PhotoScan Professional Edition, Version 1.1.0to create a DEM for each site.At each stage, the highest accuracy settings were chosen.No a priori information about camera position or orientation was recorded, so these were estimated coincidentally as part of the adjustment.In each case the initial estimates of camera position and altitude were accepted and used to generate a sparse point cloud (10 3 -10 4 points).A moderate depth filter was then used to derive a dense cloud (10 6 -10 7 points), and subsequently a mesh was constructed using the height field as the surface type.The error of the DEM was computed as the root-mean-square error based on the differences between the measured GCPs from the total station and the modeled position of the GCPs from the software.The resulting DEMs were then resampled in ArcGIS 10.3 to a resolution of 0.01 m and were clipped to remove the cones from the subsequent analyses.The DEM was then fit with an x-y plane using a method of least squares such that the DEM was flattened with a mean elevation of zero.
These processed DEMs of the four sites were analyzed to determine the surface roughness, z 0 .Lettau (1969) developed an empirical relationship to estimate z 0 : where h * is the average vertical extent or effective obstacle height, s is the silhouette area or area of the upwind face of an average element, and S is the specific area or unit ground area occupied by each obstacle.Previous studies have estimated the variables in Eq. ( 1) through a simplified standard deviation approach (which will be referred to as the Lettau-Munro method), based on the variations in elevations and the number of continuous positive groups above the mean elevation (Munro, 1989;Rees and Arnold, 2006;Brock et al., 2006).Initially, the Lettau-Munro method was applied to measure z 0 for every row and column transect of the four DEMs; however, the resulting values of z 0 did not capture the variations between sites and may have been slightly underestimated (see results).
Consequently, an alternative method was developed to estimate the effective height, silhouette area, and unit ground area of each obstacle using a similar transect approach and taking advantage of the high-resolution DEM.One problem with applying the method from Lettau (1969) is the lack of a clear definition of what constitutes an obstacle.The surface roughness will greatly vary depending on what is considered to be an obstacle, so a method must be developed that (i) objectively determines the obstacle height and (ii) yields reasonable estimates of surface roughness regardless of the resolution of the DEM.Smith (2014) states that the relationship developed by Lettau (1969) holds at low roughness densities (< 20-30 % of the surface area), beyond which the observed z 0 is less than that predicted by Lettau (1969) because the obstacles begin to aerodynamically interfere with one another.Therefore, a method was developed to select an obstacle height based on an obstacle density of 30 %.
Initially, all the relative topographic highs and lows were identified.This was done for all of the transects in each of the four cardinal directions with respect to the DEM, i.e., every east-west, north-south, west-east, and south-north transect.Every elevation change between a relative low and high was considered a potential obstacle.The depth of each obstacle was defined as the distance between two low points surrounding the obstacle's high point.In the event that an obstacle was identified, but there was no low point following the high point, i.e., the low point was outside the extent of the transect, then the depth of the obstacle could not be determined.Figure 3 shows an example of a transect from Site B, which identifies the obstacle's height and depth based on the method developed in this study.The obstacle density was then defined as the cumulative depth of all the obstacles above the obstacle threshold divided by the length of the transect.An iterative approach was then used to determine the obstacle height that causes the obstacle density to reach the 30 % threshold.
Once the obstacle height has been determined, the silhouette area and unit ground area were approximated from the height and depth of the obstacles.Specifically, the silhouette area was taken to be the height of the obstacle times a unit width and the unit ground area was estimated as the depth of the obstacle times a unit width.Based on these definitions, Eq. ( 1) may be simplified to Eq. ( 2): where d * obst is the average depth of the obstacle.The surface roughness, z 0 , was computed using the average effective obstacle height and average obstacle depth for each transect.In the event that an obstacle was identified, but did not have a depth, the obstacle's height was still used in the average.

Debris-covered glacier energy balance model
The model used in this study was a steady-state surface energy balance model for a debris-covered glacier, where where R n is the net radiation flux, H is the sensible heat flux, LE is the latent heat flux, P is the heat flux supplied by rain, and Q c is the ground heat flux (all in W m −2 ).The net radiation and sensible heat fluxes are fully described in Rounce and McKinney (2014); however, in the current study the incoming short-wave radiation was only corrected for the effects of topography as shading could not be considered due to the lack of a high-resolution DEM of the glacier.
The latent heat flux is difficult to determine without detailed knowledge of the moisture in the debris or the relative humidity at the surface.As the surface relative humidity was unknown, this study has analyzed three methods for estimating the latent heat flux: (1) assuming the debris is dry (LE = 0), (2) assuming it is dry unless the relative humidity is 100 %, at which point the surface relative humidity is assumed to also be 100 % based on the assumption that the water vapor above the surface is well mixed, and (3) assuming the surface is saturated when it is raining.These methods for modeling the latent heat flux will be referred to herein as LE Dry , LE RH100 , and LE Rain , respectively.The reservoir approach detailed by Collier et al. (2014) and the empirical relationship between debris thickness and wetness (Fujita and Sakai, 2014) were not applied to this study due to the limited amount of knowledge of moisture within the debris and how the debris properties change with respect to depth.The latent heat flux is thus estimated according to Nicholson and Benn (2006): where where ρ air is the density of air at standard sea-level pressure (1.29 kg m −3 ), P 0 is the standard air pressure at sea level (1.013 × 10 5 Pa), L e is the latent heat of evaporation of water (2.49× 10 6 J kg −1 ), A is a dimensionless transfer coefficient, u is the wind speed collected at a height of 2 m (m s −1 ), e z and e s are the vapor pressures (Pa) at height z, 2 m, and on the surface of the debris, respectively, k vk is von Kármán's constant (0.41), and z 0 is the surface roughness.
The heat flux due to precipitation was estimated following Reid and Brock (2010): where ρ w is the density of water (999.97kg m −3 ), c w is the specific heat capacity of water (4.18 × 10 3 J kg −1 K −1 ), w is the rainfall rate (m s −1 ), and T r is the temperature of rain (K), which was assumed to be equal to the air temperature.
The debris layer was broken down into layers of 0.01 m such that the nonlinear temperature profiles in the debris could be captured using a Crank-Nicolson scheme (Reid and Brock, 2010).The conductive heat flux at the surface and at the debris/ice interface were estimated following Reid and Brock (2010): where k eff is the effective thermal conductivity (W m −1 K −1 ), h is the height of each layer in the debris set at 0.01 m, and T d (1), T d (N − 1), T ice are the temperatures (K) of the first layer in the debris, the last layer before the debris/ice interface, and the temperature of the ice (273.15K), respectively.The surface temperature was computed at half-hour time steps using an iterative Newton-Raphson method approach as detailed in Reid and Brock (2010).In the event of snow, a simple snowmelt model was used (Fujita and Sakai, 2014), which applies an energy balance over the snow surface that includes net radiation, turbulent heat fluxes, and conductive heat flux with the debris layer in addition to a variable surface albedo of the snow based on the number of days since fresh snow and the air temperature.The thermal conductivity of snow was assumed to be 0.10 W m −1 K −1 (Sturm et al., 1997;Sturm et al., 2002;Rahimi and Konrad, 2012) and the surface roughness of the snow was assumed to be 0.002 m (Brock et al., 2006).If snow was on the surface, all the heat fluxes at the debris surface were assumed to be zero, with the exception of the conductive heat flux in the debris and at the debris/snow interface.If all the snow was melted on the surface, then the next time step returned to the snow-free energy balance model.
As detailed knowledge of albedo, thermal conductivity, and surface roughness was not available for the sites where temperature sensors were installed, the debris-covered glacier energy balance model was calibrated at each site from 2 June to 30 July 2014.The calibration was performed by minimizing the total sum of squares of the measured versus modeled surface temperature for each site and was done independently for the three methods used to estimate the latent heat flux.Bounds for the thermal conductivity and surface roughness were based on measured field data (see results), while the bounds for the albedo were 0.1-0.4(Inoue and Yoshida, 1980;Kayastha et al., 2000;Nicholson and Benn, 2012;Lejeune et al., 2013).A validation was then conducted at each site using data from 8 August to 12 October 2014 to assess how well the calibrated model performed.

Thermal conductivity (k)
The thermal conductivity, k, of the debris was computed using the temperature measurements from Sites 4, 11, and 13 over the time period of the study (2 June-12 October 2014) following the methods of Conway and Rasmussen (2000).The calculations used standard values for the density of rock (2700 kg m −3 ), volumetric heat capacity of rock (750 J kg −1 K −1 ± 10 %), and effective porosity (0.33) based on Nicholson and Benn (2012).Depending on the vertical spacing of temperature sensors at a site, k was computed at depths of 0.05, 0.10, and 0.20 m.The values of thermal conductivity ranged from 0.42 (±0.04) to 2.28 (±0.23)W m −1 K −1 .The average value of k for each site was 1.44 (±0.14), 1.62 (±0.16), and 0.47 (±0.04)W m −1 K −1 for Sites 4, 11, and 13, respectively.
These values agree well with other studies in the Everest region that have found the thermal conductivity to vary between 0.60 and 1.29 W m −1 K −1 (Conway and Rasmussen, 2000;Nicholson and Benn, 2012;Rounce and McKinney, 2014).In September 2013, Rounce and McKinney (2014) found the thermal conductivity to be greatly influenced by depth; however, this trend was not apparent in our current data.We believe this disparity can be explained by the time period during which the data were collected.It is likely that the temporally limited data (13-24 September 2013) presented in Rounce and McKinney (2014) represent a constantly dry surface, whereas here we observed an entire melt season, where the surface is exposed to precipitation and snow.The thermal conductivities appeared to show a trend over the monsoon season where the highest thermal conductivities were typically observed in July and August, which coincided with higher average air temperature and increased precipitation compared to the other months.As k eff is one of the parameters that is used to calibrate the model, the range of average thermal conductivity (0.47-1.62 W m −1 K −1 ) will be used to bound k eff .

Surface roughness (z 0 )
The DEMs generated using the SfM workflow had a total root-mean-square error of 0.008-0.024m.Table 2 shows that the errors in elevation (i.e., z) were smaller than in planform (i.e., x and y) with a maximum error of 0.007 m.The contrast between elevation and planimetric errors is likely a result of the identification of the GCPs in each photo during the SfM workflow, since it was easier to identify the top of the cone in each photo than it was to determine the exact point on the rim of the cone.As the error with the total station is small (maximum of 0.4 mm), this human error likely dominated the total error, although errors in estimates of both camera position and orientation will also have contributed.The DEMs were resampled to a resolution of 0.01 m such that their resolution was on the same order as their respective errors (three of the four sites had a total RMSE less than 0.01 m).The DEMs were then de-trended to account for variations in the local topography.
Initially, z 0 was estimated from Eq. (1) using the Lettau-Munro method.The average value of z 0 was 0.0037, 0.0091, 0.0022, and 0.0033 m for Sites A, B, C, and D, respectively.These values are towards the lower end of those previously reported in literature, which were estimated from wind speed profiles and range from 0.0035 to 0.060 m (Inoue and Yoshida, 1980;Takeuchi et al., 2000;Brock et al., 2010).In particular, the average value of z 0 for Sites A, C, and D was comparable or smaller to Area IV on the Khumbu Glacier (Inoue and Yoshida, 1980), which comprised small schist with bare ice.Sites A, C, and D were all debris-covered with boulders ranging up to 0.40 m, so these small estimations of z 0 are concerning.Site B also appears to be underestimated as it has similar debris characteristics to Area III on the Khumbu Glacier (Inoue and Yoshida, 1980), yet its average value of z 0 was much smaller (0.0091 m compared to 0.060 m, respectively).These apparent underestimations of z 0 led to the development of an alternative method.
The alternative method relies upon the selection of the obstacle height (or threshold) such that the obstacle density is 30 % (Smith, 2014).Figure 4 shows that as the obstacle threshold is increased, the obstacle density decreases, which www.the-cryosphere.net/9/2295/2015/The Cryosphere, 9, 2295-2310, 2015 makes intuitive sense as there will be fewer obstacles in the transect.Table 3 shows that the values of z 0 using an obstacle density of 30 % and the highest-resolution DEM (0.01 m) were 0.016, 0.043, 0.006, and 0.014 m for Sites A-D, respectively.These values agree well with the range of z 0 values previously reported (Inoue and Yoshida, 1980;Takeuchi et al., 2000;Brock et al., 2010).Furthermore, these z 0 values appear to capture the inter-site variability as Site B had the highest value of z 0 (0.043 m), which is expected since the debris cover includes large boulders up to 1 m in size (Fig. 2).Site C, which comprised the smallest grain sizes of the four sites in this study had the lowest estimation of z 0 (0.006 m).
Sites A and D had similar values of surface roughness and the standard deviations of Site A and D (0.008 and 0.012 m, respectively) appear to capture the more homogenous surface of Site A compared to the highly heterogeneous surface of Site D (Fig. 2).
The impact of DEM resolution on obstacle threshold and z 0 was analyzed to determine the robustness of this alternative method.Since the terrain is not changing and only the sampling frequency is varying, the z 0 values should remain fairly constant.Table 3 shows that the obstacle threshold increases as the resolution of the DEM becomes coarser.This occurs because the coarser resolution cannot capture the subtler changes in surface height over the debris cover.On the other hand, the estimations of z 0 remain relatively constant (±0.004 m) as the DEM resolution is reduced to 0.04 m.The consistency of this method despite variations in DEM resolution and obstacle thresholds and the objective approach for deriving the obstacle threshold using a 30 % obstacle density, lends confidence to this method.Furthermore, Nield et al. (2013) found that measures regarding surface heights are the best predictor of aerodynamic roughness, specifically for surfaces that comprise large elements or have patches of large and small elements.Therefore, it is expected that the obstacle threshold should vary for different sites as a function of their largest elements, which is consistent with the inter-site variability and the obstacle thresholds reported in this study.

Ablation stakes
Ablation stakes were installed on 18-19 May 2014 approximately 1 m into the ice at 14 sites with debris thicknesses ranging from 0.07 to 0.54 m (Table 1).The ablation stakes were measured on 9 November 2014.At 11 of the 14 sites, the ablation stakes completely melted out of the ice, indicating there was greater than 1 m of ablation.Sites 8, 13, and 15 had ablation measurements of 0.92, 0.85, and 0.89 m, respectively.These three sites had debris thicknesses of 0.20, 0.33, and 0.37 m and were oriented in the southern, northeast, and northwest directions, respectively.The lower ablation rates of Sites 13 and 15 compared to the other 12 sites is likely due to a combination of their debris thickness and aspect as they are oriented in a manner that receives less solar radiation throughout the day.Site 8 appears to be an anomaly as it has  a smaller debris thickness than eight of the sites with ablation stakes and a southerly aspect, which positions it in a manner to receive a greater amount of solar radiation throughout the day.It is possible that Site 8 had a higher albedo and/or a lower thermal conductivity, which would greatly reduce its ablation; unfortunately, these properties could not be measured in the field.Nevertheless, the ablation measurements indicate that understanding ablation rates on debris-covered glaciers is greatly influenced by slope, aspect, and properties of the debris (albedo, thermal conductivity, and surface roughness).

Model calibration
Three different methods were used to estimate the latent heat flux to determine how well each method models the measured debris temperatures.These methods are referred to as LE Rain , LE RH100 , and LE Dry .The albedo, thermal conductivity, and surface roughness for each of the three methods were optimized by minimizing the sum of squares of the surface temperature for each site (Table 4).For the LE Rain and LE RH100 model, 7 of the 10 sites had a thermal conductivity at the upper bound (1.62 W m −1 K −1 ), while for the LE Dry model 9 of the 10 sites were at the upper bound.These results indicate that the selection of the upper bound for the thermal conductivity is important and its impact on model performance is detailed in the discussion section.The albedo values ranged from 0.10 to 0.40 and had an average value around 0.32, which is consistent with albedos measured in the Khumbu (Inoue and Yoshida, 1980;Kayastha et al., 2000;Nicholson and Benn, 2012;Lejeune et al., 2013).The values of z 0 had an average around 0.014 m, which is consistent with z 0 measured in this study and those reported on other debris-covered glaciers (Inoue and Yoshida, 1980;Takeuchi et al., 2000;Brock et al., 2010).For the LE Rain and LE RH100 models, 5 of the 10 sites had a value of z 0 at its lower bound (0.006 m), which highlights the importance of measuring the surface roughness of the debris cover and will be discussed in the sensitivity analysis.
The performance of each model was assessed using the total sum of squares and the R 2 correlation coefficients.The R 2 values ranged from 0.34 to 0.92 for all three models.The average R 2 values over the calibration period for the LE Rain , LE RH100 , and LE Dry models were 0.72, 0.72, and 0.71, respectively.Figures 5c and d show the correlation between the modeled and measured surface temperature at Site 11, which had an R 2 of 0.77 and 0.75 for the calibration and validation periods, respectively.Figure 5c shows there is good agreement between the modeled and measured temperature sensors.The modeled temperatures appear to capture the daily variations in temperature well.However, there are a few days for which a positive bias in temperature can be seen during the daily high and nightly low (e.g., Fig. 5c from 16 to 18 June and 25 to 27 July).Interestingly, the overestimation of the daily high typically occurs after the nightly low has a positive bias in temperature during the previous night.The positive bias of the nightly minimum is apparent between the hours of 00:00 and 06:00 (Fig. 6).One possible explanation for the positive bias in temperature in the nightly low is an overestimation of the incoming long-wave radiation due to the poor temporal and spatial resolution of the NCEP/NCAR reanalysis data set compared to the other meteorological data from Pyramid Station.Typically, the wind speed during the night is relatively low thereby limiting the turbulent heat fluxes, which causes the incoming long-wave radiation to be a major source of energy during this time.
Nonetheless, the model performs reasonably well for all of the temperature sensors.Unfortunately, it is difficult to determine which latent heat flux model performs the best using the total sum of squares and/or the R 2 values as there was not one particular model that consistently had a lower total sum of squares and/or a higher R 2 at each site.The average R 2 value was fairly comparable for all three models.The total sum of squares of all the sites was the lowest for the LE RH100 model, followed by the LE Rain model and then the LE Dry model, but the difference between models was less than 5 %.

Model validation
Model validation was assessed from 8 August to 12 October 2014 for all three models using the R 2 values for each temperature sensor.The R 2 values for all the temperature sensors at the 10 sites ranged from 0.39 to 0.81 for all three methods.The average R 2 value for the LE Rain , LE RH100 , and LE Dry was 0.67, 0.67, and 0.68, respectively.Again, the similar performance between the three models does not provide any insight into preference for one model and is likely a result of the calibration procedure.Figure 5d shows that the LE Rain model performs well through the entire validation period.Similar to the calibration period, the LE Rain model appears to underestimate the nightly low, which causes the following daily high to be overestimated.Reid and Brock (2010) found R 2 values of 0.94 and 0.52 for temperature sensors at the surface and at a depth of 15 cm, respectively.While the R 2 value of 0.94 is higher than R 2 values found in this study, the range of R 2 is comparable.In contrast to the findings of Reid and Brock (2010), the average R 2 value for the surface temperature sensors (0.67-0.68 for all three models) was very similar to the average R 2 value of those buried in the debris (0.66-0.68 for all three models).The slightly lower R 2 values in this study may be a result of using meteorological data from an AWS located 14 km away from the glacier.Furthermore, long-wave radiation was estimated from remotely sensed data, which may also influence model performance as previously discussed.

Modeled ablation rates
Ablation rates were computed for all 15 sites that had a temperature sensor or an ablation stake.For sites that only had an ablation stake, the average calibrated parameters for that particular latent heat flux model were used.Additionally, ablation rates were estimated for the LE Rain model using the average calibrated parameters for all the sites to assess the differences between using a single set of parameters compared to optimizing the parameters at each site.The modeled ablation over the entire duration of the study period varied from 0.39 to 2.85 m among the three methods (Fig. 7).On average the LE Dry model overestimated both the LE Rain and LE RH100 models by 7.9 %.The slight variations in ablation between the models are directly related to the differences in their calibrated parameters.The slightly higher ablation rates for the LE Dry model are likely attributed to the higher values of thermal conductivities and the lack of a latent heat flux term to remove heat from the debris.Figure 7 shows there is a clear relationship between debris thickness and ablation as thin debris has higher rates of ablation compared to thicker debris, which insulates the ice to a greater Table 5. Sensitivity analysis showing percent changes relative to the total melt (m) as a function of the uncertainty associated with the calibrated parameters (α, k, z 0 ) for all sites over the study period using the LE Rain model in conjunction with the average calibrated parameters for all sites.extent thereby retarding ablation.The scatter found throughout the curve, specifically between 0.25 and 0.50 m, is due to the site-specific debris properties and the slope and aspect of each site.A comparison between the LE Rain model using the optimized parameters at each site and those using the average calibrated parameters at each site highlights the effect that site-specific properties has on ablation.Site 6, with a debris thickness of 0.08 m, is a good example as the use of average calibrated parameters increased the melt from 2.03 to 2.70 m due to an increase in thermal conductivity from 1.29 to 1.52 W m −1 K −1 .These differences in melt and the sensitivity to thermal conductivity highlight the importance of properly estimating/measuring the thermal conductivity of the debris cover.The modeled ablation rates may also be compared to the measured ablation rates.Specifically, Sites 8, 13, and 15 had measured ablation rates of 0.92, 0.85 and 0.89 m compared to their modeled ablation rates of 1.76, 0.76, and 1.22 m, respectively, for the LE Rain model.The large discrepancy between the measured and modeled ablation rates at Site 8 may be due to the lack of knowledge of the debris properties at Site 8 as previously discussed.The difference between the modeled and measured ablation rates at Site 15 may also be a result of the thermal conductivity parameter (1.61 W m −1 K −1 ), which is slightly higher than thermal conductivities previously reported in the Khumbu, which ranged from 0.60 to 1.29 W m −1 K −1 .A comparison of the daily average temperatures for Site 15 reveals there was about an hour lag between the modeled and measured tem-peratures (results not shown).Lags between temperatures are typically a result of their depth (Conway and Rasmussen, 2000), which is apparent in Fig. 5a as the 0.10 m sensor lags behind the 0.01 m sensor.It is possible that debris may have shifted over the melt season, causing the measured temperature to be at a lower depth than 0.01 m, which would greatly influence the model calibration and potentially cause the thermal conductivity to vary.Site 5 was the only other site where a slight lag was observed between the measured and modeled temperatures.The modeled ablation at Site 5 was 0.87 m using the LE Rain model, while the ablation stake melted completely out of the ice, indicating greater than 1 m of ablation.All the other model estimates of ablation were near to or greater than 1 m, which was also observed by their respective ablation stakes as they completely melted out of the ice.

Parameter
The ablation results also show strong seasonal trends with maximum melt rates occurring in June, July, and August.These ablation rates appear to taper off towards the transition seasons.Melt rates in July and August ranged from 0.3 to 2.5 cm day −1 based on the debris thickness, which is consistent with empirical relationships between mean daily ablation rate and debris thickness found on other glaciers (Nicholson and Benn, 2006).The total ablation rates may also be compared to measured surface elevation changes on Imja-Lhotse Shar Glacier derived from multiple DEMs, which were found to range from −0.82 ± 0.61 m yr −1 to −1.56 ± 0.80 m yr −1 (Bolch et al., 2011, Nuimura et al., 2012;Gardelle et al., 2013;Thakuri et al., 2015)  time periods between 1999 and 2014.It is important to note that the mass-balance estimates in these studies have been converted back to elevation changes using the ice density reported in each study.The total ablation rates are similar to the measured changes in surface elevation, which lends confidence to the results.

Sensitivity analysis
A sensitivity analysis was performed to assess how albedo, thermal conductivity, and surface roughness affect the total ablation (Table 5) based on the uncertainty with respect to each parameter.The uncertainty in thermal conductivity was ±0.40 W m −1 K −1 , which captures the approximate difference between the highest thermal conductivity measured in this study (1.62 W m −1 K −1 ) and the higher end of those previously reported (Conway and Rasmussen, 2000;Nicholson and Benn, 2012;Rounce and McKinney, 2014).The uncertainty associated with the surface roughness was ±0.010 m, which is the approximate standard deviation associated with the z 0 values for each of the three models (Table 4) and similar to the standard deviation between the four sites where z 0 was measured (±0.016 m).Lastly, the uncertainty of the albedo was estimated as ±0.10, which is the approximate standard deviation within the model calibration for each of the three models and also the difference between the mean and median albedo measured by Nicholson and Benn (2012) on Ngozumpa Glacier.The LE Rain model was used as the baseline case and the average value for each of the calibrated parameters (α, k, z 0 ) from the model optimized was used for each site.
Table 5 shows that the total ablation is most sensitive to changes in the thermal conductivity, where a ±0.40 W m −1 K −1 change causes a ±20.5 % change in total ablation on average.The uncertainty associated with the thermal conductivity is also more sensitive to thicker debris, which is consistent with the findings of Nicholson and Benn (2012).Total ablation is also moderately sensitive to changes in the albedo, where a ±0.10 change causes a ±12.0 % change in total ablation.Lastly, the total ablation is least sensitive to changes in increasing the surface roughness, as a +0.010 m increase in z 0 only caused a −7.3 % change in total ablation.However, the model was quite sensitive to a reduction in the z 0 of −0.010 m, which caused an average change in total ablation of +15.0 %.The sensitivity associated with z 0 also appears to increase with an increase in debris thickness.These results highlight the importance of properly estimating the thermal conductivity, but also show that the surface roughness and the albedo are important as well.Nicholson and Benn (2006) proposed that the temperature gradient in the debris may be assumed to be linear at a time step greater than a day, but is nonlinear for shorter time steps.This would have important implications for modeling melt on remote debris-covered glaciers where meteorological data are not available and reanalysis data sets could be used instead.The importance of temporal resolution was analyzed using 6 h and daily average data from Pyramid Station, which are consistent with the temporal resolution of NCEP/NCAR reanalysis data sets.To be consistent with this reanalysis data set such that only the effects of temporal resolution were analyzed, the wind speed and relative humidity used were instantaneous values from Pyramid Station, while all the other variables were 6 h averages.For the daily time step, all the parameters were daily averages and the temperature profile in the debris is assumed to be linear.The LE Rain model was used to model the latent heat flux.

Temporal resolution
The R 2 correlation coefficients for the sites with temperature sensors and the modeled total melt for all 15 sites were used to assess the effect of temporal resolution on model performance.The R 2 using the 6 h data ranged from 0.30 to 0.80 with an average of 0.55 over the calibration period and was significantly poorer during the validation period with R 2 values ranging from 0.15 to 0.65 with an average of 0.35. Figure 8a shows the surface temperature at Site 11 does fairly well (R 2 = 0.63) at modeling the measured surface temperatures over the calibration period.The lower R 2 values compared to the 30 min time step appear to be a result of the 6 h model underestimating the daily high, which occurs around 15:00 each day.Furthermore, Fig. 8b shows the 6 h model poorly replicates the measured data towards the transition seasons when snowfall becomes significant, which explains the poorer R 2 values for the validation period.Snowfall is problematic in the model for large time steps because the model assumes the snow is on the surface for the entire time step.Therefore, a small snow event that could melt quickly on the debris and then allow the debris to warm up during the day is perceived to remain on the snow for the 6 h time step (e.g., Fig. 8b from 27 September to 3 October).The same problem arises at the daily time step, so a snow-free model was used instead.For the daily time step, the R 2 values ranged from 0.18 to 0.63 with an average of 0.29. Figure 8c shows the daily time step is able to capture some of the temperature fluctuations over the melt season, but does not perform as well as the 30 min or 6 h models.
Since we are most interested in understanding the effects of temporal resolution, the 6 h data and daily averages from 2 June to 25 September 2014 were assessed, which is prior to the time when snowfall was recorded each day.A comparison of all the modeled and measured temperatures at the surface reveals the 6 h model underestimates the measured temperatures by an average of 1.0 (±4.3)K over the entire time period.The modeled total ablation from 2 June to 25 September reveals the ablation is consistently underestimated at all sites by an average of 11 (±5) %.The lower estimates of ablation are likely a result of the underestimation of the daily high as previous discussed.Similar to the 6 h model, the daily time step model slightly underestimates the measured temperatures at the surface on average by 0.3 (±1.9)K.The modeled total ablation is also underestimated by an average of 6 (±10) %.However, it is important to note that 5 of the 15 sites actually slightly overestimated the melt.These results suggest that if high-temporal meteorological data are not available, a first-order estimate of ablation over a melt season could be obtained using the daily time step model.It is important that the estimate is made over the entire melt season, as the daily model does not capture the daily temperature fluctuations well.Furthermore, caution should be used to avoid the transition seasons as both the daily model and the 6 h model do not have a small enough time step to properly account for snowmelt.

Thermal conductivity
One of the limitations with regards to the thermal conductivity measurements is that all the measurements were made near the surface.Therefore, the estimates of the average thermal conductivity at each site are potentially underestimated because the deeper layers that may be more compact and humid are not considered.The lack of any trends with respect to depth appears to dispute this theory; however, this is based on a limited number of measurements near the surface.Interestingly, the thermal conductivities measured at Sites 4 and 11 are similar to those estimated by Nicholson and Benn (2012) for debris cover on Ngozumpa Glacier with 10 and 20 % of the void space being filled with water (1.42 and 1.55 W m −1 K −1 , respectively).
The number of sites that reached the upper bound during the model calibration is concerning as it may indicate that the actual thermal conductivity throughout the debris is higher.To address this issue, an additional calibration was performed allowing the thermal conductivity to be unbounded.This calibration revealed that three or more out of the 10 sites for each method had thermal conductivities greater than 3.0 W m −1 K −1 with one thermal conductivity as high as 4.5 W m −1 K −1 .The lithology of the debris cover in the Everest region is predominantly granite, gneiss, and pelite (Hambrey et al., 2008).Robertson (1988) found the thermal conductivity of solid granite gneiss to be 2.87 W m −1 K −1 , so the unbounded thermal conductivities do not appear to make physical sense when one considers that the thermal conductivity of debris should be much lower than solid rock due to the pore spaces being filled with air and water.Furthermore, an optimization performed using the total sum of squares of all the surface sites reveals that increasing the thermal conductivity from 1.6 W m −1 K −1 to its minimum of 2.6 W m −1 K −1 only reduces the total sum of squares by 3 %.These results and the similar values to Nicholson and Benn (2012) lend confidence to the use of 1.62 W m −1 K −1 as the upper bound, but highlight the importance of understanding how the moisture varies within the debris and its influence on the thermal conductivity.Future work should improve measurements of the thermal conductivity by (a) accurately measuring the depth of the temperature sensors during installation and retrieval, (b) installing additional sensors (e.g., 5 cm spacing) that allow thermal conductivity within the debris to be computed at more depths, and (c) measuring moisture in the debris at various depths.

Surface roughness
The development of an alternative method for estimating z 0 was required, as the Lettau-Munro method appeared to greatly underestimate the values of z 0 .The alternative method applies the relationship developed by Lettau (1969) to a high-resolution DEM using the selection of an obstacle height (threshold) based on an obstacle density of 30 %.One of the main limitations of this study is the lack of aerodynamic roughness measurements to validate the developed methods.Previous work, e.g., Rees and Arnold (2006), has relied upon surface roughness estimates from other studies to assess the reasonableness of their results when aerodynamic data were not collected.This study relies upon the results of Inoue and Yoshida (1980), which estimated surface roughness using wind speed profiles at two sites on the Khumbu Glacier.Specifically, Sites B and C in this study have similar debris cover to Areas III and IV from Inoue and Yoshida (1980) in size.This value is similar to the higher value of 0.060 m for z 0 derived from a region on the Khumbu Glacier that consisted of large granitic boulders of 1-2 m in size lying on top of schistose rocks with a grain size varying from a few centimeters to 0.5 m (Inoue and Yoshida, 1980).The larger boulders observed by Inoue and Yoshida (1980) may explain the slightly higher value of z 0 compared to Site B. Site C, which comprised the smallest grain sizes of the four sites in this study, agrees well with the smaller value of z 0 (0.0035 m) derived by Inoue and Yoshida (1980) for an area where the supraglacial debris comprised dispersed boulders ranging in size of 0.01-0.05m.The few boulders ranging in size of up to 0.15 m may be the reason for Site C's slightly larger value of z 0 (0.006 m).Sites A and D were composed of boulders and grains that varied in size between those found in Sites B and C; therefore, we deem the value of z 0 of 0.016 and 0.014 m for Sites A and D, respectively, to be reasonable.Furthermore, these values agree fairly well with the z 0 of 0.016 m measured by Brock et al. (2010) on a debris-covered glacier in Italy that comprised a mixture of granites and schists of predominantly cobble size, with occasional boulders of < 1 m size.
Future work should seek to compare these estimates of surface roughness with roughness to determine the scale at which these two values agree.Brock et al. (2006) found there to be no significant difference between the use of a 3 and 15 m transect; however, they did state that a shorter pole would be unlikely to capture a sufficient sample of roughness elements if the vertical changes are greater than 1 m.The use of hundreds of transects over a ∼ 4 m 2 grid has the benefit of expanding the number of surface roughness elements that can be captured compared to a single transect.However, Brock et al. (2006) compared microtopographic and aerodynamic roughness over snow, slush, and ice, which is significantly different from the hummocky and heterogeneous terrain on debris-covered glaciers.Therefore, it will be important to determine the scale or fetch length at which the surface roughness agrees with the aerodynamic roughness.Nonetheless, the method developed in this paper provides an objective approach to select an obstacle height and yields consistent and reasonable estimates of z 0 for various grain sizes independent of the resolution of the DEM.

Modeled results
One of the limitations of the calibration procedure is that the LE Rain , LE RH100 , and LE Dry models all performed reasonably well.The lack of a single model clearly outperforming the others indicates that either (a) the modeling of the latent heat flux is insignificant or (b) the latent heat flux is significant, but the calibration procedure allows for changes in the latent heat flux to be compensated for via other model parameters.Brock et al. (2010) found that latent heat fluxes may be a significant energy sink when rain falls on warm debris, indicating that it is important to include the latent heat flux.
They also assessed the importance of each component of the energy balance and found that including the latent heat flux improved the correlation coefficient of their model.The average latent heat flux for both the LE Rain and LE RH100 models were comparable with values ranging from −53 to 10 W m −2 over the day.The peak instantaneous latent heat fluxes varied greatly between the two models with fluxes as high as −714 and −323 W m −2 , for the LE Rain and LE RH100 models, respectively.These values are similar to those reported by Brock et al. (2010) and support the importance of including the latent heat flux term.However, they do not yield any insight into preference between the LE Rain or LE RH100 models.These results suggest that the selection of the LE Rain or LE RH100 model should be based on data availability.Future work should seek to measure the thermal conductivity, albedo, and surface roughness, which would allow the differences between models to be evaluated.Furthermore, detailed knowledge of the debris properties, including how the thermal conductivity and water content vary with depth, would allow the performance of these models to be compared to other debris-covered glacier energy balance models (Collier et al., 2014;Fujita and Sakai, 2014).

Conclusions
Debris thickness greatly impacts ablation rates on debriscovered glaciers; however, incorporating debris cover into energy balance models is still hampered by a lack of knowledge of the debris properties.Fieldwork performed on Imja-Lhotse Shar Glacier over the 2014 melt season was used to develop new techniques to measure surface roughness, which yielded reasonable values for various grain sizes.Temperature sensors and ablation stakes installed in the debris were also used to assess the performance of a debris-covered glacier energy balance model using three different methods for estimating the latent heat flux.All three models performed well, as a result of the calibration procedure, which allowed variations in the lack of latent heat flux to be compensated for by adjusting the debris properties.However, the LE Rain and LE RH100 models yielded more reasonable values of latent heat fluxes.This suggests that in a data-scarce region either the LE Rain or LE RH100 model may be used if relative humidity or precipitation data are available.
A sensitivity analysis revealed that ablation rates were most sensitive to variations in thermal conductivity, followed closely by albedo and surface roughness.This highlights the importance of measuring the thermal conductivity and the moisture content in the debris.The effect of temporal resolution on model performance was also explored using a 6 h time step and a daily time step.The 6 h time step was found to underestimate the daily high each day, which caused the ablation rates to also be slightly underestimated.The daily time step did not model the daily average temperature as well, but yielded better estimates of ablation over the entire melt season.
Future studies should continue to work on incorporating the water content in the debris into debris-covered glacier energy balance models and determine its effect on thermal conductivity and the latent heat flux.Furthermore, an increased understanding of how the albedo may vary over the course of the day, the course of the melt season, and as a function of debris saturation, may significantly improve model performance.Lastly, the methods developed in this study have the potential to be scaled up such that maps of surface roughness on a whole glacier scale may be developed in the future, but it is imperative to determine the scale at which the surface roughness and aerodynamic roughness agree with one another.

Figure 1 .
Figure 1.Landsat 8 panchromatic image from 14 November 2014 of Imja-Lhotse Shar Glacier with the focus area of this study highlighted by the rectangular box several kilometers up-glacier from the terminus, and the site location within Nepal shown in the inset.

Figure 2 .
Figure 2. Sites A-D highlighting the variations in grain sizes that are found over the debris-covered portion of Imja-Lhotse Shar Glacier (cones 0.19 m diameter).

Figure 3 .
Figure 3. Transect from left to right of Site B showing the identification of obstacles (Obst) and their corresponding heights (h obst ) and depths (d obst ).

Figure 4 .
Figure 4.The effect of obstacle threshold (m) on both obstacle density and the estimate of surface roughness using the alternative method for Sites A-D and the DEMs with a resolution of 0.01 m.

Figure 5 .
Figure 5. Various plots for Site 11 using the LE Rain model showing (a) average daily temperatures at two depths (solid and dashed lines indicate measured and modeled temperatures, respectively), (b) average daily energy fluxes, (c) measured and modeled temperatures at a depth of 0.01 m over the calibration period, and (d) measured and modeled temperatures over the validation period.

Figure 6 .
Figure 6.Scatter plot of measured and modeled temperature for Site 11 at the surface for the LE Rain model showing the positive temperature bias overnight.

Figure 7 .
Figure 7. Modeled ablation with respect to debris thickness for all 15 sites from 18 May to 9 November 2014 for each of the three latent heat flux models, including the LE Rain model, with average values and the three measured stakes that did not exceed 1 m.

Figure 8 .
Figure 8. Modeled and measured surface temperature at Site 11 over the (a) calibration period and (a) validation period using 6 h data, and (a) entire period using daily averages.
, respectively.Site B had the highest value of z 0 (0.043 m) of the four sites in this study and consisted of larger boulders up to 1 m www.the-cryosphere.net/9

Table 1 .
Details of the debris thickness, topography, and monitoring equipment installed at each site.The use of italics notes an estimation of debris thickness.T s denotes surface temperature.

Table 2 .
Errors associated with the DEM for each site.

Table 3 .
Surface roughness (z 0 ) estimates and obstacle thresholds as a function of DEM resolution (m) at Sites A-D.

Table 4 .
Optimized values of albedo, thermal conductivity, and surface roughness for three methods of estimating latent heat flux during the calibration period.