Volume-area scaling for debris-covered glaciers

Abstract

A volume-area scaling relation is commonly used to estimate glacier volume or its future changes on a global scale. The presence of an insulating supraglacial debris cover alters the mass-balance profile of a glacier, potentially modifying the scaling relation. Here, the nature of scaling relations for extensively debris-covered glaciers is investigated. Theoretical arguments suggest that the volume-area scaling exponent for these glaciers is ~7% smaller than that for clean glaciers. This is consistent with the results from flowline-model simulations of idealised glaciers, and the available data from the Himalaya. The best-fit scale factor for debris-covered Himalayan glaciers is ~60% larger compared to that for the clean ones, implying a significantly larger stored ice volume in a debris-covered glacier compared to a clean one having the same area. These results may help improve scaling-based estimates of glacier volume and future glacier changes in regions where debris-covered glaciers are abundant.

Introduction

A volume-area scaling relation (Bahr and others, 1997, 2015) has been used extensively to estimate the total volume of glaciers on a regional to global scale (e.g. Grinsted, 2013), or to develop zero-dimensional models of glacier evolution (e.g. Radić and others, 2007). These scaling-based models are often used for numerically-efficient representation of glacier dynamics within hydrological (e.g. Zhang and others, 2015) or climate models (e.g. Kumar and others, 2019).

The volume-area scaling relation,

(1)$$V = c A^{\gamma}\comma\; $$

relates the volume of a glacier V to its area A via a dimensionless exponent γ and a scale factor c (Bahr and others, 1997, 2015). For a large set of glaciers, the glacier area typically spans several orders of magnitude. In comparison, the corresponding variation in glacier-specific c is much weaker. Therefore, a single best-fit c, along with a fixed γ, provides a reasonable statistical description of a set of glaciers (Bahr and others, 2015). In this letter, the scaling relation is interpreted in this statistical sense where all the glaciers in a given set are described with the same values of the parameters c and γ.

Equation (1) was first established empirically (e.g. Chen and Ohmura, 1990). Subsequently, Bahr and others (1997, 2015) derived it using dimensional analysis, and proved that,

(2)$$\gamma = 1 + {m + 1\over \lpar n + 2\rpar \lpar q + 1\rpar }.$$

Here, n is the Glen's law exponent which is usually taken to be 3 (e.g. Oerlemans, 2001). m and q are the exponents describing the scaling of ablation rate (b) and glacier width (w) with length (L). That is,

(3)$$b = c_m L^m\comma\; $$

(4)$$w = c_qL^q.$$

Here, c _m and c _q are the corresponding best-fit scale-factors. For clean mountain glaciers, typical values of m = 2 and q = 3/5, obtain γ = 1.375 which is consistent with observation (Bahr and others, 1997). The dependence of the scale factor c on c _q, c _m and glacier slope are known, but there is no existing analytical prescription to estimate c (Bahr and others, 2015).

On debris-covered glaciers, the insulating effect of a thickening debris layer reduces melting, and leads to a decline in ablation downglacier (e.g. Dobhal and others, 2013). This may lead to a different m for debris-covered glaciers; if so, then Eqn (2) would lead to a different characteristic γ for this type of glacier. Similarly, a possibly different characteristic q for debris-covered glaciers compared to the assumed empirical global value of 3/5 (Bahr and others, 2015) may also alter the value of γ. On the other hand, any systematic difference in the characteristic c _q or slope for these two types of glaciers may induce a corresponding change in the characteristic c value (Bahr and others, 2015). Despite the above possibilities, the nature of volume-area scaling for debris-covered glaciers has not been discussed in the literature as yet. Given the abundance of supraglacial debris in High Mountain Asia (e.g. Kraaijenbrink and others, 2017), as well as in several other glacierised regions in the world, an answer to the above problem is of interest. Establishing the appropriate scaling relation for debris-covered glaciers may help in correcting possible biases in the scaling-based estimates of global to regional-scale glacial ice volume, and the scaling model predictions of future ice-loss and global sea-level rise.

With the above motivation, I consider the following questions here.

1. (1) How does ablation rate scale with length on debris-covered glaciers?

