Distributed Global Debris Thickness Estimates Reveal Debris Significantly Impacts Glacier Mass Balance.

Supraglacial debris affects glacier mass balance as a thin layer enhances surface melting, while a thick layer reduces it. While many glaciers are debris-covered, global glacier models do not account for debris because its thickness is unknown. We provide the first globally distributed debris thickness estimates using a novel approach combining sub-debris melt and surface temperature inversion methods. Results are evaluated against observations from 22 glaciers. We find the median global debris thickness is ∼0.15 ± 0.06 m. In all regions, the net effect of accounting for debris is a reduction in sub-debris melt, on average, by 37%, which can impact regional mass balance by up to 0.40 m water equivalent (w.e.) yr-1. We also find recent observations of similar thinning rates over debris-covered and clean ice glacier tongues is primarily due to differences in ice dynamics. Our results demonstrate the importance of accounting for debris in glacier modeling efforts.

Introduction

The thickness of supraglacial debris will control how glaciers respond to climate change, since a thin layer (less than a few centimeters) will enhance melting, while a thick layer will insulate the underlying ice and suppress melting (Østrem, 1959). While this effect is well known at the point scale, recent studies in High Mountain Asia found similar elevation change rates for debris‐covered and clean‐ice glaciers (Gardelle et al., 2013; Kääb et al., 2012), a phenomenon referred to as the “debris‐cover anomaly” (Pellicciotti et al., 2015). Proposed explanations include debris‐covered glaciers experiencing enhanced melting due to ice cliffs and ponds and/or having comparatively lower emergence velocities (Brun et al., 2018). A regional study that accounted for debris on High Mountain Asia glaciers found debris delays projected mass loss (Kraaijenbrink et al., 2017). However, the effects of debris on regional and global projections of glacier mass change are largely unknown due to a lack of accurate debris thickness estimates. Given debris covers ∼4%–7% of the global glacier area (excluding the ice sheets), and up to 23% regionally (Herreid & Pellicciotti, 2020; Scherler et al., 2018), accurate debris thickness estimates are urgently needed to quantify the impact of debris on glacier mass balance from local to global scales.

Despite the importance of debris thickness for glacier evolution, in‐situ measurements of debris thickness and sub‐debris melt rates are relatively scarce and primarily limited to glaciers in Central Europe and High Mountain Asia (Tables S1 and S2). Even in‐situ measurements from well‐studied debris‐covered glaciers are typically spatially and/or temporally limited. Nonetheless, these measurements have been used to develop empirical relationships to estimate debris thickness (Kraaijenbrink et al., 2017; Mihalcea et al., 2008), or to validate debris thickness estimates from physically‐based surface temperature inversion (Foster et al., 2012; Rounce & McKinney, 2014) or sub‐debris melt inversion (Rounce et al., 2018) methods. At present, the sub‐debris melt inversion method is the only physically‐based method that can accurately estimate debris thicknesses greater than 0.5 m. However, the method is currently limited to elevation bins and requires accurate estimates of the climatic mass balance (i.e., the sum of surface accumulation, ablation, and refreezing), which must be derived from limited in‐situ measurements or remotely‐sensed elevation change data. The latter must be corrected for elevation change due to ice flow, which is challenging over steep, fast‐flowing areas. Conversely, surface temperature inversion methods have no spatial limitations but struggle to estimate thick debris and must accurately resolve the surface energy balance fluxes at the time when the satellite image is acquired.

Here we produce the first global debris thickness estimates for each of the 92,033 debris‐covered glaciers in the Randolph Glacier Inventory (RGI, version 6) (RGI Consortium, 2017; Scherler et al., 2018), excluding the ice sheets and Antarctic periphery. These estimates are then used to determine the cause of the debris cover anomaly and quantify how debris impacts sub‐debris melt and glacier mass balance on local, regional, and global scales.

Methods

