Interactive comment on “ Future high-mountain hydrology : a new parameterization of glacier retreat ” by M . Huss

This important point was raised by both reviewers and in the additional interactive comment. In the revised version, we address the issue in detail and provide an additional subfigure (see revised manuscript). The figure shows individually derived parameterizations for the 34 glaciers divided into the three size classes, as well as the empirical functions given by the equations (see Fig. 3b), as solicited particularly by Reviewer #2. Additionally, standard deviations of the different parameterizations within the size classes are displayed in order to show their range of variability and the overlap with the other functions. In general, the variability of the indivdual ∆h-parameterizations within a size class is relatively small for the ’large valley glaciers’, but high for the ’small glaciers’, as these are more importantly affected by individual characteristics of glacier geometry and small scale effects. The overlap of the three size classes is considerable, but, as discussed as response to the next comment, the use of different parameterizations relative to glacier size instead of only one, is of benefit nevertheless.


Introduction
The retreat of mountain glaciers in response to ongoing climate change is expected to have major impacts on the water resources in alpine environments all over the world (e.g., Zierl and Bugmann, 2005;Horton et al., 2006;Hagg et al., 2006;Bradley, 2006;Stahl et al., 2008).Storage of fresh water in the seasonal snow cover and glacier ice is a key element in the water cycle on different temporal and spatial scales.It importantly affects issues of irrigation and hydropower production in mountainous regions and adjacent lowlands.In the long term the anticipated strong reduction in mountain glacier extent will lead to a shortage of water in the dry summer season (Zappa et al., 2003), and therefore have high economic impacts.There is a need for simple and widely applicable methods for calculating temporal and spatial changes in glacier coverage and ice volume in order to make reliable transient projections of future changes in the hydrology of alpine catchments.The assessment of glacier geometry change requires the description of (1) the surface mass balance reflecting the climatic forcing, and (2) the ice flow dynamics.
The climatic forcing acting on the glacier can be described using mass balance models of varying complexity that relate meteorological variables to accumulation and ablation rates at the glacier surface (e.g., Arnold et al., 1996;Hock, 1999;Klok and Oerlemans, 2002;Pellicciotti et al., 2005).The ice dynamics of alpine glaciers have been assessed in many glaciological studies ranging from simple flowline models (e.g., Greuell, 1992;Oerlemans, 1997;Sugiyama et al., 2007) to complex 3-D ice flow models (Hubbard et al., 1998;Gudmundsson, 1999;Jouvet et al., 2008).Changes in Alpine glacier extent over the next decades have been Published by Copernicus Publications on behalf of the European Geosciences Union.calculated using combined models of mass balance and ice dynamics (Vieli et al., 1997;Wallinga and van de Wal, 1998;Schneeberger et al., 2003;Le Meur et al., 2007).Ice flow models, however, require considerable field data input and computational resources and are, in general, not applicable for catchment or regional scale hydrological modelling.This paper elaborates on a method to parameterize the change in ice-covered area, glacier length and ice thickness that occurs in response to global warming.Whereas widely applicable methods to calculate snow and ice melt in hydrological models exist, the simulation of the distributed glacier surface elevation change in response to changing climate forcing and ice dynamics was not yet addressed in a simple way.
The h-parameterization, initially proposed by Huss et al. (2008b), is a function relating the elevation of the glacier surface h to the surface elevation change h (equivalent to ice thickness change) occurring over a given time interval.Typically, elevation changes are small in the accumulation area and largest near the terminus of mountain glaciers (Fig. 1) as it is shown by observations (Finsterwalder and Rentsch, 1981;Arendt et al., 2002;Bauder et al., 2007), based on theoretical considerations (Jóhannesson et al., 1989) and numerical modelling (e.g., Oerlemans, 1997).For application of the h-parameterization in practice, the position along the central flowline of the glacier (as shown on the abscissa of Fig. 1b) is approximated with the surface elevation, assuming these variables to be strongly correlated.The parameterization takes into account the conservation of mass and is thus well suited for transient hydrological modelling of ice storage changes.For two glaciers in the Swiss Alps, differing in size and geometry, we assess the performance of the h-parameterization in calculating future glacier geometry over the 21st century by comparing the results to a 3-D finite element ice flow model.The changes in glacier area and consequent shifts in the discharge regime are projected for the period 2008 to 2100 based on different climate scenarios.

