Technical Note : Representing glacier dynamics in a semi-distributed hydrological model

Glaciers play an important role in high-mountain hydrology. While changing glacier areas are considered of highest importance for the understanding of future changes in runoff, glaciers are often only poorly represented in hydrological models. Most importantly, the direct coupling between the 10 simulated glacier mass balances and changing glacier areas needs feasible solutions. The use of a complex glacier model is often not possible due to data and computational limitations. The hparameterization is a simple approach to consider the spatial variation of glacier thickness and area changes. Here, we describe a conceptual implementation of the h-parameterization into the semidistributed hydrological model HBV-light, which also allows for the representation of glacier advance 15 phases, and comparison between the different versions of the implementation. The coupled glaciohydrological simulation approach, which could also be implemented in many other semi-distributed hydrological models, is illustrated based on an example application.


Introduction
Glacier melt water is an important contribution to discharge in high-mountain catchments (Köplin et al., 2013;Miller et al., 2012) and can sustain summer streamflow in many large river basins (Hagg et al., 2007;Stahl et al., 2017).When modeling the hydrology of such catchments for longer periods (>10 years), the changing glacier area has to be considered, especially when climate change is causing glacier 25 retreat.The simplest approach is to update the hydrological model with an externally simulated glacier extent, but this is unsatisfactory, as the mass balance as simulated by the hydrological model might not agree with the updated glacier extent.The use of coupled glacio-hydrological models allows the glacier extent to be linked directly to the simulated glacier mass balance and is, thus, better suited for modeling catchments with changing glacier areas (Huss et al., 2008;Stahl et al., 2008).However, modelers are 30 faced with the question of which degree of complexity is needed to represent glaciers and glacier evolution in hydrological models.The use of a fully distributed, physically-based glacier model that considers mass balance, subglacial drainage, ice flow dynamics etc. is often too data demanding and Hydrol.Earth Syst.Sci. Discuss., doi:10.5194/hess-2017-158, 2017 Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 30 March 2017 c Author(s) 2017.CC-BY 3.0 License.computationally expensive for use in a combined glacio-hydrological model where the non-glacierized part of the catchment also has to be considered.Furthermore, many semi-distributed hydrological models use simplified representations of catchment hydrology using a limited number of conceptual buckets (reservoirs) and coupling such a model with a more complex glacier model would lead to a mismatch in degree of physical and spatial representation.Hence, for hydrological modeling studies 5 there is a need for glacier models that use a similar degree of complexity and data demand as other components of the hydrological model but which are still able to represent the important glacier processes.
Recently an increasing number of hydrological models have incorporated glacier evolution models, using for example an equilibrium line altitude (ELA) shift (e.g., Linsbauer et al., 2012), volume-area 10 (V-A) scaling (e.g., Radić et al., 2008), V-A scaling and morphological image analysis (e.g., Stahl et al., 2008), other simple schemes without ice flow (e.g., Bongio et al., 2016), or more complex approaches focusing on glacier modeling (e.g., Immerzeel et al., 2012).One approach with limited glacier input data requirements, which is mass-conserving and well suited for hydrological modeling studies, is the h-parameterization, which describes the spatial distribution of the relative change in glacier surface 15 elevation in response to a change in mass balance (Huss et al., 2010).Initially, Huss et al. (2008) introduced the h-parameterization as part of their Glacier Evolution Runoff Model (GERM), while a more detailed presentation of the approach, including the derivation of generalized empirical functions applicable to unmeasured glaciers, is given in Huss et al. (2010).Since then, the h-parameterization has been applied in global scale modeling by Huss and Hock (2015) as well as in numerous studies 20 applying GERM to simulate individual glaciers or glacierized regions in the Swiss Alps (Farinotti et al., 2012;Finger et al., 2013;Gabbi et al., 2012;Huss et al., 2014;Huss and Fischer, 2016) and in Central Asia (Sorg et al., 2014).Several other glacio-hydrological models were coupled with glacier retreat simulations following the h-approach (Addor et al., 2014;Gabbi et al., 2014;Linsbauer et al., 2013;Ragettli et al., 2013;Salzmann et al., 2012;Vincent et al., 2014).However, details on its practical 25 implementation into the respective conceptual hydrological models have been provided by only few studies, for instance those by Li et al. (2015) and Duethmann et al. (2015).
As the h-parameterization is based on an empirical approximation to describe glacier retreat, it is subject to uncertainty and several limitations in terms of accurate glaciological modeling at the scale of individual glaciers (discussed in Huss et al., 2010;Linsbauer et al., 2013;Vincent et al., 2014).30 Nevertheless, for the purpose of transient hydrological modeling, particularly for regional studies covering large samples of glacierized catchments, the h-approach represents an efficient state-of-theart alternative to more complex glacier evolution models (Huss et al., 2010;Li et al., 2015).Originally, Huss et al. (2010)  with indications for a presence of recent glacier advance (Ragettli et al., 2013).Moreover, it represents a major drawback for long-term hydrological modeling covering past periods, for example with positive mass balance in the European Alps during the 1970s.A simplified scheme to incorporate short-term glacier change in case of advance as extension of the original h-approach is presented by Huss and Hock (2015).5 Here, we describe a conceptual implementation of the h-parameterization into the semi-distributed hydrological model HBV-light (Seibert and Vis, 2012), which also allows the representation of glacier advance phases, and compare different versions of the implementation.The described approach could also be implemented in many other semi-distributed hydrological models.