Our novel approach first uses a sub‐debris melt inversion method to estimate the debris thickness over binned near‐stagnant debris‐covered areas for all debris‐covered glaciers greater than 2 km2 with available geodetic mass balance and surface velocity data (Figure S1). A surface temperature inversion method is then used to estimate the distributed debris thickness over the glacier's entire debris‐covered area using the binned debris thickness estimates from the sub‐debris melt inversion method for calibration. For the remaining debris‐covered glaciers (47% of the total debris‐covered area), we extrapolate using the nearest calibrated glaciers' coefficients to estimate their distributed debris thickness.

Sub‐debris Melt Inversion Method

The sub‐debris melt inversion method relates annual sub‐debris melt (M; m w.e.) to debris thickness (h ; m) according to a second‐order reaction rate equation:

where M 0 and k are coefficients derived from fitting the modeled sub‐debris melt from a debris‐covered glacier energy balance model (Rounce et al., 2015; Text S1) for various debris thicknesses (Figure S2a).

Since debris is located in the ablation area, we can assume the sub‐debris melt is equivalent to the annual climatic mass balance, which is estimated as the elevation change plus the flux divergence (all in m w.e.). We refer to this as the observed sub‐debris melt. Elevation change data between 2000 and 2018 (Braun et al., 2019; Malz et al., 2019; Seehaus et al., 2019; Shean et al., 2020; Sommer et al., 2020) was supplemented with data from Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) DEMs, Arctic DEM v3.0 strips, and SPOT5 IPY‐SPIRIT DEMs (Text S2; Table S3). Following Rounce et al. (2018), the flux divergence is estimated from surface velocity data from 1985 to 2018 (Gardner et al., 2019) and ice thickness estimates (Farinotti et al., 2019), assuming only internal deformation and thus a ratio of column‐average ice velocity to surface velocity of 0.8. Missing velocity data are supplemented by feature tracking of Landsat and Sentinel‐2 images between 2015 and 2018 (Millan et al., 2019) (Table S3).

The inversion is only performed over near‐stagnant (<7.5 m yr−1) debris‐covered areas, and the observed sub‐debris melt data are aggregated into 10 m elevation bins, to minimize uncertainty associated with the flux divergence (since the flux divergence is near zero). Any bins where the observed sub‐debris melt exceeds our maximum modeled sub‐debris melt, or is below our minimum modeled sub‐debris melt, are discarded. The latter helps avoid estimating the debris thickness over relatively stable non‐glacierized terrain due to glacier retreat, since the debris cover and glacier extents were delineated over different time periods (RGI Consortium, 2017; Scherler et al., 2018). Uncertainty arising from the elevation change, surface velocity, and ice thickness data is incorporated into the debris thickness estimates (Text S3).

Surface Temperature Inversion Method

The distributed debris thickness is then estimated for each debris‐covered glacier grid cell using the surface temperature inversion method, based on modeled surface temperatures and Landsat‐8 surface temperature data, that is calibrated with the debris thickness estimates from the sub‐debris melt inversion method. The surface temperature data is a composite image generated in Google Earth Engine from all available Landsat‐8 images between 2013 and 2018 from the ablation season (Kraaijenbrink et al., 2017). The relationship between the debris thickness and modeled debris surface temperature (T ; °C) from a debris‐covered glacier energy balance model (Rounce et al., 2015; Text S1) (Figure S2c) is estimated using the Hill equation:

where a, b, and c are calibrated coefficients that capture the nonlinear relationship between the debris thickness and surface temperature at the time when the data were acquired. To avoid overestimating the debris thickness by only calibrating the coefficients with near‐stagnant bins at the terminus that likely comprise thicker debris, we assume the coldest debris‐covered elevation bin has a debris thickness of 0.02 m. Furthermore, the calibration minimizes the difference between the estimated debris thicknesses' sub‐debris melt rates from the two inversion methods (Equation 1). This gives more weight to thinner debris, which is appropriate since the sub‐debris melt inversion method is more sensitive to thinner debris.