2. (2) How does glacier width scale with length for debris-covered glaciers?

3. (3) What are the appropriate volume-area scaling exponent (γ), and scale factor (c) for debris-covered glaciers?

The first question is investigated theoretically. Debris-ice coupled flowline model simulations of idealised steady-state glaciers are used to validate the theoretical results obtained. Available data/estimates of volume and area of Himalayan glaciers are then utilised to obtain empirical width–length scaling relations for clean and extensively debris-covered glaciers (question 2). The answers to the first two questions would lead to the appropriate scaling relation for debris-covered glaciers (question 3) via Eqn (2). The consistency among the theoretical predictions, numerical results and empirical scaling relations for the Himalayan glaciers is checked for.

Methods

In this section, I first describe the theoretical models of idealised clean and debris-covered glaciers. These models are then utilised to establish the scaling of ablation rate with length for debris-covered glaciers. Subsequently, details about the numerical experiments performed with the idealised models, and the data analysis done for Himalayan clean and debris-covered glaciers are described.

Model for clean and debris-covered glaciers

The dynamics of clean glaciers, within a shallow-ice approximation, is described by the following equations (e.g. Oerlemans, 2001):

(5)$$\partial_t h = -{1\over w}\partial_x \lpar w u h\rpar + b\lpar z_{\rm s}\rpar \comma\; $$

(6)$$u = f_{\rm D} \vert \partial_x{z_{\rm s}}\vert ^n h^{n + 1}\comma\; $$

where h and u are local ice thickness and depth-averaged flow velocity, respectively, at a distance x along the central flowline. The mass balance b depends on surface elevation z _s. |∂_xz _s| is the absolute surface slope, and f _D is a constant controlling deformation of ice (Oerlemans, 2001). I did not consider any sliding for simplicity.

For debris-covered glaciers, the conservation of debris mass is imposed as,

(7)$$\partial_t d = -{1\over w}\partial_x \lpar w {\rm u_s} d\rpar + \alpha \vert b_{\rm d} \lpar z_{\rm s}\comma\; d\rpar \vert .$$

Here, d is the local supraglacial debris thickness, u _s = ((n + 2)/(n + 1)) u is the surface ice-flow speed, and b _d is the sub-debris ablation (hereinafter, subscript d is used to denote variables associated with debris-covered glaciers). The parameter α ≡ ν/(1 − ϕ) measures englacial debris content (Anderson and Anderson, 2016), with uniform englacial debris volume fraction ν, and supraglacial debris-layer porosity ϕ. α is assumed to be constant for a given glacier. Melt-out of englacial debris is the only source of supraglacial debris here, with sub-debris ablation taken to be (Anderson and Anderson, 2016),

(8)$$b_{\rm d}\lpar z_{\rm s}\comma\; d\rpar = {b\lpar z_{\rm s}\rpar \over 1 + d/d_0}.$$

Here, d ₀ is a characteristic length scale (~0.1 m) that is a material property (Anderson and Anderson, 2016).

Please note that this type of model has routinely been used to study steady-state and transient properties of debris-covered glaciers (Vacco and others, 2010; Anderson and Anderson, 2016; Banerjee and Wani, 2018). The limitations of the model are discussed later in the text.

Derivation of the ablation rate scaling exponent

Let us assume the following linear (i.e. m = 1) clean-glacier mass-balance profile:

(9)$$b\lpar z_{\rm s}\rpar = \beta \lpar z_{\rm s}-E\rpar \comma\; $$

where β is the balance gradient, and E is the equilibrium-line altitude (ELA). The following arguments can easily be generalised to other possible choices of b(z _s) (see Appendix A).

On an extensively debris-covered glacier, the ablation zone can be partitioned into an upper reach with thin debris cover (d(x) < d ₀), and a lower reach with a thick debris cover (d(x) > d ₀). Let the boundary between the two zones be at $x = x^\ast$ so that $d\lpar x^\ast \rpar = d_0$. If $L^\ast$ is the length of the thin-debris zone, Eqns (8) and (9) imply,

(10)$$b^\ast _{\rm d} = {1\over 2} \beta\lpar z_{\rm s}-E\rpar \approx-{1\over 2} \beta s L^\ast \comma\; $$