New glacier routine in the HBV model
The HBV model is a semi-distributed conceptual precipitation-runoff model.It has continued to be developed in Scandinavia since the 1970s (Bergström, 1976;Lindström et al., 1997) and has become a standard tool which is widely used in different model variants, particularly for modeling snowdominated catchments.Required input data are daily temperature, precipitation and potential 15 evaporation time series.Hydrological processes within a catchment are modeled by four different routines, a snow/glacier routine, a soil moisture routine, a response routine and, finally, a streamflow routing routine.Here, we describe the recent integration of a glacier evolution approach into the HBVlight software, a user-friendly and freely available version of HBV (Seibert and Vis, 2012).

Ice accumulation, melt and runoff
We first describe how ice accumulation, melt and glacial runoff are represented in HBV-light, before we describe how areal changes of the glaciers are represented.The glacier area within a catchment is conceptually simulated by two reservoirs representing glacier ice and the liquid water contained within the glacier.There can be a snowpack on top of the glacier, which also consists of a solid (snow) and a 25 liquid (water content) reservoir.The snow and glacier routine calculations are performed at each simulation time step for each elevation zone, for which elevation bands of 100 to 200 m are typically used.The elevation zones can be further subdivided according to three aspect classes (N, S, and W/E).
Depending on the temperature in relation to the threshold temperature, precipitation falls either as snow or rain.In case of rain, the precipitation is added to the water content of the snow if a snow layer is 30 present or to the water content of the glacier otherwise.If the temperature is above the threshold temperature, melt takes place in the snowpack based on a degree-day factor, and the melted snow is added to the water content of the snowpack.In the case that the water content exceeds the snow water holding capacity, the amount exceeding the snow water holding capacity flows out and is added to the For ice melt of the glacier a degree-day method is used as well, but ice melt is only simulated at times when there is no snow layer on the glacier.For temperatures above the threshold temperature, glacier melt is calculated using the degree-day factor multiplied with a glacier correction factor, which 5 represents the different albedo of ice compared to snow and typically takes values of about 1 to 2. The ice melt is added to the liquid component of the glacier, from which the outflow is computed individually for each elevation zone as suggested by Stahl et. al. (2008), extending earlier concepts by Moore (1993) to account for the enlargement of glacial conduits over the melt season.() = ()   +   •  −  •  ()  (1) 10 Q is the outflow, S the liquid water content of the glacier, S we the water equivalent of the snowpack, K min and K range are the minimum outflow coefficient and maximum range of outflow coefficient values, and A G a calibration parameter.To represent the transition from snow to ice in a simple way, at the end of each time step a user-specified fraction of the snow on top of the glacier is converted into ice and equally distributed over the whole glacier area.Typical values for this model parameter are 0.001-0.00315 which implies that the conversion of snow to ice on average takes about 1 to 3 years.
Optionally snow redistribution can be applied at the end of each time step to avoid unrealistic multiyear snow accumulation, the so-called "snow towers".During the snow redistribution, the snow layer (i.e.snowpack and snow water content) of all non-glacier areas above a certain user-specified (upper) elevation, and after reaching a certain user specified threshold, is redistributed over the non-glacier and 20 glacier areas in between a user-specified lower and upper elevation, as well as the glacier areas above the user-specified upper elevation.

