Skip to main content Accessibility help


  • Access
  • Cited by 17
  • Cited by
    This article has been cited by the following publications. This list is generated based on data provided by CrossRef.

    ADHIKARI, Surendra J. MARSHALL, Shawn and HUYBRECHTS, Philippe 2011. On Characteristic Timescales of Glacier AX010 in the Nepalese Himalaya. Bulletin of Glaciological Research, Vol. 29, Issue. , p. 19.

    Adhikari, Surendra and Marshall, Shawn J. 2011. Improvements to shear-deformational models of glacier dynamics through a longitudinal stress factor. Journal of Glaciology, Vol. 57, Issue. 206, p. 1003.

    Leclercq, Paul Willem and Oerlemans, Johannes 2012. Global and hemispheric temperature reconstruction from glacier length fluctuations. Climate Dynamics, Vol. 38, Issue. 5-6, p. 1065.

    Adhikari, S. and Marshall, S. J. 2013. Influence of high-order mechanics on simulation of glacier response to climate change: insights from Haig Glacier, Canadian Rocky Mountains. The Cryosphere, Vol. 7, Issue. 5, p. 1527.

    Banerjee, Argha and Shankar, R. 2013. On the response of Himalayan glaciers to climate change. Journal of Glaciology, Vol. 59, Issue. 215, p. 480.

    Zhang, Tong Xiao, Cunde Colgan, William Qin, Xiang Du, Wentao Sun, Weijun Liu, Yushuo and Ding, Minghu 2013. Observed and modelled ice temperature and velocity along the main flowline of East Rongbuk Glacier, Qomolangma (Mount Everest), Himalaya. Journal of Glaciology, Vol. 59, Issue. 215, p. 438.

    KATTEL, DAMBARU BALLAB and YAO, TANDONG 2013. Recent temperature trends at mountain stations on the southern slope of the central Himalayas. Journal of Earth System Science, Vol. 122, Issue. 1, p. 215.

    Singh, Vijay P. Mishra, Ashok K. Chowdhary, H. and Khedun, C. Prakash 2014. Modern Water Resources Engineering. p. 525.

    Vincent, C. Harter, M. Gilbert, A. Berthier, E. and SIX, D. 2014. Future fluctuations of Mer de Glace, French Alps, assessed using a parameterized model calibrated with past thickness changes. Annals of Glaciology, Vol. 55, Issue. 66, p. 15.

    Zhang, Tong Ju, Lili Leng, Wei Price, Stephen and Gunzburger, Max 2015. Thermomechanically coupled modelling for land-terminating glaciers: a comparison of two-dimensional, first-order and three-dimensional, full-Stokes approaches. Journal of Glaciology, Vol. 61, Issue. 228, p. 702.

    Nainwal, H.C. Banerjee, Argha Shankar, R. Semwal, Prabhat and Sharma, Tushar 2016. Shrinkage of Satopanth and Bhagirath Kharak Glaciers, India, from 1936 to 2013. Annals of Glaciology, Vol. 57, Issue. 71, p. 131.

    Banerjee, Argha and Azam, Mohd Farooq 2016. Temperature reconstruction from glacier length fluctuations in the Himalaya. Annals of Glaciology, Vol. 57, Issue. 71, p. 189.

    Gaddam, Vinay Kumar Kulkarni, Anil V and Gupta, Anil Kumar 2017. Reconstruction of specific mass balance for glaciers in Western Himalaya using seasonal sensitivity characteristic(s). Journal of Earth System Science, Vol. 126, Issue. 4,

    Gantayat, Prateek Kulkarni, Anil V. Srinivasan, J. and Schmeits, Maurice J 2017. Numerical modelling of past retreat and future evolution of Chhota Shigri glacier in Western Indian Himalaya. Annals of Glaciology, Vol. 58, Issue. 75pt2, p. 136.

    Laha, Sourav Kumari, Reshama Singh, Sunil Mishra, Aditya Sharma, Tushar Banerjee, Argha Nainwal, Harish Chandra and Shankar, R. 2017. Evaluating the contribution of avalanching to the mass balance of Himalayan glaciers. Annals of Glaciology, Vol. 58, Issue. 75pt2, p. 110.

    Rowan, Ann V. Quincey, Duncan J. Gibson, Morgan J. Glasser, Neil F. Westoby, Matthew J. Irvine-Fynn, Tristram D. L. Porter, Phillip R. and Hambrey, Michael J. 2018. The sustainability of water resources in High Mountain Asia in the context of recent and future glacier change. Geological Society, London, Special Publications, Vol. 462, Issue. 1, p. 189.

    Bolch, Tobias Shea, Joseph M. Liu, Shiyin Azam, Farooq M. Gao, Yang Gruber, Stephan Immerzeel, Walter W. Kulkarni, Anil Li, Huilin Tahir, Adnan A. Zhang, Guoqing and Zhang, Yinsheng 2019. The Hindu Kush Himalaya Assessment. p. 209.




      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Numerical modelling of historical front variations and the 21st-century evolution of glacier AX010, Nepal Himalaya
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Numerical modelling of historical front variations and the 21st-century evolution of glacier AX010, Nepal Himalaya
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Numerical modelling of historical front variations and the 21st-century evolution of glacier AX010, Nepal Himalaya
        Available formats