Study site and field data
The two study sites are located in different regions of the Swiss Alps (Fig. 2).We focus on (i) Rhonegletscher, a medium-sized valley glacier with a length of about 8 km and an area of 15.9 km 2 in 2007, and (ii) Silvrettagletscher, Hydrol.Earth Syst.Sci., 14, 815-829, 2010 www.hydrol-earth-syst-sci.net/14/815/2010/ a relatively small mountain glacier (2.8 km 2 ).For both glaciers, a variety of field data is available covering different temporal and spatial scales which makes them suitable for investigation of the proposed parameterization.
Ice volume changes for subdecadal to semicentennial intervals are available over the 20th century based on digital elevation models (DEMs) derived either from topographic maps (before 1960) or from aerial photographs (Bauder et al., 2007).Ice thickness was measured using radio-echo sounding in 2003and 2008on Rhonegletscher, and in 2007on Silvrettagletscher (Farinotti et al., 2009a).Bedrock maps were generated using an interpolation scheme accounting for the principles of ice flow dynamics (Farinotti et al., 2009a) which is applied between radio-echo sounding profiles.
Measurements of annual mass balance at stakes were carried out over several decades on both glaciers (Mercanton, 1916;Funk, 1985;Glaciological reports, 1881Glaciological reports, -2009)).Additionally, snow probings in high spatial density were performed in six accumulation seasons on each glacier (unpublished data), providing information about the snow distribution pattern over the glacier surface.Daily discharge records at a gauging station at Gletsch from the catchment of Rhonegletscher (glacierization in 2007: 47%) are available for more than five decades.
Meteorological data (daily mean air temperature and precipitation sums) are recorded at several stations close to the glaciers.For an overview of the methods used to obtain continuous meteorological time series representative for the glacier sites over the 20th century refer to Huss et al. (2008a).

h-parameterization
The h-parameterization describes the spatial distribution of the glacier surface elevation change in response to a change in mass balance.A glacier in a perfect state of equilibrium (not existing in reality) shows no ice thickness change over time because all spatial differences in surface mass balance are compensated by the ice flow dynamics.A disequilibrium of mass balance and ice flow leads to distinct spatial patterns of the glacier surface elevation change.Its distribution depends on the size, the geometry, the ice flow regime and the mass balance variability of the glacier.Consequently, the shape of the h-parameterization differs from glacier to glacier.
We discuss two possibilities for the derivation of the parameterization.The first one is suited for glaciers with observations of surface elevation change in the past; the second one is applicable to unmeasured glaciers.

Parameterizations for individual glaciers
If at least two glacier surface DEMs are available, the hparameterization is derived directly from observed surface elevation change patterns in the past.By comparing the el-evation change between two DEMs in elevation bands with an arbitrary interval (here 10 m), a distribution of surface elevation h versus the elevation change h is obtained.The hh function needs to be smoothed in order to exclude noise which would be amplified in the modelling.This applies in particular to the uppermost parts of the accumulation area, where the uncertainty in the DEMs is relatively high and the area covered by individual elevation bands is small (Bauder et al., 2007).The DEMs should be apart for a sufficiently long time period to show elevation changes that are significantly higher than the DEM uncertainty.
Given a monotonic long-term trend in glacier evolution, the quality of the h versus h relation generally increases with the time span covered and the magnitude of changes occurring during this period.The observed distribution of elevation changes is, however, affected by long-term variations in glacier mass balance that are evident in the 20th century for the European Alps (Vincent, 2002;Zemp et al., 2009;Huss et al., 2010a).For example, the distribution of h over the glacier surface is inherently different for periods of fast glacier wastage (e.g.since the 1990s), during time intervals of positive mass balance (e.g. the 1970s) or following such a period (e.g. the 1980s) due to the delayed response of glacier geometry to a change in climatic forcing.Over the next decades mostly negative mass balances are expected (e.g., Schneeberger et al., 2003).Thus, it is proposed to derive the h-h function for periods between successive DEMs characterized by persistent glacier mass loss.The h-h function is normalized in order to obtain the dimensionless h-parameterization (Fig. 3a): the abscissa is normalized with the elevation range of the glacier according to (h max − h)/(h max − h min ), where h max/min are maximum and minimum glacier surface elevation; and the ordinate is normalized with the maximum elevation change h max as h/ h max .

Parameterizations for unmeasured glaciers
As information about past glacier geometry changes is often unavailable, generalized h-parameterizations were retrieved from samples of different glacier size classes in the Swiss Alps in order to yield empirical functions applicable to unmeasured glaciers.For 34 glaciers for which repeated DEMs over the 20th century are available (Bauder et al., 2007;Huss et al., 2010b) glacier-specific hparameterizations were determined (Fig. 3c).The functions were averaged in three size classes with n members -large valley glaciers (area A > 20 km 2 , n=8), medium-sized valley glaciers (5 km 2 < A < 20 km 2 , n=12) and small glaciers (A < 5 km 2 , n=14).Most parameterizations were derived for the period between the 1930s and present; some functions only refer to the last two decades.For valley glaciers minor ice thickness changes occur over a large part of the elevation range (see Fig. 1).On small glaciers, in contrast, ice thickness changes are not restricted to regions near the   glacier terminus and the elevation dependence of h is more uniform over the glacier (Fig. 3b and c).The parameterizations for the three size classes were approximated based on least squares using empirical functions of the form where h is the normalized surface elevation change and h r the normalized elevation range.The exponent γ , prescribing the curvature of the h-function, decreases with glacier size (see numerical values in Fig. 3b).This is explained with the higher importance of ice flow on large glaciers, as well as their wider elevation range.The deviations of individual parameterizations from the mean of the size class is considerable especially for small glaciers, as these are more importantly affected by small scale characteristics of glacier geometry change.The overlap of parameterizations for different size classes decreases towards the ablation area, where the largest changes occur (Fig. 3c).
The general form of the elevation change signatures is assumed to prevail for all mountain glaciers, also including climate conditions different from Europe (e.g.Himalayas), as they are determined by the universal factors of some altitude dependence of mass balance, and gravitational ice flow.Therefore, the applicability of the h-parameterization outside of the European Alps is given, however, requires a recalibration based on repeated DEMs for very different glacier types (e.g.debris-covered glaciers, ice caps with outlet glaciers).