Glacier mass balance and area changes
To translate glacier mass balance changes into area changes in the new Glacier Area Change Routine 25 (GACR), a single-valued relation between glacier mass balance and glacier area needs to be established.This relationship, which is technically represented in the model by a lookup table, can be derived from any glaciological model.Here we suggest that the relationship (and lookup table) is computed based on the ∆h-parameterization method described in Huss et al. (2010) and shown in Figure 1.The use of a lookup table allows for periods of glacier advance (though not further than the initial glacier extent as 30 defined by the user in the "glacier profile" file).While the areal distribution of a static glacier is specified in HBV-light by means of elevation and aspect zones, for establishing the relationship between glacier area and mass balance a glacier profile is needed, which defines the initial thickness and areal distribution of the glacier at a finer resolution, typically for elevation steps of 10 m.Note that the resolution of the dynamics of the glacier routine largely depends on the number of elevation 35 intervals per elevation zone, i.e. all glacier area within each interval is either covered by a glacier or not without any gradual transition, and the percentage of glacierized area within a certain elevation zone is based on the state of the individual elevation intervals within that elevation zone.However, as described below, we also used a width scaling to allow for a reduction of glacier area of the individual elevation intervals with decreasing thickness.
For the generation of the lookup table, the total glacier volume is defined based on the definition of the 5 initial glacier profile: (2) with  the total glacier volume in mm water equivalent relative to the entire catchment area, and for each elevation band i, the (dimensionless) area a i and water equivalent h i .To generate the lookup table, the glacier is melted in steps of  corresponding to one percent of the total glacier water equivalent 10 .During each of these steps the ∆h-parameterization method of Huss et. al. (2010) is applied.For each elevation band, the corresponding elevation is normalized according to: (3) where E i,norm is the normalized elevation corresponding to the absolute elevation E i of elevation band i, and E max and E min are the maximum and minimum elevations of the glacier.The normalized water 15 equivalent change is then computed for each of the normalized elevations using the following function (Huss et al., 2010): (4) where is the normalized (dimensionless) ice thickness change of elevation band .For the initial glacier profile definition, the total glacier area in km 2 is calculated and, based on this area, one of the 20 three empirical parameterizations applicable for unmeasured glaciers from Huss et. al. (2010) are used (Figure 1).
In the next step a scaling factor  s (mm), which scales the dimensionless ∆h, is computed based on the glacier volume change of 1%, and on the area and normalized water equivalent change for each of the elevation bands: 25 (5) The new water equivalent is then computed for each elevation band as (6) If Δh i exceeds h i,old for certain elevation bands (most likely to occur at the glacier tongue), the glacier thickness is reduced to zero and the remaining portion of Δh i is included in the next iteration step (i.e. 30 the next 1% melt).Once the new water equivalent values have been computed for each elevation band, the glacier area is updated for each elevation zone.The relative glacier area for a certain elevation zone is defined as the cumulative area of the glacier covered elevation bands within that elevation zone, divided by the total area of the elevation zone.However, glaciers have an uneven distribution of ice at a particular elevation with a thinner ice layer along the edges.Therefore, a simplified representation of the 3-D glacier 5 geometry is used to scale the area within a certain elevation band (Eq.7) following the relation between glacier width and glacier thickness given in Bahr et al. (1997) as also applied by Huss and Hock (2015): The above calculations define the relationship between glacier area and mass balance and are stored in a lookup table at steps of 1% of the initial glacier mass.This means that the lookup table consists of 10 glacier areas per elevation zone for 101 different glacier mass balance situations ranging from the initial glacier mass to zero (Fig. 1).
During the actual simulation in HBV-light, the glacier extent is updated after the end of each water year (1.October).The total water equivalent of the glacier is computed.Based on the percentage of glacier 15 water equivalent in comparison to the total glacier water equivalent in the glacier profile definition, the corresponding record is extracted from the glacier lookup table and the corresponding glacier areas are applied to the different elevation zones.In the case that the glacier water equivalent exceeds its maximum, the areas corresponding to 100% are applied (i.e. the glacier can never grow larger than defined by the user in the glacier profile definition).Simulations can start, however, with a reduced 20 glacier size, by specifying the initial glacier fraction in the glacier profile file (as fraction of water equivalent).The glacier profile definition should thus contain the maximum extent of the glacier during the full simulation period.For each glacierized part of an elevation zone in HBV-light, the corresponding non-glacierized part is used to exchange the area for which the state changed from glacier to non-glacier and vice versa.In order to ensure the water balance is correct, bookkeeping is 25 done between the corresponding glacierized and non-glacierized zones.Soil moisture and snow, for example, are moved between the corresponding zones as far as it corresponds to the area exchanged.