Export citation


Due to the lack of measurements of ice velocity, mass balance, glacier geometry and other baseline data, model-based studies of glacial systems in the Nepal Himalaya are very limited. Here a numerical ice-flow model has been developed for glacier AX010 in order to study its relation to local climate and investigate the possible causes of its general retreat since the end of the Little Ice Age. First, an attempt is made to simulate the historical front variations, considering each climatic parameter separately. Good agreement between the observations and model projections can be obtained under the assumption that variations in glacier front position are a response to changes in temperature alone. The same assumption is made about future changes to explore the 21st-century evolution of the glacier. Under a no-change scenario, the glacier will retreat by another ∽600m by AD 2100, whereas it is projected to vanish completely during this century for all trends with a temperature rise larger than +2.5˚C by AD 2100 with respect to the 1980–99 mean. With constant precipitation at the 1980–99 mean, the model predicts that the glacier will cease to exist at AD 2083, 2056 or 2049 if the temperature rises linearly by 3˚C, 4.5˚C or 6˚C respectively by the end of this century. With an additional range of precipitation changes between –30% and +30%, the life expectancy of glacier AX010 varies by 18, 6 and 2 years for the respective temperature rises. Thus the role of precipitation becomes minimal for the higher trends of temperature rise.

1. Introduction

Glaciers are one of the major components of the cryosphere. They are ancient rivers of compressed snow that creep through the landscape, shaping the planet’s surface. They tend to be located in wet and cold places. Wet implies a region not too far from a moisture source; cold implies high latitude and/or altitude (Oerlemans, 2001). Glaciers are considered excellent climate-change indicators (e.g. Hall and Fagre, 2003) because their dimensions change in response to several climatic variables characterizing the local environment (Kadota and Ageta, 1992). Currently, the Earth is at its warmest since record keeping began around AD 1850 (Brohan and others, 2006). The same is true of Nepal, with the largest temperature increases found at higher altitudes (Thomas and Rai, 2005). Warming leads to the shrinkage of glacial systems, so it is not surprising that Himalayan glaciers have continuously receded over recent decades.

Glacier AX010, Shorong Himal, (27˚42′ N, 86˚34′ E; Fig. 1) is one of more than 3000 glaciers in the Nepal Himalaya. Sufficient data on this glacier exist in terms of areal extent, mass balance and ice flow (e.g. Fujita and others, 2001) to attempt a model-based study to investigate its response to climate change. Since the first detailed studies conducted in 1978/79 (e.g. Ageta and others, 1980; Ikegami and Ageta, 1991; Kadota and others, 1997), several state variables of the glacier have been recorded in 1989, 1991, 1995–99 (e.g. Fujita and others, 2001) and 2004 (Shrestha and Shrestha, 2004). In 1978, intensive observations were conducted to determine its mass balance (Ageta and others, 1980), heat balance (Ohata and Higuchi, 1980), areal extent (WGMS, 1998) and surface velocity (Ikegami and Ageta, 1991). In 1989, only the terminus position was measured (Fujita and others, 2001). In 1991, a topographic mapping of almost the whole area of the glacier was carried out (WGMS, 1998). Annual monitoring of the glacier was initiated in 1995 and continued until 1999 by means of a ground survey, the stake method, in order to obtain the mass balance, surface velocity and areal-extent data (WGMS, 2005). The ice thickness was first recorded in June 1995 (Kadota and others, 1997). Radio-echo sounding was conducted at three different locations on the glacier. Its areal extent was also monitored in 2004 (Shrestha and Shrestha, 2004).

Fig. 1. Location of Shorong Himal and topographical map of glacier AX010 as surveyed in 1979 (Ikegami and Ageta, 1991, modified). Solid circles on the glacier are positions of mass-balance stakes. The dotted line along the glacier is the assumed flowline. The solid lines across the glacier are 20m interval contour lines. The grey spot labelled ‘pond’ next to the glacier terminus is a glacier-fed lake.

Besides these observational records, there have been a few model-based studies of glacier AX010. The mass-balance model by Ageta and Kadota (1992) showed that its mass balance is highly sensitive to summer temperature, not only because of summer melting but also because the glacier is in the monsoon climate regime. It therefore receives almost all of its precipitation in summer, with temperature determining the fraction accumulating as snow. Furthermore, Kayastha and others (1999), based on their energy-balance model, demonstrated the high sensitivity of the mass balance to global radiation and the surface albedo of the glacier. Kadota and others (1997) projected the surface profile for 2005 and 2015 for unchanged climatic conditions from 1995 onwards (summer temperature of 2.9˚C at the snout, and summer precipitation of 1.4 m). Shrestha and Shrestha (2004) reported that the glacier surface measured in summer 2004 was remarkably close to Kadota and others’ (1997) prediction for 2005.

