Relative e ff ect of slope and equilibrium line altitude on the retreat of Himalayan glaciers

Introduction Conclusions References


Introduction
The Hind-Kush Himalayan region possesses one of the largest concentrations of mountain glaciers and melt water from these glaciers is an important source for many rivers in the region.Numerous investigations have been carried out to understand changes in glaciers in the Himalayas (Kulkarni et al., 2007;Kulkarni et al., 2011;Bolch et al. , 2010;Yong et al., 2010).These investigations show that a majority of glaciers in the Himalayas are retreating.However recent investigations in Karakoram mountain range indicate that some glaciers are advancing (Scherler et al., 2011b;Hewitt, 2005).In addition, rate of retreat is different for individual glaciers depending upon numerous geomorphological parameters like area-altitude distribution, length, size, slope, debris cover etc. (Deota et al., 2011;Kulkarni et al., 2005).Field observations carried out in a large glacier like Siachen have shown that the glacier is almost stationary or has shown little retreat since 1995.This has lead to the erroneous conclusion that glaciers in the North-West Himalayas are not affected by global warming (Ganjoo and Kaul, 2009).
In our opinion, the different rates of retreat/advance of glaciers within a region, over which the climatic conditions do not change significantly, is due to the important role played by the dynamics of ice movement, which in turn is controlled by the mean slope and length of the glacier.In this paper, we provide an explanation for this apparent contradiction of advancing glaciers in a global warming scenario, by using a simple model to understand relative importance of slope, length, and Equilibrium Line Altitude (ELA).ELA is considered as a good indicator of glacier mass balance (Benn and Lehmkul, 2000).The model shows that rate of retreat of glaciers can be different even if environmental changes are similar in a given region.Locations of the Himalayan glaciers considered are shown in Fig. 1.This is done without explicitly using the concepts of inertia and response time of a glacier.The approach is different from the ones followed in previous studies such as Pelto and Hedlund (2001).
In what follows, mean slope s is defined as follows where h max is the altitude at the top of the glacier, h min is the altitude at the snout and L is the length of the glacier.Schematic view is shown in Fig. 2.
Published by Copernicus Publications on behalf of the European Geosciences Union.