Sensitivity test of different model variants
To illustrate the new Glacier Area Change Routine (GACR) and its different components on the simulation results, we calibrated the HBV-light model for one example catchment, the Alpbach 30 catchment in the Swiss Alps.This catchment is one of the glacierized headwater catchments in the Rhine River basin, located in Central Switzerland.The catchment area is about 21 km² and elevations range from 1022 m up to 3192 m a.s.l. with a mean elevation of 2194 m.The catchment contains the glacier Glatt Firn.According to the glacier inventory for the year 2010 the glacierized area amounts to about 4 km², corresponding to a catchment glacier coverage of about 20%.35 Only the static part of the glacier routine is used, i.e. the complete dynamic part of the glacier 5 routine is disabled; the glacier area is not adjusted but stays exactly as defined by the user in the model setup during the whole simulation.

2) Full new GACR (GACR)
The full version of the model as described in section 2.1, with the static and dynamic part of the glacier routine included.10

3) GACR without glacier width scaling (GACR-w)
The application of glacier width scaling (Eq.7) by elevation band is disabled.In practice, this corresponds to a 2-D representation of glacier area change.A change in glacier area is realized only when the mean glacier water equivalent of an elevation band (Eq.6) reaches zero.As a result, the area change will mostly occur at the glacier terminus.15

4) GACR without glacier advanceglacier retreat only (GACR-a)
The original method described by Huss et. al. (2010) considers only the parameterization of glacier retreat and not glacier advance.In the new GACR, glacier advance up to the initial state is enabled by means of the lookup table generation.To demonstrate the effect of neglecting temporary glacier advance, we used a version that applies only glacier retreat.In periods with a positive annual glacier 20 mass balance the glacier area is kept constant.
For each of these four versions, we calibrated the model 10 times, using a genetic algorithm (Seibert, 2000) with 3500 model runs per calibration trial.The 10 independent calibration trials allowed parameter uncertainty effects to be considered by taking several optimized parameter sets into account.25 The simulation period was January 1 st , 1901 to December 31 st , 2006, and was preceded by a three year spin-up period.As an objective function, the average of the Nash-Sutcliffe efficiency for daily discharge, the relative volume error of the total discharge, root mean squared error of the snow cover simulations and absolute mean relative error of the glacier volume estimates were used.The estimates of observed glacier volume were based on different glacier cover datasets for three particular years 30 during the simulation period as described below.
The simulation period  is a period in which glaciers of the European Alps retreated considerably, yet this period also covers diverse climate conditions including a period between the 1960s and the 1980s that was characterized by rather balanced conditions or temporarily by glacier 35 advance.For the setup and the calibration of the model in terms of glacier conditions, several To setup the HBV-light model for the Alpbach catchment, the spatial modeling units were discretized as follows: Firstly, the glacierized and non-glacierized catchment area fractions for the state at simulation 15 start 1901 were distinguished.Therefore all areas within the Alpbach catchment that were glacier covered according to the underlying map or glacier inventory for a specific year were summed up as one model glacier.Both the non-glacierized and the glacierized model areas were further divided into area fractions per elevation zones, and then further differentiated within each elevation zone into area fractions for three aspect classes.20 For the application of the h-parametrisation, in addition to the main model setup an initial distribution of glacier thickness/volume over surface elevation was defined (Figure 1b).In the present implementation of HBV-light this thickness distribution is referred to as "glacier profile" and glacier thickness is expressed in mm water equivalent per elevation bands, with an interval that is smaller than the elevation zone interval of the overall catchment model.Here, an elevation band interval of 10 m was 25 applied for the initial glacier profile, whereas an interval of 100 m had been used for the setup of the whole catchment model.As no data on glacier thickness for the state at 1901 (start of simulation) was available, an initial glacier profile was reconstructed based on the different glacier area datasets, the glacier thickness data available for later years, and physically-based glacier scaling relationships (Bahr et al., 1997).For the final conversion of glacier ice thickness to water equivalents used in the modeling, 30 an ice density of 0.9 kg•l -1 was applied.