In this study, the relation between glacier AX010 and local climate change is investigated under a series of reasonable assumptions. This provides an opportunity to understand how a small valley glacier in the Nepal Himalaya responds to climate change. In addition, we investigate which climatic parameter played the decisive role in driving the history of the glacier terminus fluctuations during the past two centuries.

2. The Numerical Ice-Flow Model

The numerical ice-flow model used in this study describes ice flow along a central flowline. The prognostic equation of this glacier model is the continuity equation, describing how the change in ice thickness is related to the flux divergence and the net surface mass balance. Along a flowline, conservation of mass can be approximated as (e.g. Oerlemans, 2001):


where H is ice thickness, w is glacier width, M is mean annual net surface balance, t is time, x is distance along the flowline and U is vertically averaged velocity parallel to the bed. This velocity is comprised of two components, one originating from deformation U d and the other from basal sliding U s:


Oerlemans (2001) reported that the flow parameter values f d = 1.9×10–24 Pa–3 s–1 and f s = 5.7×10–20 Pa–3m2 s–1 provided good results when applied to several glaciers. However, it is not obvious a priori that these are also the best values for glacier AX010. These parameters depend on the specifics of a glacier, primarily its geometric setting, and may also vary temporally at a particular glacier, thereby serving a tuning purpose. We therefore adopt the above-mentioned values only as a reasonable starting point. A detailed description of how the flow parameters were tuned for our glacier is given in section 3. Equation (2) results from a generalized flow law and relates the total volume flux through a vertical section to the driving stress S d:


where ρ is ice density (850 kgm –3), g is gravitational acceleration (9.81ms–2) and h is surface elevation.

Although the model is basically one-dimensional, the two-dimensional glacier geometry is accounted for by a varying width factor for each gridpoint (see Oerlemans, 1997). Assuming a trapezoidal-shaped valley, the width factor (w in Equation (1)) depends on the local ice thickness at the central flowline and the side slope of the valley cross-section:


Here w 0 is the bottom width of the glacier. The second term on the right-hand side accounts for the width contribution of the valley’s side walls inclined at a slope of 45˚.

A finite-difference scheme was used to numerically integrate the continuity equation (Equation (1)) for ice thickness forward in time, under the additional condition that ice thickness cannot be negative (i.e. H i≥0). The initial condition at t = 0 year is zero ice thickness everywhere (i.e. Hi(t=0) = 0m). The time integration is done with a forward explicit scheme. A time-step of 0.01 years was found to be sufficient to obtain a stable and accurate solution. This was checked by considering a steady-state glacier with a simple geometry (i.e. a constant bed slope and a unit glacier width with full lateral symmetry together with a linear mass-balance profile with elevation). Considering a linear stability criterion, the model time-step (Δt) should satisfy the condition that Δt≤(Δx)2/4D. Here Δx is the horizontal grid spacing and D is the average diffusivity (of the order of 2000 m2 a–1 on our glacier). A horizontal grid spacing of 10m was employed throughout this study, which satisfies the above criterion. Further steady-state experiments were performed to confirm whether the model conserves mass, and whether the ice flux, ice thickness and ice velocity along the central flowline were in agreement with theoretical considerations. As expected for our model set-up, all results exhibited appropriate symmetry with respect to the midway point. Furthermore, the ice flux anywhere in the glacier was equal to the integrated upstream mass balance, and the mean surface mass balance over the entire glacier was zero to within an acceptable accuracy.

3. Basic Input Data for Glacier AX010

In this section, we discuss how basic input data to run the model on glacier AX010 were collected. The input data related to the glacier geometry were extracted from a topographic map corresponding to the year 1996. Due to the fairly simple geometry of the glacier, a central flowline was easily constructed, ensuring that it remained perpendicular to all the contours. This line runs from the head of the glacier for 3 km downstream. A digital elevation model (DEM) was prepared based on the topographic map. From the DEM, we extracted the surface profile of the glacier along the assumed flowline. In addition, we accounted for the ice thicknesses measured in 1995 (Kadota and others, 1997). The expression that exhibits the relation between the ice thickness and the local driving stress (Equation (3)) was used to reconstruct the glacier bed profile in places where direct ice-thickness measurements were lacking (Oerlemans, 2001). We also recorded the surface width of the glacier from the DEM by simply measuring the length of a contour (within the glacier domain) that passes through every gridpoint. The bottom width of the glacier at each gridpoint was then calculated using Equation (4).

Besides the extraction of the geometric input, we plotted the mass-balance data recorded during the years 1995–99 against the difference between their elevation and the corresponding equilibrium-line altitude (ELA). Before plotting, the 1998 dataset was smoothed by taking a two-point average. A few outliers were removed to ensure a good fit around the ELA. The measurements were well fitted linearly, giving a value for the balance gradient of 0.1 mw.e. a–1m–1, with a correlation coefficient of 0.91 (Fig. 2). This value for the balance gradient is corroborated by the value of 10mmm–1 found by Harper and Humphrey (2003). At any time t, the mass-balance profile at the glacier surface, M(h), can now be described as:


Fig. 2. Extraction of the balance gradient: mass-balance data recorded during 1995–99 (WGMS, 2005) are plotted against the elevation with respect to the ELA. This gives a balance gradient of 0.01 mw.e. a–1m–1 (R 2 = 0.91). Triangles represent data outliers that were not used for the analysis.

Here the first term on the right-hand side is the reference balance profile, independent of time, and the second term represents the balance perturbation, assumed independent of altitude h.

With these geometric data and the associated mass-balance profile, the model was initially run with the flow parameters reported in section 2. ΔM was set to 0 mw.e. a–1. By slightly adjusting the ELA, a steady-state glacier was produced, whose terminus position was at the same location as on the 1996 topographic map. The corresponding value of the ELA was 5174ma.s.l. Comparing the modelled with the observed ice surface gave a poor fit for several parts of the glacier. Although some of these mismatches may merely reflect the non-steady-state conditions of the real glacier, we chose to fine-tune the model by adjusting the flow parameters until a good fit was obtained. We could have obtained a similarly good match between observed and simulated surface elevation by adjusting the bed topography (Oerlemans, 1997) but opted to respect the known ice-thickness measurements. A best fit between the modelled and observed ice surface was obtained by dividing both flow parameters by a factor of 30. We conserved the ratio between the sliding and deformation parameters f s and f d, mainly because we had no field evidence to support a different partitioning between these two flow mechanisms.

Finally a few small adjustments were made to the bed profile in zones without thickness measurements. The final adjusted bed is shown in Figure 3. Here the ice surface results from a time-dependent simulation with historic mass-balance forcing as discussed further below (section 5). Down from x = 1300 m, the bottom width was taken as a constant value of 80 m. These adjusted parameters were used for all the subsequent experiments.

Fig. 3. Longitudinal profile of the glacier along the central flowline for AD 1996. The black bold curve is the bed elevation reconstructed from the observed surface elevation (thin curve) using Equation (3). The dotted curve is the simulated surface elevation. Open circles along the bed profile represent the locations where ice thickness was measured.

4. Basic Sensitivity Experiments

In the model, the climate forcing is introduced by variations in the local net mass balance. Starting from the net balance curve, mass-balance variations can be imposed by assuming either a balance shift or an ELA shift. Since the mass balance versus elevation with respect to the ELA exhibits a linear relationship (Fig. 2), the two approaches are equivalent. Here we express mass-balance changes as a balance shift as was done by Huybrechts and others (1989).

Basic model sensitivity versus stepwise perturbations in the net mass balance is shown in Figure 4a. It is clear that the model exhibits a similar behaviour for both positive and negative shifts in the mass balance except for a higher perturbation in negative mass balance (<–0.20mw.e. a–1). An abrupt increase in sensitivity can be seen for a mass-balance perturbation between –0.20 and –0.25ma–1. The total length of the glacier suddenly reduces from 1050 to 650 m for the respective perturbations. Closer inspection of the ice-thickness profile reveals a geometric origin. There exists a thin ice layer downstream from x∽900m for a balance perturbation of –0.20ma–1 that is associated with a bedrock sill (Fig. 3). Increasing the negative mass-balance perturbation from –0.20ma–1 to –0.25ma–1 is apparently enough to melt the entire tongue over a distance of almost 500 m.

Fig. 4. (a) Steady-state glacier length versus mass-balance perturbations; and (b) reaction of the glacier front position to a stepwise change in surface mass balance of given magnitude (mw.e. a–1). The corresponding e-folding length response time (τ) is given in years.

Next we explore the length response time of glacier AX010. This timescale is defined as the e-folding time period required to reach (1 – e–1), i.e. 63% of the final response after a stepwise perturbation in climate conditions. This intrinsic timescale yields information on how quickly a glacier adjusts to a sudden change in mass balance. This knowledge is very useful to determine how far back in time the model needs to be spun up (e.g. to simulate the observed historical front variations, as discussed in section 5).

To estimate the length response time, stepwise changes in the mass balance in the range –0.5 to +0.5mw.e. a–1 were applied to a steady-state reference glacier corresponding to the year 1996. Changes in glacier length and their corresponding response times are shown in Figure 4b. The presence of the icefall between x∽900 and x∽1200m along the longitudinal glacier profile causes a rapid retreat of the glacier starting from ∽1200m for the larger negative mass-balance changes. This can be clearly seen in the figure for a mass-balance perturbation of –0.50ma–1. The length response time of glacier AX010 was found to be in the range 50–85 years and similar for both positive and negative perturbations. Although this estimate matches well with the length response time of other glaciers as estimated using the same method for Glacier d’Argentière, France, (27–45 years; Huybrechts and others, 1989) Nigardsbreen, Norway, (63– 73 years; Oerlemans, 1997) and Sofiyskiy glacier, Russian Altai, (73–84 years; De Smedt and Pattyn, 2003), we believe that our result is quite high for such a small subtropical valley glacier. Such a high value results mainly from the adjustments made in the flow parameters (reduction by a factor of 30): the smaller these parameters, the stiffer the ice becomes, thereby making the glacier respond more slowly to a climate forcing. In fact, reducing the flow parameters is equivalent to making the ice colder. This might, however, be true for high-altitude Himalayan valley glaciers such as glacier AX010.

