Impact of debris cover on glacier ablation and atmosphere–glacier feedbacks in the Karakoram

The Karakoram range of the Hindu-Kush Himalaya is characterized by both extensive glaciation and a widespread prevalence of surficial debris cover on the glaciers. Surface debris exerts a strong control on glacier surface-energy and mass fluxes and, by modifying surface boundary conditions, has the potential to alter atmosphere– glacier feedbacks. To date, the influence of debris on Karakoram glaciers has only been directly assessed by a small number of glaciological measurements over short periods. Here, we include supraglacial debris in a high-resolution, interactively coupled atmosphere–glacier modeling system. To investigate glaciological and meteorological changes that arise due to the presence of debris, we perform two simulations using the coupled model from 1 May to 1 October 2004: one that treats all glacier surfaces as debris-free and one that introduces a simplified specification for the debris thickness. The basin-averaged impact of debris is a reduction in ablation of ∼ 14 %, although the difference exceeds 5 mw.e. on the lowest-altitude glacier tongues. The relatively modest reduction in basin-mean mass loss results in part from non-negligible sub-debris melt rates under thicker covers and from compensating increases in melt under thinner debris, and may help to explain the lack of distinct differences in recent elevation changes between clean and debriscovered ice. The presence of debris also strongly alters the surface boundary condition and thus heat exchanges with the atmosphere; near-surface meteorological fields at lower elevations and their vertical gradients; and the atmospheric boundary layer development. These findings are relevant for glacio-hydrological studies on debris-covered glaciers and contribute towards an improved understanding of glacier behavior in the Karakoram.