Implementation
Once the h-parameterization is established, its application for the calculation of future glacier extent requires the following field data input, which can be potentially obtained for a large number of glaciers, also in remote regions of the world, from readily available data sets (e.g.DEMs from the Shuttle Radar Topography Mission (SRTM); glacier inventories): 1. a DEM and glacier outlines for model initialization, 2. the spatial distribution of surface mass balance calculated using a model driven by climate variables (e.g., Klok and Oerlemans, 2002;Machguth et al., 2006;Huss et al., 2008a), and 3. an estimate of the glacier bedrock topography.This information is required to update the glacier extent and determines the size of the ice reservoir.Ice thickness distribution and total ice volume can be calculated based on an inversion of surface topography (Farinotti et al., 2009a,b).Simpler methods (Haeberli and Hoelzle, 1995;Bahr et al., 1997) can be applied to estimate the overall ice volume and, thus, to obtain a first guess of the mean ice thickness.
The h-parameterization is applied to update the 3-D glacier geometry at the end of each hydrological year.The lost or gained ice volume, calculated using a mass balance model, is converted into a distributed ice thickness change according to the parameterization.It is assumed that the redistribution of the annual surface accumulation by ice flow is instantaneous.Integration of the dimensionless h-function over the normalized elevation range h r of the glacier, multiplied with the area A (in m 2 ) of each elevation band i and the ice density ρ ice must equal the total annual change in glacier mass B a : (2) B a (in kg) is given by the mass balance computation.f s (in m) is a factor that scales the magnitude of the dimensionless ice thickness change (ordinate in Fig. 3) and is chosen for each year such that Eq. ( 2) is satisfied.h r refers to the annually updated glacier extent.An updated glacier surface elevation h 1 for the initial altitude h 0 of elevation band i is calculated as Surface lowering is restricted not to exceed that year's surface mass balance in ice equivalent at the glacier terminus.This is important for years with strongly negative mass balance, when the scaling of the h-parameterization (see Eq. 2) would result in unrealistically high surface lowering on the glacier tongue.This corresponds to non-dynamic downwasting of the glacier tongue, currently observed on many alpine valley glaciers (Paul et al., 2004).In this case, the parameterization is rescaled (Eq.2) so that its integration yields the total annual volume change.The hparameterization is not applied near the border of the glacier (i.e.where the ice thickness is smaller than 10 m).It is assumed that in these regions thinning due to ice dynamics is small.Ice thickness change is calculated based on the local mass balance in ice equivalent.Glacier extent is determined by updating the ice thickness distribution; the glacier disappears where ice thickness drops to zero.

Glacier Evolution Runoff Model (GERM)
In this paper, the framework of glacier surface updating based on the h-parameterization is the glacio-hydrological model GERM.The model is designed to calculate the runoff from highly glacierized drainage basins and includes components for computing snow accumulation distribution, snow and ice melt based on the temperature-index approach, evaporation, and runoff routing.All components of GERM are described in Huss et al. (2008b).
The model is forced by daily temperature and precipitation, and is calibrated and validated using different field data.For Rhonegletscher and Silvrettagletscher the glacier module of GERM is tuned so that the results match the observed (i) changes in ice volume, (ii) in-situ point-based annual mass balance, and (iii) in-situ winter accumulation.We use the parameter set calibrated in Huss et al. (2008a) for both glaciers.Runoff calculated by GERM for the Rhonegletscher catchment is validated against long-term daily discharge records.
The glacio-hydrological model reproduces observed annual runoff volumes with an average deviation of 1.4% over the period 1957-2006 (Fig. 4a).The Nash and Sutcliffe (1970)criterion is R 2 =0.58 for annual, and R 2 =0.96 for monthly runoff.The model thus simulates the seasonality of runoff reasonably well (Fig. 4b) and is consistent with long-term ice volume changes and seasonal glacier mass balance variability.