5. Simulation of Historical Front Variations

The fairly short record of the terminus positions clearly reflects a continuous receding trend of glacier AX010 since 1978 (Fujita and other, 2001; Shrestha and Shrestha, 2004). In addition, it is believed that the glacial systems in the Nepal Himalaya were much longer during the Little Ice Age (LIA) than today (Thomas and Rai, 2005). With this information, an attempt was made to simulate the historic record of the glacier’s front positions. A comparison between modelled and observed front positions is expected to explain more about the type of climatic variables decisively affecting the glacier’s mass balance (Huybrechts and others, 1989). Since the mass-balance record of the glacier is limited to a few years, a proper forcing function has to be established. Here we implemented mass-balance perturbations as a function of temperature and precipitation separately. As glacier AX010 is a summer-accumulation-type glacier, summer balance is considered representative for the annual balance. The winter precipitation, which accounts for <20% of the annual value, does not contribute much to accumulation, as it is blown away by high-speed gusts that can exceed 140 kmh–1, especially in January (Tartari and others, 1998), and ablation in winter is negligible.

A first experiment was conducted under the assumption that the mass-balance perturbation is a sole function of summer temperature T s. We combined T s data from Darjeeling and Kathmandu stations (Fig. 5a), which were selected as the nearest stations to the glacier (∽150km away, to the east and west respectively) that have a long series of data. As the corresponding climatic anomalies of these stations match well for overlapping periods, one can assume that these datasets are likely to represent climatic variations around the glacier well. These data only extend back to 1882, so proxy data have to be used for the preceding period. We have used the general trend of the temperature anomaly derived from the Dasuopu (Qinghai– Tibetan Plateau) ice core for the period 1600–1881 (Feng and Hu, 2005).

Fig. 5. (a)Darjeeling and Kathmandu summer mean temperature anomalies with respect to the 1951–80 mean. (b) Mass-balance perturbation solely based upon the temperature anomaly, constructed by combining a general trend of T s anomaly derived from the Dasuopu ice core (Feng and Hu, 2005) for the period before 1882, the Darjeeling anomaly for 1882–1980, and the Kathmandu anomaly thereafter. (c) Simulated and observed front positions.

Assuming that the same linear balance model also applies to the past, the mass balance is perturbed, ΔM (mw.e. a–1), equally for all altitudes in the experiments as follows:


Here ΔT (˚C) is the temperature anomaly, C 1 (mw.e. a–1 ˚C–1) controls the terminus response amplitude, and C 2 (˚C) is a parameter to be optimized in order to obtain the right glacier length (Huybrechts and others, 1989). Parameter C 1 gives the sensitivity of the mass balance to the summer temperature. This glacier property depends on the geometric and climatic setting, and is unique for each glacier. Parameter C 2, on the other hand, serves a tuning purpose. It controls the absolute length of the glacier but hardly affects the amplitude of length variations. Both parameters are considered invariant in time.

The value of C 1 for glacier AX010 has been derived from the mass-balance data recorded during the years 1995–99 assuming an atmospheric lapse rate of –6.5˚Ckm–1. Excluding the data of 1998 when the ELA was higher than the highest altitude of the glacier, the maximum difference in ELA is 23 m. The corresponding difference in mass balance is 0.111 m (WGMS, 2005). Consequently, a 1˚C change in temperature causes the mass balance to change by 0.7 mw.e. a–1. Our value of C 1 = 0.7 mw.e. a–1 ˚C–1 agrees well with the range of mass-balance sensitivities to summer temperature (0.7–1.0mw.e. a–1 ˚C–1) as estimated by Kayastha and others (1999).

Our simulation of historic front variations started at AD 1200 (t = 0 year) with zero ice thickness. For the mass balance, the ELA was set at 5210ma.s.l. (the 1995–99 mean) with a value of ΔM = +0.32mw.e. a–1 until AD 1600. This period to spin the model up is long enough to have negligible effect on the reconstructed history of the terminus position since the LIA. Normally, a glacier carries in its memory the mass-balance history over a period equivalent to one or two times the response time (Oerlemans, 1997). After AD 1600, the model was forced with ΔM values (Equation (6)) derived from climate reconstructions. The best fit between the modelled and observed front positions (Fig. 5c) was achieved for C 2 = +0.05˚C. The corresponding evolution of mass-balance forcing is shown in Figure 5b.

As shown in Figure 5c, the model projects that the glacier was slightly advancing up to a maximum length of around 2500 m until ∽AD 1850, with a general retreat thereafter. The simulated and the observed history of retreat match remarkably well (root-mean-square error (rmse) ∽27 m; linear correlation coefficient 0.95). This seems to suggest that the glacier length history can be defined well by the temperature anomaly alone. In Figure 3, a comparison between the modelled and observed AD 1996 glacier stand can be seen in terms of ice thickness (rmse ∽11 m). The model underestimates the ice thickness in the accumulation zone, but this disagreement will probably have little effect on the results of the historical front simulations, as variations in front position are mostly governed by changes in ablation area (De Smedt and Pattyn, 2003).