Motivation and hypothesis
The role of the mean-slope in determining equilibrium lengths of glaciers is well known.Using simple arguments, one can derive an expression for the equilibrium length, in terms of the slope, equilibrium line altitude and mean thickness (Oerlemans, 2008).Change in any of these parameters would result in advance/retreat.Empirical evidence of the role of slope is also available.In Fig. 3, we show the variation of retreat rates of glaciers in the Parbati and Baspa basins with mean slope.The retreat rates were derived over a 11 year time period from satellite images.One can see that there is a trend of decreasing retreat rate with increasing slope.This set had only retreating glaciers and can be assumed to have nearly the same change in environmental condition (like change in ELA).In addition, it can be seen that the variability decreases with increasing slope, suggesting that the slope has a major role to play.
Using the more extensive satellite data of Scherler et al. (2011), the distribution of retreat rates for low slopes (s < 0.15) and high slopes (s > 0.25) are plotted in Fig. 14a.To quantify the differences in the probability distribution, the Kolmogorov-Smirnov test was performed.The value of the statistic obtained was 0.2857, indicating that the distributions are similar, but distinct.For identical distributions, the K-S statistic is zero.A limitation of this dataset is that it is available for eight years, in which inter-annual and shortterm variation would also be seen.For low slopes, the maxima is at −10 m yr −1 while for large slopes the maxima is a 0 m yr −1 with more advancing glaciers.This suggests that slope in addition to the climate sensitivity term 1/s, there is another part proportional to s contributing to the advance.
The slope, in addition to determining the sensitivity to changes in ELA can also influence the advance due to dynamics.A commonly used equation for glacier simulations (Oerlemans, 1988;Adhikari and Huybrechts, 2009) is of the following form.
where, H (x) represents the thickness at a point x along the flow-line, U (x) the mean velocity of ice along the flow-line and Ḃl is the mass balance with dimensions [L]/[T ].
The mass balance term (with underbrace) is a function of the altitude and climate forcing and is prescribed as a function of x and t.The other term (underlined) depends on the velocity U which is a function of the shear-stress, which depends mainly on the bottom and surface slopes and parameters which represent basal slip and sliding.The first term represents the thermodynamics (snowfall, melting) and the second the dynamics (gravity, slope effects).The net movement (advance/retreat) of the glacier front is due to the integrated effect of both the processes.Given that gravity is always present, one can consider the following scenarios 1.The mass balance term is significant and unchanging: The system would evolve to a steady state, wherein the accumulation, ablation and flow are in balance.This is the usual equilibrium scenario.
2. The mass balance term is negligible (this is a hypothetical situation): In this case, the ice would flow down the slope due to gravity.This can also be thought of as the limiting case where the thermodynamic processes are unimportant and only ice dynamics, controlled by the slope is the relevant process.
In the first case changes to the glacier length are possible only if the external conditions (mass balance) change with time.
In the second case, the glacier would continue to advance, with the thickness reducing (since volume is conserved).Our hypothesis is that the advance/retreat of a glacier can be understood in terms of theses two tendencies and that they can be added linearly.The first term is the tendency of the glacier front to advance due to gravity, which includes both ice-deformation and sliding.This is the new idea proposed in this paper.The second term, is equivalent to the climate sensitivity of glacier length to ambient air temperature as derived by Oerlemans (2008), multiplied by the rate of change of temperature with time.
The effect of the two processes is assumed to be linear and the net advance/retreat is the sum of the two.These terms can be further expressed as and F 1 is a function of L,H and s and is obtained from ice-flow simulations.dL/dh e is estimated from equilibrium simulations.dh e /dt is related to the environmental change and is expected to be related to the rate of increase of mean global temperature.
Given the lengths, slopes and retreat rates for a set of real glaciers, one can then find a least-squares, best-fit estimate of α and dh e /dt.These values can then be used to predict the retreat/advance for other glaciers.A schematic view of the modelling approach is shown in Fig. 5.One can also use it to understand the contribution of different terms to the advance/retreat.This way of splitting the change in glacier length is different from the manner it is done in simple models (Cuffey and Paterson, 2010), for the accumulation zone, ablation zone or at the glacier terminus.In those models, the mass balance component and ice-flow term are both present, leading to a linear system, which responds to an abrupt change in environmental changes.In our model, the thermodynamics and dynamics processes are split.To the best of our knowledge this particular way of decomposing the change in length of the glacier has not been done before.
We are aware that the model proposed above is highly simplified and therefore has many limitations.A number of factors such as debris cover, orientation of the glacier, area-altitude distribution, variation of slope, number of snow avalanches (Hewitt, 2005) play a role in the glacier dynamics and all these factors are not explicitly included.The attempt is to show that, even with these limitations, some aspects of glacier dynamics can be explained.Due to the coupling of the slope and other factors, established by studies such as those by Scherler et al. (2011a), to a first order, the role of these factors could be indirectly represented by the slope.
The climatic conditions would also vary across the glaciers considered.The objective is to show that if ELA changes are similar, glacier retreat will be different, depending upon geomorphological parameters.Our premise is that on a climate time-scale, to the first order, the ELA changes resulting from global climate change are of similar magnitude.

Numerical model
For simulating glaciers, numerical models of varying degrees of sophistication are possible (Kotlarski et al., 2010).Among them, Adhikari and Huybrechts (2009), used a simple model to simulate the variation of glacier AX010, and study scenarios for its future evolution.This simple model, based on a formulation due to Oerlemans (1988), seems to be quite effective in simulating observed data over the past fifty years.We have developed a FORTRAN code based on the same formulation and that code has been used for all the simulations presented here.
The model consists of a partial differential equation for the variation of mean thickness H along the flow-line (x).
where, w b is the width at the bottom of the glacier, U the mean velocity of ice along the flow-line, λ the side-slope and Ḃl is the mass balance with dimensions The velocity of ice U is split in to two components, the sliding velocity U s and the deformation velocity U d , which are modelled as follows where τ is the shear stress and f s ,f d are parameters to be derived from measurements and then tuned numerically, g is the acceleration due to gravity and h b is the elevation at the bottom/base of the glacier.
The prognostic equation for H can be recast in the following form. where The mass balance, which is a function of x and t is modelled as follows where h EL is the equilibrium line altitude, β is estimated from historical data and B l−hist (t) has to be specified either from observations or using some proxy data.
Values of the parameters used were: Adhikari and Huybrechts (2009) mention that f s and f d were multiplied by a factor (1/γ ) to match the results with observations.The values of γ used by them range from 2.5 to 30.We find that a factor of 3.25 gives the best match.
The code was validated by simulating the glacier AX010 with the mass balance data present in appendices of Adhikari (2007) and the results were found to be nearly identical to Adhikari and Huybrechts (2009).While Adhikari and Huybrechts (2009) have used their numerical model to simulate the historical variation of particular glaciers and project the future scenarios, we have used the model in a different way.We use it primarily to simulate idealized glacier flow at two extreme conditions: (a) with zero mass balance and (b) equilibrium conditions.