To capture the nonlinear relationship between the debris thickness and surface temperature, we identify all snow‐free days, according to the energy balance model, within one month of the Landsat‐8 acquisition date. Since the debris is warming rapidly in the morning when the satellite image is acquired (∼10:30 a.m. local), we linearly interpolate the hourly modeled surface temperatures to estimate the surface temperatures that coincide with the acquisition time. We then fit coefficients a, b, and c using the interpolated modeled surface temperatures; thus, producing up to 60 potential sets of coefficients. We estimate the debris thickness from the Landsat‐8 surface temperature data for each set and select the coefficients that minimize the difference between the estimated sub‐debris melt rates and those from the sub‐debris melt inversion method (Equation 1). To ensure the difference between the sub‐debris melt and surface temperature inversion methods is minimized, we then fix coefficients b and c, and recalibrate a. The calibrated coefficients are then used to estimate the distributed debris thickness over the entire debris‐covered glacier area (Figure S1).

Extrapolated Debris Thickness Estimates

For glaciers lacking calibration data (47% of the total debris‐covered area), that is, those that are missing data, too small, or not stagnant at their terminus, we use coefficients (M 0, k, a, b, and c) from the nearest calibrated glaciers to estimate the debris thickness. To ensure the extrapolated debris thicknesses are physically reasonable, each uncalibrated glacier loops through the parameter sets of its nearest neighbors and accepts the debris thickness estimates if they meet two conditions: (i) the debris thickness estimate for the coldest elevation bin, that is, the thinnest debris, is less than 0.2 m, and (ii) the debris thickness estimate for the warmest elevation bin, that is, the thickest debris, is less than 3.0 m. The uncalibrated glaciers loop through the parameter sets until 10 reasonable estimates are found. The debris thickness is calculated to be the median of these 10 reasonable estimates for each grid cell. A leave‐one‐glacier‐out cross‐validation, that is, applying the extrapolation to calibrated glaciers, reveals the median (±normalized median absolute deviation) difference is 0.04 ± 0.20 m. The median distance to the nearest calibrated glacier is 12 km (Figure S3).

Sub‐debris Melt Enhancement Factors

Sub‐debris melt enhancement factors are estimated from the debris thickness using the debris thickness‐melt curves (Equation 1) normalized by the clean‐ice melt (Figure S2b). For debris thicknesses less than 0.02 m, the enhancement factor is linearly interpolated between the value of one for clean ice and the enhancement factor associated with a debris thickness of 0.02 m. For each climate grid cell containing debris‐covered glaciers, clean‐ice melt is estimated using an energy balance model (Text S1).

Modeled Regional Mass Balance

The Python Glacier Evolution Model (PyGEM) (Rounce et al., 2020) is used to quantify the impact of debris thickness on regional mass balance for Central Europe, Caucasus/Middle East, and High Mountain Asia (Central Asia, South Asia West, South Asia East). PyGEM estimates the climatic mass balance of every glacier using 10 m elevation bins and a monthly time step. The model is forced with monthly ERA5 air temperature and precipitation data, and computes glacier melt using a degree‐day model, accumulation using a temperature threshold, and refreezing based on the annual air temperature. The sub‐debris melt enhancement factors for each elevation bin are used to account for debris. The glacier area is assumed constant.

We calibrate each glacier using the calibration scheme from Huss and Hock (2015) and regional specific mass balance data from 2000 to 2018 for High Mountain Asia (Shean et al., 2020) and 2006–2015 for the other regions (Hock et al., 2019). We then run the calibrated model from 2000 to 2018 without applying the enhancement factors for debris‐covered elevation bins and compare the results to quantify the impact of accounting for debris. We perform an additional simulation in High Mountain Asia, where we calibrate the model using data from every glacier to determine how the use of regional mass balance data for calibration impacts the results. The regional specific mass balances are within ±0.01 m w.e. yr−1 compared to the simulations run using regional specific mass balance for calibration, thus justifying the use of regional data for the other regions. Uncertainties associated with the debris thickness estimates are evaluated from additional simulations using enhancement factors derived from the 16% and 84% debris thickness estimates (Text S3).