Another, similar experiment was conducted under the assumption that the balance perturbation is a sole function of summer precipitation P s. For this purpose, we used the precipitation anomalies from Kathmandu and Darjeeling stations (Fig. 6a), as well as proxy data derived from the Dasuopu ice core that have a correlation coefficient of 0.57 with Nepal rainfall for the period 1901–95 (Duan and others, 2004). Here also a linear mass-balance model is assumed to be applicable to the past:


Fig. 6. (a)Darjeeling and Kathmandu summer precipitation anomalies (by ratio) with respect to the 1951–80 mean. (b) Mass-balance perturbation solely based upon the precipitation anomaly, constructed by combining a general trend of P s anomaly derived from the Dasuopu ice core (Duan and others, 2004) for the period before 1870, the Darjeeling anomaly for 1870–1972 and the Kathmandu anomaly thereafter. (c) Simulated and observed terminus positions.

where ΔP is the dimensionless precipitation anomaly expressed as a ratio. The value of C 3 (mw.e. a–1) describes the amplitude of the terminus response. Likewise, C 4 (mw.e. a–1) is a tuning parameter controlling mean glacier length.

Parameter C 3 was estimated from the empirical equations defined by Kadota and Ageta (1992) based on records of mass balance and climatic parameters in 1978/79. During that season, T s and P s at the glacier terminus (4958ma.s.l.) were recorded as 2.4˚C and 1.45 m, respectively. Using the same atmospheric lapse rate of –6.5˚Ckm–1, we found that the mean T s at the glacier head (5360ma.s.l.) would be about – 0.2˚C. These equations were used to plot the mass balance against temperature on glacier AX010 for a varying precipitation. The plot revealed that a 10% change (both positive and negative) in precipitation, i.e. 0.145 m, causes a variation in mass balance between 0.145 m (for T s≤–3.2˚C) and 0 m (for T s≥3.4˚C). For the same change in P s, the average value of mass-balance change for a temperature range –0.2˚C (at the glacier head: 0.123 mw.e. a–1) to +2.4˚C (at the terminus: 0.036 mw.e. a–1) was found to be 0.08 mw.e. a–1. This means that, assuming the linear balance model, a 1 m change in precipitation results in a 0.55 mw.e. a–1 change in mass balance, which is used as the value of C 3. The physical meaning is that on average 55% of the summer precipitation falls as snow and nourishes the glacier.

Further experimentation revealed a suitable value of C 4 = 0.40mw.e. a–1. The resulting mass-balance perturbations and front positions are shown in Figure 6b and c. The figure reveals that the length history of glacier AX010 as solely defined by the precipitation forcing does not match the recorded length history appreciably (rmse ∽110 m). This suggests that precipitation by itself is not a good predictor to explain the glacier history. This result could have been anticipated since there is a little correspondence between the longer-term precipitation trend and glacier length. To investigate the future evolution of glacier AX010, it may thus suffice to assume that temperature can solely define the glacier’s response.

6. Future Evolution of the Glacier

Shrestha and Shrestha (2004) mentioned that the observed trends in temperature and precipitation in Nepal are in good agreement with most climate-model predictions. We therefore opted to choose projected scenarios rather than prescribe fully arbitrary ones to force the model in the future. All the scenarios considered in this study are based upon the Intergovernmental Panel on Climate Change (IPCC) projections for the high-altitude Himalaya and the Tibetan Plateau (Christensen and others, 2007). We selected three different scenarios of temperature change by the end of the 21st century: +3˚C, +4.5˚C and +6˚C with respect to the 1980–99 mean. Likewise, five different scenarios of precipitation change were chosen in the range –30% to +30% with respect to the 1980–99 mean. The future evolution simulations follow on from the experiment of the historical front variations as solely defined by the temperature forcing.

In the first series of experiments we keep the precipitation constant at the 1980–99 mean. Only the temperature anomaly is varied, linearly by +3˚C, +4.5˚C and +6˚C by AD 2100, respectively. Smaller temperature rises are excluded because +3˚C and +6˚C are the lower and upper limits on IPCC projections for the region considered. Figure 7a displays the applied mass-balance perturbation throughout this century. The resulting evolution of the glacier length is shown in Figure 7b. As shown by the solid curve in this figure, the glacier length is reduced to 940 m by AD 2100 if the temperature remains constant at the 1980–99 mean. However, the glacier ceases to exist at AD 2083, 2056 or 2049 if the temperature increases linearly by the end of the 21st century by +3˚C, +4.5˚C or +6˚C, respectively. During its decay, the glacier retains ice mass in a number of segments depending upon the basal topography. The possible number of segments for different scenarios of the temperature rise is shown in Figure 7c. It is seen that if the temperature remains constant, the glacier will not disintegrate into segments within this century. For all of the rising-temperature scenarios considered here, however, the glacier breaks into a number of pieces before complete drawdown. Fragmentation starts at around AD 2020, with the maximum possible number of segments equal to three.

