Modelling glacier change in the Everest region, Nepal Himalaya

Abstract. In this study, we apply a glacier mass balance and ice redistribution model to examine the sensitivity of glaciers in the Everest region of Nepal to climate change. High-resolution temperature and precipitation fields derived from gridded station data, and bias-corrected with independent station observations, are used to drive the historical model from 1961 to 2007. The model is calibrated against geodetically derived estimates of net glacier mass change from 1992 to 2008, termini position of four large glaciers at the end of the calibration period, average velocities observed on selected debris-covered glaciers, and total glacierized area. We integrate field-based observations of glacier mass balance and ice thickness with remotely sensed observations of decadal glacier change to validate the model. Between 1961 and 2007, the mean modelled volume change over the Dudh Koshi basin is −6.4 ± 1.5 km3, a decrease of 15.6% from the original estimated ice volume in 1961. Modelled glacier area change between 1961 and 2007 is −101.0 ± 11.4 km2, a decrease of approximately 20% from the initial extent. The modelled glacier sensitivity to future climate change is high. Application of temperature and precipitation anomalies from warm/dry and wet/cold end-members of the CMIP5 RCP4.5 and RCP8.5 ensemble results in sustained mass loss from glaciers in the Everest region through the 21st century.


Introduction
High-elevation snow and ice cover play pivotal roles in Himalayan hydrologic systems (e.g. Viviroli et al., 2007;Immerzeel et al., 2010;Racoviteanu et al., 2013). In the monsoon-affected portions of the Himalayas, meltwater from seasonal snowpacks and glaciers provides an important source of streamflow during pre-and post-monsoon seasons, while rainfall-induced runoff during the monsoon dominates the overall hydrologic cycle . Against this backdrop, changes in glacier area and volume are expected to have large impacts on the availability of water during the dry seasons (Immerzeel et al., 2010), which will impact agriculture, hydropower generation, and local water resources availability. In the current study, our main objectives are to calibrate and test a model of glacier mass balance and redistribution, and to present scenarios of catchment-scale future glacier evolution in the Everest region.

Study area and climate
The ICIMOD (2011) inventory indicates that the Dudh Koshi basin in central Nepal contains a total glacierized area of approximately 410 km 2 (Fig. 1). The region contains some of the world's highest mountain peaks, including Sagarmatha (Mount Everest), Cho Oyu, Makalu, Lhotse, and Nuptse. The Dudh Koshi River is a major contributor to the Koshi River, which contains nearly one-quarter of Nepal's exploitable hydroelectric potential. Approximately 110 km 2 , or 25 % of the total glacierized area, is classified as debris-covered (Fig. 2), with surface melt rates that are typically lower than those observed on clean glaciers due to the insulating effect of the debris (Reid and Brock, 2010;Lejeune et al., 2013).
The climate of the region is characterized by pronounced seasonality of both temperature and precipitation. At 5000 m (see analysis below), mean daily temperatures range between atures ranging between −25 and +10 • C. During the monsoon period (June-September), temperatures at 5000 m are greater than 0 • C and variability is low. The majority of annual precipitation (approximately 77 %, derived from gridded climate fields, see below) falls between 1 June and 30 September during the summer monsoon (Wagnon et al., 2013). An additional 14 % of precipitation occurs during the pre-monsoon period (March-May), with little or no precipitation during the post-monsoon and winter seasons. The interaction between moisture advected from the Indian Ocean during the monsoon and the two-step topography of the Dudh Koshi region (foothills, main ranges) results in two spatial maxima of precipitation (Bookhagen and Burbank, 2006).

Himalayan glaciology
The current status of glaciers varies across the Hindu Kush Himalayan (HKH) region. Most areas have seen pronounced glacier retreat and downwasting in recent years Kääb et al., 2012;Yao et al., 2012), though some areas, such as the Karakoram and Pamir ranges, have experienced equilibrium or even slight mass gain (Gardelle et al., , 2013Jacob et al., 2012). In the Everest region ( Fig. 1), Gardelle et al. (2013) find an average annual rate of mass loss of −0.26 ± 0.13 m w.e. yr −1 between 2000 and 2011, while Nuimura et al. (2012) estimate mass loss rates of calculated for glaciers in the Khumbu region (Bolch et al., 2008a, b). Areal extents of glaciers in Sagarmatha National Park decreased 5 % during the second half of the 20th century (Bolch et al., 2008b;Salerno et al., 2008;Thakuri et al., 2014). These estimates do not distinguish between debriscovered and clean-ice glaciers.
One consequence of glacier retreat in the Himalayas is the formation of proglacial lakes, which may pose a risk to downstream communities. Terminus retreat at Lumding and Imja glaciers, measured at −42 and −34 m yr −1 , respectively, between 1976 and 2000 increased to −74 m yr −1 at both glaciers between 2000 and 2007 (Bajracharya and Mool, 2010). Rapid terminus retreat results in the growth of proglacial lakes which are dammed by lateral and terminal moraines (Bolch et al., 2008b;Benn et al., 2012;Thompson et al., 2012). The failure of moraine dams in the Koshi River basin has led to 15 recorded glacier lake outburst flood (GLOF) events since 1965, with flows up to 100 times greater than average annual flow (Chen et al., 2013), and the frequency of GLOFs in the Himalayas is believed to have increased since the 1940s (Richardson and Reynolds, 2000). Changes in glacier extents and volumes in response to climate change thus have important impacts not only on water resources availability but also on geophysical hazards.
The climate sensitivity of a glacier depends primarily on its mass balance amplitude. Glaciers in wetter climates typi- cally extend to lower elevations, and are thus more sensitive to temperature changes than those in dry climates (Oerlemans and Fortuin, 1992). Himalayan glaciers, and glaciers of the Dudh Koshi in particular, present a unique challenge as observations of temperature and precipitation at high elevations are scarce. Regionally, the climate varies from monsoon-dominated southern slopes to relatively dry leeward high-elevation regions. Accordingly, equilibrium line altitudes (ELAs) in the region vary both spatially and temporally but generally range from 5200 m in the south to 5800 m in northern portions of the basin (Williams, 1983;Asahi, 2010;Wagnon et al., 2013). Nearly 80 % of the glacierized area in the Dudh Koshi basin lies between 5000 and 6000 m (Fig. 2), and the region is expected to be sensitive to climatic changes.