where $b_{\rm d}^\ast \equiv b_{\rm d}\lpar x^\ast \rpar$, and s is the bedrock slope. Here, I have used that $\lpar 1 + d\lpar x^\ast \rpar /d_0\rpar = 2$, and assumed the ice surface to be approximately parallel to the bedrock near ELA.

Now, in a steady state Eqn (7) implies,

(11)$$\partial_x \lpar w u_{\rm s} d\rpar = w \alpha \vert b_{\rm d} \lpar z_{\rm s}\comma\; d\rpar \vert \comma\; $$

At $x = x^\ast$, the left-hand side can be approximated by ${w^\ast u_{\rm s}^\ast d_0 }/{L^\ast }$, and the right-hand side is given by $\alpha w^\ast b_{\rm d}^\ast$, where $u_{\rm s}^\ast \equiv \bar {u}_s\lpar x^\ast \rpar$ and $w^\ast \equiv w\lpar x^\ast \rpar$. Thus,

(12)$${w^\ast u_{\rm s}^\ast d_0} \approx \alpha w^\ast L^\ast b_{\rm d}^\ast $$

The above relation simply states that to maintain a steady state the supraglacial debris flux at $x^\ast$ must equal the rate of total melt-out debris production over the stretch $L^\ast$ at the top of the ablation zone.

Combining Eqns (10) and (12), one obtains,

(13)$$b^\ast _{\rm d}\approx \sqrt{{ s \beta d_0 \over 2 \alpha} u^\ast _{\rm s}}.$$

Now, using $u^\ast _{\rm s}\sim h^{n + 1}$ (from Eqn (6)), and h = V/A ~ $L^{\lpar {m_{\rm d} + 1}\rpar /\lpar {n + 2}\rpar }$ (from Eqns. (2 and 4)), one obtains that $b^\ast _{\rm d} \sim L^{\lpar {\lpar n + 1\rpar \lpar m_{\rm d} + 1\rpar \rpar }/{2\lpar n + 2\rpar }}$. Or equivalently, m _d = ((n + 1)(m _d + 1))/2(n + 2). This equation is solved for m _d to obtain,

(14)$$m_{\rm d} = {n + 1\over n + 3} = {2\over 3}\comma\; $$

where the second equation is for n = 3. Apart from its dependence on $u_{\rm s}^\ast$, $b_{\rm d}^\ast$ also scales with slope as s ^1/2. This slope dependence gets absorbed into the corresponding scale factor c _m (Bahr and others, 2015). The above discussions also imply that $L^\ast \sim L^{ 2/ 3}$.

Please note that the above relations are valid only in the limit of a uniform englacial debris concentration, a linear clean-glacier mass-balance profile (i.e. m = 1), and a steady glacier.

Repeating the above arguments for a quadratic (m = 2) clean-glacier mass-balance profile (see Appendix A) obtains,

(15)$$m_{\rm d} = {2n + 3\over n + 3} = {3\over 2}\comma\; $$

where the second equality assumes n = 3.

Numerical methods

Numerical flowline-model simulations of idealised glaciers having a constant width (i.e. q = 0) and a linear bedrock, were used to verify the theoretical results derived above.

The clean glaciers were simulated with an explicit finite-difference implementation (Oerlemans, 2001) of Eqns (5 and 6), with the mass-balance profile as given by Eqn (9). The spatial and temporal discretisation steps were 25 m and 0.01 year, respectively. A linear bedrock with a maximum elevation of 6000 m and bedrock slope in the range 0.1–0.4 was used. β was taken to be 0.007 a⁻¹. These values were chosen to mimic typical Himalayan conditions (Laha and others, 2017). Each simulation was run with a fixed ELA over 4000 years to produce the corresponding steady-state. The length of the modelled steady-state clean glaciers for ELA values of 5400, 5500,…, 5900 and 5950 m, was in the range of 3–17 km.