Ice flow model
A three-dimensional finite element ice flow model coupled to a surface mass balance model is used to simulate the flow dynamics of Rhonegletscher.The model is based on a Volume of Fluid (VOF) formulation (Scardovelli and Zaleski, 1999) and solves the nonlinear Stokes equations.It has been successfully applied to different types of fluids (Maronnier et al., 2003) and was first employed by Jouvet et al. (2008) to calculate glacier dynamics.Model application to Rhonegletscher is described in detail by Jouvet et al. (2009).
Mass conservation along the ice-air interface yields a transport equation which can be used to determine the evolution of glacier geometry.Given the initial shape of the glacier, the ice flow velocity u is obtained by solving a 3-D nonlinear Stokes problem derived from a regularized version of Glen's law with a temperature independent rate factor (Jouvet et al., 2008(Jouvet et al., , 2009)).A nonlinear sliding law is prescribed along a given part of the ice-bedrock interface.
Both rate factor and sliding coefficient minimize the rootmean-square error of the difference between observed (see Mercanton, 1916) and simulated surface velocities.Then, glacier geometry is updated each year by computing the volume fraction of ice ϕ, which satisfies the transport equation where bδ A is a source term (i.e.surface mass balance) acting on the ice-air interface A .A decoupling algorithm allows the equations to be solved using different numerical techniques (Jouvet et al., 2009).The nonlinear Stokes problem is solved on a fixed, unstructured finite element mesh consisting of tetrahedrons, whereas the transport equation is solved using a fixed, regular grid of smaller cells.
The ice flow model is implemented with the same scheme for calculation of glacier surface mass balance as GERM, using the parameters calibrated with field measurements.GERM and the ice flow model are driven using identical meteorological time series in daily resolution.
For Rhone-and Silvrettagletscher results of the combined ice flow and mass balance model are validated against repeated observations of glacier geometry.The model is initial-  ized employing the glacier surface DEM for 1874 (Rhone), and 1938 (Silvretta), respectively.Over the 20th century, changes in glacier area and surface elevation along a longitudinal profile are reproduced by the model for both glaciers (Fig. 5).Model performance is assessed by evaluating the root-mean-square (rms) error of the difference between simulated and observed glacier surface elevation at all grid cells for six DEMs for each glacier.We find an average rms error of 18 m for Rhonegletscher, and 11 m for Silvrettagletscher.The errors remain similar over the entire modelling period.These values are regarded as the a priori uncertainty in the 3-D ice flow model in the present configuration.Disagreement of modelled and measured surface elevation is attributed to (1) the description of ice dynamics in the model, (2) the uncertainty in the bedrock elevation (the deviations are largest in regions not covered by radio-echo sounding measurements), and (3) uncertainties in the calculated mass balance distribution.

Future climate forcing
Our assumptions on future climate are based on Frei (2007), who statistically evaluated the results of 16 regional climate models within the PRUDENCE project (Christensen and Christensen, 2007).Three scenarios for the seasonal change in air temperature and precipitation are defined according to Huss et al. (2007).Scenario 2 represents the median (most probable), Scenario 1 a glacier-friendly evolution (cold-wet) and Scenario 3 is most unfavourable to glacier existence (warm-dry) (Fig. 6).The three scenarios envelope a 95% confidence interval of expected future climate change in the Alps and thus cover the range of uncertainty in climate system evolution as given by different regional climate models.Expected changes in climate variables in 2050 are given for the three scenarios in Table 1.Daily meteorological time series for the period 2008 to 2100 are constructed by applying seasonal changes in air temperature and precipitation (Table 1) to detrended series recorded in the past.

h-parameterization versus ice flow model
As the h-parameterization is derived from observed changes in glacier surface elevation over the 20th century, a validation is not possible for this period.The parameterization aims at simplifying glacier retreat calculations for the next decades.Therefore, we compare the performance of the h-parameterization to the coupled surface mass balance and 3-D ice flow model that is run into the future using realistic climate change scenarios.Our model comparison does not imply that the predicted future glacier extent is "correct", but it allows the investigation of the appropriateness of the glacier retreat parameterization versus the ice flow model.The latter is assumed to provide the best possible estimate of future 3-D evolution of glacier geometry.
The ice flow model and GERM, which is implemented with the h-parameterization for calculating the distributed ice thickness change, are forced using the same climatic input for the period 2008 to 2100.For both glaciers an upto-date DEM (glacier geometry and outline) is available for 2007 serving as initialization for the future model runs.We apply the specific h-parameterization for Rhonegletscher (Fig. 3a) and simulate the evolution of Silvrettagletscher using the approximation (Eq. 1) derived for the sample of small alpine glaciers (Fig. 3b).Mass balance is calculated in daily resolution for every grid cell applying the parameters of the glacier module calibrated for the past.Model runs were performed for the Scenarios 1 to 3. Because the identical scheme for mass balance calculation is used to drive both the ice flow model and GERM, differences in simulated surface elevation and glacier extent are attributable to the description of ice flow dynamics only; the ice flow model is thus regarded as the reference result.
www.hydrol-earth-syst-sci.net/14/815/2010/ Hydrol.Earth Syst.Sci., 14, 815-829, 2010  Table 2. Comparison of surface elevation calculated using the h-parameterization with the ice flow model for three 30-year (Rhonegletscher) and 20-year periods (Silvrettagletscher), respectively, and three climate scenarios (Sc1-3).rms d represents rootmean-square errors of the elevation difference at all grid cells averaged over the selected period.According to our results, large alpine valley glaciers, such as Rhonegletscher, will be likely to lose up to 90% of their ice volume until 2100.The spatial distribution of the icecovered area as obtained from the ice flow model agrees with results of the h-parameterization and the distributed surface elevation change is captured over almost the entire 21st century (Fig. 7).Towards 2100 surface elevation near the glacier terminus is underestimated by the parameterization.This may partly be related to a mass balance-altitude feed-back that tends to continuously amplify an initial difference in calculated surface elevation.