Introduction
The Karakoram region of the greater Himalaya (∼ 74-77 • E,34-37 • N; Fig. 1) is extensively glacierized, with an ice-covered area of ∼ 18 000 km 2 (Bolch et al., 2012).Supraglacial debris is widespread, and covers an estimated ∼ 18-22 % of the glacierized area (Scherler et al., 2011;Hewitt, 2011), a fraction that is approximately twice as large as the Himalayan average of ∼ 10 % (Bolch et al., 2012).The region has received a great deal of public and scientific attention in recent years due to evidence of stable or even slightly positive mass balances in the 2000s (Hewitt, 2005;Scherler et al., 2011;Gardelle et al., 2012Gardelle et al., , 2013;;Kääb et al., 2012) that are in contrast with predominantly negative balances of glaciers in the rest of the Hindu-Kush Himalaya (HKH; Cogley, 2011) Knowledge of the hydrological response of Karakoram glaciers to climate change is critical, since their meltwater contributes to freshwater resources in this highly populated region of South Asia (Kaser et al., 2010;Lutz et al., 2014).However, due to logistical constraints and political instability, field observations of glaciological and meteorological conditions in the Karakoram are sparse in space and time, in particular at high altitudes (Mihalcea et al., 2006(Mihalcea et al., , 2008a;;Mayer et al., 2014).Although observational records have been supplemented in recent decades by remote-sensing data (e.g., Gardelle et al., 2012Gardelle et al., , 2013;;Kääb et al., 2012), large gaps remain in our understanding of the important drivers of glacier change in this region, including regional atmospheric conditions, local topography and glacier debris cover, as well as interactions between them.Physically based numerical modeling has the potential to supplement observations and provide additional insight into contemporary glacier dynamics as well as to provide a methodology for predictions of future glacier response.
The prevalence of debris cover has a strong potential influence on glacier behavior in the Karakoram, as field studies have shown that debris cover can significantly alter the ice ablation rate compared to that of clean ice (e.g., Østrem, 1959;Fujii, 1977;Inoue and Yoshida, 1980).Ice melt is enhanced beneath debris cover less than a few centimeters thick, due to increased absorption of solar radiation.Conversely, ice ablation decreases exponentially as the thicknesses increases above this depth, due to insulation of the ice from atmospheric energy sources.Surficial debris also drastically alters glacier surface conditions, by permitting the surface temperature to exceed the melting point and by modifying the surface roughness and saturation conditions, which impacts the surface-energy fluxes (Inoue and Yoshida, 1980;Takeuchi et al., 2001;Brock et al., 2010) and the atmospheric boundary layer (Granger et al., 2002).Therefore, there is a strong potential for debris-covered ice to affect atmosphere-glacier feedbacks in this region.
Two main issues arise in attempting to include the influence of debris cover in simulations of Karakoram glaciers.First, the debris thickness, extent and thermal properties are largely unknown and their specification is highly uncertain.Second, the spatial distribution of meteorological forcing data is complicated by highly heterogeneous surface conditions in the ablation zone (e.g., Nicholson and Benn, 2012) and the complex topography, with current approaches that use elevation-based extrapolation appearing to be inadequate (Reid et al., 2012).Here, we investigate the influence of debris cover on Karakoram glacier surface-energy and mass exchanges and feedbacks with the atmosphere over an ablation season, using an interactively coupled atmosphere and glacier climatic mass balance (CMB) model that includes debris cover.By comparing a debris-free simulation to a simulation where we include debris cover with a simple specification of thickness, we first quantify differences in the surface energy balance and mass fluxes.We then assess feedbacks between the atmosphere and glacier surfaces using the coupled model and differences in boundary layer development.

Methods
The modeling tool employed in this study is the interactively coupled high-resolution atmosphere and glacier climatic mass balance model WRF-CMB, which explicitly resolves the surface-energy and CMB processes of alpine glaciers at the regional scale (Collier et al., 2013).The coupled model has been previously applied to the study region neglecting debris cover and was capable of reproducing the magnitudes of the few available observations of glacier CMB in this region.The changes introduced to the atmospheric and glacier-CMB model components for this study are described in Sect.2.1 and 2.2, respectively.We compare two WRF-CMB simulations for the period of 1 May to 1 October 2004: the first treated all glacier surfaces as debris-free (CLN) and the second introduced a simplified debris thickness specification (DEB), which is described in Sect.2.3.

Regional atmospheric model
The atmospheric component of WRF-CMB is the Advanced Research version of the Weather Research and Forecasting (WRF) model version 3.6.1 (Skamarock and Klemp, 2008).In this study, WRF was configured with three nested domains, of 30, 10 and 2 km resolution, which were centered over the Karakoram region (Fig. 1).The domains had 40 vertical levels, with the model top located at 50 hPa.For these simulations, debris cover is introduced in the 2 km domain only, since it provides the best representation of both the complex topography and glacier extents.
The model configuration was based on the previous application of WRF-CMB over this region (Collier et al., 2013, Table 1).However, for this study, the land surface model was updated to the Noah-MP (multiparameterization) scheme (Niu et al., 2011), which provides an improved treatment of snow physics in non-glacierized grid cells compared with the Noah scheme (Chen and Dudhia, 2001) that was previously used, by prognosing the energy balance and skin temperature of the vegetation canopy and snowpack separately, introducing multiple layers in the snowpack, and providing an improved treatment of frozen soils.Note that the prognosis of surface and subsurface conditions for glacierized grid cells is performed by the CMB model, which is discussed in the next section.The adaptive time stepping scheme was used, which greatly increased the execution speed of the simulations.Horizontal diffusion was also changed to be computed in physical space rather than along model levels, whereby diffusion acts on horizontal gradients computed using a verti-cal correction term rather than on the gradients on coordinate surfaces.We adopt this approach because it may be more accurate in complex terrain where the vertical levels are significantly sloped and because it provided a clear improvement in simulated precipitation in recent applications of WRF-CMB.Finally, for the finest-resolution domain (hereafter WRF D3), slope effects on radiation and topographic shading were accounted for and a cumulus parameterization was neglected, since at 2 km resolution, it is assumed to be convectionpermitting (e.g., Weisman et al., 1997).
The USGS land-cover data used by WRF were updated to incorporate more recent glacier inventories.Over the Himalayan region, we used the glacier outlines from the Randolph Glacier Inventory v. 3.2 (RGI; Pfeffer et al., 2014).For the Karakoram itself, we used the inventory of Rankl et al. (2014), which was obtained by updating the RGI manually on the basis of Landsat scenes.To determine which grid cells in each WRF domain were glacierized, the outlines were rasterized on a grid with a resolution that was 50 times higher than the original grid spacing of the domain.The fractional glacier coverage of grid cells was calculated on this finer grid, and a threshold of 40 % coverage was used to classify a grid cell as "glacier".The soil categories and vegetation parameters were also updated to be consistent with the glacier outlines.
The atmospheric model was forced at the boundaries of the coarse-resolution domain with the ERA-Interim reanalysis from the European Centre for Medium-Range Weather Forecasts (ECMWF; Dee et al., 2011).The spatial and temporal resolution of the ERA-Interim data are approximately 80 km (T255 spectral resolution) and 6 hourly, respectively.Snow depths in ERA-Interim over the Karakoram are unrealistic (more than 20 m; Collier et al., 2013), therefore an alternative initial snow condition was provided by the Global EASE-Grid 8 day blended SSM/I and MODIS snow cover data set for snow water equivalent (Brodzik et al., 2007), by assuming a snow density of 300 kg m −3 and specifying a depth of 0.5 m for areas with missing data, such as over large glaciers.This assumption affected 0.7, 5 and 40 % of grid points in D1-D3, respectively.We note that analysis of summer (June-July-August-September) mean fields over the Karakoram in ERA-Interim indicate that near-surface air temperatures were close to the 1979-2014 mean in 2004, while precipitation was significantly below average.

Glacier CMB model with debris treatment
The original basis of the glacier CMB model is the processbased model of (Mölg et al., 2008(Mölg et al., , 2009)).The model solves the full energy balance equation to determine the energy for snow and ice ablation.The computation of the specific column mass balance accounts for: surface and subsurface melt, refreezing and changes in liquid water storage in the snowpack, surface vapor fluxes and solid precipitation.The CMB model was adapted for interactive coupling with WRF by Collier et al. (2013) and modified to include supraglacial debris by Collier et al. (2014).For the version employed in this study, a time-varying snowpack is introduced on top of a static debris layer, both of which overly a column of ice resolved down to a depth of 7.0 m.The vertical levels in the subsurface used for these simulations are presented in Table 2.
A full description of the debris modifications is given by Collier et al. (2014), however we provide a brief summary here.The debris layer is resolved into 1 cm layers and has an assumed porosity function that decreases linearly with depth.The properties of each layer in the debris are computed as weighted functions of whole-rock values and the contents of the pore space (air, water or ice) using values presented in Table 3.For the whole-rock values, the albedo was based on 50 spot measurements on a debris-covered glacier in Nepal (Nicholson and Benn, 2012); the density and thermal conductivity were selected as representative values spanning major rock types taken from Daly et al. (1966); Clark (1966), respectively; and, the specific heat capacity was taken from Conway and Rasmussen (2000).Moisture in the debris and its phase are modeled using a simple reservoir parameterization.When debris is exposed at the surface, the surface vapor pressure is parameterized as a linear function of the distance between the surface and the saturated horizon.
Surface temperature is predicted using an iterative approach to determine the value that yields zero net flux in the surface energy balance equation.Initial test simulations with WRF-CMB over the Karakoram gave unrealistically low surface temperatures as a result of excessive nighttime damp-

Snow
variable Debris every 0.01 m Ice 0.1, 0.2, 0.3, 0.4, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 5.0, 7.0 m ing of the turbulent fluxes, in particular the sensible heat flux (QS) over debris-free glacier surfaces at high elevations.The stability corrections are based on the bulk Richardson number (specifically, those provided in Braithwaite, 1995) and have been used previously in glacier CMB modeling (e.g., Mölg et al., 2008Mölg et al., , 2009;;Reid et al., 2012).In the most stable conditions, the turbulent fluxes are fully damped, which resulted in decoupling of the surface and the atmosphere and excessive radiative cooling.Even in less stable conditions, the damping of modeled turbulent fluxes has been found to be excessive compared with eddy covariance measurements over glaciers (Conway and Cullen , 2013).Congruent with previous modeling studies of glacier surface-energy fluxes, we therefore limit the maximum amount of damping in stable conditions to 30 % (Martin and Lejeune, 1998;Giesen et al., 2009).In addition, we adopt a minimum wind speed of 1 m s −1 to be consistent with neighboring non-glacierized grid cells simulated by the Noah-MP LSM (land surface model; Niu et al., 2011).However, test simulations in early April indicate that the second correction has a minimal impact on wind speeds and turbulent fluxes in glacier grid cells and, thus, may be unnecessary.
To prevent errors arising from blended snow and debris layers, such as constraints on possible temperature solutions or excess melting, an adaptive vertical grid in the snowpack was introduced.For snow depths of up to 1 m, the nearest integer number of 10 cm layers are assigned, while areas of the snowpack that exceed one meter are resolved into the nearest integer number of 50 cm layers.Snow depths between 1 and 10 cm are assigned a single computational layer, and depths less than 1 cm are not treated with a unique layer.Normalized linear interpolation is performed to calculate temperature changes over regions of the snowpack where the layer depths have changed.This procedure conserves the bulk heat content of the snowpack, except when the depth crosses the minimum threshold of 1 cm.In both simulations, timestep changes in the bulk heat content of the snowpack in WRF D3 were small (less than 0.01 K).The CMB model is not designed for detailed snowpack studies and therefore only prognoses a bulk snow density.Since the total snow depth is not modified by the interpolation scheme, snow mass is conserved.
The debris-free version of the CMB model normally has levels located at fixed depths in the subsurface, with the thermal and physical properties of each layer computed as a weighted average of the snow and ice content.However, to isolate the influence of debris on glacier energy and mass  (Brock et al., 2010) fluxes, the CLN simulation also employs the adaptive vertical grid in the snowpack in this study.A test simulation from 1 April to 1 May 2004 was performed to compare the two adaptive and non-adaptive grids, with reasonable agreement in daily mean simulated snow depth (R 2 = 0.99; mean deviation, MD = −1.9cm) and snow melt (R 2 = 0.87; MD = 6.8 × 10 −4 kg m −2 ).

Specification of debris extent and thickness in WRF D3
The RGI and the inventory of Rankl et al. (2014) provide glacier outlines that include debris-covered glacier areas when detected, but they do not delineate these areas.To define debris-covered areas in WRF D3, the clean ice/firn/snow mask of Kääb et al. (2012) was rasterized on the same highresolution (40 m) grid used to compute glacierized grid cells (cf.Sect.2.1) For each WRF pixel in D3, the percent coverage of debris was determined and the same threshold of 40 % was used to classify a glacier pixel as debris-covered.Figure 2a provides an example of the delineation for the Baltoro glacier 76 • 26 E, 35 • 45 N).We note that any debris-covered glacier areas that are not detected during the generation of the glacier outlines are missed.Specifying the debris thickness was more complex, since this field varies strongly over small spatial scales.For example, Nicholson and Benn (2012) reported very heterogeneous debris thicknesses on the Ngozumpa glacier, Nepal, varying between 0.5 and 2.0 m over distances of less than 100 m.Spatial variability arises from many factors, including hillslope fluxes to the glacier; surface and subsurface transport and, the presence of ice cliffs, melt ponds and crevasses (e.g., Brock et al., 2010;Zhang et al., 2011).The few available field measurements do not support a relationship between debris thickness and elevation (e.g., Mihalcea et al., 2006;Reid et al., 2012).However, measurements on the Tibetan Plateau (Zhang et al., 2011) in Nepal (Nicholson and Benn, 2012), and in the Karakoram (Mihalcea et al., 2008a) indicate that thicker values are more prevalent near glacier termini, while thinner ones are more ubiquitous up-glacier.
In this study, we adopt a simple linear approach that was informed by this observed relationship to specify debris thickness over the areas identified as debris-covered in WRF D3.For this method, the distance down-glacier was com- puted, starting from the top of the debris-covered area of each glacier and moving along its centerline (Fig. 2b, d).Centerline data were provided by Rankl et al. (2014) for both main glacier trunks and their tributaries.We then assumed a fixed gradient to distribute debris over areas identified as being debris-covered as a function of distance down-glacier, with a single thickness specified in each 2 km grid cell.We tested two gradients, 1.0 and 0.75 cm km −1 , which gave thicknesses exceeding 40 and 30 cm, respectively, at the termini of the longest glaciers in the Karakoram (thicknesses derived using the 0.75 cm km −1 gradient are summarized in Fig. 2c).Where centerline information was unavailable (i.e., outside of the black contour in Fig. 2d), a constant thickness of 10 cm was assigned to each debris-covered pixel.For clarity, these data are not included in Fig. 2c.
Both gradients are consistent with the ASTER-derived debris-thickness data for the Baltoro glacier of Mihalcea et al. (2008a) after averaging the data onto the WRF-D3 grid.However, these data show a nonlinear increase near the terminus, and indicate that the 1 cm km −1 gradient distributes too much debris in the middle ablation zone, while the 0.75 cm km −1 value distributes too little near the terminus.Here, we focus our discussion on the 0.75 cm km −1 gra-dient simulation and suggest that our analysis thus represents a conservative estimate of the impact of debris.However, since the nonlinear increase is located close to the terminus, we assume the lower gradient is most valid at the regional scale.For comparison, the 1 cm km −1 gradient decreases the basin-mean mass loss between 1 July and 1 October by a further ∼ 4 % compared with the lower value.
This approach underestimates peak thicknesses at the termini of the Baltoro, which exceed 1 m (e.g., Mihalcea et al., 2008a).However, it is well established that ablation decreases exponentially with debris thicknesses above a few centimeters (e.g., Østrem, 1959;Loomis, 1970;Mattson et al., 1993).As the debris layer is resolved into 1 cm layers, including debris depths of up to 1 m would therefore greatly increase the computational expense of the CMB model, with likely only a small change to the amount of sub-debris ice melt.In addition, features such as meltwater ponds and ice cliffs in the ablation zone absorb significantly more energy than adjacent debris-covered surfaces.These features may give compensatory high melt rates (e.g., Inoue and Yoshida, 1980;Sakai et al., 1998Sakai et al., , 2000;;Pellicciotti et al., 2014;Immerzeel et al., 2014a) that support using a thinner average or "effective" debris thickness when assigning an average value to each 2 km grid cell in WRF D3.
After applying this method, WRF D3 contains a total of 5273 glacierized grid cells, 821 of which are debris-covered glacier cells, which gives a proportion of debris-covered glacier area in WRF D3 of ∼ 16 %.
In the following section, we evaluate and compare the DEB and CLN simulations, often focusing on altitudinal profiles where variables are averaged in 250 m elevational bins.Note that for these profiles, there are only 3 (9) glacierized grid cells present below 3250 m a.s.l.(above 7000 m), compared with at least 17 and up to 1100 grid cells in between these altitudes.In addition, when computing basin-averaged quantities, we excluded the bordering 10 grid points in WRF D3 (5 of which are specified at the boundaries).