Historical and projected climate trends
Analyses of climate trends in the region are limited, primarily due to the lack of long-term records (Shrestha and Aryal, 2011). Available studies indicate that the mean annual temperatures have increased in the region, and particularly at high elevations (Shrestha et al., 1999;Rangwala et al., 2009;Ohmura, 2012;Rangwala and Miller, 2012). Reported mean annual temperature trends range between 0.025 and 0.06 • C yr −1 for the periods 1971 to and 1977to 1994, respectively (Shrestha and Aryal, 2011Qi et al., 2013). Changes in temperature are particularly important for monsoon-type glaciers, which are sensitive to the elevation of the rain/snow threshold during the monsoon season . Results from the CMIP5 (Climate Modelling Intercomparison Project) ensemble suggest that temperatures in the region will increase between 1.3 and 2.4 • C over the period 1961-1990to 2021-2050(Lutz et al., 2013, which correspond to rates of 0.021 to 0.040 • C yr −1 . Precipitation amounts, timing, and phase will affect glacier responses on both annual and decadal timescales. In the greater Himalayas, trends in precipitation totals appear to be mixed and relatively weak (Mirza et al., 1998;Gautam et al., 2010;Dimri and Dash, 2012;Qi et al., 2013), though the observational network is composed mostly of low-elevation valley stations that may not reflect changes in snowfall amounts at higher elevations. General circulation model projections suggest both increased monsoon precipitation (Kripalani et al., 2007) and delayed monsoon onset www.the-cryosphere.net/9/1105/2015/ The Cryosphere, 9, 1105-1128, 2015   temperature fields and resampled SRTM data, with period mean in grey, (b) daily standard deviation (σ ) of γ T , and (c) mean daily coefficient of determination (R 2 ) calculated from the linear regression of resampled SRTM elevations and APRHODITE cell temperatures. All temperature/elevation regressions are significant. (Ashfaq et al., 2009;Mölg et al., 2012) in the 21st century, while the change in total annual precipitation is mixed. In the Himalayas, the CMIP5 ensemble shows projected changes in precipitation between −8 to +15 % (Lutz et al., 2013;Palazzi et al., 2013).

Models of glacier change
In spite of the recent observed changes in glaciers in the Everest region, the reported climatic trends, the expected glacier sensitivity to climatic change, and the importance of glacier water resources in the region, few studies have attempted to model the historical or future response of these glaciers to climate change . Empirical mass balance and snowmelt and ice melt models have been developed from field observations (Ageta and Higuchi, 1984;Ageta and Kadota, 1992;Nakawo et al., 1999) and reanalysis products (Fujita and Nuimura, 2011;Rasmussen, 2013), and such approaches have been used to quantify glacier contributions to streamflow (Racoviteanu et al., 2013;Nepal et al., 2013). Projections of higher ELAs in the region (Fujita and Nuimura, 2011) and volume areascaling approaches (Shi and Liu, 2000;Cogley, 2011) indicate continued mass wastage in the future, yet impact studies on the response of glaciers to climate change require models that link mass balance processes with representations of glacier dynamics.
One-and two-dimensional models of glacier dynamics have been applied previously to the Khumbu Glacier  and the East Rongbuk Glacier , respectively. However, these and higher-order models of glacier dynamics are severely limited by input data availability (e.g. bed topography, ice temperatures, basal water pressure) and uncertainties in key model parameters, and have not been applied at catchment scales in the region. Debris-covered glaciers, which compose 25 % of total glacierized area, present additional modelling challenges, and validation is also limited by the availability of data. Relatively coarse methods of simulating future glacier change (e.g. Stahl et al., 2008) can be improved by applying models that can reasonably simulate key glaciological parameters (thickness, velocity, and mass redistribution).
The main objective of this study is to apply a glacier mass balance and redistribution model to the Dudh Koshi River basin, Nepal. To accomplish this, we (1) develop downscaling routines for temperature and precipitation; (2) calibrate and test the model with available field and remotely sensed observations; and (3) explore the modelled sensitivity of glaciers in the Everest region to future climate change with a suite of temperature and precipitation changes from the CMIP5 ensemble.

Daily climate fields
There are few observations of temperature and precipitation in the basin, and no temperature records longer than 15 years are available. To generate high-resolution fields of temperature (T ) and precipitation (P ) as inputs to the model, we use data from the APHRODITE (Asian Precipitation -Highly-Resolved Observational Data Integration Towards Evaluation of Water Resources) project (Yatagai et al., 2009(Yatagai et al., , 2012. APHRODITE products have been previously used to test regional climate model simulations in northern India  and the western Himalaya , and to compare precipitation data sets in the Himalayan region (Palazzi et al., 2013). For this study, we use APHRODITE T fields (V1204R1) that are based on daily station anomalies from climatological means, interpolated on 0.05 • grids and then resampled to 0.25 • fields, and we refer to Yatagai et al. (2012) for more details. The APHRODITE P fields (V1101) are based on a similar technique using precipitation ratios but incorporate a weighted interpolation scheme based on topographical considerations (Yatagai et al., 2012).
To generate high-resolution fields of T and P for the glacier mass balance model, we extract a 196 (14 × 14) grid cell subset of the daily APHRODITE T and P fields that covers the Koshi basin (Fig. 1). Approximate elevations for each 0.25 • grid cell are extracted from a resampled gap-filled Shuttle Radar Topography Mission (SRTM V4; Farr et al., 2007) digital elevation model (DEM). Based on this subset we derive relations between elevation and temperature and precipitation respectively at coarse resolution. We then use these relations in combination with the 90 m SRTM DEM to produce high-resolution daily climate fields.

Temperature
Downscaled temperature fields at daily 90 m resolution are computed as where γ T is the daily vertical temperature gradient (Fig. 3) derived from the 0.25 • APRHODITE temperatures and SRTM elevations, T 0 is the daily temperature intercept, and C DOY is a bias correction based on the day of year (Fig. 4). The bias-correction factor is computed from the mean daily temperature difference between observed and estimated mean daily temperatures at four stations operated by the Italian Everest-K2-National Research Centre (EVK2CNR; Fig. 1, Table 1), and it ranges from 3 to 6 • C. The EVK2CNR stations are independent of the APHRODITE product.

Precipitation
To calculate high-resolution daily precipitation fields from the APHRODITE subset, we prescribe daily precipitationelevation functions from the 0.25 • APHRODITE precipitation fields and resampled SRTM data. For each day, we calculate the mean precipitation in 500 m elevation bins (P 500 ) and prescribe a fitted linear interpolation function to estimate precipitation on the 90 m SRTM DEM (Fig. 5).
As APHRODITE fields are based on interpolated station data (Yatagai et al., 2012), there is a large uncertainty in the precipitation at high elevations. Independent tests of the precipitation downscaling approach were conducted by comparing precipitation observations from the EVK2CNR stations with precipitation estimated using the station elevation and the daily precipitation-elevation functions (Fig. 6). As EVK2CNR stations are not capable of measuring solid precipitation (Wagnon et al., 2013), we only examine days where only liquid precipitation (T > 0) is expected.
While orographic forcing of moist air masses typically produces increased precipitation with elevation, in very highelevation regions (i.e. those greater than 4000 m) both observations and models indicate that precipitation totals will decrease above a certain elevation (Harper and Humphrey, 2003;Mölg et al., 2009). This is due in part to the drying effect from upwind orographic forcing but is also related to the low column-averaged water vapour content indicated by the Clausius-Clapeyron relation. Given that there are no precipitation observations at elevations above 5300 m, and available evidence suggests that precipitation will likely decrease at high elevations, we scale estimated precipitation using a correction factor p cor : where p cor decreases from 1 at the height of a calibrated threshold elevation (Z c ; Table 2 Above 7500 m, we assume that precipitation amounts minus wind erosion and sublimation (Wagnon et al., 2013) are likely to be negligible. The total area above 7500 m represents only 1.2 % of the total basin area.

Glacier mass balance and redistribution
Following the methods of Immerzeel et al. (2012) and Immerzeel et al. (2013), daily accumulation and ablation between 1961 and 2007 are estimated from the gridded T and P fields. All calculations are based on the 90 m SRTM DEM. Daily accumulation is equal to the total precipitation when T < 0 • C, which is a conservative threshold with respect to other studies that have used values of 1.5 or 2 • C (Oerlemans and Fortuin, 1992), but this value has been used in previous Himalayan models . Daily ablation is estimated using a modified degree-day factor (ddf M ) that varies with DEM-derived aspect (θ ) and surface type: where ddf is the initial melt factor (in mm • C −1 d −1 ), and R exp is a factor which quantifies the aspect (or exposure) dependence of ddf . Initial values for melt factors for snow, ice, and debris-covered glaciers (Azam et al., 2014) are given in Table 2. The extent of debriscovered glaciers was extracted from the ICIMOD (2011) glacier inventory. To redistribute mass from accumulation to ablation areas, we use a simplified flow model which assumes that basal sliding is the principal process for glacier movement and neglects deformational flow. While cold-based glaciers have been observed on the Tibetan Plateau (Liu et al., 2009), warm-based glaciers and polythermal regimes have been identified on the monsoon-influenced southern slopes of the Himalayas (Mae et al., 1975;Ageta and Higuchi, 1984;Kääb, 2005;Hewitt, 2007). Our assumption in this case is a necessary simplification of the sliding and deformational components of ice flow, which have not yet been modelled at the basin scale in the Himalayas.
Glacier motion is modelled as slow, viscous flow using Weertman's sliding law (Weertman, 1957), which describes glacier movement as a combination of both pressure melting and ice creep near the glacier bed. Glacier flow is assumed to be proportional to the basal shear stress (τ b , Pa): Material roughness coefficient N m −2 s 1/3 1.80 × 10 9 ±5.00 × 10 8 1.51 × 10 8 ddf C Clean ice melt factor Here, v (unitless) is a measure of bedrock roughness, R (Pa m −2 s) is a material roughness coefficient, u is the sliding speed (m s −1 ) and n (unitless) is the creep constant of Glen's flow law, here assumed to equal 3 (Glen, 1955). The roughness of the bedrock (v) is defined as the dimension of objects on the bedrock divided by the distance between them. Smaller values for v indicate more effective regelation. R is a material roughness coefficient that controls the viscous shearing (Fowler, 2010). Basal shear stress (τ b ) is defined as where ρ is ice density (kg m −3 ), g is gravitational acceleration (m s −2 ), H is ice thickness (m), and β is surface slope ( • ). We assume that motion occurs only when basal shear stress exceeds the equilibrium shear stress (τ 0 = 80 000 N m −2 ; Immerzeel et al., 2012), and combine Eqs. (5) and (6) to derive the glacier velocity: For each time step, glacier movement in each cell is thus modelled as a function of slope, ice thickness, and assumed bedrock roughness. The total outgoing ice flux at each time step is then determined by the glacier velocity, the horizontal resolution, and the estimated ice depth. Ice transported out of a specific cell is distributed to all neighbouring downstream cells based on slope, with steeper cells receiving a greater share of the ice flux.
As avalanches can contribute significantly to glacier accumulation in steep mountainous terrain (Inoue, 1977;Scherler et al., 2011b), the model incorporates an avalanching component which redistributes accumulated snowfall (Bernhardt and Schulz, 2010). The approach assumes that all snow in a given cell is transported to the downstream cell with the steepest slope whenever snow-holding depth and a minimum slope angle is exceeded. The snow-holding depth is deep in flat areas and shallow in steep areas and decreases exponentially with increasing slope angle.
Based on field observations and an analysis of the slopes of glacierized pixels in the catchment (Fig. 7), we assign a threshold avalanching angle (β TH ) of 50 • . Change in ice thickness at each time step is thus the net result of ice flow through the cell, ablation, and accumulation from both precipitation and avalanching. Changes in glacier area and volume are calculated at daily time steps, and pixels with a snow water equivalent greater than 0.2 m w.e. are classified as glacier. The model does not assume steady-state conditions, and reported changes in volume and area thus represent transient states within the model.

Model initialization
Initial ice thickness for each glacierized grid cell is derived from Eq. (6): with a minimum prescribed slope of 1.5 • . We use τ 0 here, as the actual basal shear stress depends on the ice thickness. In the Dudh Koshi basin, Eq. (8) produces a total estimated glacier volume of 32.9 km 3 , based on the ICIMOD (2011) glacier inventory and SRTM DEM. While volumearea scaling relations are uncertain (Frey et al., 2013), empirical relations from Huss and Farinotti (2012) and Radić and Hock (2010) applied to individual glaciers generate basin- The Cryosphere, 9, 1105-1128, 2015 www.the-cryosphere.net/9/1105/2015/ wide volume estimates of 31.9 and 27.5 km 3 , respectively, which lends some support to the approach used here. From the initial ice thicknesses we estimate glacier thicknesses and extents in 1961 by driving the glacier mass balance and redistribution model with modified APHRODITE temperature fields. To simulate the observed climate in the region prior to 1961, temperatures in the initialization run are decreased by −0.025 • C yr −1 (Shrestha and Aryal, 2011), for a total decrease of −1.2 • C over the 47-year initialization period. Precipitation is left unchanged in the model initialization, and we use uncalibrated model parameters (Table 2).
Mass change at the end of the 47-year initialization period is close to zero, indicating that near-equilibrium conditions have been realized. Additional runs of the initialization period, with temperatures fixed at −1.2 • C, yield relatively small changes in glacier thickness (Fig. 8). However, it is possible that there are significant uncertainties in our estimates of initial (1961) thicknesses and extents, given the forcings and parameter set used, and the lag in glacier geometry responses to climate forcings.

Model calibration
From the modelled 1961 ice thicknesses and extents, the model is calibrated with six parameters: degree-day factors for clean ice (ddf C ), debris-covered ice (ddf I ), snow (ddf S ), and debris covered ice on the Khumbu Glacier (ddf K ), material roughness coefficient R, and elevation of the precipitation maximum Z C (Table 2). Initial simulations showed anomalous flow velocities of the Khumbu Glacier, which may be due to the assumption that basal sliding is the main process of movement. This may not hold given the steep icefall above the glacier tongue and the large high-altitude accumulation area. We have corrected for this by calibrating a specific melt factor for this glacier, though improved representation of the glacier dynamics should reduce the need for a separate ddf K . Twenty parameter sets (Table 3) were developed by varying the six calibration factors within specified ranges (Table 2). Initial values for each parameter were selected from published studies.
For each of the 20 runs (Table 4), we quantify the model skill by scoring (a) modelled and observed glacier extents at the termini of four large glaciers in the catchment (ICI-MOD, 2011), (b) the geodetically derived mean basin-wide glacier mass balance of −0.40 m w.e. yr −1 over the period 1992-2008 (Nuimura et al., 2012), (c) a mean velocity of 10 m yr −1 for debris-covered glaciers (Nakawo et al., 1999;Quincey et al., 2009), and (d) the total glacierized area in 2007 (410 km 2 ; ICIMOD, 2011). These tests gauge the ability of the model to accurately reproduce key glacier parameters: extent, mass change, and velocity. Scores are derived from the difference between modelled and observed quantities, with a score of zero indicating a perfect match. Scores for all four metrics are added to obtain an overall ranking of the 20 parameter sets and are weighted equally.
The glacier extent score denotes the relative deviation from a perfect match of the four large glacier termini at the end of the calibration period (Fig. 1). There are eight test polygons in total that include ice-covered and adjacent icefree areas. For example, if only 20 % of the glacier polygons in Fig. 1 are ice covered then the score equals 0.8. The mass balance score is based on the relative offset from the catchment mean mass balance of −0.40 m w.e. yr −1 over the period 1992-2008: If the modelled mean mass balance (B m ) equals −0.20 m w.e. yr −1 , then the mass balance score (S MB ) is 0.5. The total ice area score is based on the departure from the total glacierized area at the end of the simulation (410 km 2 , ICIMOD, 2011). If the simulated ice extent is 300 km 2 , then the score is 0.27 ((410-300)/410). Finally the flow velocity score quantifies the deviation from a mean glacier velocity of debris-covered tongues from 1992 to 2008 (10 m yr −1 ). For example, if the average simulated flow velocity is 2 m yr −1 , then the score is 0.8. The final score used to select the optimal parameter set is a simple addition of the four scores.

Model validation
Temperature and precipitation fields developed for this study were tested independently using point observations of mean daily temperature and total daily precipitation at the four EVK2NCR sites. We calculate mean bias error (MBE) and root mean square error (RMSE) to evaluate the skill of the elevation-based downscaling.

Glacier sensitivity to future climate change
To examine the sensitivity of modelled glaciers to future climate change, we drive the calibrated model with temperature and precipitation anomalies prescribed from eight CMIP5  climate simulations that represent cold/warm and dry/wet end-members (Table 5; Immerzeel et al., 2013). Decadal T and P anomalies relative to 1961-1990 are extracted from the CMIP5 end-members. Temperature trends are strong in all CMIP5 simulations, with ensemble mean temperature increases to 2100 as great as +8 • C in late winter and early spring (January-April). Precipitation anomalies do not show any significant trends and vary between 0.4 and 1.8 times the baseline period. Uncertainties in our scenarios of future climate change are examined through the mean and standard deviation of modelled ice areas and volumes derived from the eight CMIP5 models. As the model is empirically based and we assume only changes in T and P (all other state and input variables remain unchanged), we stress that the resulting glacier change realizations are a reflection of the modelled sensitivity to climate change, as opposed to physically based projections.

APHRODITE downscaling
Daily vertical temperature gradients calculated from the APHRODITE temperature fields and resampled SRTM range from −0.010 to −0.004 • C m −1 and are highly significant (Fig. 3). Calculated γ T values are most negative in the premonsoon (mid-April) and least negative during the active phase of the summer monsoon (mid-June to late August). This is likely a function of the increased moisture advection in the monsoon and pre-monsoon periods, which re-sults in a less negative moist adiabatic lapse rate. These findings are consistent with temperature gradient observations between −0.0046 • C m −1 (monsoon) and −0.0064 • C m −1 (pre-monsoon) in a nearby Himalayan catchment (Immerzeel et al., 2014b). The standard deviation in calculated γ T is lowest during the monsoon and greatest in the winter. At all four EVK2CNR stations, daily temperatures estimated from APHRODITE vertical gradients are greater than observed, with mean daily differences ranging from −1 to +8 • C (Fig. 4). Micro-meteorological conditions may contribute to the larger biases observed at Pyramid (winter) and Pheriche (summer). During the summer monsoon pewww.the-cryosphere.net/9/1105/2015/ The Cryosphere, 9, 1105-1128, 2015 riod (mid-June to mid-September), the mean difference for all stations is approximately 5 • C. We develop a bias correction for the day of year (DOY) based on the mean temperature bias from the four stations, which ranges from 3.22 to 6.00 • C. The largest bias coincides with the approximate onset of the summer monsoon (DOY 150, or 31 May). A possible mechanism for this is the pre-monsoon increase in humidity at lower elevations, which would be well-represented in the gridded APHRODITE data but not at the higher elevation EVK2CNR stations. The increased humidity would result in a less negative derived temperature gradient, and thus greater errors at the high-elevation stations. The variability in calculated temperature gradients is sharply reduced at onset of the monsoon, which supports this hypothesis. Biascorrected estimates of daily temperature (Fig. 9) have root mean squared errors (RMSE) of 1.21 to 2.07 • C and mean bias errors (MBE) of −0.87 to 0.63 • C. Based on the calculated daily temperature gradients, intercepts, and the bias correction, we estimate the height of the 0 • C isotherm (Z T =0 ) for the period 1961-2007 to examine melt potential and snow-line elevations. Mean monthly values of Z T =0 range from 3200 m (January) to 5800 m (July), though it can reach elevations of over 6500 m on occasion. This corresponds to meteorological observations from Langtang Valley, Nepal , and from the Khumbu Valley (http://www.the-cryosphere-discuss.net/ 7/C1879/2013/tcd-7-C1879-2013.pdf).
Daily precipitation-elevation functions (Fig. 5) exhibit strong decreases in precipitation above 4000 m, particularly in the monsoon and pre-monsoon periods. Absolute precipitation totals are greatest during the monsoon period, but large precipitation events can still occur in the post-monsoon period (October-November). As often observed in highelevation environments, daily precipitation totals observed at the EVK2CNR stations are not well captured by the downscaling process (Fig. 6). This is likely due to the difficulties in estimating precipitation in complex terrain Pellicciotti et al., 2012) and to errors in the precipitation measurements. For daily liquid precipitation (T > 0 • C), RMSEs range between 2.05 and 8.21 mm, while MBEs range from −0.85 to 1.77 mm. However, accumulated precipitation totals (Fig. 6) and mean monthly precipitation values show greater coherence, which lends some support for the downscaling approach used. At Pyramid (5035 m), the highest station with precipitation observations, the fit between cumulative predicted and observed precipitation is quite close. However, at Pheriche (4260 m), predicted precipitation is nearly double that observed over the period of record, which suggests that further refinements to the precipitation downscaling method are needed.

Model results and validation
For the calibration runs, we report here volume and area values averaged between 1 November and 31 January. Reported uncertainties are the standard deviation in modelled values from the 20 simulations. Modelled ice volumes from the 20 calibration runs (Fig. 10) decrease from 41.0 km 3 in 1961 to between 31.6 and 37.1 km 3 in 2007, with a 20-member mean of 34.5 ± 1.5 km 3 at the end of the simulation period. The ensemble mean modelled glacierized area in the calibration runs decreases from 499 km 2 to 392 ± 11 km 2 , with a final range of 374 to 397 km 2 .
Parameters for the calibrated model were chosen from Run 5, which had the lowest additive score of the 20 parameter sets (Table 4). Run 5 generates glacier volume and area totals that are lower but within 1 standard deviation of the model mean (Fig. 10). The selected parameter set contains degree-day factors ( Table 2) that are all slightly higher than those observed by Azam et al. (2014) at Chhota Shigri Glacier but are similar to values obtained for snow and ice by Singh et al. (2000) at Dokriani Glacier, Garhwal Himalaya. The value of the material roughness coefficient in the selected parameter set lies between the values used previously in Baltoro (Pakistan) and Langtang (Nepal, Fig. 1) catchments (Immerzeel et al., 2013, Supplementary Information).
Spatially distributed output from the calibrated model (Run 5), 1961-2007, is summarized in Fig. 11. Mean annual ablation (Fig. 11a) ranges from 0 to 4.00 m w.e. yr −1 , though most modelled values are less than 1.80 m w.e. yr −1 . Debris-covered termini, despite having lower degree-day factors, are nevertheless subjected to large melt rates due to their relatively low elevation and consequently higher temperatures. Our model generates maximum melt rates at the transition between debris-covered and clean glacier ice, at elevations of approximately 5000 m (Fig. 2). This is consistent with geodetic observations of mass change in the catchment (e.g. Bolch et al., 2008b). Maximum mean annual snowfall (Fig. 11b) amounts of up to 0.50 m w.e. yr −1 are observed at 6268 m (the calibrated value of Z C , Table 2), but due to the precipitation scaling function (Eq. 2) the highest peaks receive zero snowfall amounts. The calibrated height of Z C (6268 m) is similar to the elevation of maximum snowfall (between 6200 and 6300 m) estimated for the Annapurna range in mid-Nepal ( Fig. 1; Harper and Humphrey, 2003).
Modelled glacier velocities during the calibration period are less than 10 m yr −1 over debris-covered glacier termini and between 30 and 100 m yr −1 between the accumulation and ablation zones. While there are differences in both the spatial pattern and magnitude of modelled and observed velocities (e.g. Quincey et al., 2009), we feel that our simplification of glacier dynamics is unavoidable in the current study, and the development of higher-order physically based models will lead to improved representations of glacier flow.

Mass balance
Over the entire domain, modelled mean annual mass balances (b a ; Fig. 11c   range between −1.45 and +0.11 m w.e. (Fig. 13), with low variability amongst the different pa-rameter sets. Surface mass balance observations at the same site from 2007 to 2012 range between −0.67 and +0.46 m w.e. (Wagnon et al., 2013). As model and observation periods do not overlap, direct comparisons between modelled and observed mass balances are not possible. However, the mean mass balance observed at Mera Glacier between 2007 and 2012 is −0.08 m w.e., whereas the mean modelled mass balance between 2000 and 2006 is −0.16 m w.e. We note that our reconstructed mass balance series at Mera Glacier shows strong similarities to the reconstructed mass balance at Chhota Shigri Glacier (Azam et al., 2014), with balanced conditions in the late 1980s and early 1990s. Standard deviations of observed and modelled mass balance are 0.51 and 0.29 m w.e., respectively, and the greater variability in observed b a is likely linked to the short observation period (5 years) and to enhanced local variability which cannot be captured with downscaled climate fields. The mass balance model, although it may underestimate the The Cryosphere, 9, 1105-1128, 2015 www  1962 1967 1972 1977 1982 1987 1992 1997 2002  inter-annual variability, is able to simulate a surface mass balance that is in a plausible and realistic range.

Modelled and observed glacier thickness
At the end of the calibrated run , modelled ice thicknesses range between 0 and 620 m, though 98 % of these are less than 205 m (Fig. 11d). Similar ice thicknesses have been estimated for the large debris-covered Gangotri Glacier, Indian Himalaya, using slope, surface velocities, and simple flow laws (Gantayat et al., 2014). Due to the model formulation, low-angle slopes on glacier termini may result in unrealistic estimates of ice depth, and a minimum surface slope of 1.5 • is prescribed in the model. Radio-echo surveys in 1999 indicated that centerline ice thicknesses on the Khumbu Glacier decreased from approximately 400 m at Everest Base Camp to less than 100 m near the terminus (Gades et al., 2000). Our model accurately captures this decrease in the upper portions but overestimates ice thickness in the relatively flat terminus. Recent observations of ice thickness obtained from ground penetrating radar (GPR) surveys in the basin are examined in detail below. Estimates of glacier thickness extracted from the calibrated model and are compared with depth profiles found with GPR surveys conducted at Mera Glacier (Wagnon et al., 2013) and Changri Nup Glacier (C. Vincent, unpublished data). To facilitate the comparison, we obtained surface elevations and bedrock depths from the GPR surveys, and we matched these to the modelled ice thicknesses of the corresponding pixels (Fig. 14). At the lower elevation profile on Mera Glacier (5350 m), the shape of the bedrock profile is  (Fig. 1). Ice depth estimates for all 20 calibration runs are given in grey, and the results for Run 5 are shown as a dashed black line. similar to the model, but ice thicknesses are approximately half what is observed or less. This may be due in part to the surface slope extracted from the DEM, which controls the modelled ice thickness. The transect at 5350 m was collected in a flat section between two steeper slopes, which would likely be mapped as a steep slope in the DEM. For the profile at 5520 m both the shape and the depths of the bedrock profile are generally well-captured by the model. At the Changri Nup cross section, which lies on a relatively flat section of the main glacier body, modelled ice depths are approximately two-thirds of the observed. Modelled ice depths do not appear to be highly sensitive to the range of model parameters used in the 20 calibration runs, though variability is higher for Mera Glacier than for Changri Nup.

Modelled and observed glacier extents and shrinkage
Modelled historical changes in glacier area (Fig. 10) exhibit greater variability than modelled ice volumes. This is largely due to the sensitivity of the modelled glacier area to large snowfall events, as snowfall amounts greater than the 0.2 m w.e. threshold are classified as glacier. To compare modelled and observed extents we use the mean extent at the end of the ablation season (1 November-31 January). Using semi-automated classifications of Landsat imagery, glacier extents in the Dudh Koshi basin were constructed for 1990constructed for , 2000constructed for , and 2010constructed for (ICIMOD, 2011Bajracharya et al., 2014a, available at rds.icimod.org). As the glacier change signal is greatest at lower elevations, and errors in glacier delineation due to persistent snow cover are possible at higher elevations, we consider the change in glacier area below 5500 m, which roughly equals the equilibrium line altitude in the catchment.

Glacier sensitivity to future climate change
Decadal temperature and precipitation anomalies extracted from members of the CMIP5 ensemble that capture a range of climate scenarios (Table 5) are applied to the historical APHRODITE T and P fields. The calibrated glacier mass and redistribution model is then used to explore the sensitivity of modelled glaciers to future climate change in the Dudh Koshi basin. From initial glacier volumes and extents (Eq. 8), the mean projected changes in total ice volume at 2050 are −39.3 and −52.4 % for RCP4.5 and RCP8.5 emissions scenarios, respectively ( Table 6). The minimum projected volume change at 2050 is −26 % (cold/wet), and the maximum is −70 % (warm/dry). At 2100 the projected mean total volume loss is estimated at −83.7 % for RCP4.5 scenarios, and −94.7 % for RCP8.5, with a range between −70 and −99 %. Radić et al. (2014) and Marzeion et al. (2012), respectively, estimate mean glacier volume changes in south-east Asia of −50 and −60 % for RCP4.5 scenarios and −75 and −70 % for RCP8.5 by 2100. In all scenarios presented here, the rate of ice loss decreases towards the end of the simulation period (Fig. 16), which indicates a shift towards equilibrium mass balance conditions. Increased precipitation may slow the rate of future mass loss, but it is not sufficient to offset the increases in glacier melt due to increased temperatures. Changes in the timing and magnitude of monsoon precipitation may thus be less important than previously believed (Mölg et al., 2012;Bolch et al., 2012). The main difference between the RCP4.5 and RCP8.5 scenarios is the magnitude of the temperature increase, which leads to greater losses of ice volume in the RCP8.5 scenarios. This is due in part to the increased melt but also to the expansion of the ablation area and the change in precipitation phase from solid to liquid. Based on the daily temperature gradients and projected monthly temperature increases, the elevation of the 0 • C isotherm may increase by 800 to 1200 m by 2100. A potential snow-line elevation of 7000 m in August would expose 90 % of the current glacierized area to melt and severely restrict snow accumulation during the monsoon.
With a distributed model we can examine the possible impact of future climate change on Everest-region glacier area and thickness with respect to elevation. The patterns of decreases in ice area (Fig. 17) and ice thickness (Fig. 18) with elevation illustrate the combined effects of increased melt rates due to warmer temperatures and the insulating effect of debris cover. The greatest losses in glacier area, both relative and absolute, are expected at elevations close to the current ELA (approx. 5500 m), where the greatest amount of debris-free ice area currently exists. At lower elevations, where glaciers are exclusively debris-covered (Fig. 2), modelled glacier thicknesses are greater (Fig. 11), melt rates are lower, and modelled changes in glacier area and volume will be less than those near the ELA.
Wet and cool scenarios for both the RCP4.5 and RCP8.5 scenarios show the possible survival of debris-covered glaciers between 4000 and 4500 m, albeit with greatly reduced thicknesses (Fig. 18). In both warm and dry scenarios, glaciers below 5500 m could be eliminated, and in the RCP8.5 scenario, glacier thicknesses between 6000 and 6500 m could experience reductions by the year 2100. According to these scenarios, no changes are expected in the glacier volumes at elevations above 7000 m.

Discussion
Through a multi-parameter calibration and validation with independent data sets, we model the mass balance and mass redistribution of glaciers in the Dudh Koshi basin over the period 1961-2007. Temperature and precipitation changes specified from end-members of the CMIP5 ensemble are applied to historical climate fields to examine the sensitivity of glaciers in the region to future climate change. Expected increases in temperature will result in sustained mass losses that are only partially offset by increases in precipitation. We can identify three main sources of uncertainty in our approach: parametric, structural, and climate inputs. These are discussed below. Although considerable progress is made in this study by the systematic integration of field-based observations into our modelling approach, there are still a number of key challenges to be addressed in the future.

Structural uncertainty
The glacier mass balance and redistribution model used in this study has precedents in other studies  and has been calibrated here with observational data. While the model is a simplification of complex ice flow and dynamical processes, it is an important tool that can be used to explore the sensitivity of glaciers in the region Year Ice volume (% of initial) RCP4.5 RCP8.5 Figure 16. Sensitivity of modelled glacier volumes to decadal T and P anomalies from four RCP4.5 (blue) and four RCP8.5 (red) ensemble members (see Table 5 for details). Realizations are given as thin lines, and ensemble means are thick lines. All realizations are smoothed with a loess filter (span = 0.05) to minimize interannual variations. to future climate change. Given the forcings (−1.2 • C over 47 years) and parameter set (uncalibrated) used in the initialization, and the lag in actual glacier geometry response to climate change, it is possible that there are additional uncertainties in our estimates of initial ice volumes.
Our assumption of stationary debris cover may also be incorrect in the long-term, as glacier wastage typically leads to increased debris concentrations and the development of a debris cover. However, the median glacier slope above 5500 m is greater than 20 • (Fig. 7), and the development of debris cover on such slopes is unlikely (cf. Fig. 3b, Scherler et al., 2011a) as de-glaciation proceeds. Until higher-order models of glacier dynamics (e.g. Adhikari and Huybrechts, 2009;Clarke et al., 2015) are sufficiently advanced and explicitly include the effects of debris cover, and the additional input data (bedrock topography, ice temperatures) are wellconstrained, simple modelling approaches will still be required for basin-scale analyses of glacier change scenarios.

Parametric uncertainty
Our calibration approach relies on 20 sets of six different parameters with values taken randomly from pre-assigned initial values and ranges (Table 3). Model results from the 20 parameter sets (Figs. 12,13,14) suggest that the parametric uncertainty is well-constrained. The selected set of calibrated parameters is similar to those used in other regions , but a much larger and more computationally expensive Monte Carlo-type simulation must be undertaken to reduce the parametric uncertainty. Additional calibration data sets would also be beneficial, and these could include a greater number of ice depth measurements from debris-covered and clean-ice glaciers, remotely sensed snow cover, and glacier mass balance.

Input climate data uncertainty
The lack of high-elevation temperature and precipitation data to force the mass balance model is one of the key challenges that nearly all Himalayan modelling studies face. In this study, we derive temperature gradients and precipitationelevation functions from the 0.25 • gridded APHRODITE data, which in turn is based primarily on low-elevation stations. The downscaling approach is then tested with semiindependent station data from the EVK2CNR network of stations in the Dudh Koshi basin. While temperatures can be skillfully modelled after applying a bias correction based on the day of year, our ability to predict precipitation ranges from very good (at Pyramid) to very poor (at Pheriche). Difficulty in quantifying precipitation and precipitation gradients in high-mountain areas is likely one of the largest sources of uncertainty in mountain hydrology Nepal et al., 2013). Further investigations into high-elevation precipitation gradients, through field studies, remote sensing derivatives, and/or the use of high-resolution numerical weather models, will help to increase our understanding of glacier nourishment in the region. An analysis of the sensitivity of modelled glacier change to the rain/snow threshold temperature is also recommended.

Response times
Glaciers in the region are highly sensitive to temperature changes. Precipitation increases of 15 % (mostly during the monsoon season) will be unable to counter the loss of glacier mass due to increased melt rates. For intense warming scenarios, our ensemble mean volume change is more negative than regional estimates given by both Marzeion et al. (2012) and Radić et al. (2014). The potential loss of lower-elevation www.the-cryosphere.net/9/1105/2015/ The Cryosphere, 9, 1105-1128, 2015 glaciers in the study area raises the question of glacier response times. The actual response times of glaciers in the region can be approximated from modelled thicknesses and mass balance rates near the glacier terminus, following the methods of Jóhannesson et al. (1989): where H is a representative glacier thickness andḃ a (ḃ a > 0) is the mean annual mass balance near the terminus. Given our modelled ice thicknesses and mean annual mass balances at the termini of glaciers throughout the catchment, Eq. (10) suggests that the smaller glaciers in the southern portions of the basin have total glacier response times on the order of 20-50 years, while the large debris-covered glaciers have response times of 200-500 years. These first-order estimates reflect the time it takes the glaciers to reach a new equilibrium state in response to a step change in climate  and are in agreement with the modelled persistence of debris-covered termini and loss of smaller, lowelevation glaciers. Our scenarios suggest that future reductions in glacier area will occur mainly in clean ice regions between accumulation areas and debris-covered termini. We anticipate that the hypsometric distribution of ice will become bi-modal as glacier mass loss proceeds: debris-covered tongues will continue to exist (in reduced states) at low elevations but will become separated from their high-elevation accumulation zones (Kääb, 2005). Current examples of this type of glacier change can be found at Chorabari Glacier, Garhwal Himalaya (Dobhal et al., 2013), and at Lirung Glacier (central Nepal) in nearby Langtang Valley (Immerzeel et al., 2014a), where glacier wastage above the debris-covered termini has left stagnant debris-covered ice below and small high-elevation ice masses above. Model scenarios from this study are thus consistent with field observations and suggest that this will become a familiar picture in the coming decades.

Conclusions
In the mountains of high Asia, changes in glacier volumes will impact the timing and magnitude of streamflows, particularly in the pre-monsoon period . Our study advances the current understanding of Himalayan glacier evolution under climate change and examines the basin-scale evolution of glaciers in the Dudh Koshi basin of central Nepal using a distributed glacier mass balance and redistribution model. We constrain the glacier model parameters with observations where possible and calibrate against observations of net glacier mass change, velocities on debriscovered termini, and glacier extents. Our work represents a first-order estimate of future glacier change and is subject to considerable uncertainty from a number of sources.
Temperature and precipitation anomalies from endmember scenarios extracted from the CMIP5 RCP4.5 and RCP8.5 ensemble  are applied to historical downscaled climate fields, and the model is used to explore the sensitivity of glaciers in the Dudh Koshi basin to future climate change. Modelled glacier sensitivity to temperature change is high, with large decreases in ice thicknesses and extents for even the most conservative climate change scenario. Future climate scenarios with increased precipitation and reduced warming result in decreased mass losses, though increases in precipitation are insufficient to offset the dramatic increase in mass loss through increased melting.
Glaciers in the region appear to be highly sensitive to changes in temperature, and projected increases in precipitation are insufficient to offset the increased glacier melt. While we have identified numerous sources of uncertainty in the model, the signal of future glacier change in the region is clear and compelling. Advancements in the representation of ice dynamics (Clarke et al., 2015) and understanding of highaltitude precipitation will result in improved catchment-scale estimates of glacier sensitivity to future climate change in high mountain Asia.