Results

Debris Thickness Estimates

The distributed debris thickness estimates of some of the most well‐studied debris‐covered glaciers in Alaska, Europe, and High Mountain Asia show that debris thickness is less than 0.02 m in the uppermost debris‐covered areas and increases toward the termini of the glaciers (Figure 1) as expected (Anderson & Anderson, 2018). The maximum debris thickness of these glaciers ranges from ∼0.6 m on Kennicott and Lirung to more than 2 m on Koxkar and Ngozumpa. The debris thickness of most glaciers also increases toward the lateral margins, where unstable valley walls can provide an additional source of debris. Glaciers with prominent medial moraines may have thicker debris in the center of the glacier (e.g., Suldenferner; Figure 1b).

Figure 1

Modeled debris thickness (h ) for well‐studied debris‐covered glaciers located in Alaska (a), Central Europe (b, c), Central Asia (d), South Asia West (f), and South Asia East (e, g, h). Note the saturated, nonlinear colorbar. Accumulation areas, where debris cover is absent, are not shown. Arrows indicate glacier flow direction toward the terminus.

To evaluate our debris thickness estimates, we compare individual measurements with geolocation information to the nearest pixel and also aggregate measurements/pixels into 50 m elevation bins comparing those with more than five observations (Table S1; Figures 2a and 2b). Both comparisons face issues with measurements being limited in space and/or quantity, varying considerably over small scales (<30 m), having limited geolocation information, and changing over time due to glacier flow and/or local debris redistribution. Our estimates agree reasonably well with these in‐situ measurements (Figures 2a and 2b; Table S4) considering the method's global application. The scatter is expected considering the small‐scale variability with point measurements and that its unlikely the limited in‐situ measurements are representative of the entire elevation bin which includes thicker debris near the lateral moraines (Text S4; Figure S4). Interestingly, our debris thickness uncertainty estimates (Text S3; Figure S5) are comparable to the root mean square error of in‐situ versus modeled debris thicknesses, regardless of debris thickness magnitude (Table S4), which lends confidence to our uncertainty estimates.

Figure 2

Comparison of modeled and observed (a, b) debris thickness (h ) for elevation bins and individual points/pixels and (c) sub‐debris melt rates using the observed debris thickness. Binned debris thicknesses are scaled and colored by the number of observations (n ). Sub‐debris melt rates are colored by debris thickness. Error bars represent spatial variability within the corresponding 50 m elevation bin according to the median and normalized median absolute deviation. Details of the observations are provided in Tables S1 and S2.

The method's ability to reproduce debris thicknesses exceeding 1 m (e.g., Koxkar (Wu & Liu, 2012), Ngozumpa (Nicholson & Benn, 2012); Figure 1) suggests the method can reasonably estimate the total volume of debris. Similarly, despite a nonlinear debris thickness‐melt relationship, the method's ability to reproduce debris thicknesses less than 0.2 m (e.g., Kennicott in Alaska (Anderson et al., 2019), Suldenferner in Central Europe (Nicholson & Boxall, 2020), Khumbu in High Mountain Asia (Rowan & Gibson, 2020); Figure 1) suggests sub‐debris melt rates are well estimated.

The global debris thickness estimates show considerable regional‐scale variability (Figure 3). The global median (16%–84% confidence interval) debris thickness of all debris‐covered areas is 0.15 m (0.09–0.21 m) and the total volume of debris is 7.3 km3 (4.3–12.8 km3). The global mean debris thickness of 0.35 m (0.21–0.62 m) is higher than the median since the debris thickness histograms are right‐skewed (Figure 3). Regionally, the median debris thickness varies from 0.04 m (0.02–0.07 m) in North Asia to 0.19 m (0.12–0.27 m) in South Asia East. Other regions with thicker debris include Alaska and Caucasus/Middle East. Regions with thinner debris include Svalbard and Southern Andes.