Results
Figure 2a shows the observed glacier profile for the initial state at simulation start in the year 1901 as well as for the three different years for which data was available from which the glacier profile could be Hydrol.Earth Syst. Sci. Discuss., doi:10.5194/hess-2017-158, 2017 Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 30 March 2017 c Author(s) 2017.CC-BY 3.0 License.derived.The observed decrease of glacier area occurred at all elevations.Figure 2b shows the glacier profile for the simulation with the full new GACR model version.With this version, glacier retreat also occured at all elevations.This is due to the h-approach, which distributes the change in glacier mass balance over the different elevation zones, in combination with the implemented width scaling, which relates a decrease in glacier thickness to a reduction of the glacier area.In order to compare the 5 simulated and observed glacier profiles, Figure 2c shows the differences between simulated and observed glacier area for the different elevation zones.The h-approach by definition results in zero change in glacier thickness at the very top of the glacier.The lower the position on the glacier, the larger the change in thickness can be.This pattern can be seen in Figure 2b, where there is, contrary to the observed data of Figure 2a, hardly any change in glacier area in the higher elevation zones.As a 10 result, the difference between simulated and observed areas in Figure 2c is positive for the higher elevation zones (for the years 1973 and 2003).This is compensated for by a negative difference between simulated and observed glacier areas for the lower elevation zones.Overall, the new GACR is able to depict the major pattern of long-term glacier area change over the elevation zones in the Alpbach catchment.15 The simulations using the full GACR also correspond well with the observations in terms of total catchment glacier area (Fig. 3).In addition, Figure 3 illustrates the differences in the four different model versions to simulate the changes in total catchment glacier area (Fig. 3a) and the resulting effects on the change in glacier water equivalent and cumulative runoff (Fig. 3b and c), which are relevant within the scope of hydrological modeling.Among all model versions the new full GACR is best in 20 representing the pattern of change in total glacier area based on the comparison with available observation based data (Fig. 3a).Whereas there is a considerable mismatch of the simulated and observed glacier area around the year 1940, for the later years the simulated and observed glacier areas are in good agreement.The model version that does not incorporate glacier advance is just as effective in meeting the final state of the glacier area in the year 2003 as the full version.In terms of glacier area 25 the results of both versions, GACR-a and GACR, are only different during phases dominated by positive glacier mass balance.As soon as the annual glacier water equivalent (glacier volume) decreases to its previous minimum again, the reduction in glacier area continues.For glaciers with a net negative mass balance over time, differences can therefore be rather small.However, if there are more and longer periods of glacier advance, differences might become more apparent.However, in case of overall net 30 positive glacier mass balance, the fact that the maximum glacier extent cannot exceed what is specified in the glacier profile file becomes an obvious limitation of the new GACR routine.For the version GACR-w the glacier stays at its maximum size a bit longer than for the full new GACR and the version GACR-a, since elevation bands need to be melted completely before the glacier area starts to reduce.In contrast, in the full new GACR and GACR-a simulations width scaling is applied as soon as the glacier 35 mass balance becomes slightly negative, and therefore a reduction in glacier area can be observed immediately.For simulations with only the static glacier routine the glacier area stays constant (horizontal grey line in Figure 3a).
The constant area allows for (much) higher melt rates in comparison to the other model versions once the glacier has partly melted, since a larger area is contributing to the overall melt than in the version including a GACR.This can also be clearly observed in Figure 3b, where the model version with a 5 stationary glacier area shows much stronger dynamics in glacier water equivalent change.As a result the cumulative glacier runoff (Fig. 3c) is highest for simulations with stationary glacier area, especially during the second half of the simulation period when the difference in glacier area in comparison to the other versions is more notable.Generally, the larger the glacier area, the more discharge is generated by the glacier.The stationary glacier area model (no GACR) results in the potentially largest amount of 10 glacier runoff, followed by the simulations without width scaling, for which the 10 different model calibrations resulted in the largest spread.The difference between the versions GACR and the GACR-a is minor, with the latter likely resulting in an underestimation of generated glacier runoff, due to the smaller area during phases of glacier growth.