Fig. 7. (a)Applied mass-balance perturbations, dM; (b) the corresponding evolution of glacier length; and (c) the number of glacier segments for the first series of scenarios with constant precipitation (at the 1980–99 mean) and temperature varied in the range +3 to +6˚C by AD 2100.

Before considering another series of experiments with a varying precipitation rate, we explore the threshold temperature for the complete drawdown of the glacier within this century, by imposing temperature rises of 0–6˚C in 0.1 ˚C increments. We determine +2.5˚C as the threshold value above which the glacier disappears completely within this century. The year of glacier drawdown for each 0.1˚C increment in temperature (above the threshold value) is plotted in Figure 8. The general nature of the curve is exponential. If the temperature increases by less than 2.5˚C with respect to the 1980–99 mean, glacier AX010 is expected to survive until at least the end of this century. On the other hand, if the temperature rises linearly by more than 6˚C by AD 2100, the glacier will only last until the 2040s.

Fig. 8. Time (year AD) of complete glacier drawdown as a function of a linear temperature rise until AD 2100. For complete melting of glacier AX010 within this century, the temperature rise by AD 2100 must be higher than a threshold of 2.5˚C with respect to the 1980–99 mean.

In the next series of experiments, we maintain a constant rise in temperature but additionally vary the precipitation in the range –30 to +30%. For a constant rise in temperature of 3˚C by AD 2100, the model projects that the glacier will vanish during the 18 year period AD 2075–93, the timing being determined by the percentage of precipitation change: AD 2075, 2079, 2083, 2088 or 2093 if the precipitation linearly changes by –30%, –15%, 0%, +15% or +30% respectively by AD 2100. In all these scenarios, the glacier is likely to break into two pieces between ∽AD 2030 and ∽AD 2065, thereafter retaining ice only in a single segment.

If the temperature rises linearly by 4.5˚C by AD 2100, for all scenarios of precipitation change the glacier will vanish within the 6 year period AD 2054–60. The model projects fragmentation into a maximum of three pieces between AD 2026 and 2029 for all scenarios of precipitation change, and then ice retention in two segments until the year of glacier drawdown.

If the temperature rises linearly by 6˚C by AD 2100, for all scenarios of precipitation change the glacier will vanish within the 2 year period AD 2048–50. It is likely to fragment into a maximum of three segments between AD 2025 and 2027, then retain the ice mass in two patches until AD 2038–42 depending on the percentage of precipitation change. The remaining ice is then retained in a single segment until the year of glacier drawdown.

These results demonstrate that the given range of precipitation change (–30 to +30%) produces a narrower range of year of glacier drawdown (2 years) for a higher temperature rise (+6˚C) than for a smaller temperature rise (18 years at +3˚C). This implies that for a higher temperature change, precipitation will have less effect on mass balance, and thus affect the glacier’s behaviour minimally. This could be mainly due to the overriding effects of temperature over precipitation not only by enhancing the rate of ablation but also by reducing the fraction of snow accumulation.

7. Conclusions and Recommendations

The simulation of historical front variations of glacier AX010 shows that temperature, but not precipitation, can solely define the history of glacier fluctuations. This does not imply that precipitation plays no role in defining the evolution of Himalayan glaciers such as glacier AX010 that are mainly fed by the Asian summer monsoon; but it indicates the overriding effect of summer temperature over summer precipitation in defining ice melting and the snow accumulation fraction. Had we had better quantitative information regarding 19th-century glacier front positions, we could have refined the ratio of the contribution of temperature to that of precipitation in defining the historical fluctuations of the glacier. Investigation of the future behaviour of the glacier confirms that it will shrink continuously throughout this century, even if the climatic situation remains constant at the 1980–99 mean. The glacier is projected to vanish completely within this century for all trends with a temperature rise of ≥2.5˚C by AD 2100 with respect to the 1980–99 mean and precipitation changes tested in this study. The fidelity of these projections is restricted by the assumptions made in this study. An important constraint of the model is probably the choice of rather low flow parameters. A more extensive set of ice-thickness measurements on more locations of the glacier would also have helped to minimize possible errors made in defining the bed profile. In addition, the collection of any information left by end moraines would have been useful to better validate the model during the LIA glacier maximum.


We thank the Vlaamse Interuniversitaire Raad (VLIR), Belgium, for awarding a scholarship to S.A. to perform this research at the Vrije Universiteit Brussel. Thanks are also due to two reviewers, H. Blatter and M. Kessler, for thoughtful suggestions. We acknowledge support and suggestions from S. Marshall, G.K.C. Clarke, B. De Smedt and R.B. Kayastha.