Period
The rms error of surface elevation obtained using the parameterization in comparison to the ice flow model is around 20 m throughout the entire modelling period (Rhonegletscher) and is comparable for all three climate scenarios (Table 2).Good agreement is also found for Silvrettagletscher, for which the generalized h-parameterization was used (Fig. 3b).Average rms errors are between 10 and 20 m and are consistent as well for all scenarios (Table 2).Lower rms errors in comparison to Rhonegletscher are explained by smaller ice thicknesses.These results indicate that the h-parameterization is able to approximate glacier surface elevation changes as predicted by a sophisticated 3-D ice flow model for a wide range of possible future climate evolutions and for different glacier types.Misfits are in the same order as the a priori uncertainty in the ice flow model inferred for the 20th century.
In Fig. 8 the comparison of the glacier surface elevation calculated using the h-parameterization to the ice flow model is shown on a longitudinal profile.For Rhonegletscher the parameterization yields particularly satisfying results until 2050.The distribution of the surface elevation change along the central flow line is almost perfectly reproduced for all scenarios.Towards the end of the modelling period larger misfits occur in the ablation area.These are most pronounced for Scenario 3 (warm-dry).Scenario 1 (cold-wet) shows better agreement for the entire 21st century (Fig. 8).
For this scenario ice flow dynamics are similar as in the period for which the h-parameterization was derived.This might explain the good performance.Future changes in area and length of Rhonegletscher are captured (Fig. 9).Ice-covered areas agree with each other within 3% throughout the entire modelling period.For Scenario 2 (median) the mean misfit in glacier area is only 1.5%.Also glacier length is reproduced within a difference of 3% (Fig. 9).Towards 2100 larger errors are evident for Scenario 3.
The parameterization performs equally well for the smaller Silvrettagletscher.In general, it predicts slightly too fast glacier surface lowering in the ablation area and too little thinning in the accumulation area (Fig. 10).This is related to the fact that the h-parameterization used was not derived specifically for Silvrettagletscher and thus ignores its particular characteristics.Nevertheless, ice-covered area and glacier length calculated either using the parameterization or the ice flow model agree well (Fig. 11).The mean misfit of glacier area is <5% for Scenarios 1 and 3 and <1% for Scenario 2. The time period with simulated ice coverage of <20% of the initial glacier area was excluded in this evaluation.Glacier length is reproduced within 10%.The complete disappearance of Silvrettagletscher during the 21st century is very likely.Whereas according to Scenario 2 no ice is left around 2070, the glacier survives in strongly reduced form only according to Scenario 1 (Fig. 11).

Future runoff
Runoff from the drainage basin of Rhonegletscher is simulated for the 21st century using GERM implemented with the h-parameterization.In order to take into account yearto-year climate variability we perform multiple model runs using randomly disturbed meteorological series (Huss et al., 2008b).Due to the representation of glaciers as a dynamically changing storage component in the water cycle, the shift from an ice melt to a snow melt dominated runoff regime can be simulated transiently.Our results are in line with previous studies (Braun et al., 2000;Zierl and Bugmann, 2005;Stahl et al., 2008;Huss et al., 2008b) and indicate that this transition will take place until the end of the 21st century.During this time frame, water resource management will have to adapt to the major changes in the hydrological regime of high alpine catchments.Over the next decades annual runoff from the Rhonegletscher catchment is expected to increase considerably (Fig. 12).For 2040 calculated runoff is between 23% and 36% higher than in 1961-1990, depending on the scenario.This is due to a rapid release of water from long-term glacial storage, because of increasing temperatures and consequently strong glacier retreat.Whereas annual runoff volume reaches a tipping point before 2100 in Scenario 2 (median) and 3 (warm-dry), the model results indicate steadily increasing runoff for Scenario 1 (cold-wet) due to relatively constant reduction of the glacier ice volume and enhanced precipitation at the same time (Table 1 and Fig. 12).
The application of the ice flow model leads to a slightly different glacier evolution relative to the hparameterization (see e.g.Fig. 8).We evaluate the resulting differences in annual runoff from the catchment of Rhonegletscher in decadal intervals by comparing simulated glacier ice volume.Averaged over the century, the misfits are negligible (0.1% of annual runoff).Higher differences with both positive and negative sign occur for individual decades.The rms error of annual runoff obtained using the h-parameterization or based on the ice flow model over all scenarios and the entire study period is 3.2%.
Our results show the gradual transition from a glacier melt to a snow melt dominated runoff hydrograph (Fig. 13).The transient increase in runoff over the next four decades is mainly concentrated in the summer months.Strong ice melt in combination with precipitation events also leads to an increased potential for destructive floods in this period.This is due to the significantly reduced storage capacity of glaciers that have only marginal snow and firn coverage.
After complete disappearance of the ice coverage, lowflow conditions in the summer months will prevail.They will be even intensified by increased potential evaporation in a warmer climate.On the scale of the European Alps this is expected to have regional (Horton et al., 2006) to transnational impacts (Bradley, 2006;Junghans et al., 2010) on irrigation and agriculture, the biosphere and shipping traffic on major streams.The glacier discharge peak is shifted from early August to end of May.For the Rhonegletscher catchment, we find a decrease in August runoff by 55% for the year 2100 according to Scenario 2 compared to the period 1961-1990 (Fig. 13), and even by 94% for Scenario 3. Conversely, our results show a doubling of winter runoff (Scenario 2), and a more than fivefold increase for Scenario 3. Results of Scenario 1 indicate significantly smaller changes in the hydrograph.