The numerical model for debris-covered glaciers was the same as above, except that the sub-debris ablation was computed using Eqn (8) with d ₀ = 0.1 m, and that the supraglacial debris-thickness profile was evolved by solving Eqn (7) with an upwind finite-difference scheme. α was varied between 0.0002 and 0.0012 such that the maximum modelled debris thickness was in the range of 0.1–3 m. This is comparable to Himalayan conditions (Kraaijenbrink and others, 2017; Banerjee and Wani, 2018). The length of the modelled steady-state debris-covered glaciers was in the range of 3–28 km, which is again typical of Himalayan debris-covered glaciers (Banerjee and Wani, 2018).

Empirical analysis of data from the Himalaya

For all the clean and debris-covered glaciers in the Himalaya that are larger than 1 km², area (RGI, 2017) and estimated volume (Kraaijenbrink and others, 2017) were fitted to Eqn (1) to obtain the corresponding volume-area scaling relation. The volume estimates used here (Kraaijenbrink and others, 2017) were obtained from a variant of plastic-ice approximation (Frey and others, 2014). Available debris-cover extents (Kraaijenbrink and others, 2017) were utilised to select the clean (<5% debris-cover), and extensively debris-covered (> 50% debris cover) Himalayan glaciers. Similarly, the data of length and area (RGI, 2017) of these two types of glaciers were used to obtain the corresponding width-scaling relation (Eqn (4)).

Results and discussions

In this section, I first discuss some general features of the simulated steady-state glaciers. Then, the simulation results are used to validate the theoretical predictions for scaling of ablation rate with length, and that of volume with area for the idealised glaciers with uniform width. Subsequently, the empirical width and volume scaling relations for Himalayan glaciers are presented. In the end, I discuss the implications of the differences between the scaling relations for debris-covered and clean glaciers.

The steady-state simulated glaciers

As shown in Figure 1, the simulations reproduced various well-known characteristics of debris-covered glaciers, e.g. a down-glacier increase in debris-thickness, a non-monotonic mass-balance profile with low melt rate near the terminus, a relatively longer ablation zone, and a smaller accumulation-area ratio compared to the corresponding clean glaciers (Banerjee and Shankar, 2013). Glacier length, ice thickness, typical ablation rate and debris-thickness of the modelled glaciers were in the reasonable range as well (Kraaijenbrink and others, 2017; Laha and others, 2017; Banerjee and Wani, 2018). I confirm that the debris-flux increased monotonically down-glacier, and that divergence of the debris flux matched the melt-out rates of debris at all sites.

Fig. 1. (a) Ice-thickness profiles of the simulated debris-covered/clean glaciers are plotted with thick/thin solid lines. Line colours represent the ELA. Bedrock slope is 0.1, and englacial debris volume fraction α =0.0008. (b) The corresponding mass-balance and debris-thickness (inset) profiles of simulated debris-covered glaciers are shown.

It can be seen in Figure 1 that the maxima of the debris-thickness profiles were located somewhat upstream of the terminus. This non-monotonicity of the modelled debris-thickness profile may be an artefact of the simple boundary condition used at the terminus, where the debris flux out of the last site is assumed to leave the system. For real debris-covered glaciers or for model glaciers with different boundary conditions (e.g. Anderson and Anderson, 2016) such a feature may or may not be present. The present boundary condition also led to some grid-size dependence of the steady state.

Scaling of ablation rate for the simulated glaciers

Figure 2 shows that for the modelled constant-width idealised glaciers, the characteristic sub-debris ablation rate $b^\ast _{\rm d}$ scales as ~ L ^0.65±0.00, i.e. m _d ≈ 2/3. This is just as predicted by Eqn (14) for a linear clean glacier mass-balance profile (Eqn (9)). The same figure also shows that the mean ablation rate on clean glaciers scales with L almost linearly with m = 0.90 ± 0.00. The small departure from linearity may be related to the corrections due to the ice-elevation feedback effects.

Fig. 2. The characteristic ablation rates $b^\ast _{\rm d}$ on the simulated debris-covered glaciers (red circles) scales as ~ L ^0.67, and mean ablation on the simulated clean glaciers (blue circles) scales as ~ L. This is consistent with the theoretical analysis (see text). The best-fit curves (solid lines) yield exponents that are somewhat smaller with m = 0.90 ± 0.0 and m _d = 0.65 ± 0.00.

