the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Possible biases in scalingbased estimates of glacier change: a case study in the Himalaya
Argha Banerjee
Disha Patil
Ajinkya Jadhav
Approximate glacier models are routinely used to compute the future evolution of mountain glaciers under any given climatechange scenario. A majority of these models are based on statistical scaling relations between glacier volume, area, and/or length. In this paper, longterm predictions from scalingbased models are compared with those from a twodimensional shallowice approximation (SIA) model. We derive expressions for climate sensitivity and response time of glaciers assuming a timeindependent volume–area scaling. These expressions are validated using a scalingmodel simulation of the response of 703 synthetic glaciers from the central Himalaya to a step change in climate. The same experiment repeated with the SIA model yields about 2 times larger climate sensitivity and response time than those predicted by the scaling model. In addition, the SIA model obtains area response time that is about 1.5 times larger than the corresponding volume response time, whereas scaling models implicitly assume the two response times to be equal to each other. These results indicate the possibility of a low bias in the scaling model estimates of the longterm loss of glacier area and volume. The SIA model outputs are used to obtain parameterisations, climate sensitivity, and response time of glaciers as functions of ablation rate near the terminus, massbalance gradient, and mean thickness. Using a linearresponse model based on these parameterisations, we find that the linearresponse model outperforms the scaling model in reproducing the glacier response simulated by the SIA model. This linearresponse model may be useful for predicting the evolution of mountain glaciers on a global scale.
 Article
(1297 KB) 
Supplement
(1103 KB)  BibTeX
 EndNote