Discussion and conclusion 15
The glaciological part of the coupled model as described above is a simple representation of glacier processes, but it allows glacier dynamics to be considered at a level of complexity which is similar to the hydrological model.Compared to an infinitely thick glacier for which the area does not change, the approach described here, which allows for area changes as a result of simulated mass balance changes, is certainly a more realistic representation and the changing area clearly affects variables such as the 20 simulated runoff.A major simplification is that only one glacier is considered in each subcatchment, which means that if there are several glaciers these are simulated as one virtually aggregated glacier.
The ∆h-parameterization approach of Huss et al. (2010) was found to be suitable and utilizes existing empirical functions to distribute changes along the glacier.This reduces the need for new calibration parameters.While the ∆h-parametrization could also be based on data for specific glacier(s) (as done, 25 for instance, by Duethmann et al., 2015), which would better represent local conditions, in reality such data is rarely available.A re-evaluation of the empirical ∆h-parameterizations would be most useful for applications beyond the European Alps, where the data for the empirical ∆h -parameterizations by Huss et al. (2010) came from.
Several adaptions were needed for the implementation in a semi-distributed model.Most importantly 30 the use of a lookup table to represent the mass-area relationship allows for advancing glaciers.
Furthermore, the geometric width scaling for individual elevation bands allows for the representation of a decreasing glacier area with decreasing thickness.The example simulations shown in this technical note illustrate the effect of these modifications, which maintain the conceptual model approach.The simulations demonstrate that the new glacier evolution routine is, in general, capable of reproducing 35 Hydrol.Earth Syst.Sci. Discuss., doi:10.5194/hess-2017-158, 2017 Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 30 March 2017 c Author(s) 2017.CC-BY 3.0 License.observed area changes.However, given the limited data this should not be taken as a proof that the model is correct, even if the simulations seem glacio-hydrologically reasonable.Allowing for advancing glaciers and changing areas due to glacier thinning makes a difference in the simulations (Fig. 2).Both these aspects are also important as they enable a comparison between simulated and observed glacier area (or also glacier length) (see Fig. 2a and Fig. 1).This is crucial for model calibration and validation 5 as glacier areas and glacier lengths are much more frequently available than other glacier observations.Besides its simplicity, the presented approach also has other limitations.As in all dynamic glacier modeling approaches, one challenge is to obtain initial thickness distributions along the glacier.
Furthermore, glacier advance is only possible up to the initial state.In most cases this is not a major limitation as long as suitable information on early glacier extents are available as most climate data and 10 scenarios lead to retreating glaciers.If needed, one could provide a larger initial glacier extent (with some thickness profile) to establish the mass-area relation to create the lookup table.In this case the actual simulations would start at a certain fraction of this hypothetical maximum situation.
We present a conceptually stringent approach, which allows changing glacier areas to be considered in an approximate but realistic way.This approach could in principle also be implemented into other 15 hydrological models.A full validation is challenging due to limited data in the considered catchment, as it is the case in most catchments, but the simulations of the full GACR agree reasonably well with available observations.Nevertheless, a validation based on the comparison with observations alone does not clearly determine how much better the full GACR is compared to the other tested versions.Still, based on glacio-hydrological reasoning, all the different components of the GACR are important.20 In many hydrological model applications of partially glacierized catchments that do not specifically target the contributions of glaciers to runoff, glacier areas are not directly updated and studies with a coupled glacio-hydrological approach often describe little details of the glacier routine, especially when it comes to the question on whether simulated glacial mass balances are translated into glacier area changes and, if so, how this is done.In a recent review on hydrological modeling of glacierized 25 catchments in central Asia (Chen et al., 2016), for instance, this issue is not discussed at all.The main advantage of the coupled glacio-hydrological approach as described in this technical note is that glacier mass balance and area changes are consistent with the hydrological model.This also allows the model to be used to simulate future scenarios.While the approach described in this technical note is a rather simple representation of glacier processes, it enables this important representation of basic behavior of 30 glaciers in high mountain catchments including changing glacier areas.