Discussion
In order to investigate the suitability of alternative simple approaches to calculate the change in glacier surface area in response to future climate change we compare two supplementary parameterizations to the ice flow model for Rhonegletscher.We consider (1) the "AAR-method", and (2) "nondynamic downwasting".
In previous studies assessing the glacier retreat in the 21st century a simple method based on the accumulation area ratio (AAR) was applied (Horton et al., 2006;Schaefli et al., 2007;Paul et al., 2007).It is assumed that the glacierized area A gl is proportional to the multi-year mean accumulation area A acc , which is calculated based on climate variables.A gl is derived using the relation where AAR s is the accumulation area ratio required to yield a balanced mass budget.We assume AAR s =0.6 (Schaefli et al., 2007).The accumulation area A acc is calculated as the mean of the previous 10 years in order to average out yearto-year meteorological variability.An updated glacier area A gl is determined annually using Eq. ( 5).In this approach it is implicitly assumed that the glacier is always in balance with the prevailing climate conditions, which is not justified for large ice masses in a changing climate.A major drawback of the "AAR-method" is that it does not account for the conservation of mass.For transient hydrological modelling, in particular, this represents a serious problem as ice volume can disappear without contributing to discharge.Given the mass balance and the ice thickness distribution, the change in glacier area can also be calculated by rigorously assuming non-dynamic downwasting.Ice thickness change at each grid cell equals the annual mass balance in ice equivalent.This approach completely neglects a transport of mass from the accumulation to the ablation area.Thus, the lowering in the ablation area is expected to be too fast and surface elevation is supposed to increase at high altitudes.An advantage of this approach compared to the "AAR-method" is that ice volume is transiently modelled, i.e. mass is conserved.
Model runs using both supplementary parameterizations were performed for climate Scenarios 1-3 based on the identical meteorological input and the same scheme to calculate glacier mass balance as for the ice flow model.Table 3 shows that the rms error of surface elevation obtained with both supplementary parameterizations is significantly higher than for the h-parameterization (see Table 2 for comparison).The "AAR-method" yields a considerable underestimation of glacier area for the entire modelling period.rms errors of surface elevation compared to the ice flow model are in the range of 60 to 120 m for Rhonegletscher (Table 3).For predicting glacier area over a time step that is much longer than the glacier response time (Hoelzle et al., 2003) the "AARmethod" may yield suitable results.
Assuming no redistribution of snow and ice accumulation by ice flow leads to a monotonically increasing misfit relative to the ice flow model (Table 3).Whereas rms errors are lower than for the "AAR-method" in the first decades of the 21st century, "non-dynamic downwasting" performs worse with increasing length of the modelling period.This is explained by unhindered and continuous growth of the ice thickness in the accumulation area and too fast glacier recession in the ablation area.
Although misfits in calculated glacier surface elevation are comparable between the "AAR-method" and "non-dynamic downwasting" (Table 3), the differences in reproducing glacier area and length as simulated by the ice flow model are significant (Fig. 14).According to the "AAR-method" Rhonegletscher almost immediately loses its entire valley glacier tongue; averaged over the modelling period glacier area is underestimated by 44%, and calculated glacier length is lower by 66% compared to the ice flow model (Scenario 2).Most Alpine glaciers are presently out of equilibrium, i.e. their size is too large for the current climate conditions.The "AAR-method" does not account for the time scale glaciers need to lose their excess ice mass and thus predicts too fast glacier retreat.The misfits are largest around 2050 when the glacier is far from equilibrium and become smaller as the ice-covered area decreases (Fig. 14).
Assuming "non-dynamic downwasting", a relatively good match with glacier area and length given by the ice flow model is achieved (Fig. 14), being only slightly worse than for the h-parameterization."Non-dynamic downwasting" yields the best results for Scenario 3 that is based on extreme climate change leading to rapid deceleration of ice flow and consequent disintegration of the glacier.
The testing of supplementary parameterizations for calculating future glacier retreat shows that they yield results inferior to those of the h-parameterization, which is able to reproduce distributed glacier surface elevation change, as well as glacier area and length accurately based on the same field data input.The use of the "AAR-method" for transient projections of future glacier extent is discouraged as volume change time scales are neglected; the violation of mass conservation results in strong underestimation of runoff over the next decades (see Huss et al., 2008b)."Non-dynamic downwasting" leads to an unrealistic distribution of glacier surface elevation, however, might be a suitable alternative to the h-parameterization for small glaciers or in the case of very rapid climate change.
The major source of uncertainty for the prediction of future glacier area and consequent changes in the water balance and the runoff regime of high alpine drainage basins are the changes in the climate system, and how they act on the snow coverage and the glacier surface.State-of-the-art climate models indicate a wide range of possible climate trends in the future (IPCC, 2007).The spread is mainly due to unknown future emission of anthropogenic greenhouse gases.However, different climate models also exhibit diverging results when forcing them with the same CO 2 input (e.g., Christensen and Christensen, 2007).The use of ensembles, allowing a probability distribution of the climate model results to be determined is therefore highly important.
With climate change several feedback mechanisms that are only partly incorporated in our mass balance model could contribute to a bias in the projections of glacier wastage and future runoff.Whereas we take into account positive feedbacks of increasing fraction of precipitation occurring in liquid form, and the prolongation of the melting season, additional positive backcoupling mechanisms as darker glacier tongues due to dust deposition (Oerlemans et al., 2009) and lakes forming in front of glacier termini, e.g. at Rhonegletscher (Farinotti et al., 2009a), are not considered.We do not account for feedbacks that locally reduce ice melt rates as well, such as increasingly debris-covered glacier tongues (Kellerer-Pirklbauer et al., 2008), and decreased direct radiation as the glacier lowers into deeply incised and shaded valley bottoms.Because of the various poorly understood processes that may compensate for each other, it is difficult to provide a joint uncertainty estimate of errors in future projections due to shortcomings of the glacier surface mass balance modelling procedure.
Investigation of the changes in climatic forcing at the elevation of the long-term equilibrium line of Alpine glaciers based on long-term measurement has shown that the degreeday factors of temperature-index models have changed significantly over the 20th century (Huss et al., 2009).Changes over short time periods are also evident from the analysis of energy balance measurements (Braithwaite, 1995;Carenzo et al., 2009).The use of constant parameters in empirical models calibrated in the past therefore adds to the uncertainty in simulated glacier extent over the 21st century.This uncertainty does, however, not compromise the comparison of the h-parameterization with the ice flow model as both are driven by the same mass balance calculation scheme.
Our results show that the h-parameterization reproduces results given by the application of a complex ice flow model, and that the parameterization is thus suitable for the calculation of glacier retreat in response to future climate change.However, as the h-parameterization represents a rigorous simplification of three-dimensional glacier geometry changes, some restrictions are discussed.
The h-parameterization does not perform well in simulating glacier advance.In theory, the pattern of the altitudinal distribution of the ice thickness changes only differs by sign for glacier retreat and advance (Jóhannesson et al., 1989).However, this is only true after reaching a steady state.Alpine glaciers take one to several decades to reach a new equilibrium geometry after a step like change in climate (Hoelzle et al., 2003).The parameterization is applied at the end of every year, assuming an immediate transformation of the local mass change into a distributed surface elevation change over the entire glacier.As the transient changes in surface elevation are more complex for glacier advance than for retreat (see e.g., Finsterwalder and Rentsch, 1981), and the h-parameterization is explicitly derived for periods of persistently negative mass balances, the short-term glacier geometry change in the case of advance cannot be simulated accurately.For this task we recommend the use of an ice flow model.
Changes in the shape of the h-parameterization over time are likely.With climate change glaciers diminish in area which might lead to a lower curvature of the optimal h-parameterization, i.e. to a lower exponent γ in Eq. (1).However, given the good performance of the h-parameterizations over a time interval of 90 years (see e.g.Figs. 8 and 10) the uncertainty induced by these factors is assumed to be acceptable.
In order to test the spatial transferability of the generalized h-parameterizations (see Fig. 3), we performed two crossvalidation experiments based on the functions derived for the 34 glaciers individually.(1) For calculating an average h-function for each size class, one glacier is omitted; the mean rms error of the omitted function and the size class parameterization is evaluated.This is repeated for all items in the size class.(2) As in (1), one glacier is omitted, but the averaged function used for comparison is based on glaciers in arbitrary size classes (same number of functions as in (1), randomly chosen).We find that, overall, the rms errors are significantly higher (T-test at the 95% level) in experiment (2) compared to (1), indicating that the use of the size class specific parameterization is of benefit.For one fifth of the analyzed glaciers, however, the division into size classes is not an advantage.