Impact of slope
For medium to large glaciers (length >1.5 km), the meanslope and length are expected to play a major role, therefore we ignore variations of the bed topography and perform idealized simulations with the base varying linearly and consider different slopes.
The first set of ice-flow simulations were performed with zero mass balance.This was done to simulate gravity effects on a mass of ice and quantify the part of motion of glaciers which is due to just flow down the incline, in the absence of snowfall/melting.Such a flow does not generally occur in nature, since some mass balance is usually present.However, these simulations provide an indication of the tendency of a mass of ice to flow and we use the initial trend (after one year) to estimate F 1 as a function of L, H and s.
In the second set, mass balance varying as a linear function of altitude was imposed and simulations performed varying the equilibrium line altitude, keeping the origin of the glacier constant.An assumption here is that for glaciers considered, it is only difference h max − h ELA which matters.

Ice-flow simulations
Results of ice-flow simulations, with no snowfall or melting, which were performed to estimate F 1 as a function of L, H and s, are presented in this section.The simulation was started with a block of ice of length L 0 and uniform thickness H 0 and with the mass balance term set to zero.Schematic view of the process is shown in Fig. 4. L 0 was varied from 2 to 6 km. and H 0 varied from 50 to 250 m.The model was integrated for upto 2000 years to understand the qualitative behavior.However, only the velocity at initial stages, i.e. after one year was used for further calculations.It should be noted that the functional form of F 1 , i.e. the dependence on L, H and s is not affected by the time at which the front velocity is chosen.Only the constant of multiplication changes.Initially, the thickness has a top-hat profile and soon diffusive processes make the distribution smooth, with a maxima towards the lower end.Gravitational force causes the ice mass to stretch and flow down the incline.Initially, length increases sharply with time and later there is nearly linear growth.The rate of change increases with the slope.The maximum thickness decreases with time and falls more rapidly with increasing slope.The effect of the mean slope on the dynamics is to increase the average ice-velocity.The velocity of the front increases with slope and decreases gradually with time.
For application to real glaciers, the value of dL/dt after one year, given the length and thickness at the end of the previous year is used.On varying L 0 , the velocity was found to increase linearly The dependence on H 0 was found to follow a three-fourths power law.The variation with s was found to be linear (Figs. 6 and 7).The following expression was found to be a good fit for the simulated values of dL/dt.Using the above and additional simulations for various values of H 0 and L 0 , the following formula for ice-front velocity, i.e. the advance/retreat rate of the block, was determined where L is in kilometres, H is in metres and dL/dt in metres/year.

