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

Abstract

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 $\sim 14$ %, although the difference exceeds 5 $\mathrm{m}\mathrm{w}.\mathrm{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 debris-covered 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.

Topographic height shaded in units of $\mathrm{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 omitted. See PDF]

Introduction

The Karakoram region of the greater Himalaya ($\sim 74$–77${}^{\circ }$ E,34–37${}^{\circ }$ N; Fig. ) is extensively glacierized, with an ice-covered area of $\sim \mathrm{18\hspace{0.17em}000}$ ${\mathrm{km}}^{\mathrm{2}}$ . Supraglacial debris is widespread, and covers an estimated $\sim 18$–22 % of the glacierized area , a fraction that is approximately twice as large as the Himalayan average of $\sim 10$ % . 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 that are in contrast with predominantly negative balances of glaciers in the rest of the Hindu-Kush Himalaya HKH; 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 . 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 . Although observational records have been supplemented in recent decades by remote-sensing data e.g.,, 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.,. 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 and the atmospheric boundary layer . 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., and the complex topography, with current approaches that use elevation-based extrapolation appearing to be inadequate . 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.

WRF configuration.

Model configuration
Horizontal grid spacing	30, 10, 2 $\mathrm{km}$ (domains 1–3)
Min/max time step	30/200, 10/60, 2/13 s
Vertical levels	40
Model top pressure	50 $\mathrm{hPa}$
Model physics
Radiation	CAM
Microphysics	Thompson
Cumulus	Kain–Fritsch (none in D3)
Atmospheric boundary layer	Yonsei University
Surface layer	Monin–Obukhov (revised MM5)
Land surface	Noah-MP
Dynamics
Top boundary condition	Rayleigh damping
Horizontal diffusion	Computed in physical space
Lateral boundaries
Forcing	ERA-Interim, T255 spectral resolution
	updated 6 hourly

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 . 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.  and  , 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. .

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 . 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. ). 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 Table . However, for this study, the land surface model was updated to the Noah-MP (multiparameterization) scheme , which provides an improved treatment of snow physics in non-glacierized grid cells compared with the Noah scheme 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 vertical 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 convection-permitting e.g.,.

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; ). For the Karakoram itself, we used the inventory of , 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;. 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 $\mathrm{m}$;, 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 , by assuming a snow density of 300 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$ and specifying a depth of 0.5 $\mathrm{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 process-based model of . 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 and modified to include supraglacial debris by . 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 $\mathrm{m}$. The vertical levels in the subsurface used for these simulations are presented in Table .

A full description of the debris modifications is given by , however we provide a brief summary here. The debris layer is resolved into 1 $\mathrm{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 . For the whole-rock values, the albedo was based on 50 spot measurements on a debris-covered glacier in Nepal ; the density and thermal conductivity were selected as representative values spanning major rock types taken from , respectively; and, the specific heat capacity was taken from . 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.

Subsurface layer depths.

Snow	variable
Debris	every 0.01 $\mathrm{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

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 damping 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 ) and have been used previously in glacier CMB modeling e.g.,. 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 . Congruent with previous modeling studies of glacier surface-energy fluxes, we therefore limit the maximum amount of damping in stable conditions to 30 % . In addition, we adopt a minimum wind speed of 1 $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ to be consistent with neighboring non-glacierized grid cells simulated by the Noah-MP LSM (land surface model; ). 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.

Physical properties in the CMB model.

Density ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$)
ice	915	–
whole rock	2700
water	1000	–
Specific heat capacity ($\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$)
air	1005	–
ice	2106	–
whole rock	750
water	4181	–
Thermal conductivity ($\mathrm{W}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{K}}^{-\mathrm{1}}$)
air	0.024	–
ice	2.51	–
whole rock	2.50
water	0.58	–
Surface roughness length ($\mathrm{m}$)
ice	0.001
debris	0.016
Albedo
ice	0.30
firn	0.55
fresh snow	0.85
debris	0.20
Emissivity
ice/snow	0.98	–
debris	0.94

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 $\mathrm{cm}$ layers are assigned, while areas of the snowpack that exceed one meter are resolved into the nearest integer number of 50 $\mathrm{cm}$ layers. Snow depths between 1 and 10 $\mathrm{cm}$ are assigned a single computational layer, and depths less than 1 $\mathrm{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 $\mathrm{cm}$. In both simulations, timestep changes in the bulk heat content of the snowpack in WRF D3 were small (less than 0.01 $\mathrm{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 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^{\mathrm{2}}$= 0.99; mean deviation, MD $=-$1.9 $\mathrm{cm}$) and snow melt ($R^{\mathrm{2}}$= 0.87; MD $=6.8\times {10}^{-\mathrm{4}}$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$).

Specification of debris extent and thickness in WRF D3

The RGI and the inventory of 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 was rasterized on the same high-resolution (40 m) grid used to compute glacierized grid cells (cf. Sect. ) 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 a provides an example of the delineation for the Baltoro glacier 76${}^{\circ }$ 26${}^{\prime }$ E, 35${}^{\circ }$ 45${}^{\prime }$ 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, reported very heterogeneous debris thicknesses on the Ngozumpa glacier, Nepal, varying between 0.5 and 2.0 $\mathrm{m}$ over distances of less than 100 $\mathrm{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.,. The few available field measurements do not support a relationship between debris thickness and elevation e.g.,. However, measurements on the Tibetan Plateau in Nepal , and in the Karakoram indicate that thicker values are more prevalent near glacier termini, while thinner ones are more ubiquitous up-glacier.

(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 . 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 $\mathrm{cm}{\mathrm{km}}^{-\mathrm{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 omitted. See PDF]

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 computed, starting from the top of the debris-covered area of each glacier and moving along its centerline (Fig. b, d). Centerline data were provided by 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 $\mathrm{cm}{\mathrm{km}}^{-\mathrm{1}}$, which gave thicknesses exceeding 40 and 30 $\mathrm{cm}$, respectively, at the termini of the longest glaciers in the Karakoram (thicknesses derived using the 0.75 $\mathrm{cm}{\mathrm{km}}^{-\mathrm{1}}$ gradient are summarized in Fig. c). Where centerline information was unavailable (i.e., outside of the black contour in Fig. d), a constant thickness of 10 $\mathrm{cm}$ was assigned to each debris-covered pixel. For clarity, these data are not included in Fig. c.

(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 ${}^{\circ }$C.

[Figure omitted. See PDF]

Both gradients are consistent with the ASTER-derived debris-thickness data for the Baltoro glacier of after averaging the data onto the WRF-D3 grid. However, these data show a nonlinear increase near the terminus, and indicate that the 1 $\mathrm{cm}{\mathrm{km}}^{-\mathrm{1}}$ gradient distributes too much debris in the middle ablation zone, while the 0.75 $\mathrm{cm}{\mathrm{km}}^{-\mathrm{1}}$ value distributes too little near the terminus. Here, we focus our discussion on the 0.75 $\mathrm{cm}{\mathrm{km}}^{-\mathrm{1}}$ gradient 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 $\mathrm{cm}{\mathrm{km}}^{-\mathrm{1}}$ gradient decreases the basin-mean mass loss between 1 July and 1 October by a further $\sim \mathrm{4}$ % compared with the lower value.

This approach underestimates peak thicknesses at the termini of the Baltoro, which exceed 1 $\mathrm{m}$ e.g.,. However, it is well established that ablation decreases exponentially with debris thicknesses above a few centimeters e.g.,. As the debris layer is resolved into 1 $\mathrm{cm}$ layers, including debris depths of up to 1 $\mathrm{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., 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 $\sim 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 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (above 7000 $\mathrm{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).

Results

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 $\mathrm{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 a 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 $\sim 5100$ $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{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 non-glacierized areas on the 1 $\mathrm{km}$ resolution grid. For example, Figure b 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. b). 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. c 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. b. 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 ). The variability in LST at this point is well captured in DEB, including days with maxima exceeding $\sim 30$ ${}^{\circ }$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 a. The month of May is characterized by basin-mean accumulation (Fig. b), consistent with the findings of 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 $\sim 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 , which is primarily due to increased precipitation as a result of changing the atmospheric diffusion scheme (Table ) 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.

Time series of (a) basin-mean cumulative glacier CMB in $\mathrm{kg}{\mathrm{m}}^{-\mathrm{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 omitted. See PDF]

The basin-mean vertical balance profile indicates that between 1 July and 15 September 2004, the zero-balance altitude is located at $\sim 5700$ $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{l}.$ (Fig. ). For comparison, annual equilibrium line altitudes (ELAs) in the Karakoram are estimated to range from 4200 to 4800 $\mathrm{m}$ . 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.,. Below $\sim 5700$ $\mathrm{m}$, there is a $\sim 18$ % reduction in total ablation in DEB compared with CLN (of 5.3 $\mathrm{m}\mathrm{w}.\mathrm{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 ). Elevational profiles reveal even stronger impacts in the ablation areas, as the number of grid cells with exposed debris increases towards lower elevations (Fig. a,b; cf. Fig c). Net shortwave radiation (SWnet) increases 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.

The cumulative vertical balance profile, averaged in 250 m elevation bins between 3000 and 7500 $\mathrm{m}\mathrm{a}.\mathrm{s}.\mathrm{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 omitted. See PDF]

As a result of these changes to the surface-energy dynamics, total-column melt decreases by $\sim 18$ % below 5000 $\mathrm{m}$ (Fig. c), 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. c). The other mass fluxes are not strongly affected (Table ; Fig. c). While surface vapor fluxes are small when spatially and temporally averaged, they represent a non-negligible mass flux in total, with $\sim 1.8\times {10}^{\mathrm{5}}$ $\mathrm{kg}$ of sublimation and $2.0\times {10}^{\mathrm{4}}$ $\mathrm{kg}$ of deposition at snow and ice surfaces. Vapor exchange between the debris and the atmosphere also totals $-1.0\times {10}^{\mathrm{4}}$ $\mathrm{kg}$ over the simulation period.

(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 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$ and $\mathrm{kg}{\mathrm{m}}^{-\mathrm{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 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{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 omitted. See PDF]

Simulated daily mean ablation (corresponding to sub-debris-ice and total-column values in DEB and CLN, respectively) shows a general decrease with both topographic height and increase in debris thickness (Fig. ). Although melt rates below 3500 $\mathrm{m}$ have been estimated to be small due to insulation by thick debris cover , our results suggest that appreciable rates, of up to $\sim \mathrm{2}$ $\mathrm{cm}\mathrm{w}.\mathrm{e}.{\mathrm{day}}^{-\mathrm{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, reported rates of $\sim \mathrm{2}$ to 14 $\mathrm{cm}\mathrm{w}.\mathrm{e}.{\mathrm{day}}^{-\mathrm{1}}$ under debris covers of $\sim \mathrm{1}$ to 38 $\mathrm{cm}$ on the Hinarche glacier (74${}^{\circ }$ 43${}^{\prime }$ E, 36${}^{\circ }$ 5${}^{\prime }$ N) in 2008. reported rates of 1–6 $\mathrm{cm}\mathrm{w}.\mathrm{e}.{\mathrm{d}}^{-\mathrm{1}}$ on the Baltoro glacier in 2004 over elevations of $\sim 4000$–4700 $\mathrm{m}$ and thicknesses of 0 to 18 $\mathrm{cm}$, and the modeled melt rates over a similar period compare well with their Østrem curve (cf. their Fig. 7).

Daily mean ablation rate versus debris thickness for DEB (circle markers) and CLN (horizontal line markers), with the range of topographic height value of each data point indicated by the color of the marker. 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 $\mathrm{cm}$ thickness results from the specification of debris where centerline information was unavailable (cf. Sect. )

[Figure omitted. See PDF]

A spatial plot of the total cumulative mass balance in DEB delineates regions of glacier mass gain and loss in the Karakoram (Fig. a). 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 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ (Fig. b). The strong decrease in mass loss in these areas changes the cumulative basin-mean mass balance on 15 September from $-$919 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ in CLN to $-$831 in DEB. Considering the whole simulation period, the basin-mean values are $-$856 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$ in CLN and $-$737 in DEB (a reduction of $\sim 14$ %) on 1 October 2004, with differences exceeding 5 $\mathrm{m}\mathrm{w}.\mathrm{e}.$ on the lowest debris-covered tongues.

Mean glacier surface-energy and climatic-mass fluxes.

Surface energy fluxes ($\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$)	DEB	CLN
net shortwave (SWnet)	153.2	149.9
net longwave (LWnet)	$-$85.6	$-$83.1
sensible heat (QS)	4.3	8.3
latent heat (QL)	$-$13.6	$-$13.8
conduction (QC)	3.0	20.7
penetrating SW (QPS)	$-$21.4	$-$28.5
precipitation (QPRC)	$\sim \mathrm{0}$	$\sim \mathrm{0}$
residual energy (QM)	38.9	53.6
Mass fluxes ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{2}}$)	DEB	CLN
total-column melt	$-$0.67	$-$0.71
snow refreeze	0.11	0.11
sublimation	$-$0.02	$-$0.02
deposition	$\sim \mathrm{0}$	$\sim \mathrm{0}$
evaporation	$\sim \mathrm{0}$	–
condensation	$\sim \mathrm{0}$	–
surface accumulation	0.07	0.07

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 $\sim 6400$ $\mathrm{m}$ (Fig. a). The presence of debris results in up to 700 additional hours with surface temperatures above 273.15 $\mathrm{K}$ in DEB compared with CLN (Fig. b), which provide a strong heat flux to the atmosphere. Considering all hours between 1 July and 15 September, an extra $3.5\times {10}^{\mathrm{7}}$ $\mathrm{W}$ of energy is transferred to the atmosphere in DEB by the sensible heat flux.

Cumulative CMB in [$\mathrm{kg}{\mathrm{m}}^{-\mathrm{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 .

[Figure omitted. See PDF]

The change in surface boundary conditions produces higher basin-mean near-surface air temperatures, of up to 2–3 $\mathrm{K}$ at the lowest glacierized elevations (Fig. a), consistent with observations of higher air temperatures over debris-covered glacier areas during the ablation season . The vertical gradient in 2 m air temperature below 5000 $\mathrm{m}$ is more than one degree higher in DEB than CLN ($-$0.0074 compared with $-$0.0062 $\mathrm{K}{\mathrm{m}}^{-\mathrm{1}}$; $\sim -0.0073$ for both simulations above this elevation). Basin-mean accumulated precipitation ranges from 50 to 175 $\mathrm{mm}\mathrm{w}.\mathrm{e}.$ below 5000 $\mathrm{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 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 $\mathrm{m}$ to more than 90 % above $\sim 5500$ $\mathrm{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 . Finally, higher surface roughness values over debris result in a decrease of near-surface horizontal wind speeds at lower elevations (Fig. c). 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. , ).

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 .

[Figure omitted. See PDF]

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. a) which reduces the stability of the glacier surface layer and enhances turbulent mixing. Peak negative QS values in DEB exceed $-$200 $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$, consistent with eddy covariance measurements of this flux over supraglacial debris . 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 $\mathrm{km}$ in the afternoon compared with only a couple of hundred meters in CLN (Fig. b). 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. c), 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.,. The absence of daytime katabatic flows over debris-covered areas is consistent with the findings of .

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 atmosphere–glacier interactions in an explicitly resolved framework. The first-order impact of debris was estimated, with thickness determined using a fixed gradient of 0.75 $\mathrm{cm}{\mathrm{km}}^{-\mathrm{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 glacio-hydrological 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 $\sim 14$ % over the region and exceeding 5 $\mathrm{m}\mathrm{w}.\mathrm{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.,. 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 and in the eastern Himalaya, where the monsoon circulation is more dominant . Finally, we showed that debris induces significant alterations to the atmospheric boundary layer development and along-glacier winds, through changes in the turbulent heat flux.

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

[Figure omitted. See PDF]

A comparison of the simulated diurnal cycle of (a) the turbulent flux of sensible heat, QS [$\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$]; (b) the planetary boundary layer (PBL) height [km]; and, (c) the along-glacier wind speed [from the lowest model level; $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$], which is positive for up-glacier flow. The data are averaged over exposed debris in DEB (black curve) and the corresponding grid cells in CLN (gray curve).

[Figure omitted. See PDF]

Simulated ice-ablation rates in DEB under thicker debris $\sim ($O(10 $\mathrm{cm}$)) at lower elevations are consistent with the findings of of non-negligible melt energy under debris covers exceeding 1 $\mathrm{m}$ using a degree-day modeling approach on the Baltoro glacier, and with the measured rates reported by . 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 ice-melt rates for pixels with debris thickness exceeding 20 $\mathrm{cm}$ show some correlation with the debris moisture content (water: $R^{\mathrm{2}}$ = 0.3; ice: $R^{\mathrm{2}}$ = $-$0.69). However, our results suggest near-surface air temperature ($R^{\mathrm{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 ${}^{\circ }$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 , 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 .

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 , 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 ice-melt rates 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 except close to the terminus. However, our approach results in peak thicknesses of less than $\sim 15$ $\mathrm{cm}$ on glaciers less than 20 $\mathrm{km}$ in length, while other studies in the Himalaya and elsewhere have reported much higher depths on glaciers of similar lengths e.g.,. 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 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.,, 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 . These methods are therefore best suited for debris-covered 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;. Nonetheless, by providing an estimate of the controlling influence of debris, these simulations contribute to a greater understanding of glacier behavior in the Karakoram.

Acknowledgements

The authors thank M. Rankl and C. Kienholz for providing glacier centerline information from . They also thank G. Diolaiuti for providing the debris thickness data of , whose field and laboratory work was completed as a part of the SHARE (Station at High Altitude for Research on the Environment) project and supported by EvK2CNR. E. Collier was supported in part by a Queen Elizabeth II graduate scholarship. L. Nicholson is funded by an Austrian Science Fund (FWF grant V309). F. Maussion acknowledges support by the Austrian Science Fund (FWF project P22443-N21). A. B. G. Bush acknowledges support from the Natural Sciences and Engineering Research Council, and the Canadian Institute for Advanced Research (Earth System Evolution Program). This work was partly carried out under the Collaborative Adaptation Research Initiative in Africa and Asia (CARIAA) with financial support from the UK Government's Department for International Development and the International Development Research Centre, Ottawa, Canada. The views expressed in this work are those of the creators and do not necessarily represent those of the UK Government's Department for International Development, the International Development Research Centre, Canada or its Board of Governors. The supercomputing resources for this study were provided by Compute/Calcul Canada (https://computecanada.ca) and WestGrid (https://www.westgrid.ca) on the Bugaboo cluster. We thank three anonymous reviewers for their feedback, which greatly improved the manuscript, and V. Radic for acting as editor. Edited by: V. Radic