Conclusions
In this study a simple parameterization for calculating the future change in ice-covered area, glacier volume and 3-D geometry is proposed and compared to results obtained from a state-of-the-art ice flow model coupled to a surface mass balance model.The evolution of Rhonegletscher and Silvrettagletscher was simulated for the period 2008-2100 using three different climate change scenarios based on ensemble evaluation of 16 regional climate models.We find a dramatic retreat for both glaciers.Whereas the smaller Silvrettagletscher is likely to disappear around 2070, Rhonegletscher most probably still exists at the end of the 21st century, however, with a projected ice volume loss of more than 90% compared to today.Transient simulation of runoff from a highly glacierized catchment quantitatively shows that glacier wastage will importantly impact on the hydrological regime.
The h-parameterization is able to reproduce glacier retreat over the 21st century given by a complex 3-D ice flow model and shows satisfying results for a wide range of scenarios of future climate change and different glacier types.Input data requirements are limited to distributed glacier surface mass balance in annual resolution, and the initial ice thickness distribution.These data can be obtained from often readily available data sets, such as glacier inventory data and land surface DEMs with a wide spatial coverage (e.g.SRTM).The parameterization proposed in this paper is thus applicable to any mountain glacier.The hparameterization is computationally cheap and therefore represents a viable alternative to complex ice flow modelling.It offers the possibility to calculate the transient evolution of large glacier samples in regional impact studies.
This study underlines that climate change will have major impacts on alpine hydrology and the water resource management of mountainous regions.Glacier retreat and the release of freshwater from long-term glacial storage is expected to be a key element in projections of high alpine runoff over the next decades, and requires a description in hydrological modelling that meets the current level of glaciological research.As the application of ice flow models is normally not in the scope of regionally distributed hydrological models, this necessitates simple alternatives.Based on the good performance of the h-parameterization in reproducing (i) distributed glacier surface elevation, (ii) glacier area and (iii) length given by a 3-D ice flow model, we recommend this parameterization for the approximation of future glacier change for projections of stream-flow runoff from glacierized drainage basins over the 21st century.

)Fig. 1 .
Fig. 1.(a) Observed long-term surface elevation changes in a longitudinal glacier profile for a large alpine valley glacier (Grosser Aletschgletscher).(b) Distribution of ice thickness change along the glacier in 1957-1999.Changes in ice thickness are small in the accumulation area and largest near the glacier terminus.The variability in the equilibrium line altitude (ELA) over the 20th century is shown.