As shown in Appendix Figure A1, the simulation results corresponding to a quadratic clean-glacier mass-balance profile produced relatively smaller values of m = 1.5 and m _d = 1.2 as compared to the expected/predicted values of m = 2 and m _d = 3/2.

Volume-area scaling for the simulated glaciers

For the simulated constant-width clean glaciers with linear mass-balance profile (q = 0 and m = 1), and corresponding debris-covered glaciers (q = 0 and m _d = 2/3), Eqn (2) predicts γ = 1.4 and γ_d = 1.33, respectively. Simulation results (Fig. 3) confirmed that these theoretically predicted exponents provided accurate descriptions of the variation of V with A for the constant-width clean and debris-covered glaciers. A similar match between the simulated and predicted γ and γ_d values was obtained for the quadratic clean-glacier mass-balance profile as well (Fig. A1).

Fig. 3. Volume of simulated steady glaciers with bed slope of 0.1 is plotted with solid circles as a function of area. In (a) and (b), symbol colours denote englacial debris volume fraction α and debris-covered fraction, respectively. All the glaciers here have a constant width of 1 km. Power-law behaviour with the predicted exponents for clean/debris-covered glaciers are shown with solid/dashed lines as guide to the eye. The inset in (a) shows volume-area scaling for glaciers with more than 50% debris cover for bed slope of 0.1 (black symbols), 0.2 (purple symbols) and 0.3 (cyan symbols).

For debris-covered glaciers, the scale factor increased weakly with increasing α (Fig. 3), and decreased with increasing bedrock slope (Fig. 3, inset). For the smallest englacial debris concentration used, namely, α = 0.0002, there were systematic deviations from the predicted power-law scaling, particularly for smaller glaciers with a relatively low debris-covered fraction (Fig. 3). This suggested a smooth transition from the clean glacier value of γ = 1.4 to γ_d = 1.33 as α and/or debris-covered fraction increased.

Empirical width-scaling exponents for Himalayan glaciers

An empirical analysis of the data for clean and extensively debris-covered Himalayan glaciers obtained w = (0.43 ± 0.02) L ^0.60±0.02 and w = (0.19 ± 0.04) L ^0.80±0.09, respectively. While the above best-fit q on clean Himalayan glaciers matched the recommended global value (Bahr, 1997), the best-fit q _d on debris-covered Himalayan glaciers was relatively larger. Also, the best-fit scale factor (c _q) was significantly smaller for the debris-covered ones compared to that for the clean ones.

Given the above empirical values of q = 0.6 and q _d = 0.8 for the Himalayan glaciers, and assuming m = 2 and m _d = 1.5, Eqn (2) obtained γ = 1.375 and γ_d = 1.28 for the clean and debris-covered Himalayan glaciers, respectively.

Empirical volume-area scaling for Himalayan glaciers

For the above clean and debris-covered Himalayan glaciers, the best-fit volume-area scaling relations were, V = (0.024 ± 0.001) A ^1.38±0.02 and V = (0.040 ± 0.002) A ^1.30±0.02, respectively. These best-fit forms are consistent with the theoretical predictions for the volume-area scaling exponents presented above. Thus, theoretical calculations, simulation results for idealised glaciers and empirical data from the Himalaya were all consistent with each other as far as the values of the volume-area scaling exponent is concerned. These pieces of evidence indicate that the volume-area scaling relations for extensively debris-covered and clean glaciers are different, with γ larger than γ_d by about 7%. Moreover, the best-fit volume-area scaling relations for Himalayan glaciers suggest that the scale factor c for debris-covered glaciers is about 1.6 times that for the clean glaciers. This can be ascribed to a significantly smaller best-fit c _q, and a gentler slope of the debris-covered glaciers (Bahr and others, 2015). The median values of slope for the debris-covered and clean Himalayan glaciers are 16 and 21^°, respectively.

Notably, the alternative expression for the exponent given by Bahr and others (2015) γ = (2q + 1)/(q + 1), is inconsistent with the above values of q _d = 0.8 and γ_d = 1.3 obtained from empirical data of debris-covered Himalayan glaciers. This suggests that some of the assumptions involved in deriving this alternative expression may not be valid for debris-covered glaciers.