In the coming decades, shrinking mountain glaciers will contribute significantly to global eustatic sealevel rise (e.g. Radić et al., 2014; Hock et al., 2019; Marzeion et al., 2020) and impact the hydrology of glacierised basins worldwide (e.g. Huss and Hock, 2018; Immerzeel et al., 2020). The reliability of the predicted changes in global sea level and those in the regional hydrology of various river basins is, thus, intimately tied to the accuracy of the predicted total ice loss from mountain glaciers for any given climate scenario.
Instantaneous (annual) glacier surface mass balance can be calculated readily using climate model outputs. In contrast, any prediction of the longterm evolution of a glacier requires simulating the slow (decadal) changes in glacier area and geometry. Ideally, this is to be done by solving the dynamical iceflow equations (e.g. Hutter, 1983). However, the numerical cost of such a computation on a global scale is high, even if simplified approximate descriptions of the iceflow equations, like shallowice approximation (SIA) (Hutter, 1983) or its higherorder variants, were to be used (Egholm et al., 2011; Clarke et al., 2015). Onedimensional SIAbased modelling tools are promising developments in this regard (Maussion et al., 2019; Zekollari et al., 2019; Rounce et al., 2020). The uncertainties associated with various input parameters, e.g. an uncertain glacier bedrock, limit the benefit of using the physically based iceflow models (Farinotti et al., 2016). Consequently, a majority of the recent estimates of the global to regional scale evolution of mountain glaciers relies on lowdimensional approximate parameterisations of glacier dynamics (e.g. Radić et al., 2014). The results from these simplified models have provided critical inputs for multimodel ensembleaveraged estimates of future sealevel rise (Hock et al., 2019; Marzeion et al., 2020), assessments of regional to global vulnerability to sealevel rise (e.g. Kulp and Strauss, 2019), and understanding the coevolution of glaciers, river runoff, and climate in glacierised regions like highmountain Asia (e.g. Zhao et al., 2014; Zhang et al., 2015; Kumar et al., 2019).
While some of the approximate parameterisations of glacier dynamics are empirical prescriptions for adjusting the hypsometry of the transient glaciers (Raper and Braithwaite, 2006; Huss et al., 2010; Huss and Hock, 2015), a majority of them are primarily based on a statistical volume–area (or volume–area–length) scaling relation. This volume–area scaling equation relates glacier volume V to glacier area A as
where γ is a dimensionless scaling exponent, and c is a scale factor (Bahr et al., 2015). This relation was established empirically (e.g. Chen and Ohmura, 1990) and subsequently proved using dimensional analysis (Bahr et al., 1997, 2015). The derivation utilised the empirical sublinear scaling of glacier width and ablation rate with the glacier length (Bahr, 1997).
Theoretically, the scaling exponent γ is timeindependent and can be expressed as $\mathit{\gamma}=\mathrm{1}+\frac{m+\mathrm{1}}{m+n+\mathrm{3}}$ (Bahr et al., 2015). Here, n is the powerlaw exponent of Glen's rheology of ice (Glen, 1955), and m is the scaling exponent of ablation rate with glacier length (Bahr, 1997). For an individual glacier, the scalefactor c captures the control of all the glacierspecific factors (except area) on its volume (Bahr et al., 2015). There is no available theoretical prescription for obtaining the value of c for an arbitrary glacier. c may be calibrated for a particular glacier based on available independent measurements of area and volume over an epoch, but its time dependence can be accessed only with a detailed model simulation (Bahr et al., 2015). For a large enough ensemble, glacier area typically spans a few orders of magnitude. However, the corresponding c values vary over a relatively restricted range (Bahr et al., 2015). This allows an approximate statistical description of any set of glaciers using Eq. (1), where a single bestfit c and a fixed γ are used (Bahr et al., 2015). Such a bestfit scaling relation provides a fairly accurate estimate of the total ice volume of a large set of glaciers, but the corresponding predictions for the individual glaciers have relatively large uncertainties (Bahr et al., 2015). Note that there is no theoretical constraint for c to be timeindependent for a given set of nonsteady glaciers (Bahr et al., 2015).
It is the above statistical interpretation of the scaling relation, where a bestfit timeinvariant c and a constant γ are used to describe an ensemble of glaciers, that is exploited in the scalingbased approximate models of glacier dynamics (e.g. Radić et al., 2007). Hereinafter, we refer to the models based on such an approach (e.g. Radić et al., 2007) as “scaling models”. As the present study investigates the possibility of biases in scaling model predictions of glacier evolution, we restrict ourselves to the above statistical interpretation of the scaling relation.
The performance of scaling models in simulating the transient glacier response has previously been tested against various dynamical iceflow models (e.g. SIA, higherorder approximations, or Stokes' model) in one to three dimensions using both idealised (Radić et al., 2007; Adhikari and Marshall, 2012) and realistic geometries (Radić et al., 2008; Farinotti and Huss, 2013). The uncertainties introduced by a scalingmodel parameterisation of the evolution of glaciers with realistic geometries were considered by Farinotti and Huss (2013). The spirit of the present study is quite similar to that of Farinotti and Huss (2013), except that we are investigating the possible intrinsic biases of scaling models in a situation where the parameters (c and γ) are known accurately. The specific objectives of the present study are

to obtain analytical predictions for climate sensitivity and response time of glaciers in a scaling model;

to compare the climate sensitivity and response time of a large number of synthetic glaciers with realistic geometries, as obtained from a scaling model and a 2D SIA model;

to investigate the possibility of longterm biases in scaling model estimates of changes in glacier area and volume with respect to corresponding SIA results;

to find convenient parameterisations of glacierresponse properties obtained from the SIA simulations and develop an accurate linearresponse model.
Note that a linearresponse model introduced in the last objective is a lowcomplexity model obtained in the limit of a relatively small deviation around a steady state (e.g. Oerlemans, 2001). To apply this model to a large number of glaciers, the response time and climate sensitivity need to be specified for each of them. A lack of accurate and numerically convenient parameterisations of these dynamical properties may have limited its application (Harrison et al., 2001; Lüthi, 2009; Bach et al., 2018). Here, we aim to obtain parameterisations of the glacierresponse properties as functions of a few easily accessible properties of the glaciers, using results from 2D SIA simulations of a large ensemble of synthetic glaciers with realistic geometries.
The paper is organised as follows. First, we theoretically derive the glacierresponse properties within a timeinvariant scaling assumption (Sects. 2.1 and 3.1). Then, we compare the performance of a representative scaling model (Radić et al., 2007) with that of a twodimensional SIA model, in simulating the response of 703 idealised Himalayan glaciers in the Ganga basin to a hypothetical step rise in equilibrium line altitude (ELA) (Sects. 2.2 and 3.2). We use the response properties obtained from the scaling model to test the above analytical expressions for glacierresponse properties. The corresponding SIA results are used to obtain parameterisations for the linearresponse properties of realistic glaciers. The accuracy of the scaling model and a linearresponse model in reproducing the SIAderived longterm loss of total glacier area and volume is assessed for the above 703 glaciers. The performance of the linearresponse model is also tested for an independent set of 164 glaciers in the western Himalaya without any further calibration. We also discuss the applicability of the linearresponse model for actual computation of future glacier loss for a set of transient glaciers forced by any arbitrary timevariation ELA (Sect. 3.3).
2.1 Theoretical methods
For a theoretical analysis of the glacierresponse properties implied by a scaling model, we consider a set of hypothetical glaciers that respond to a warming climate such that the volume–area scaling relation (Eq. 1) is valid, and c is a given timeinvariant constant. Then, the fractional changes in area and volume of these glaciers, in the limit of small changes, are related as follows.
where ΔV and ΔA are the changes in area and volume, and the mean ice thickness is $h=V/A$. The above equation is the basis of the scaling models of glacier evolution (e.g. Radić et al., 2007). We have derived analytical expressions for glacierresponse time and climate sensitivity starting from this equation, essentially following the line of arguments by Harrison et al. (2001).
2.2 Numerical methods
We simulated the response of an ensemble of synthetic clean glaciers with realistic geometries to a hypothetical step change in ELA using three different methods (scaling, SIA, and linearresponse models). For this exercise, we considered all the 814 glaciers larger than 2 km^{2} in the Ganga basin, the central Himalaya (Fig. S1 in the Supplement). The icefree bedrock for each glacier was obtained using available icethickness estimates (Kraaijenbrink et al., 2017) and surfaceelevation data (NASA et al., 2019). The following idealised elevationdependent linear massbalance profile was used,
Here β is the balance gradient, z is the surface elevation, and E is the equilibriumline altitude (ELA). b_{0} is a cutoff on maximum accumulation taken to be 1.0 m yr^{−1}. The choice of β is described later. In our massbalance model, we neglected complicating factors like supraglacial debris cover and its effects on ablation, and the avalanche contribution to accumulation (Laha et al., 2017). The debris effects are expected to modify the scaling relations as well (Banerjee, 2020). Overall, the simulated glaciers cannot be considered faithful copies of the actual Himalayan glaciers. Rather, they constituted an ensemble of synthetic clean glaciers with realistic geometries (e.g. Farinotti and Huss, 2013) to be used here for a comparative study of the performances of the three models.
2.2.1 A 2D SIA model
The iceflow dynamics was implemented within a twodimensional SIA (Hutter, 1983) as a numerically efficient nonlinear diffusion problem (Oerlemans, 2001). While SIA may not be the best method for simulating valley glaciers due to its limitation in describing ice flow influenced by longitudinal stresses and/or steep bedrock slopes (Le Meur et al., 2004), there is enough evidence in the literature that SIA does a reasonable job of describing both the steady and transient dynamics of valley glaciers (e.g. Leysinger Vieli and Gudmundsson, 2004; Le Meur et al., 2004; Radić et al., 2008). The contribution of sliding to the flow was neglected here for simplicity.
The value of Glen's flowlaw exponent was assumed to be 3 (e.g. Oerlemans, 2001). For the sake of simplicity, we did not tune any of the model parameters to match the observed ice thickness and/or flow velocity on any of these glaciers. The only exception was ELA which was tuned to obtain the initial steady state as described below. In order to avoid possible dependence of the results on any specific choice of parameters, we picked the parameters related to mass balance and flow from random distributions. The rate constant of Glen's law was picked randomly from the set $\mathit{\{}\mathrm{0.5},\mathrm{0.6},\mathrm{\dots},\mathrm{1.4},\mathrm{1.5}\mathit{\}}\times {\mathrm{10}}^{\mathrm{24}}$ Pa^{−3} s^{−1} for each of the glaciers. This range of values is comparable to those used to model mountain glaciers previously (Radić et al., 2008). The balance gradient β was also picked randomly from the set of values $\mathit{\{}\mathrm{0.005},\mathrm{0.006},\mathrm{\dots},\mathrm{0.009},\mathrm{0.010}\mathit{\}}$ yr^{−1} for each glacier. This range of β values is comparable to the observed massbalance gradients in the Himalaya (e.g. Wagnon et al., 2013).
The model was integrated using a linearised implicit finitedifference scheme (Hindmarsh and Payne, 1996), with a noslip boundary condition at the ice–bedrock interface and a noflux boundary condition at the domain boundary. An iterative conjugategradient method was employed within the implicit scheme, with a spatial grid size of 100 m × 100 m and time steps of 0.01 years. To avoid the known problem of a possible violation of mass conservation in SIA on steep terrains (Jarosch et al., 2013), we smoothed the bedrock with a centrally weighted 3×3 movingwindow averaging. In addition, the conservation of ice was explicitly monitored by tracking the total accumulation and ablation on the glacier surface and the ice flux out of the glacier boundary in the ablation zone. The cumulative net gain of ice matched the total ice in the domain to within one part per 10^{9} at any time t. Only on three glaciers (out of the total of 814) was a violation of conservation due to steep bedrock observed, and these three were not considered in our analysis (Fig. S2). One more glacier had to be removed where an erroneously mapped truncated tributary led to an unrealistic piling up of ice (Fig. S2).
The SIA simulation was run starting with an empty bedrock, with the initial E being the median elevation. The simulation was continued until an approximate steady state was reached such that the absolute value of the net specific balance was less than 10^{−4} m yr^{−1}. Subsequently, E was moved up or down, and the simulation was repeated until the extent of the steady state was similar to the present glacier extent (RGI Consortium, 2017) (Fig. S2). Once the desired steady state was found (see Fig. S3 for a few examples), the glaciers were perturbed by a 50 m step rise in ELA. Subsequently, the annual values of area and volume were recorded for the next 1000 years (Fig. S4). The mean and standard deviation of the modelled ELA for these 810 glaciers were 5480 and 445 m, respectively.
Out of the total 810 simulated glaciers from the Ganga basin, on 98 glaciers the fractional change in glacier area at t=1000 was more than 50 %, and these were excluded from the analysis. This was necessary as a linearresponse model can only be applied to glaciers with small relative changes (Oerlemans, 2001). We confirmed that the nature of our results does not depend on the precise value of this cutoff (Fig. S6). An additional nine glaciers had response time larger than 500 years and they were removed. This was done to avoid a possible overestimation of the response time whenever its magnitude was comparable to or larger than the total simulation period of 1000 years (Fig. S7). The removal of these nine glaciers led to a reduction in the number (total area) of simulated glaciers by only ∼1 % (∼2 %).
Finally, we were left with an ensemble of 703 synthetic Himalayan glaciers (Fig. S1), with area in the range of 2.2–156.0 km^{2} (a median value 5.5 km^{2}). The steady glaciers modelled with SIA had, on average, 1.25 times larger area and 1.66 times larger ice thickness (Figs. S3, S8) compared to the corresponding estimates of Kraaijenbrink et al. (2017). The higher thickness of the modelled glaciers can be ascribed to a larger modelled area, a steady mass balance, and an uncalibrated SIA model. The total area and volume of these 703 synthetic glaciers were 6865 and 847 km^{3}, respectively. This set covered 86 % of the total 810 glaciers numberwise and 89 % areawise. The distributions of glacier area and mean slope for the two sets of 810 and 703 synthetic glaciers are shown in Fig. S8.
2.2.2 Scaling model
The response of the above set of 703 steadystate glaciers to a 50 m instantaneous rise in ELA was also computed with a scaling model (Radić et al., 2007). The SIAderived initial steadystate volume, area, and hypsometry (with the bin size of 25 m) for each of the glaciers were used as the starting point. For any of the modelled glaciers, the scaling and SIA models used the same massbalance parameters. At any time t during the evolution, the massbalance function (Eq. 3) was summed over the instantaneous glacier hypsometry to obtain the net volume loss for that particular time step. The corresponding area loss was then obtained using Eq. (2). The reduction in the area was assumed to have taken place in the lowest elevation band/s of each glacier (Radić et al., 2007). The scaling exponent was fixed at γ=1.286 because of the assumed linear massbalance profiles of the simulated glaciers (i.e., m=1). The annualresolution time series of area and volume were recorded for 1000 years for each of the glaciers.
2.2.3 Glacierresponse properties
For each of the 703 glaciers, the time series of volume and area as obtained using the SIA and scaling models were separately fitted to linearresponse forms (e.g. Eq. 9 below) to obtain the corresponding bestfit values of the four linearresponse parameters (the climate sensitivities and the response times of area and volume) for each of them (Fig. S4).
Note that applying a step change in ELA of a steadystate glacier to obtain the stepresponse function is a standard prescription for obtaining glacierresponse properties (Oerlemans, 2001; Leysinger Vieli and Gudmundsson, 2004; Harrison et al., 2001; Bach et al., 2018). Within a linearresponse assumption, the step responses of volume and area have an exponential form (e.g. Eq. 9 below). The asymptotic exponential decay time is the response time of the glacier, and the asymptotic magnitude of the decay is the climate sensitivity. Because of the deviations of the simulated response from a pure exponential decay (Fig. S4), the bestfit response time may be slightly different from the efolding time, which has been used in some of the previous studies (e.g. Leysinger Vieli and Gudmundsson, 2004; Bach et al., 2018). However, we take the bestfit asymptotic decay time to be the response time. By definition, it minimises the deviation between the predictions of the SIA and linearresponse models and thus improves the performance of the latter in reproducing SIA results to some extent. We confirm that the difference between the above two definitions of the response time is small.
The bestfit linearresponse properties obtained from the scaling model results for the 703 glaciers were used to verify the corresponding theoretical expressions obtained from scaling theory (Eqs. 8, 11, 12, 13 below). The bestfit response times and climate sensitivities obtained from the SIA simulations of the 703 glaciers were used to fit for empirical relations that are motivated by the corresponding expressions derived from the scaling theory. These fitted forms would allow estimation of the response properties of any given glaciers as functions of properties like mean thickness, mass balance gradient, and so on. All the above fits were performed on a loglog scale, and R^{2} values of the fits were noted.
2.2.4 A linearresponse model
The bestfit empirical parameterisations for climate sensitivity and response time obtained by fitting the SIA results as described above (given in Eqs. 14–17 later) were used to run a linearresponse model simulation for any given glacier. This model was applied to simulate the response of the above 703 synthetic Himalayan glaciers to a 50 m step change in ELA at t=0. We emphasise that for the linearresponse model, we did not use the bestfit response properties of the individual glacier derived from the SIA simulations. Rather, the parameterisations of the same obtained by fitting the SIAderived response properties (given in Eqs. 14–17 later) were utilised. These parameterisations thus allow the model to be applied to any other set of Himalayan glaciers without the need for simulating them with SIA first.
To assess the uncertainty of the linearresponse model output, the uncertainty of each of the fit parameters was set equal to the corresponding standard error, and the 95 % uncertainty band for the linearresponse model outputs was generated using a Monte Carlo method.
To test the applicability of the above linearresponse model that was calibrated using SIA results for the 703 central Himalayan glaciers, the same model was applied to a different set of 204 glaciers from the western Himalaya. The parameterisations developed for the central Himalayan glaciers as discussed above (given in Eqs. 14–17 later) were used to estimate the response properties of each of these western Himalayan glaciers using input values of corresponding massbalance gradient, mean thickness, and ablation rate near the terminus. For these western Himalayan glaciers, SIA and scaling model simulations were also performed following the procedures as detailed above. The glaciers showing more than 50 % change at the 500year mark in the corresponding SIA simulations were left out as before, and the time series of total area and total volume of 164 western Himalayan glaciers obtained using the three different models were compared.
3.1 Theoretical results
Below, we derive some relevant consequences of the timeinvariant scaling assumption, including expressions for the climate sensitivity and response time of glacier area and volume. These results are expected to be generally valid for all scaling models that are based on Eq. (2).
3.1.1 The rates of area and volume change
Equation (2), which was derived from Eq. (1) assuming a timeindependent c, implies
Here $\dot{V}$ and $\dot{A}$ denote the corresponding rates of change of glacier volume and area, respectively. If the net specific balance is δb (m yr^{−1}), then the annual rate of volume loss $\dot{V}=\mathit{\delta}bA$. This, together with Eq. (4), implies
Thus, in the scaling models, the rate of change of area scales with the glacier area with an exponent (2−γ). This is consistent with empirical observations for real glaciers as well (Banerjee and Kumari, 2019). As the scale factor $\frac{\mathit{\delta}b}{\mathit{\gamma}c}$ on the righthand side (RHS) of Eq. (5) is proportional to the net specific mass balance, this may be a convenient way of obtaining mean regional thinning rates from relatively straightforward remotesensing measurements of the rate of area change. However, the accuracy of this relation is contingent on the validity of the assumption of a timeindependent c.
3.1.2 Area response time
To compute the area response time, let us consider a constant perturbation, i.e., a step change in ELA applied to a steady glacier for time t≥0 (e.g. Oerlemans, 2001). Let's denote the corresponding instantaneous net negative balance at t=0 by δb_{0}A, the asymptotic (t→∞) shrinkage of glacier area by $\mathrm{\Delta}{A}_{\mathrm{\infty}}\equiv A\left(\mathrm{0}\right)A(t\to \mathrm{\infty})$, and that of ice volume by ΔV_{∞}. Then, we have (Harrison et al., 2001)
Here b_{t} is the ablation rate near the terminus. The area response time of the glacier can be expressed as ${\mathit{\tau}}_{A}\approx \mathrm{\Delta}{A}_{\mathrm{\infty}}/\dot{A}$. Therefore, using the above expressions for $\dot{A}$ (Eq. 5) and ΔA_{∞} (Eq. 7), we obtain
Here the symbol τ^{∗} is a convenient shorthand notation for the timescale ${\left(\frac{{b}_{\mathrm{t}}}{\mathit{\gamma}h}+\mathit{\beta}\right)}^{\mathrm{1}}$. In the above derivation, ΔV_{∞} that appears in Eq. (7) is eliminated with the help of Eq. (2). Equation (8) is comparable with the expression of area response time as given by Harrison et al. (2001) or Lüthi (2009).
3.1.3 Volume response time
The instantaneous change in volume (ΔV(t)) for a steady glacier perturbed by a small step change in ELA at t=0 is given by
where τ_{v} is the volume response time and ΔV_{∞} is the volume sensitivity (e.g. Lüthi, 2009). Now, V(t),V(0), and V(t→∞) appearing in Eq. (9) can be expressed in terms of A(t),A(0), and A(t→∞), respectively, with the help of corresponding scaling relations (Eq. 1). This, in the limit of small fractional changes in area, yields
Comparing the above two equations and using Eq. (8), one obtains
Thus, all scaling models implicitly assume the area and volume response times of a glacier to be equal to each other. However, it is known that for mountain glaciers area response time is larger than the volume response time within an SIA model (Oerlemans, 2001; Leysinger Vieli and Gudmundsson, 2004). Therefore, the assumed equality of the two response times in scaling models (Eq. 11) contradicts the existing SIA results. This is an intrinsic bias that is present in any scaling model.
After a step change in ELA, as the ablation zone shrinks, the initial net negative balance of a glacier gradually decays to zero over a period determined by the corresponding response time. A longer area response time in SIA implies that this reduction in the ablation zone is slower here than that in a scaling model. A corresponding feedback of a larger ablation zone on the net mass balance should then lead to a higher longterm volume loss in an SIA model than that in a scaling model. This indicates the possibility of a low bias in scaling model estimates of the climate sensitivity of volume, or equivalently, that in the longterm changes in glacier volume due to any rise in ELA.
3.1.4 Climate sensitivity of area and volume
An expression for the climate sensitivity of glacier area (ΔA_{∞}), which is the asymptotic change in area due to a change in ELA by δE, is obtained by eliminating ΔV_{∞} from Eq. (7) using Eq. (2),
Here we have used the definition of τ^{∗} (Eq. 8), and that of δb_{0}≈βδE for a step change in ELA by δE. The RHS of the above equation is denoted by α^{∗} for convenience.
The corresponding expression for $\frac{\mathrm{\Delta}{V}_{\mathrm{\infty}}}{V}$ is then obtained using Eq. (2),
Again, Eq. (13) is comparable to the expression of volume sensitivity as derived by Harrison et al. (2001), where the authors used an arbitrary thickness scale H, instead of the denominator of γh appearing in the definition of α^{∗} above.
Strictly speaking, the climate sensitivity of area and volume with respect to a change in ELA should be defined as $\frac{\mathrm{\Delta}{A}_{\mathrm{\infty}}}{\mathit{\delta}E}$ and $\frac{\mathrm{\Delta}{V}_{\mathrm{\infty}}}{\mathit{\delta}E}$, respectively. However, in this paper, we use ΔA_{∞} and ΔV_{∞} as the corresponding sensitivities to simplify the notation.
3.2 Numerical results
3.2.1 Volume–area scaling and a timedependent scale factor in the SIA model
Following Eq. (1), a powerlaw relation between the area and volume of the 703 glaciers with an exponent $\mathit{\gamma}=\mathrm{1}+\frac{m+\mathrm{1}}{m+n+\mathrm{3}}=\mathrm{1.286}$ is expected as m=1 and n=3. The ensemble of glaciers modelled with SIA did conform to the above powerlaw relation V=cA^{1.286} at any time t with a single bestfit c. The scale factor slowly decreased with time. For example, Fig. 1a shows the powerlaw fits at t=0 and t=500 years (R^{2}=0.9), where the bestfit c values were 0.053±0.001 and 0.47±0.001 km^{3−2γ}, respectively. This implies a ∼11 % reduction in c for the ensemble over the period of 500 years after the step change in ELA was applied. A timedependent c is consistent with the theoretical arguments of Bahr et al. (2015).
The slow and systematic decline in c for the ensemble of shrinking glaciers simulated with the SIA model contradicts the basic assumption of scaling models that c is timeinvariant. A decreasing c would mean Eq. (2) is violated, with $\frac{\mathrm{\Delta}V}{V}=\mathit{\gamma}\frac{\mathrm{\Delta}A}{A}+\frac{\mathrm{\Delta}c}{c}$. Note that all three fractional changes involved in this relation are negative. Therefore, for any given $\left\mathrm{\Delta}A\right$, the corresponding $\left\mathrm{\Delta}V\right$ is going to be larger in the SIA model than that in a scaling model where $\frac{\mathrm{\Delta}c}{c}$ is assumed to be zero (Eq. 2). Even though the decline in c is only about 11 %, it may be associated with a stronger low bias in the longterm change predicted by scaling models. This is because a larger volume change in SIA would lead to a thinner glacier, and a corresponding surfaceelevation feedback to mass balance is likely to amplify the corresponding longterm mass loss over time.
The dependence of the glacierspecific scale factor on the mean slope is known (Bahr et al., 2015) and has been incorporated in modified scaling relations where volume is a powerlaw function of both area and slope (e.g. Grinsted, 2013; Zekollari and Huybrechts, 2015). For the simulated 703 glaciers, the mean slope increases with time as the area is lost preferentially from the gently sloping lower ablation zone. For example, the median slope of the 703 simulated glaciers decreased from 0.41 at t=0 to 0.37 at t=500 years. This ∼10 % reduction in slope is expected to lead to a ∼5 % decline in c (Bahr et al., 2015). So, at least part of the time dependence of c for transient glaciers in SIA simulation is explained by the slopedependence of c. However, there may be other factors contributing to the decline in c for the transient glaciers as discussed below.
3.2.2 Area and volume response times
The theoretical predictions for glacier area and volume response time (Eq. 11) worked rather well for the scaling model results (Fig. 2c and d), with bestfit relations of ${\mathit{\tau}}_{V}=(\mathrm{0.914}\pm \mathrm{0.002}){\mathit{\tau}}_{A}$ with R^{2}=0.99 and ${\mathit{\tau}}_{A}=(\mathrm{1.066}\pm \mathrm{0.008}){\mathit{\tau}}^{\ast}$ with R^{2}=0.80.
For SIA simulations, the data showed that τ_{A}>τ_{V} and that the two response times were still proportional to each other (Fig. 3c: ${\mathit{\tau}}_{V}=(\mathrm{0.687}\pm \mathrm{0.004}){\mathit{\tau}}_{A}$, with R^{2}=0.94). Also, τ_{A} was proportional to τ^{∗} to a good approximation (Fig. 3d: ${\mathit{\tau}}_{A}=(\mathrm{2.56}\pm \mathrm{0.04}){\mathit{\tau}}^{\ast}$, with R^{2}=0.53). Interestingly, the value of the proportionality constant in the latter relation as obtained from SIA was about 2.4 times larger than the corresponding value obtained in the scaling model. This underlines the relatively large underestimation of area response time by the scaling model. Similarly, the volume response time was about 1.8 times larger in the SIA simulation than the corresponding scaling model value. This implies that for a given ELA perturbation, the glacier response is much faster in the scaling model compared to that in the SIA model for the ensemble of 703 synthetic glaciers.
Apart from the overall underestimation of area and volume response times by the scaling model, another serious limitation of scaling models that emerges from the above analysis is that here the area and volume response times are equal to each other (Eq. 11 and Fig. 3c). In contrast, the SIA model predicted τ_{A}≈1.5τ_{V}. The ratio of the two response times obtained from the 2D SIA model here is generally consistent with earlier results based on 1D flowline models (Oerlemans, 2001; Leysinger Vieli and Gudmundsson, 2004). The equality of the two response times in the scaling model led to a linear trajectory in the V−A plane for the transient glaciers (Fig. 1b). In comparison, a relatively larger area response time, together with slow initial changes in area (Figs. S4, S10), led to curved V−A plane trajectories for individual transient glaciers in the SIA model simulations. In particular, a slowly changing area means the V−A trajectories bend downward, causing c to decrease for the transient ensemble (Fig. 1). Moreover, at the early stages of response, glaciers simulated by a scaling model lose area much quicker than those simulated by an SIA model (Fig. 1b). The associated net massbalance feedbacks then lead to a subdued longterm volume response in the scaling model and a comparatively stronger volume response in the SIA model, just as predicted in Sect. 3.1.3.
3.2.3 The climate sensitivity of glacier area and volume
For the 703 glaciers simulated by the scaling model, the fitted asymptotic fractional changes in area and volume, or equivalently the corresponding (fractional) climate sensitivities, were proportional to each other (Fig. 2a: $\frac{\mathrm{\Delta}{V}_{\mathrm{\infty}}}{V}=(\mathrm{1.383}\pm \mathrm{0.003})\frac{\mathrm{\Delta}{A}_{\mathrm{\infty}}}{A}$, with R^{2}=0.99). Here, the bestfit constant of proportionality is close to but about 8 % larger than γ=1.286 as predicted by Eq. (2).
In contrast, the SIA simulations obtained $\frac{\mathrm{\Delta}{V}_{\mathrm{\infty}}}{V}=(\mathrm{1.93}\pm \mathrm{0.02})\frac{\mathrm{\Delta}{A}_{\mathrm{\infty}}}{A}$, with R^{2}=0.85 (Fig. 3a). In this case, the constant of proportionality was ∼1.5γ, compared to the corresponding value of ∼γ in the scaling model. This larger value of the ratio of the two climate sensitivities in the SIA model is consistent with the observed decline in c for the transient glaciers simulated with this model (Fig. 1). Note that no theoretical prediction is available for the ratio of asymptotic fractional changes in volume and area in an SIA model.
Figure 2b shows that in the scaling model, climate sensitivity of glacier volume is proportional to α^{∗} $\left(\frac{\mathrm{\Delta}{V}_{\mathrm{\infty}}}{V}=(\mathrm{0.655}\pm \mathrm{0.008}){\mathit{\alpha}}^{\ast},\text{with}{R}^{\mathrm{2}}=\mathrm{0.67}\right)$. This is in line with Eq. (13), except that the constant of proportionality is significantly less than γ. A similar proportionality between the SIAderived bestfit $\frac{\mathrm{\Delta}{V}_{\mathrm{\infty}}}{V}$ and α^{∗} is shown in Fig. 3b, with $\frac{\mathrm{\Delta}{V}_{\mathrm{\infty}}}{V}=(\mathrm{1.71}\pm \mathrm{0.03}){\mathit{\alpha}}^{\ast}$. However, in this case the fit is relatively noisy with R^{2}=0.48.
The above relations suggest that the climate sensitivity of volume in the SIA simulation was about 2.6 times larger than that in the scaling model. Similarly, the climate sensitivity of glacier area obtained from the SIA model was also about 1.9 times larger than that obtained from the scaling model. This trend of a relatively large (by about a factor of 2) underestimation of climate sensitivity of glacier volume and area by the scaling model is consistent with the effects of a relatively faster shrinkage of the ablation zone in the early stages of the response as discussed in Sect. 3.1.3 and 3.2.2.
3.2.4 The total glacier loss estimated using the three models
Starting with an initial volume (area) of 847 km^{3} (6865 km^{2}), the 703 glaciers simulated by SIA lost a total of 194 km^{3} (726 km^{2}) of volume (area) in 500 years due to the step rise in ELA by 50 m. As shown in Fig. 4, both the scaling and the linearresponse models underestimated the longterm change in total area in this experiment, with estimated area changes of 334 and 623 km^{2}, respectively. The scalingmodel prediction for area change was only 46 % of the corresponding SIA estimate, while the linearresponse model estimate was 86 % of that of SIA. Similar trends were seen for the magnitudes of estimated volume change as well, with the respective scaling and linearresponse model estimates being ∼31 % and ∼75 % of the corresponding SIA prediction (Fig. 4). We confirmed that the nature of the above results does not depend on the chosen cutoff of 50 % change that was used to select the 703 glaciers (Fig. S6). In fact, with a smaller cutoff, the linearresponse model estimates were even closer to the corresponding SIA estimates (Fig. S6). This is expected as linearresponse models are derived in the limit of small fractional changes (Oerlemans, 2001).
The low bias in the longterm changes of glacier area and volume computed with the scaling model is consistent with the underestimation of corresponding climate sensitivities by this model (Sect. 3.2.3). On shorter timescales of multiple decades, an underestimation of response times by about a factor of 2 (Sect. 3.2.2) partly compensates for a corresponding underestimation of the climate sensitivities (Sect. 3.2.3), and the deviations between the SIA and the scaling model are not that prominent (Fig. 4). The biases in the scaling model become clearer over multiple centuries (Fig. 4).
Note that, depending on the details of the scaling and SIA models compared, or the set of glaciers simulated, the actual magnitude of the biases in scalingmodelderived climate sensitivity, response time, and longterm glacier change could be different from these here. However, based on the theoretical arguments and numerical evidence presented, similar qualitative trends are expected if the above exercise were to be repeated with a more detailed model and/or for a more realistic set of glaciers.
The above results indicate the possibility of a negative bias in scaling model estimates of future changes in mountain glaciers and the corresponding contribution to sealevel rise. As an example, let us consider a recent comparison (Hock et al., 2019) of projected endofthecentury sealevel rise contribution of glaciers from six different models, with five of them being based on some form of scaling. In that intercomparison study, the hypsometricadjustmentbased model (Huss and Hock, 2015) consistently predicted the largest fractional change of global glacier volume and area under various climate scenarios (Table 3 of Hock et al., 2019). In another recent comparison, similar trends are seen as far as globalscale fractional volume loss by 2100 are concerned (Figs. S17–S20 of Marzeion et al., 2020), although on a regional scale there are differences. However, it is difficult to draw a definite conclusion about any potential bias in scaling models from the abovementioned studies as there are wide differences among the model runs in terms of the initial conditions, climate forcing, and massbalance parameterisations used. An intercomparison of the models where the same set of glaciers, with the same initial geometry and volume, are simulated under the same massbalance forcing – similar to the strategy used in the present study – is necessary to identify possible biases in the existing scaling models.
The above results show that the linearresponse model outperformed the scaling model, producing a closer match with the SIA results for the 703 synthetic glaciers from the Gangetic Himalaya. However, this linearresponse model was calibrated using the SIA results for the same set of glaciers. Therefore, this match is not enough to establish the effectiveness of the linearresponse model. To confirm the improved performance of the linearresponse model compared to that of the scaling model, we applied both the models to simulate a different set of 164 glaciers in the western Himalaya (Fig. S1). The bestfit linearresponse properties obtained from SIA simulation of the 703 central Himalayan glaciers were first fitted to obtain four equations (Eqs. 14–17) that relate the response properties to $\mathit{\beta},\mathit{\gamma},h$, and b_{t} as described before. The same equations were used to estimate the response properties of each of the 164 western Himalayan glaciers as required for the linearresponse model simulations. In this independent experiment, the linearresponse model again outperformed the scaling model in reproducing the corresponding SIA results (Fig. S9). This confirms that the linearresponse model, along with Eqs. (14)–(17), can be used for computing longterm glacier changes accurately.
3.3 The effects of glacier geometry
Can the biases in the scaling model described above be artefacts arising out of some peculiarities of the geometry of the specific set of glaciers being simulated and not relevant in general for scaling model computations of globalscale mass loss of mountain glaciers? To rule out this possibility, we simulated the response of a set of highly idealised synthetic glaciers using both a flowline model (Banerjee, 2017) and the above scaling model (Radić et al., 2007). Note that this flowline model included sliding as well. All of these synthetic glaciers have the same constant width, the same linear bedrock with a constant slope, and the same linear massbalance profile. Only the ELA was varied between glaciers. Even for this highly idealised set of glaciers, the scaling model estimates for the evolution of total area and volume showed biases compared to those obtained from the flowline model (Fig. S9), and these biases were qualitatively very similar to those depicted in Figs. 1 and 4. Again, the scaling model predicted relatively smaller climate sensitivities, a relatively faster area response, and a low bias in the longterm changes, compared to corresponding flowlinemodel estimates (Fig. S9).
The above flowlinemodel experiment provides an additional piece of evidence that the scalingmodel biases discussed in this paper are, in general, expected to be present in scaling model simulations of any set of glaciers. We reemphasise that even though the biases are expected to be qualitatively similar to those presented here, the magnitudes of the biases are likely to depend on the detailed characteristics (related to geometry, flow, and massbalance processes) of the glaciers studied and the models used.
3.4 The linearresponse model and its application to real glaciers
As described above, we have used results from the 2D SIA model simulations of the response of 703 synthetic Himalayan glaciers to a 50 m step change in ELA, to obtain the following bestfit parameterisations of the glacierresponse properties (i.e., $\frac{\mathrm{\Delta}{V}_{\mathrm{\infty}}}{V},\frac{\mathrm{\Delta}{A}_{\mathrm{\infty}}}{A},{\mathit{\tau}}_{A}$, and τ_{V}).
Here as defined before, ${\mathit{\tau}}^{\ast}\equiv {\left(\frac{{b}_{\mathrm{t}}}{\mathit{\gamma}h}+\mathit{\beta}\right)}^{\mathrm{1}}$, ${\mathit{\alpha}}^{\ast}\equiv \frac{\mathit{\beta}\mathit{\delta}E{\mathit{\tau}}^{\ast}}{\mathit{\gamma}h}$, and δE=50 m.
With the estimated glacierspecific response properties obtained from Eqs. (14)–(17), it is possible to compute the evolution glacier volume and area accurately for any glacier and for any arbitrary ELA forcing function. For this, the following general solution of the linearresponse equation is used.
Here ΔE(t) is the given (arbitrary) ELA forcing function. This equation simply states that any continuous ELA change can be interpreted as the sum total of a series of discrete impulses, and the corresponding net response is given by a superposition of suitably delayed responses due to each of the impulses. An analogous expression can be obtained for the area evolution just by replacing all the V values in the above equation with A values.
Note that the above formulation does not require the initial state to be steady. As long as the glacier is close to a steady state, a linearresponse theory will be a good approximation (Oerlemans, 2001). However, an additional initial condition, i.e., the value of ΔV(0), is needed to apply the linearresponse model to transient glaciers. ΔV(0) is the initial departure from a steady state and can be obtained from the observed rate of volume loss ($\dot{V}$) simply as $\mathrm{\Delta}V\left(\mathrm{0}\right)={\mathit{\tau}}_{V}\dot{V}$. Thus, the linearresponse model can be used to evolve the area and volume of a real set of glaciers for any arbitrary timedependent ELA forcing given the initial rates of change of volume and area, initial thickness, massbalance gradient, and melt rate near glacier terminus.
Since the above parameterisation of linearresponse properties (Eqs. 14–17) is derived from SIA simulations of an ensemble of Himalayan glaciers, when applying them to any other glacierised region in the world, it may be necessary to simulate a few tens of glaciers (having a representative range of area and slope) from that region using SIA first and confirming the accuracy of the above parameterisations.
Due to the noise present in the fits (Fig. 3), the linearresponse model predictions for an individual glacier would have significant uncertainties. However, for a large set of glaciers, the linearresponse model provides accurate estimates of the total area and volume evolution (Fig. 4, S6 and S9).
3.5 Limitations of the present study
Because of the idealised descriptions of ice flow and the massbalance profile (as discussed in Sect. 2.2), and the absence of model calibration to match the available observed data of surface velocity, ice thickness, recent mass balance, etc., the glaciers simulated here are not faithful copies of the Himalayan ones. For a set of more realistic glaciers, the magnitudes of the corresponding biases in scalingmodelderived climate sensitivity and response time could be different from those obtained here. However, based on the theoretical arguments and numerical evidence presented, similar qualitative trends are expected if the above exercise were to be repeated for a more realistic model that includes higherorder mechanics, a more realistic massbalance model, and so on. Similarly, the parameterisations for the linearresponse properties given here are obtained from 2D simulations of 703 synthetic Himalayan glaciers with some idealisations (Sect. 2.2) and without any tuning of model parameters. The fit parameters in Eqs. (14)–(17) may be different for a different set of glaciers. The parameterisations may also change if a more detailed and calibrated model of the same glaciers is used. However, the protocol used here to obtain the parameterisation for linearresponse properties can be directly applied without any change for any set of glaciers and for any iceflow and/or massbalance model. While applying the linearresponse model to any other region, it may be useful to obtain the response properties of a few tens of representative glaciers using flowmodel simulations and check if any recalibration of the parameterisation as given in Eqs. (14)–(17) is necessary.
We performed a theoretical analysis of the response of mountain glaciers within a timeindependent scaling assumption. In addition, the step response of 703 steadystate synthetic Himalayan glaciers with realistic geometries and idealised massbalance profiles were simulated with three different models: a scaling model, a 2D SIA model, and a linearresponse model. The results obtained are as follows.

Analytical expressions for climate sensitivity and response time of glacier area and volume are derived within a timeindependent scaling assumption. These expressions are validated using results from the scaling model simulation of the ensemble of 703 glaciers.

The response of the glaciers simulated with the 2D SIA model reveals that the initial steady states and the transient states follow the volume–area scaling relation, with the bestfit scale factor reducing slowly with time.

For the ensemble of glaciers studied, the scaling model obtains relatively smaller climate sensitivities of glacier area and volume by a factor of about 1.9 and 2.6, respectively, compared to those obtained from the SIA model. This results in a low bias in the longterm changes predicted by the scaling model.

For the ensemble of glaciers studied, the scaling model underestimates volume (area) response time by a factor of ∼1.8 (2.4) compared to the corresponding SIA estimates.

For the scaling model, τ_{A}≈τ_{V}, and $\frac{\mathrm{\Delta}{V}_{\mathrm{\infty}}}{V}\approx \mathit{\gamma}\frac{\mathrm{\Delta}{A}_{\mathrm{\infty}}}{A}$. In contrast, for the SIA simulations, τ_{A}≈1.5τ_{V} and $\frac{\mathrm{\Delta}{V}_{\mathrm{\infty}}}{V}\approx \mathrm{1.5}\mathit{\gamma}\frac{\mathrm{\Delta}{A}_{\mathrm{\infty}}}{A}$.

The relatively larger ratio of the two response times in the SIA simulations, along with an initial slow change in the area, leads to curved V−A trajectories, a decreasing c, and a relatively larger longterm volume loss for the transient glaciers due to a corresponding massbalance feedback.

A linearresponse model based on the parameterisations of SIAderived response properties helps reduce the biases in the longterm glacier changes predicted by the scaling model for the idealised central Himalayan glaciers. The improved performance of this model is validated on an independent set of 164 glaciers in the western Himalaya.
Based on the theoretical arguments and numerical evidence presented here, it is possible that qualitatively similar biases may generally be present in the longterm glacier changes computed with scaling models. However, the actual magnitudes of such biases in scaling models may be different from those obtained here for a set of synthetic Himalayan glaciers with idealised massbalance profiles. Possible biases in scaling models may, in turn, lead to a low bias in the corresponding estimates of the longterm sealevel rise contribution from shrinking mountain glaciers. On a multidecadal scale, a faster response due to shorter response times in the scaling model can compensate for the effects of smaller climate sensitivities to some extent. However, the low biases in scalingmodelderived changes in glacier area and volume are likely to become apparent over longer timescales of multiple centuries. The linearresponse model presented above could potentially be useful in predicting the longterm global glacier change and sealevel rise due to its accuracy and numerical efficiency.
The glacier model codes are available in the repository: https://github.com/DishaPatil/glacier_models (last access: 28 September 2020) (also available at https://doi.org/10.5281/zenodo.4050370, Patil, 2020).
The supplement related to this article is available online at: https://doi.org/10.5194/tc1432352020supplement.
AB designed the study, did the theoretical analysis, and wrote the paper. AJ and DP wrote the codes. AJ, DP, and AB ran the simulations. All three authors contributed to the analysis of the simulated data and discussions.
The authors declare that they have no conflict of interest.
The authors acknowledge valuable inputs from reviewer Eviatar Bach, the anonymous reviewer, and editor Valentina Radić. Deepak Suryavanshi has contributed to the initial development of the SIA code.
This research has been supported by the Ministry of Earth Sciences, Government of India (grant nos. MoES/PAMC/H&C/80/2016PCII and MoES/PAMC/H&C/79/2016PCII).
This paper was edited by Valentina Radic and reviewed by Eviatar Bach and one anonymous referee.
Adhikari, S. and Marshall, S. J.: Glacier volumearea relation for highorder mechanics and transient glacier states, Geophys. Res. Lett., 39, L16505, https://doi.org/10.1029/2012GL052712, 2012. a
Bach, E., Radić, V., and Schoof, C.: How sensitive are mountain glaciers to climate change? Insights from a block model, J. Glaciol., 64, 247–258, https://doi.org/10.1017/jog.2018.15, 2018. a, b, c
Bahr, D. B.: Width and length scaling of glaciers, J. Glaciol., 43, 557–562, https://doi.org/10.3189/S0022143000035164, 1997. a, b
Bahr, D. B., Meier, M. F., and Peckham, S. D.: The physical basis of glacier volumearea scaling, J. Geophys. Res.Sol. Ea., 102, 20355–20362, https://doi.org/10.1029/97JB01696, 1997. a
Bahr, D. B., Pfeffer, W. T., and Kaser, G.: A review of volumearea scaling of glaciers, Rev. Geophys., 53, 95–140, https://doi.org/10.1002/2014RG000470, 2015. a, b, c, d, e, f, g, h, i, j, k, l
Banerjee, A.: Brief communication: Thinning of debriscovered and debrisfree glaciers in a warming climate, The Cryosphere, 11, 133–138, https://doi.org/10.5194/tc111332017, 2017. a
Banerjee, A.: Volumearea scaling for debriscovered glaciers, J. Glaciol., https://doi.org/10.1017/jog.2020.69, in press, 2020. a
Banerjee, A. and Kumari, R.: Glacier area and the variability of glacier change, eartharxiv[preprint], https://doi.org/10.31223/osf.io/y2vs6, 2019. a
Chen, J. and Ohmura, A.: Estimation of Alpine glacier water resources and their change since the 1870s, IAHSAISH P, 193, 127–135, 1990. a
Clarke, G. K., Jarosch, A. H., Anslow, F. S., Radić, V., and Menounos, B.: Projected deglaciation of western Canada in the twentyfirst century, Nat. Geosci., 8, 372–377, https://doi.org/10.1038/NGEO2407, 2015. a
Egholm, D. L., Knudsen, M. F., Clark, C. D., and Lesemann, J. E.: Modeling the flow of glaciers in steep terrains: The integrated secondorder shallow ice approximation (iSOSIA), J. Geophys. Res.Earth, 116, F02012, https://doi.org/10.1029/2010JF001900, 2011. a
Farinotti, D. and Huss, M.: An upperbound estimate for the accuracy of glacier volume–area scaling, The Cryosphere, 7, 1707–1720, https://doi.org/10.5194/tc717072013, 2013. a, b, c, d
Farinotti, D., Brinkerhoff, D. J., Clarke, G. K. C., Fürst, J. J., Frey, H., Gantayat, P., GilletChaulet, F., Girard, C., Huss, M., Leclercq, P. W., Linsbauer, A., Machguth, H., Martin, C., Maussion, F., Morlighem, M., Mosbeux, C., Pandit, A., Portmann, A., Rabatel, A., Ramsankaran, R., Reerink, T. J., Sanchez, O., Stentoft, P. A., Singh Kumari, S., van Pelt, W. J. J., Anderson, B., Benham, T., Binder, D., Dowdeswell, J. A., Fischer, A., Helfricht, K., Kutuzov, S., Lavrentiev, I., McNabb, R., Gudmundsson, G. H., Li, H., and Andreassen, L. M.: How accurate are estimates of glacier ice thickness? Results from ITMIX, the Ice Thickness Models Intercomparison eXperiment, The Cryosphere, 11, 949–970, https://doi.org/10.5194/tc119492017, 2017. a
Glen, J. W.: The Creep of Polycrystalline Ice, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 228, 519–538, https://doi.org/10.1098/rspa.1955.0066, 1955. a
Grinsted, A.: An estimate of global glacier volume, The Cryosphere, 7, 141–151, https://doi.org/10.5194/tc71412013, 2013. a
Harrison, W. D., Elsberg, D. H., Echelmeyer, K. A., and Krimmel, R. M.: On the characterization of glacier response by a single timescale, J. Glaciol., 47, 659–664, https://doi.org/10.3189/172756501781831837, 2001. a, b, c, d, e, f
Hindmarsh, R. C. and Payne, A. J.: Timestep limits for stable solutions of the icesheet equation, Ann. Glaciol., 23, 74–85, https://doi.org/10.3189/S0260305500013288, 1996. a
Hock, R., Bliss, A., Marzeion, B., Giesen, R. H., Hirabayashi, Y., Huss, M., Radić, V., and Slangen, A. B.: GlacierMIP–A model intercomparison of globalscale glacier massbalance models and projections, J. Glaciol., 65, 453–467, https://doi.org/10.1017/jog.2019.22, 2019. a, b, c, d
Huss M. and Hock, R.: A new model for global glacier change and sealevel rise, Front. Earth Sci., 3, 54, https://doi.org/10.3389/feart.2015.00054, 2015. a, b
Huss, M. and Hock, R.: Globalscale hydrological response to future glacier mass loss, Nat. Clim. Change., 8, 135–140, https://doi.org/10.1038/s415580170049x, 2018. a
Huss, M., Jouvet, G., Farinotti, D., and Bauder, A.: Future highmountain hydrology: a new parameterization of glacier retreat, Hydrol. Earth Syst. Sci., 14, 815–829, https://doi.org/10.5194/hess148152010, 2010. a
Hutter, K.: Theoretical Glaciology, D. Reidel Publ. Co., Dordrecht, the Netherlands, 1983. a, b, c
Immerzeel, W. W., Lutz, A. F., Andrade, M., Bahl, A., Biemans, H., Bolch, T., Hyde, S., Brumby, S., Davies, B. J., Elmore, A. C., Emmer, A., Feng, M., Fernández, A., Haritashya, U., Kargel, J. S., Koppes, M., Kraaijenbrink, P. D. A., Kulkarni, A. V., Mayewski, P. A., Nepal, S., Pacheco, P., Painter, T. H., Pellicciotti, F., Rajaram, H., Rupper, S., Sinisalo, A., Shrestha, A. B., Viviroli, D., Wada, Y., Xiao, C., Yao, T., and Baillie, J. E. M.: Importance and vulnerability of the world's water towers, Nature, 577, 364–369, https://doi.org/10.1038/s415860191822y, 2020. a
Jarosch, A. H., Schoof, C. G., and Anslow, F. S.: Restoring mass conservation to shallow ice flow models over complex terrain, The Cryosphere, 7, 229–240, https://doi.org/10.5194/tc72292013, 2013. a
Jóhannesson, T., Raymond, C., and Waddington, E. D.: Timescale for adjustment of glaciers to changes in mass balance, J. Glaciol., 35, 355–369, https://doi.org/10.3189/S002214300000928X, 1989.
Kulp, S. and Strauss, B. H.: Global DEM errors underpredict coastal vulnerability to sea level rise and flooding, Front. Earth Sci., 4, 36, https://doi.org/10.3389/feart.2016.00036, 2016. a
Kumar, P., Saharwardi, M. S., Banerjee, A., Azam, M. F., Dubey, A. K., and Murtugudde, R.: Snowfall Variability Dictates Glacier Mass Balance Variability in HimalayaKarakoram, Sci. Rep.UK, 9, 18192, https://doi.org/10.1038/s41598019545539, 2019. a
Kraaijenbrink, P. D. A., Bierkens, M. F. P., Lutz, A. F., and Immerzeel, W. W.: Impact of a global temperature rise of 1.5 degrees Celsius on Asia's glaciers, Nature, 549, 257–260, https://doi.org/10.1038/nature23878, 2017. a, b
Laha, S., Kumari, R., Singh, S., Mishra, A., Sharma, T., Banerjee, A., Nainwal, H. C., and Shankar, R.: Evaluating the contribution of avalanching to the mass balance of Himalayan glaciers, Ann. Glaciol., 58, 110–118, https://doi.org/10.1017/aog.2017.27, 2017. a
Le Meur, E., Gagliardini, O., Zwinger, T., and Ruokolainen, J.: Glacier flow modelling: a comparison of the Shallow Ice Approximation and the fullStokes solution, C. R. Phys., 5, 709–722, https://doi.org/10.1016/j.crhy.2004.10.001, 2004. a, b
Leysinger Vieli, G. J.M. C. and Gudmundsson, G. H.: On estimating length fluctuations of glaciers caused by changes in climatic forcing, J. Geophys. Res., 109, F01007, https://doi.org/10.1029/2003JF000027, 2004. a, b, c, d, e
Lüthi, M. P.: Transient response of idealized glaciers to climate variations, J. Glaciol., 55, 918–930, https://doi.org/10.3189/002214309790152519, 2009. a, b, c
Marzeion, B., Hock, R., Anderson, B., Bliss, A., Champollion, N., Fujita, K., Huss, M., Immerzeel, W., Kraaijenbrink, P., Malles, J.H., and Maussion, F.: Partitioning the Uncertainty of Ensemble Projections of Global Glacier Mass Change, Earth's Future, 8, e2019EF001470, https://doi.org/10.1029/2019EF001470, 2020. a, b, c
Maussion, F., Butenko, A., Champollion, N., Dusch, M., Eis, J., Fourteau, K., Gregor, P., Jarosch, A. H., Landmann, J., Oesterle, F., Recinos, B., Rothenpieler, T., Vlug, A., Wild, C. T., and Marzeion, B.: The Open Global Glacier Model (OGGM) v1.1, Geosci. Model Dev., 12, 909–931, https://doi.org/10.5194/gmd129092019, 2019. a
NASA/METI/AIST/Japan Spacesystems and U.S./Japan ASTER Science Team: ASTER Global Digital Elevation Model V003 [Data set], NASA EOSDIS Land Processes DAAC, https://doi.org/10.5067/ASTER/ASTGTM.003, 2019. a
Oerlemans, J.: Glaciers and climate change, A. A. Balkema Publishers, Lisse, 2001. a, b, c, d, e, f, g, h, i, j
Patil, D.: DishaPatil/glacier_models: Glacier Models code released for publication (Version v1.0.0), Zenodo, https://doi.org/10.5281/zenodo.4050370, 2020. a
Radić, V., Hock, R., and Oerlemans, J.: Volumearea scaling vs flowline modelling in glacier volume projections, Ann. Glaciol., 46, 234–240, https://doi.org/10.3189/172756407782871288, 2007. a, b, c, d, e, f, g, h
Radić, V., Hock, R., and Oerlemans, J.: Analysis of scaling methods in deriving future volume evolutions of valley glaciers, J. Glaciol., 54, 601–612, https://doi.org/10.3189/002214308786570809, 2008. a, b, c
Radić, V., Bliss, A., Beedlow, A. C., Hock, R., Miles, E., and Cogley, J. G.: Regional and global projections of twentyfirst century glacier mass changes in response to climate scenarios from global climate models, Clim. Dynam., 42, 37–58, https://doi.org/10.1007/s0038201317197, 2014. a, b
Raper, S. C. and Braithwaite, R. J.: Low sea level rise projections from mountain glaciers and icecaps under global warming, Nature, 439, 311–313, https://doi.org/10.1038/nature04448, 2006. a
RGI Consortium: Randolph Glacier Inventory – A Dataset of Global Glacier Outlines: Version 6.0: Technical Report, Global Land Ice Measurements from Space, Colorado, USA, Digital Media, https://doi.org/10.7265/N5RGI60, 2017. a
Rounce, D. R., Khurana, T., Short, M. B., Hock, R., Shean, D. E., and Brinkerhoff, D. J.: Quantifying parameter uncertainty in a largescale glacier evolution model using Bayesian inference: application to High Mountain Asia, J. Glaciol., 66, 175–187, https://doi.org/10.1017/jog.2019.91, 2020. a
Wagnon, P., Vincent, C., Arnaud, Y., Berthier, E., Vuillermoz, E., Gruber, S., Ménégoz, M., Gilbert, A., Dumont, M., Shea, J. M., Stumm, D., and Pokhrel, B. K.: Seasonal and annual mass balances of Mera and Pokalde glaciers (Nepal Himalaya) since 2007, The Cryosphere, 7, 1769–1786, https://doi.org/10.5194/tc717692013, 2013. a
Zekollari, H. and Huybrechts, P.: On the climategeometry imbalance, response time and volumearea scaling of an alpine glacier: insights from a 3D flow model applied to Vadret da Morteratsch, Switzerland, Ann. Glaciol., 56, 51–62, https://doi.org/10.3189/2015AoG70A921, 2015. a
Zekollari, H., Huss, M., and Farinotti, D.: Modelling the future evolution of glaciers in the European Alps under the EUROCORDEX RCM ensemble, The Cryosphere, 13, 1125–1146, https://doi.org/10.5194/tc1311252019, 2019. a
Zhang, Y., Hirabayashi, Y., Liu, Q., and Liu, S.: Glacier runoff and its impact in a highly glacierized catchment in the southeastern Tibetan Plateau: past and future trends, J. Glaciol., 61, 713–730, https://doi.org/10.3189/2015JoG14J188, 2015. a
Zhao, L., Ding, R., and Moore, J. C.: Glacier volume and area change by 2050 in high mountain Asia, Global Planet. Change, 122, 197–207, https://doi.org/10.1016/j.gloplacha.2014.08.006, 2014. a