Land surface temperature
For model evaluation, we compared simulated daytime land surface temperature (LST) with daily fields from the MODIS Terra MOD11A1 and Aqua MYD11A1 data sets, which have spatial resolutions of 1 km.Only MODIS data with the highest quality flag were used for the comparison and WRF-CMB data were taken from the closest available time step in local solar time.We focussed on daytime LST, because this field had a higher number of valid pixels at lower elevations over the simulation period than nighttime LST. Figure 3a shows mean elevational profiles of LST over glacierized pixels for composite MODIS data and for the CLN and DEB simulations.Although both modeled profiles are lower than in MODIS, the simulated profile in DEB is in much closer agreement than CLN, as mean LST exceeds the melting point below ∼ 5100 m a.s.l.
Examination of the MODIS LST data suggests that they may contain a positive bias, as a result of blending of different glacier surface types as well as glacierized and nonglacierized areas on the 1 km resolution grid.For example, Figure 3b shows an example of MODIS Terra LST on 5 August 2004 around the Baltoro glacier, a time slice that was selected for the low number of missing values in this region.MODIS exceeds the melting point over most of the glacier, including over smaller, largely debris-free tributary glaciers, due to blending with valley rock walls.The data are also higher over glacier areas with debris-covered fractions that fall below the threshold of 40 % used to define a WRF pixel as a debris-covered (cf.Fig. 2b).Therefore, the binary definition of debris-free and debris-covered glacier surface types, as well as inaccuracies in the glacier mask, also likely contribute to lower LST in WRF-CMB.
To examine temporal variations, Fig. 3c shows a time series of LST for July and August 2004 from all three data sets at a pixel on the Baltoro glacier tongue, which is denoted by a black circle in Fig. 3b.This pixel was selected since it falls within the glacier outline on the MODIS grid and because the debris coverage in 2004 appears to be 100 % (cf.Fig. 2 in Mihalcea et al., 2008a).The variability in LST at this point is well captured in DEB, including days with maxima exceeding ∼ 30 • C and higher or lower periods, while as expected, CLN greatly underpredicts LST and its variability.