Implications for scaling-based models

A relatively smaller γ_d (by ~6%) and a larger c (by ~60%) for extensively debris-covered Himalayan glaciers compared to their clean counterparts as established above, suggest a possible underestimation of the volume of debris-covered glaciers when derived using the clean-glacier scaling relation. For a given value of glacier area, this fractional bias is given by, ΔV/V ≈ Δc/c + Δγ lnA, and the negative first term on the right-hand side is likely to have a larger magnitude than the positive second term (Δc/c ~ −0.6, and Δγ ~ 0.1).

The above underestimation of volume implies a corresponding negative bias in the estimate mean ice-thickness. A low bias in estimated thickness may lead to an overestimation of instantaneous area loss on debris-covered glaciers when computed using the relation ΔA = ΔV/γh, if the fractional negative bias in h is larger than the fractional positive bias in γ. Any such bias in the estimated area change for debris-covered glaciers may affect the asymptotic response due to a feedback between glacier area and net balance. These issues may compromise the accuracy of the predictions from scaling-based glacier models in regions where extensively debris-covered glaciers are abundant.

Limitations

The theoretical results on scaling of ablation with glacier length, and that of volume with area on steady debris-covered glaciers presented above, are based on assumptions of a uniform englacial debris-concentration, and specific forms of debris-thickness-dependent ablation rate. There are several processes in a typical Himalayan debris-covered glacier that violate these assumptions to some extent. The englacial debris concentration is not constant, and varies spatially (Nakawo and others, 1986) even within a single glacier. It also varies with time due to stochastic and systematic changes in headwall-erosion rates (Banerjee and Wani, 2018). Additional sources of debris due to rockfall, debris avalanche and degradation of lateral moraines (Nakawo and others, 1986) add further noise to the debris-thickness distribution (Shah and others, 2019). Complex mass-balance processes like avalanche activity (Laha and others, 2017), and stochastic increase in local ablation due to thermokarst features like supraglacial ponds and ice-cliffs (Sakai and others, 2000), may modify the mass-balance profile. Despite ignoring these complicating factors, the close match between theoretically predicted and empirically observed scaling exponents for Himalayan glaciers may indicate that the idealised model employed here was able to capture the dominant processes on debris-covered glaciers rather well.

Conclusions

To summarise the results from the above analysis of scaling relations for debris-covered glaciers:

1. (1) The characteristic ablation rate on idealised debris-covered glaciers scale with glaciers length with the exponent m _d ≈ 1.5. m _d is smaller than the corresponding clean-glacier exponent m = 2.

2. (2) For the extensively debris-covered glaciers in the Himalaya, the glacier width scales with length with an exponent q _d ≈ 0.8. The corresponding exponent for the Himalayan clean glaciers is q = 0.6.

3. (3) The above values of m _d and q _d predict a volume-area scaling exponent γ_d = 1.28 for debris-covered glaciers. γ_d is about 7% smaller than the corresponding clean-glacier exponent γ = 1.375.

4. (4) The empirical volume-area scaling exponents for clean and debris-covered Himalayan glaciers (γ = 1.38 ± 0.00 and γ_d = 1.30 ± 0.02) are consistent with the above predictions.

5. (5) The Himalayan debris-covered glaciers are characterised by a ~60% larger value of the scale factor c.

The difference in the scaling relations for the two types of glaciers – particularly, a larger scale factor c for the debris-covered glaciers – when not taken into account, may cause a significant negative bias in corresponding scaling-derived estimates of volume and ice thickness. Related biases could also be present in the results of scaling-based simple models of glacier dynamics, with an overestimation of the instantaneous area change of debris-covered glaciers for a given mass-balance forcing. A larger scale factor for the debris-covered Himalayan glaciers implies that the stored ice volume in a debris-covered glacier is larger than that in a clean glacier having similar area. The results presented here may be useful in quantifying and correcting possible biases in existing scaling-based estimates of the volume of debris-covered glaciers or that of the future ice loss at these glaciers.