Figure 3

Regional debris thickness (h ) and sub‐debris melt enhancement factors (E ). Bars show median values with the tick showing values of 0.1 m and 0.5, respectively. Circles show relative debris‐covered area (Scherler et al., 2018). Histograms are normalized by their total area. Enhancement factors greater than one are darker red. Numbers refer to the total debris‐covered area (km2). Debris‐covered glaciers are shown in blue and RGI regions are outlined in black.

Debris Thickness Impact on Sub‐debris Melt

The impact of debris thickness on sub‐debris melt is quantified using sub‐debris melt enhancement factors (hereafter referred to as enhancement factors), which equal the ratio between sub‐debris melt and clean‐ice melt. Enhanced sub‐debris melt beneath thin debris (typically 0.03–0.05 m) coincides with enhancement factors greater than 1, while thick debris can have enhancement factors less than 0.1, corresponding to a greater than 90% reduction in sub‐debris melt. The agreement and lack of systematic bias between modeled sub‐debris melt rates and in‐situ measurements (Figure 2; r = 0.85, RMSE = 0.08 mm w.e. d−1) provides confidence in our enhancement factors. Note the model underestimates high melt rates (>50 mm w.e. d−1), which may be due to local topography or differences in debris properties.

For individual glaciers, spatial variations in modeled debris thickness results in enhanced sub‐debris melting in the uppermost debris‐covered areas and limited melting near the glacier termini (Figure S6). Some debris‐covered glaciers (e.g., Koxkar, Ngozumpa) have enhancement factors less than 0.5 over up to 84% of the debris‐covered area. Others (e.g., Kennicott, Miage, Baltoro) have values greater than one over up to 45% of the debris‐covered area.

Regionally, the nonlinear debris thickness‐melt relationship means that the spatial distribution of debris thickness, and corresponding enhancement factors, are important (Figure 3). The global mean enhancement factor is 0.63 (0.51–0.73), varying from a minimum of 0.39 (0.28–0.5) in Iceland to a maximum of 0.80 (0.68–0.88) in North Asia. In all regions, the net effect of accounting for debris is a reduction in sub‐debris glacier melt.

Debris Thickness Impact on Regional Mass Balance

The impact of debris on regional glacier mass balance is exemplified for some regions with significant debris cover: Central Europe, Caucasus/Middle East, Central Asia, South Asia West, and South Asia East. The net impact depends on the spatial distribution of debris thickness, relative amount of debris‐covered area to total glacier area, and mass balance gradient. The latter is essential since the glacier termini, where debris is thickest, coincides with where the climatic mass balance is most negative for clean‐ice glaciers. We estimate how much accounting for debris thickness affects regional specific mass balance (i.e., mass change per unit area) by running simulations that (1) account for debris thickness using enhancement factors and (2) assume no debris cover.

Accounting for debris changes the regional specific mass balance by +0.40 m w.e. yr−1 (0.25–0.64 m w.e. yr−1) for Caucasus/Middle East and +0.18 m w.e. yr−1 (0.06–0.38 m w.e. yr−1) for Central Europe (Figure 4). The change in regions in High Mountain Asia varies from +0.26 m w.e. yr−1 (0.19–0.36 m w.e. yr−1) in South Asia East to +0.07 m w.e. yr−1 (0.03–0.12 m w.e. yr−1) in South Asia West; interestingly, the changes as a percentage of regional mass balance are similar, varying from +30 to 39%. These results show that debris substantially reduces regional mass loss despite covering a small fraction of the total glacier area (5%–23%).

Figure 4

Climatic mass balance (b ), area‐altitude distributions (AAD), and sub‐debris melt enhancement factors (E ) as a function of elevation for various regions. Climatic mass balances from PyGEM are shown accounting for debris (green) and assuming clean ice (blue) with colored text showing the regional glacier‐wide mass balance (m w.e. yr−1). Area‐altitude distributions are shown for the entire glacier (blue) and debris‐covered area (AAD ; green). Shaded areas represent uncertainty (16% and 84% confidence interval) from debris thickness estimates. Dotted horizontal line shows the modeled mean annual equilibrium line altitude.