Ageta, Y. and Kadota, T.. 1992. Predictions of changes of glacier mass balance in the Nepal Himalaya and Tibetan Plateau: a case study of air temperature increase for three glaciers. Ann. Glaciol., 16, 89–94.
Ageta, Y., Ohata, T., Tanaka, Y., Ikegami, K. and Higuchi, K.. 1980. Mass balance of glacier AX010 in Shorong Himal, east Nepal during the summer monsoon season. Seppyo, J. Jpn. Soc. Snow Ice, 41, 34–41.
Brohan, P., Kennedy, J.J., Harris, I., Tett, S.F.B. and Jones, P.D.. 2006. Uncertainty estimates in regional and global observed temperature changes: a new data set from 1850. J. Geophys. Res., 111(D12), D12106. (10.1029/2005JD006548.)
Christensen, J.H. and 16 others. 2007. Regional climate projections. In Solomon, S. and 7 others, eds. Climate change 2007: the physical science basis. Contribution of Working Group I to the Fourth Assessment Report of the Intergovernmental Panel on Climate Change. Cambridge, etc., Cambridge University Press, 847–940.
De Smedt, B. and Pattyn, F.. 2003. Numerical modelling of historical front variations and dynamic response of Sofiyskiy glacier, Altai mountains, Russia. Ann. Glaciol., 37, 143–149.
Duan, K., Yao, T. and Thompson, L.G.. 2004. Low-frequency of southern Asian monsoon variability using a 295-year record from the Dasuopu ice core in the central Himalayas. Geophys. Res. Lett., 31(16), L16209. (10.1029/2004GL020015.)
Feng, S. and Hu, Q.. 2005. Regulation of Tibetan Plateau heating on variation of Indian summer monsoon in the last two millennia. Geophys. Res. Lett., 32(2), L02702. (10.1029/2004GL021246.)
Fujita, K., Kadota, T., Rana, B., Kayastha, R.B. and Ageta, Y.. 2001. Shrinkage of Glacier AX010 in Shorong region, Nepal Himalayas in the 1990s. Bull. Glaciol. Res. 18, 51–54.
Hall, M.H.P. and Fagre, D.B.. 2003. Modeled climate-induced glacier change in Glacier National Park, 1850–2100. BioScience, 53(2), 131–140.
Harper, J.T. and Humphrey, N.F.. 2003. High altitude Himalayan climate inferred from glacial ice flux. Geophys. Res. Lett., 30(14), 1764. (10.1029/2003GL017329.)
Huybrechts, P., de Nooze, P. and Decleir, H.. 1989. Numerical modelling of Glacier d’Argentière and its historic front variations. In Oerlemans, J., ed. Glacier fluctuations and climatic change. Dordrecht, etc., Kluwer Academic Publishers, 373–389.
Ikegami, K. and Ageta, Y.. 1991. Ice flow of glacier AX010 in the Nepal Himalaya. Bull. Glacier Res. 9, 17–22.
Kadota, T. and Ageta, Y.. 1992. On the relation between climate and retreat of Glacier AX010 in the Nepal Himalaya from 1978 to 1989. Bull. Glacier Res. 10, 1–10.
Kadota, T., Fujita, K., Seko, K., Kayastha, R.B. and Ageta, Y.. 1997. Monitoring and prediction of shrinkage of a small glacier in the Nepal Himalaya. Ann. Glaciol., 24, 90–94.
Kayastha, R.B., Ohata, T. and Ageta, Y.. 1999. Application of a mass-balance model to a Himalayan glacier. J. Glaciol., 45(151), 559–567.
Oerlemans, J. 1997. A flowline model for Nigardsbreen, Norway: projection of future glacier length based on dynamic calibration with the historic record. Ann. Glaciol., 24, 382–389.
Oerlemans, J. 2001. Glaciers and climate change. Lisse, etc., A.A. Balkema.
Ohata, T. and Higuchi, K.. 1980. Heat balance study on glacier AX010 in Shorong Himal, east Nepal. Seppyo, J. Jpn. Soc. Snow Ice, 41, 42–47.
Shrestha, M.L. and Shrestha, A.B.. 2004. Recent trends and potential climate change impacts on glacier retreat/glacier lakes in Nepal and potential adaptation measures. In Organisation for Economic Cooperation and Development, OECD Global Forum on Sustainable Development: Development and Climate Change. Paris, Organisation for Economic Cooperation and Development, 5–14.
Tartari, G., Verza, G.P. and Bertolami, L.. 1998. Meteorological data at the Pyramid Observatory Laboratory (Khumbu Valley, Sagarmatha National Park, Nepal). Mem. Ist. Ital. Idrobiol., 57, 23–40.
Thomas, K. and Rai, S.C. 2005. An overview of glaciers, glacier retreat, and subsequent impacts in Nepal, India and China. Kathmandu, World Wildlife Fund Nepal Programme.
World Glacier Monitoring Service (WGMS). 1998. Fluctuations of glaciers 1990–1995 with addendas from earlier years (Vol. VII), ed. Haeberli, W., Hoelzle, M., Suter, S. and Frauenfelder, R.. IAHS/UNEP/UNESCO, World Glacier Monitoring Service, Zürich.
World Glacier Monitoring Service (WGMS). 2005. Fluctuations of glaciers 1995–2000 (Vol. VIII), ed. Haeberli, W., Zemp, M., Frauenfelder, R., Hoelzle, M. and Kääb, A.. IAHS/UNEP/UNESCO, World Glacier Monitoring Service, Zürich.