Fig. 2 .Fig. 2 .
Fig. 2. Study site overview.(a) Location of the glaciers in Switzerland.(b) Hydrological drainage Rhonegletscher (dashed).The gauging station at Gletsch is marked.The contour line interval is 100 colours indicate ice thickness.(c) Map of Silvrettagletscher.Note that the scale differs from (b).
-parameterization derived for Rhonegletscher from observed ice thickness changes in the 20th The variability in the equilibrium line altitude (ELA) is shown.(b) Empirical ∆h-parameterizations glacier size classes applicable to unmeasured glaciers derived from digital elevation model comparison aciers.The equations refer to a numerical approximation (according to Eq. 1) of the line for each size ) Approximations (see (b) for equations) including individually derived parameterizations for all 34 (thin lines).The colour indicates the size class.Error bars show the uncertainty of the approximation; calculated as the standard deviation of the individual parameterizations within the size class. 23

Fig. 3 .
Fig. 3. (a) h-parameterization derived for Rhonegletscher from observed ice thickness changes in the 20th century.The variability in the equilibrium line altitude (ELA) is shown.(b) Empirical hparameterizations for three glacier size classes applicable to unmeasured glaciers derived from digital elevation model comparison for 34 glaciers.The equations refer to a numerical approximation (according to Eq. 1) of the line for each size class.(c) Approximations (see (b) for equations) including individually derived parameterizations for all 34 glaciers (thin lines).The colour indicates the size class.Error bars show the uncertainty of the approximation; they are calculated as the standard deviation of the individual parameterizations within the size class.

Fig. 5 .
Fig. 5. Validation of the coupled ice flow and surface mass balance model over the 20th century using field data.Comparison of observed and calculated glacier surface on a longitudinal profile in 1959 and 2007.Model runs were initialized in 1874 (Rhone) and in 1938 (Silvretta).

Fig. 6 .Fig. 7 .Fig. 6 .
Fig. 6.Deviations of annual (a) mean air temperature and (b) precipitation sums from the period 198 Measured data for the 20th century are displayed by bars, annual trends as assumed by Scenarios 1-3 ( are shown until 2100.

Fig. 8 .Fig. 9 .
Fig. 8. Evolution of glacier surface elevation along a longitudinal profile calculated for Rhonegletscher using the 3D ice flow model and the ∆h-parameterization.Results for different climate change scenarios are shown.

Fig. 8 .
Fig. 8. Evolution of glacier surface elevation along a longitudinal profile calculated for Rhonegletscher using the 3-D ice flow model and the h-parameterization.Results for different climate change scenarios are shown.

Fig. 9 .
Fig. 9. Calculated (a) glacier area and (b) length based on the ice flow model and the h-parameterization according to the three climate scenarios (Rhonegletscher).

Fig. 12 .
Fig. 12. Calculated annual runoff volume from the Rhonegletscher catchment based on three climate scenarios.Shaded areas indicate the range of variability (±2 standard deviations σ ) obtained by performing multiple model runs with random meteorological variability around the expected average of each scenario (solid lines).The dot shows the measured mean annual runoff in the climatic normal period, ±2σ are indicated.

Fig. 13 .Fig. 13 .
Fig. 13.Projected change in the hydrological regime of the Rhonegletscher catchment according to S 2 (median).Measured discharge in the past refers to the 1961-1990 average.Numbers in brackets ind glacierization of the drainage basin.

Fig. 14 .Fig. 14 .
Fig. 14.Comparison of simple alternative approaches for calculating future glacier change with the i model.(a) Glacier area, and (b) length (Rhonegletscher).

Table 3 .
Comparison of surface elevation calculated using supplementary simple parameterizations to calculate glacier retreat with the ice flow model for three 30-year periods and scenarios in the future (Rhonegletscher).The root-mean-square errors rms d of the difference in elevation at all grid cells averaged over the period considered.