Equilibrium mass balance simulations
Results of simulations with mass balance are presented in this section.Objective was to estimate dL/dh e and average thickness H as functions of L and s.The simulations were started with zero ice and integrated with a specified gradient balance β.Once steady state was reached, the equilibrium length and thickness were determined.A schematic view of the process is shown in Fig. 8. Observed equilibrium line altitudes, lengths and mean slopes for a few Himalayan glaciers are listed in Table 1.Therefore simulations were performed for slopes varying from 0.075 to 0.2 and h max − h ELA from 0 to 1600.The value of β used by Adhikari and Huybrechts is 0.01.For Chhota Shigri, the value was estimated to be 0.009.Since these values are quite close, simulations were performed with β = 0.009.For convenience, we define With linear variation of mass balance with altitude (gradient β = 0.009) and zero net mass balance, equilibrium shapes on varying base slope and h e are as follows.
Equilibrium values of length increase with h e (Fig. 9).Also, one can see that for the same value of h e , length decreases with increasing slope.Oerlemans (2008) derived the following estimate for equilibrium length.
where H m is the average thickness.For large L the term 2H m /(sL) is small since it varies as 1/ √ L and the slope is approximately 2/s.On comparing the two expressions and values from the graph, for L > 1 km, the following expression serves as a good approximation The variation of equilibrium value of average thickness H is similar to that of L. It increases with h e and for the same value of h e , it decreases with increasing slope.A curve fit for the average thickness, as a function of L and s is as follows where L and H are expressed in metres.This fit is similar in form to the one used by Oerlemans (2008).
To check the usefulness of the curve fit, it is compared with observed values of thickness and length for a large number of glaciers from the World Glacier Inventory (WGI) database (http://nsidc.org/data/g01130.html).One can see from Fig. 10, that the fit is within the range of observed values.

Application to real glaciers
As per our hypothesis, This simplifies to where L is expressed in km and H is determined from Eq. ( 13).
Given the lengths, slopes and retreat rates of real glaciers, one can find a least-squares, best-fit estimate of α and dh e /dt.For a few glaciers, the control set, these quantities have been listed in Table 2. Locations of these glaciers are given in Fig. 1.Retreat rates considered have been for a period of 25 years.The retreat rate for Khumbu is taken from Rai et al. (2005).
The best fit values are: α = 0.04053 and dh e /dt = −0.6659.Using the values of α and dh e /dt, the computed and observed values of retreat for the fitted set (glaciers listed in Table 2) is shown in Table 3 and the predicted set (AX010, Zemu and Gangotri) in Table 4. Observed retreat rates are from Basnett et al. (2011) for Zemu and Kumar et al. (2008) for Gangotri.

Discussion
As one can see from Tables 3 and 4, the computed values of retreat are close to what is observed.RMS error for the first set is 1.61 m yr −1 and the second set is 3.82 m yr −1 .One should note that these values are comparable to errors of measurement using field data.However, these differences are not significant, since our main aim is to explain the balance between the two processes determining advance/retreat and not to match the exact value for any particular glacier.
While the role of debris cover is not explicitly modelled, since the set of glaciers used for estimating the constants include glaciers with different amounts of debris cover, its effect is implicitly present.
For most of the glaciers, the dL/dt dyn term is small and their behaviour is dominated by the climate term, leading to a net retreat.Relative roles played by the length and slope are brought out in Fig. 11.In the first bar-chart, the glaciers are arranged in increasing order of slope.The retreat, is inversely proportional to slope.The tendency to advance depends on both the slope and length.
Although the lengths of Zemu and Gangotri glaciers are similar the rate of retreat is quite different.Zemu glacier is almost stationary while Gangotri glacier is retreating at the rate of 19 m per year.The proposed model suggests that the large difference in rate of retreat between these glaciers is on account of the difference in slope.Slope of Zemu glacier is almost double that of Gangotri glacier.The higher slope of Zemu glacier causes the advance due to gravity to be comparable to the ablation term leading to an almost zero rate of advance/retreat.The retreat of Gangotri glacier is more sensitive to changes in mass balance because the slope of the glacier is much smaller.For the AX010 glacier, even though the slope is high, its short length, causes the advance term to be negligible and the climate term to dominate, leading to a net retreat.
Net retreat for both the sets is compared in Fig. 12, indicating a reasonably good match between computed and observed values.

Application to larger datasets: Parbati and other basins
In the previous section, our model was applied to a small set of glaciers for which long term records of retreat are available.Application for a larger set is preferable to validate the model.However for larger sets which are available from satellites, the time period is less, which would not be appropriate for our model.Over a time period of 8 to 10 years, the inter-annual variations of snowfall/melting would dominate the observed retreat and the climate trends are not very clear.This is evident in the wide range of retreats from satellite data (maximum around 60 m yr −1 ) than those from long term onsite observations (maximum around 20 m yr −1 ).This has to be kept in mind while comparing these results.First, we show application to a larger set of 38 glaciers in the Parbati basin, for which retreat rates over a 11 year period are available.A subset of 15 glaciers was used to determine the coefficients in the model.The best fit values α and dh e /dt were found to be 0.05 and −1.64, respectively.These coefficients were used to compute retreat rates for the complete set of 38 glaciers.The scatter plot of observed versus computed retreat rates is shown in Fig. 13.One can see that the comparison is good for low values of retreat and reasonable for glaciers with high retreat (above 40 m yr −1 ).There is scatter, is likely to be due to the shorter duration of data.The scatter could also be due to local variations in ELA change and boundary effects, which are not accounted for in our empirical model.
Now we apply the model to the set of glaciers of Scherler et al. (2011b).Since retreat rates are over a 8 year period, we do not compare for individual glaciers, but only the distribution.Distribution of retreat rates for a simulated set of glaciers, with length distribution similar to that of Scherler etal and slopes centered around two values: 0.14 (low slope) and 0.28 (high slope) is shown in Fig. 14b.One can see that qualitatively, the distributions are similar to that for observed data (Fig. 14a).In particular, the shift of 10 m yr −1 in the peak from low to high slopes is captured.

Conclusions
Using simulations with a numerical ice-flow model and simple hypotheses, we have demonstrated the relative effect of slope (dynamics) and equilibrium line altitude (thermodynamics) on the retreat of Himalayan glaciers.We have shown that dynamics, as determined by the length and mean-slope can explain major differences in the behavior of glaciers when subject to similar environmental changes.Decomposition of the glacier front velocity in terms of slope and ELA is a novel approach and as far as we know, quantification in these terms has not been done before.The drastically different responses of Gangotri and Zemu, glaciers of nearly the same length, is explained well by this model.In the case of Zemu, advance due to slope is around 11.7 m yr −1 which is balanced by retreat due to climate change, while for Gangotri, the climate term dominates.Therefore, using only the observed retreat as an indicator of climate change, leads to erroneous conclusions.
The model has also been applied to a larger set of glaciers in the Parbati basin, and other regions.For these glaciers, with retreat data over a shorter time-period, the distributions are well simulated.While we have concentrated on mainly Himalayan glaciers, the concept is quite general and hence should be applicable to glaciers in other regions of the world.

Fig. 1 .Fig. 2 .
Fig. 1.Map showing the positions of the glaciers considered for the study.International borders are shown.

Fig. 3 .
Fig. 3. Observed variation of retreat rates with slope, for glaciers in the Parbati and Baspa basins.The solid circle represents the mean value and bars, the standard deviation.The set consists of 57 glaciers.

Fig. 6 .
Fig. 6.Variation of velocity of the front as a function of slope for two different values of L 0 .

Fig. 7 .
Fig. 7. Top panel: Variation of velocity of the front as a function of thickness H 0 with slope and L 0 kept constant.Bottom panel: Variation of velocity of the front as a function of thickness L 0 with slope and H 0 kept constant.

Fig. 9 .
Fig. 9. Variation of equilibrium length as function of h e and s.

Fig. 10 .
Fig. 10.Observed variation of mean thickness (filled squares), with error bars as function of L from WGI database (2009).Solid line represents the curve fit for a value of s of 0.15.

Fig. 11 .
Fig. 11.Bar chart showing the balance of the opposing tendencies of advance due to dynamics (red boxes), retreat due to thermodynamics (blue boxes) and the net movement (black bars).The glaciers are arranged in increasing order of slope.Figures on the top are for the fitted set and those on the bottom are for the predicted set.

Fig. 12 .
Fig. 12. Scatter-plot of observed retreat versus predicted retreat for the fitted set (top) and the predicted set (bottom).

Fig. 13 .
Fig. 13.Scatter-plot of observed retreat versus predicted retreat for the Parbati basin.The thinner lines represent the band of uncertainty (10 m yr −1 ).

Fig. 14 .
Fig. 14.Distribution of retreat rates.Top: from observed data of Scherler et al. (2011) and Bottom: for a simulated set of glaciers with low slope (0.14) and high slope (0.28).The bin size is 10.

Table 1 .
Observed values of slope, h e = h max − h ELA and length for a few Himalayan glaciers.

Table 3 .
Computed and observed retreats for the fitted set.RMS error is 1.61 m yr −1 .

Table 4 .
Computed and observed retreats for the predicted set.RMS error is 3.82 m yr −1 .