Glacier surface energy and climatic-mass-balance dynamics
The basin-mean cumulative glacier CMB for both simulations is shown in Figure 4a.The month of May is characterized by basin-mean accumulation (Fig. 4b), consistent with the findings of Maussion et al. (2014) of the importance of spring precipitation in this region.On average, the melt season lasts from approximately mid-June until mid-September, over which period more than ∼ 90 % of grid cells categorized as debris-covered are exposed.As a result, there is a significant decrease in net ablation, as is discussed at the end of this section.Note that the basin-averaged CMB during summer is less negative than in a previous debris-free model run (Collier et al., 2013), which is primarily due to increased precipitation as a result of changing the atmospheric diffusion scheme (Table 1) through the albedo effect.The decrease in ablation is likely an improvement, since the previous estimate showed a negative bias in comparison with in situ glaciological measurements.To isolate the impacts of debris, we focus our analysis on the period of 1 July to 15 September 2004, when more than 35 % of debris pixels are exposed on average over the Karakoram.The basin-mean vertical balance profile indicates that between 1 July and 15 September 2004, the zero-balance altitude is located at ∼ 5700 m a.s.l.(Fig. 5).For comparison, annual equilibrium line altitudes (ELAs) in the Karakoram are estimated to range from 4200 to 4800 m (Young and Hewitt, 1993).We note that the absence of avalanche accumulation in our simulations may contribute to an overestimate of the zero-balance altitude, as this process is regionally important and produces ELAs that are often located hundreds of meters below the climatic snowline (e.g., Benn and Lehmkuhl, 2000;Hewitt, 2005Hewitt, , 2011)).Below ∼ 5700 m, there is a ∼ 18 % reduction in total ablation in DEB compared with CLN (of 5.3 m w.e.), which we anticipate represents an underestimate due to the nonlinear debris thickness observed near the Baltoro terminus.
The presence of surface debris has a noticeable impact on basin-mean surface-energy fluxes between 1 July and 15 September 2004 (Table 4).Elevational profiles reveal even stronger impacts in the ablation areas, as the number of grid cells with exposed debris increases towards lower elevations (Fig. 6a,b creases due to the lower surface albedo, while net longwave radiation (LWnet) becomes more negative due to stronger emission by warmer debris surfaces.The turbulent flux of sensible heat becomes a smaller energy source or even sink, while that of latent heat (QL) becomes slightly more negative.The conductive heat flux (QC) transitions from a small energy gain in CLN to a strong sink in DEB, due to solar heating of the debris, and extracts nearly twice as much energy from the surface as LWnet at the lowest glacierized elevations.Finally, both penetrating shortwave radiation (QPS) and the energy available for surface melt (the residual of the surface-energy budget; QM) decrease strongly towards lower elevations in DEB, as the overlying snow cover goes to zero, while in CLN theses fluxes provide strong energy sinks.
As a result of these changes to the surface-energy dynamics, total-column melt decreases by ∼ 18 % below 5000 m (Fig. 6c), with the small difference above this elevation reflecting overlying snow cover and some compensating increases in melt under thinner debris, which are prevalent (cf.Fig. 2c).The other mass fluxes are not strongly affected (Table 4; Fig. 6c).While surface vapor fluxes are small when spatially and temporally averaged, they represent a nonnegligible mass flux in total, with ∼ 1.8 × 10 5 kg of sublimation and 2.0 × 10 4 kg of deposition at snow and ice surfaces.Vapor exchange between the debris and the atmosphere also totals −1.0 × 10 4 kg over the simulation period.
Simulated daily mean ablation (corresponding to subdebris-ice and total-column values in DEB and CLN, respectively) shows a general decrease with both topographic height and increase in debris thickness (Fig. 7).Although melt rates below 3500 m have been estimated to be small due to insulation by thick debris cover (Hewitt, 2005), our results suggest that appreciable rates, of up to ∼ 2 cm w.e.day −1 occur under the thickest layers at lower elevations.For the thinnest debris layers (of a few centimeters), ablation is enhanced in DEB compared with CLN.Simulated values are consistent with the few available field measurements of glacier ablation in this region.For example, Mayer et al. (2010) reported rates of ∼ 2 to 14 cm w.e.day −1 under debris covers of ∼ 1 to 38 cm on the Hinarche glacier ( 74• 43 E, 36 • 5 N) in 2008.Mihalcea et al. (2006) reported rates of 1-6 cm w.e.d −1 on the Baltoro glacier in 2004 over elevations of ∼ 4000-4700 m and thicknesses of 0 to 18 cm, and the modeled melt rates over a similar period compare well with their Østrem curve (cf.their Fig. 7).
A spatial plot of the total cumulative mass balance in DEB delineates regions of glacier mass gain and loss in the Karakoram (Fig. 8a).Accumulation is higher in the western part of the domain, where more precipitation falls over the simulation period (not shown).Differences between DEB and CLN are small over most of the domain, with the exception of lower-altitude glacier tongues where differences exceed 2.5 m w.e. (Fig. 8b).The strong decrease in mass loss in these areas changes the cumulative basin-mean mass balance on 15 September from −919 kg m −2 in CLN to −831 in DEB.Considering the whole simulation period, the basinmean values are −856 kg m −2 in CLN and −737 in DEB (a reduction of ∼ 14 %) on 1 October 2004, with differences exceeding 5 m w.e. on the lowest debris-covered tongues.

Atmosphere-glacier feedbacks
The total number of hours for which the surface temperature reaches or exceeds the melting point ranges from more than 1500 at low-altitude glacier termini to less than 50 above ∼ 6400 m (Fig. 9a).The presence of debris results in up to 700 additional hours with surface temperatures above 273.15K in DEB compared with CLN (Fig. 9b), which provide a strong heat flux to the atmosphere.Considering all hours between 1 July and 15 September, an extra 3.5 × 10 7 W of energy is transferred to the atmosphere in DEB by the sensible heat flux.
The change in surface boundary conditions produces higher basin-mean near-surface air temperatures, of up to 2-3 K at the lowest glacierized elevations (Fig. 10a), consistent with observations of higher air temperatures over debriscovered glacier areas during the ablation season (Takeuchi et al., 2000(Takeuchi et al., , 2001;;Reid et al., 2012).The vertical gradient in 2 m air temperature below 5000 m is more than one degree higher in DEB than CLN (−0.0074 compared with −0.0062K m −1 ; ∼ −0.0073 for both simulations above this elevation).Basin-mean accumulated precipitation ranges from 50 to 175 mm w.e.below 5000 m and increases approximately linearly with elevation above this level.The area-averaged differences between CLN and DEB are very small, with a slight decrease (increase) at the lowest (highest) elevations in www.the-cryosphere.net/9/1617/2015/The Cryosphere, 9, 1617-1632, 2015 Here, "ablation" refers to sub-debris ice melt in DEB (i.e., only snow-free pixels are selected) and total column melt (surface and englacial) in CLN for the same pixels and time periods.The concentration of data points at 10 cm thickness results from the specification of debris where centerline information was unavailable (cf.Sect.2.3) DEB, consistent with warmer and thus less humid conditions that contribute to slower cooling and saturation of air moving upslope and a shift of surface precipitation up-glacier.The simulated frozen fraction increases approximately linearly from 0 % below 3250 m to more than 90 % above ∼ 5500 m (not shown).These results are consistent with estimates of annual precipitation, which indicate that valleys are drier and precipitation increases up towards accumulation areas, and with previously reported frozen fractions (Winiger et al., 2005;Hewitt, 2005).Finally, higher surface roughness values over debris result in a decrease of near-surface horizontal wind speeds at lower elevations (Fig. 10c).It is noteworthy that changes in atmosphere-glacier feedbacks due to the presence of surface debris also help to drive the differences in observed ablation (cf. Figs. 6,7).
Figure 11 illustrates alterations to the diurnal cycles of the turbulent flux of sensible heat (QS), the planetary boundary layer (PBL) depth, and the along-glacier component of the near-surface winds over exposed debris pixels in DEB and their equivalents in CLN.Solar heating of the debris surface drives a strongly negative daytime QS in DEB (Fig. 11a) which reduces the stability of the glacier surface layer and enhances turbulent mixing.Peak negative QS values in DEB exceed −200 W m −2 , consistent with eddy covariance measurements of this flux over supraglacial debris (Collier et al., 2014).In comparison, QS in CLN is approximately 1 order since of magnitude smaller and positive.As a result of energy transfer by QS, a deep convective mixed layer develops in DEB, with the mean PBL height reaching nearly 1.5 km in the afternoon compared with only a couple of hundred meters in CLN (Fig. 11b).Finally, near-surface along-glacier winds in DEB are primarily anabatic during the day (directed up-glacier, which is defined here as positive) and katabatic during the evening and early morning (down-glacier and negative; Fig. 11c), compared with sustained katabatic flows (glacier winds) in CLN, resulting from cooling of the air near the ice surface, which is constrained at the melting point (e.g., van den Broeke, 1996).The absence of daytime katabatic flows over debris-covered areas is consistent with the findings of Brock et al. (2010).

Discussion and conclusions
In this study, surficial debris was introduced to the coupled atmosphere-glacier modeling system, WRF-CMB.The model provides a unique tool for investigating the influence of debris cover on both Karakoram glaciers and atmosphereglacier interactions in an explicitly resolved framework.The first-order impact of debris was estimated, with thickness determined using a fixed gradient of 0.75 cm km −1 with distance down-glacier in debris-covered areas, focusing on the period of 1 July to 15 September 2004 when more than 35 % of debris-covered pixels were exposed.The findings presented in this study have important implications for glaciohydrological studies in the Karakoram, as they confirm that neglecting supraglacial debris will result in an overestimation of glacier mass loss during the ablation season, of ∼ 14 % over the region and exceeding 5 m w.e. at the lowest elevations.In addition, exposed debris alters near-surface meteorological fields and their elevational gradients, which are often key modeling parameters used to extrapolate forcing data from a point location (e.g., an automatic weather station) over the rest of the glacier surface (e.g., Marshall et al., 2007;Gardner et al., 2009;Reid et al., 2012).The lapse rate in air temperature at lower elevations is more than one degree steeper in DEB, as a result of surface temperatures exceeding the melting point and a higher net turbulent transfer of sensible heat to the atmosphere that produces higher near-surface air temperatures, and is higher than values reported for smaller debris-covered glaciers (Reid et al., 2012) and in the eastern Himalaya, where the monsoon circulation is more dominant (Immerzeel et al., 2014b).Finally, we showed that debris induces significant alterations to the eling approach on the Baltoro glacier, and with the measured rates reported by Mayer et al. (2010).The authors of the latter study suggest the mechanism is more efficient heat transfer in the debris in the presence of moisture during the ablation season despite its thickness.In this study, mean icemelt rates for pixels with debris thickness exceeding 20 cm show some correlation with the debris moisture content (water: R 2 = 0.3; ice: R 2 = −0.69).However, our results suggest near-surface air temperature (R 2 = 0.91) is a stronger driver on average of the melt rates simulated below thick debris.The interactive nature of the simulation may permit a positive feedback mechanism, in which higher surface temperatures over thicker debris transfer energy to the atmosphere, in turn promoting higher air temperatures and further melt.Even when the air temperature is below 0 • C, energy conduction when the debris surface temperature exceeds this threshold likely also contributes to sub-debris ice melt, which is supported by our simulations.In combination with surface heterogeneity in the ablation zone (e.g., the presence of meltwater ponds and ice cliffs) and recent changes in ice-flow velocities (Quincey et al., 2009;Scherler and Strecker , 2012), both the simulated melt rates under thicker debris and enhanced melt under thinner debris help to explain the lack of significant differences in recent elevation changes between debris-free and debris-covered glacier surfaces in the Karakoram (Gardelle et al., 2013).
In surface energy balance studies of supraglacial debris, the latent heat flux is often neglected where measurements of surface humidity are unavailable, due to the complexity of treating the moist physics of debris.In the DEB simulation, the latent heat flus over exposed debris was non-negligible and primarily negative; furthermore, it contributed to a vapor loss that comprised 5.5 % of the total considering all glacierized pixels.Thus, our study suggests that neglecting QL and surface vapor exchange may be inappropriate assumptions, even for basin-scale studies.We further note that the simple parameterization developed for QL tended to underestimate the vapor-pressure gradient in the surface layer (Collier et al., 2014), suggesting that the importance of QL is underestimated in this study.However, the treatments of QL and the debris moisture content represent key sources of uncertainty in our simulations, since (i) they were developed in a different region and (ii) these fields impact sub-debris icemelt rates (Collier et al., 2014) but are not well measured or studied.
The alterations to the glacier energy and mass fluxes and to atmosphere-glacier interactions presented in this study are based on the ablation season of 2004 only and are sensitive to the debris thickness field, with small adjustments to the thickness gradient, resulting in significant changes in basin-mean glacier CMB.The gradient was consistent with ASTER-derived thickness data on the Baltoro glacier (Mihalcea et al., 2008a) except close to the terminus.However, our approach results in peak thicknesses of less than ∼ 15 cm on glaciers less than 20 km in length, while other studies in the Himalaya and elsewhere have reported much higher depths on glaciers of similar lengths (e.g., Mihalcea et al., 2008b;Rounce and McKinney, 2014).Thus, the impact on glacier ablation that we reported likely represents an underestimate, due to nonlinear effects near termini and the likely presence of steeper thickness gradients on shorter glaciers.Additional sources of uncertainty in our results include (i) the temporal discrepancy between our study period and the clean snow/ice mask of Kääb et al. (2012) used to delineate debris-covered areas, which was generated using Landsat data from the year 2000; and, (ii) our binary assignment of surface types as "debris-covered" or "debris-free" using a 40 % threshold.
There have been numerous recent efforts to more precisely determine debris thickness fields using satellite-derived surface temperature fields (e.g., Suzuki et al., 2007;Mihalcea et al., 2008a;Foster et al., 2012;Brenning et al., 2012), which is an appealing solution due to the wide spatial and temporal coverage of remote-sensing data.However, none of these studies have successfully reproduced field measurements without using empirically determined relationships or calibration factors (Mihalcea et al., 2008a;Foster et al., 2012).These methods are therefore best suited for debriscovered glaciers for which the necessary measurements to compute the relationships or factors are available, and their applicability for regional-scale studies such as this one is uncertain.Thus, important future steps for glacier CMB studies in the Karakoram include increasing the accuracy and spatial detail of the debris thickness field and its physical properties; improving our understanding of moisture fluxes between the debris and the atmosphere and accounting for subgrid-scale surface heterogeneity (e.g., by introducing a treatment of ice cliffs; Reid and Brock, 2014).Nonetheless, by providing an estimate of the controlling influence of debris, these simulations contribute to a greater understanding of glacier behavior in the Karakoram.

Figure 1 .
Figure 1.Topographic height shaded in units of km for (a) all three model domains in WRF-CMB, which are centered over the Karakoram and configured with grid spacings of 30, 10 and 2 km, and (b) a zoom-in of the finest resolution domain, WRF D3.

Figure 2 .
Figure 2. (a) Debris-covered (gray) and debris-free (blue) glacier areas, calculated on a 40 m grid for the Baltoro glacier and surrounding areas.The distance down-glacier over debris-covered areas, which is multiplied by a fixed gradient to map debris, is shown for (b) the Baltoro glacier and (d) the entire WRF-D3 region.In (d), the black contour delineates the region where centerline information was available from (Rankl et al., 2014).Outside of this region, distance down-glacier was not computed and debris-covered areas are shaded in gray to indicate these data are missing.(c) A box plot of debris thickness values, assuming a fixed gradient of 0.75 cm km −1 and averaging in 250 m elevation bins over WRF D3.The thick blue and thin black lines indicate the mean and median thicknesses in each bin.The total number of debris-covered pixels is given as a text string at the upper end of the range.

Figure 3 .
Figure 3. (a) Mean elevational profiles of daytime land surface temperature (LST) from DEB, CLN and composite MODIS Terra MOD11A1/Aqua MYD11A1 data sets, averaged from 1 June to 1 September 2004 and in 250 m elevation bins over glacierized pixels in WRF.(b) A sample time slice of MODIS Terra LST from 5 August 2004 on its native grid, overlaid on the Baltoro glacier outline and debris-covered area.(c) Time series of LST from 1 July to 1 September 2004 from the same data sets as in panel (a), taken from a pixel on the Baltoro tongue, which is denoted by a black circle in panel (b).The unit for all plots is • C.

Figure 4 .
Figure 4. Time series of (a) basin-mean cumulative glacier CMB in kg m −2 and (b) the daily maximum percentage of debris pixels that are exposed in DEB, for DEB (black curve) and CLN (gray) over the whole simulation period of 1 May to 1 October 2004.

Figure 5 .
Figure5.The cumulative vertical balance profile, averaged in 250 m elevation bins between 3000 and 7500 m a.s.l., over all glacier pixels and from 1 July to 15 September 2004.Solid gray circle markers denote results from the CLN simulation, while those from DEB are plotted with black markers.The shape of the black marker indicates the range of the mean debris thickness in that elevational band.

Figure 6 .Figure 7 .
Figure 6.(a)The percentage of debris-covered pixels in each 250 m elevation bin that are exposed on average between 1 July and 15 September 2004.Minimum and maximum values over the same period are indicated by gray shading.Elevational profiles of mean glacier (b) surface-energy and (c) mass fluxes, in units of W m −2 and kg m −2 respectively, with the solid (dashed) lines denoting data from DEB (CLN).For (c), evaporation and condensation in DEB are not shown as their profiles are approximately zero (less than 0.02 kg m −2 for all elevational bands).Note that these profiles correspond to an amalgamation of all glacierized grid cells, rather than the mean elevational profile along glacier.

Figure 8 .
Figure 8. Cumulative CMB in [kg m −2 ] between 1 July and 15 September 2004 for (a) the DEB simulation and (b) the difference between DEB and CLN.The black contour delineates the region where centerline information was available from Rankl et al. (2014).

Figure 9 .
Figure 9.The number of hours where the surface temperatures reaches or exceeds the melting point between 1 July and 15 September for (a) the DEB simulation and (b) the difference between CLN and DEB.The black contour delineates the region where centerline information was available from Rankl et al. (2014).

Figure 10 .Figure 11 .
Figure 10.Elevational profiles of near-surface (a) air temperature [K] and (b) accumulated precipitation [mm w.e.], and (c) wind speed at a height of 10 m [m s −1 ], from the DEB (black circle markers) and CLN (gray squares) simulations between 1 July and 15 September 2004.

Table 3 .
Physical properties in the CMB model.

Table 4 .
Mean glacier surface-energy and climatic-mass fluxes.