Introduction
A volume-area scaling relation (Bahr and others, Reference Bahr, Meier and Peckham1997, Reference Bahr, Pfeffer and Kaser2015) has been used extensively to estimate the total volume of glaciers on a regional to global scale (e.g. Grinsted, Reference Grinsted2013), or to develop zero-dimensional models of glacier evolution (e.g. Radić and others, Reference Radić, Hock and Oerlemans2007). These scaling-based models are often used for numerically-efficient representation of glacier dynamics within hydrological (e.g. Zhang and others, Reference Zhang, Hirabayashi, Liu and Liu2015) or climate models (e.g. Kumar and others, Reference Kumar, Saharwardi, Banerjee, Azam and Dubey2019).
The volume-area scaling relation,
relates the volume of a glacier V to its area A via a dimensionless exponent γ and a scale factor c (Bahr and others, Reference Bahr, Meier and Peckham1997, Reference Bahr, Pfeffer and Kaser2015). 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, Reference Bahr, Pfeffer and Kaser2015). 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, Reference Chen and Ohmura1990). Subsequently, Bahr and others (Reference Bahr, Meier and Peckham1997, Reference Bahr, Pfeffer and Kaser2015) derived it using dimensional analysis, and proved that,
Here, n is the Glen's law exponent which is usually taken to be 3 (e.g. Oerlemans, Reference Oerlemans2001). m and q are the exponents describing the scaling of ablation rate (b) and glacier width (w) with length (L). That is,
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, Reference Bahr, Meier and Peckham1997). 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, Reference Bahr, Pfeffer and Kaser2015).
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, Reference Dobhal, Mehta and Srivastava2013). 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, Reference Bahr, Pfeffer and Kaser2015) 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, Reference Bahr, Pfeffer and Kaser2015). 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, Reference Kraaijenbrink and Bierkens2017), 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) How does ablation rate scale with length on debris-covered glaciers?
(2) How does glacier width scale with length for debris-covered glaciers?
(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, Reference Oerlemans2001):
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, Reference Oerlemans2001). I did not consider any sliding for simplicity.
For debris-covered glaciers, the conservation of debris mass is imposed as,
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, Reference Anderson and Anderson2016), 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, Reference Anderson and Anderson2016),
Here, d 0 is a characteristic length scale (~0.1 m) that is a material property (Anderson and Anderson, Reference Anderson and Anderson2016).
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, Reference Vacco, Alley and Pollard2010; Anderson and Anderson, Reference Anderson and Anderson2016; Banerjee and Wani, Reference Banerjee and Wani2018). 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:
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 0), and a lower reach with a thick debris cover (d(x) > d 0). 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,
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,
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,
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,
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,
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, Reference Bahr, Pfeffer and Kaser2015). 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,
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, Reference Oerlemans2001) 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−1. These values were chosen to mimic typical Himalayan conditions (Laha and others, Reference Laha, Kumari, Singh, Mishra and Sharma2017). 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 = 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, Reference Kraaijenbrink and Bierkens2017; Banerjee and Wani, Reference Banerjee and Wani2018). 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, Reference Banerjee and Wani2018).
Empirical analysis of data from the Himalaya
For all the clean and debris-covered glaciers in the Himalaya that are larger than 1 km2, area (RGI, 2017) and estimated volume (Kraaijenbrink and others, Reference Kraaijenbrink and Bierkens2017) were fitted to Eqn (1) to obtain the corresponding volume-area scaling relation. The volume estimates used here (Kraaijenbrink and others, Reference Kraaijenbrink and Bierkens2017) were obtained from a variant of plastic-ice approximation (Frey and others, Reference Frey, Macguth and Huss2014). Available debris-cover extents (Kraaijenbrink and others, Reference Kraaijenbrink and Bierkens2017) 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, Reference Banerjee and Shankar2013). Glacier length, ice thickness, typical ablation rate and debris-thickness of the modelled glaciers were in the reasonable range as well (Kraaijenbrink and others, Reference Kraaijenbrink and Bierkens2017; Laha and others, Reference Laha, Kumari, Singh, Mishra and Sharma2017; Banerjee and Wani, Reference Banerjee and Wani2018). 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.
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, Reference Anderson and Anderson2016) 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.
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).
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, Reference Bahr1997), 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, Reference Bahr, Pfeffer and Kaser2015). 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 (Reference Bahr, Pfeffer and Kaser2015) γ = (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, Reference Nakawo, Iwata, Watanabe and Yoshida1986) even within a single glacier. It also varies with time due to stochastic and systematic changes in headwall-erosion rates (Banerjee and Wani, Reference Banerjee and Wani2018). Additional sources of debris due to rockfall, debris avalanche and degradation of lateral moraines (Nakawo and others, Reference Nakawo, Iwata, Watanabe and Yoshida1986) add further noise to the debris-thickness distribution (Shah and others, Reference Shah, Banerjee, Nainwal and Shankar2019). Complex mass-balance processes like avalanche activity (Laha and others, Reference Laha, Kumari, Singh, Mishra and Sharma2017), and stochastic increase in local ablation due to thermokarst features like supraglacial ponds and ice-cliffs (Sakai and others, Reference Sakai, Takeuchi, Fujita and Nakawo2000), 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) 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) 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) 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) 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) 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.
Acknowledgements
I acknowledge useful discussions with James Ferguson and R. Shankar. Critical comments from the anonymous reviewers and editor Hester Jiskoot led to significant improvements of the content and the presentation.
Appendix A
Scaling of ablation rate for glaciers with quadratic clean-glacier mass-balance profile
The general form of a quadratic mass-balance profile for clean glaciers is given by,
Here, z 0 is a reference elevation higher than the top of the bedrock, E is the ELA and β1 is a mass-balance parameter having unit of m−1 a−1. The corresponding sub-debris ablation is taken to be, b d(z, d) = b(z)/(1 + d/d 0) as before. Then, the sub-debris ablation at $x^\ast$ is,
Here, I have assumed that $\lpar z_0-E\rpar \gg s L^\ast$, and linearised Eqn (A1) around z s = E. Now, taking (z 0 − E) ~ s L,
Now, just as for the m = 1 case discussed in the main text, a steady-state condition implies, ${d_0 u^\ast _{s} }/{L^\ast }\sim {\alpha b^\ast _{\rm d}}$. Together the last two relations imply,
Or, m d = 1/2 + ((n + 1)(m d + 1))/2(n + 2), which can be simplified to obtain,
Where in the last equation n = 3 is used.
As shown in Appendix Figure A1, for idealised glaciers with the quadratic form of b(z s) (Eqn (A1)), m d = 1.16 ± 0.03, which is lower than the above predicted value of 1.5. However, with m d = 1.16 and q = 0, Eqn (2) predicts γd = 1.43, which matches the best-fit γd = 1.43 ± 0.00 obtained for the simulated glaciers. For the above quadratic mass-balance profile (Eqn (A1)), scaling of mean ablation in the simulated clean glaciers yielded m = 1.50 ± 0.03 (Fig. A1), which is lower than the naive expectation of m = 2. However, the best-fit γ = 1.52 ± 0.00 is again consistent the corresponding predicted value of γ = 1.5 (from Eqn (2), setting m = 1.5 and q = 0).
In the above numerical experiments, the best-fit ablation rate scaling exponents are about 25% smaller than the predicted values for clean and debris-covered glaciers. This is possibly related to the effects of ice-elevation feedback. However, given these best-fit m and m d, the predicted γ and γd match with the corresponding modelled values.