Figure captions
derived the h-parameterization for periods dominated by negative mass balances and glacier retreat.The missing representation of glacier advance is related to uncertainties in regions 35 Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2017-158,2017 Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 30 March 2017 c Author(s) 2017.CC-BY 3.0 License.
Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2017-158,2017   Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 30 March 2017 c Author(s) 2017.CC-BY 3.0 License.liquid water reservoir of the glacier.If the temperature is below the threshold temperature, part of the water content in the snow layer refreezes.
Syst.Sci.Discuss., doi:10.5194/hess-2017-158,2017   Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 30 March 2017 c Author(s) 2017.CC-BY 3.0 License.To demonstrate the effect of the different parts of the GACR, four different versions of the GACR were used, where for each version certain components of the new glacier routine were disabled:1) Stationary glacier area (no GACR) Hydrol.Earth Syst.Sci.Discuss., doi:10.5194/hess-2017-158,2017   Manuscript under review for journal Hydrol.Earth Syst.Sci. Discussion started: 30 March 2017 c Author(s) 2017.CC-BY 3.0 License.observation based datasets from diverse sources were used: the glacier area for the state around the years of 1901 (start of simulation period) and 1940 was based on digitized historical topographic maps, in Switzerland known as "Siegfriedkarte".For both years, 1901 and 1940, the area of the Alpbach catchment is reconstructed from two adjacent map sheets.To describe the glacier area around the simulation start in 1901, maps from the years 1894 and 1899 were used, and to describe the glacier area 5 around 1940, maps from the years 1933 and 1942 were used.Glacier areas for the years 1973, 2003, and 2010 were extracted from the glacier inventories by Müller et al. (1976), Paul et al. (2011) and Fischer et al.(2014), respectively.For the years 1973 and 2010 gridded datasets of estimated glacier thickness based on the method presented in Huss and Farinotti (2012) were also used (unpublished data provided by Matthias Huss).In addition, streamflow observations (station Erstfeld, Bodenberg, period 10 1960-2006) from the Swiss Federal Agency of the Environment (FOEN) and a gridded snow water equivalent (SWE) climatology product from the Institute of Snow and Avalanche Research (WSL-SLF, covering Nov-May for the period 1972-2006) were used to calibrate the model.

Figure 1 :
Figure 1: ∆h-parameterization and its implementation in HBV-light: a) ∆h-parameterization functions for three glacier size classes from Huss et al. (2010), b/c) pre-simulation application of the ∆hparameterization for a medium glacier size to the example glacier profile data of the Alpbach catchment by melt in steps of ∆M =1%M to generate the lookup table.b) absolute glacier volume per elevation 5 band, c) relative glacier area per elevation band as relative fraction of the initial glacier area a initial of the elevation interval.For each elevation interval, the resulting glacier water equivalent/glacier volume (b) and glacier area (c) are shown as grey lines, for visibility, only results of steps of ∆M =5% are shown here.The initial profile (100% M) and profiles for a glacier volume reduction by 20, 40, 60, and 80% are highlighted by 10 colored labeled lines.

Figure 2 : 20 Figure 3 :FiguresFigure 1
Figure 2: Observed and simulated glacier areas per elevation zone for years when observations are available: a) Glacier areas for the different elevations derived from maps or remote sensing, b) Glacier 15 areas for the different elevations as simulated with the full new GACR model version, c) Diffrences between observed and simulated glacier areas.

Figure 2 :
Figure 2: Observed and simulated glacier areas per elevation zone for years when observations are available: a) Glacier areas for the different elevations derived from maps or remote sensing, b) Glacier areas for the different elevations as simulated with the full new GACR model version, c) Diffrences between observed and simulated glacier areas.5

Figure 3 :
Figure 3: Comparison of the simulations by the different versions of the glacier routine: a) Total glacier area, b) Change in glacier storage, c) Cumulative glacier runoff.The range of simulation results 5 represents the results from 10 model (equally suitable) parameterizations for each of the different versions of the glacier area change routine (GACR) and for the version without glacier area change routine.