Debris Cover Anomaly

Our debris thicknesses are derived from elevation change and surface temperature data that include surface features like ice cliffs and ponds, which can melt 3–17 times faster than the surrounding debris‐covered ice (Brun et al., 2018; Miles et al., 2018). Consequently, the mixed‐pixel effects reduce modeled debris thickness estimates (Rounce et al., 2018) and increase enhancement factors (Text S4). Despite their higher melt rates, the enhancement factors for cliffs and ponds on debris‐covered glaciers in previous studies vary from 0.81 to 1.85 (Table S5), which are comparable to our enhancement factors for thin debris (Figure S6).

In High Mountain Asia, the mean enhancement factor of 0.64 indicates melt beneath debris‐covered areas is, on average, reduced by 36% compared to clean ice. If we assume ice cliffs and ponds have enhancement factors of 0.81–1.85 and comprise 6%–14% of the debris‐covered areas (Brun et al., 2018; Herreid & Pellicciotti, 2018; Miles et al., 2018), then enhancement factors for the remaining debris‐covered area are 0.44–0.63. We can therefore estimate that enhanced melting due to cliffs and ponds accounts for 3%–35% (0.01–0.20 of the 0.37–0.56 reduction) of the debris cover anomaly and attribute the remaining 65%–97% to reduced emergence velocities over debris‐covered areas compared to clean ice. Hence, less ice flux into the slow‐moving, debris‐covered portion of the glacier is the dominant cause of the debris cover anomaly in High Mountain Asia.

Conclusions

Our global debris thickness estimates provide the first assessment of how debris affects the sub‐debris melt of debris‐covered glaciers at local, regional, and global scales. At the glacier scale, the spatial variations in debris thickness can cause glacier tongues to have an inverted ablation gradient meaning the glacier may experience higher rates of mass loss further upglacier than at its terminus, which is consistent with previous studies (e.g., Benn & Lehmkuhl, 2000). This can cause glaciers to stagnate, develop supraglacial and proglacial lakes, and become a flood hazard for downstream communities (Benn et al., 2012; Rounce et al., 2017; Shugar et al., 2020). These proglacial lakes also substantially increase glacier‐wide mass loss rates (King et al., 2019). Our enhancement factors will therefore enable glacier evolution models to account for these inverted ablation gradients and help predict the development of glacial lakes and their impact on glacier‐wide mass change.

At the regional scale, enhancement factors reveal the net effect of accounting for debris is a reduction in sub‐debris melt in all regions, which can significantly impact regional mass balance. Furthermore, we find the primary cause of the debris cover anomaly is a reduction in ice flux into the lowermost portions of debris‐covered glaciers. Regional glacier models should therefore seek to more accurately account for glacier dynamics. These enhancement factors will also improve ice thickness estimates for debris‐covered glaciers, which currently do not account for debris (Farinotti et al., 2017). Future work should prioritize systematically measuring debris thickness to enable proper validation and error assessments and measuring sub‐debris melt beneath thin debris to determine local and regional variability of the critical debris thickness (e.g., Fyffe et al., 2020). Furthermore, advancing models' ability to estimate melt beneath thin debris, improving the meteorological forcing, incorporating site‐specific measurements of debris properties and their spatial variation, and developing ways to for account for local topography may significantly improve and reduce uncertainty associated with the debris thicknesses and enhancement factors.

Overall, these global debris thickness estimates fundamentally advance our ability to account for debris thickness in glacier reconstructions, debris transport models, hazard assessments, and landscape evolution models. Future work should seek to account for debris in global glacier evolution models to determine how debris cover will affect projections of glacier mass change and runoff at regional and global scales.

Supporting Information S1

Click here for additional data file.