{'fraction_non_alphanumeric': 0.02873201438848921, 'fraction_numerical': 0.018300359712230214, 'mean_word_length': 4.7567730802415875, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Representing glaciers in a regional climate model

A glacier parameterization scheme has been developed and implemented into the regional climate model REMO. The new scheme interactively simulates the mass balance as well as changes of the areal extent of glaciers on a subgrid scale. The temporal evolution and the general magnitude of the simulated glacier mass balance in the European Alps are in good accordance with observations for the period 1958–1980, but the strong mass loss towards the end of the twentieth century is systematically underestimated. The simulated decrease of glacier area in the Alps between 1958 and 2003 ranges from −17.1 to −23.6%. The results indicate that observed glacier mass balances can be approximately reproduced within a regional climate model based on simplified concepts of glacier-climate interaction. However, realistic results can only be achieved by explicitly accounting for the subgrid variability of atmospheric parameters within a climate model grid box.


Introduction
The twentieth century has been a period of worldwide glacier retreat. This phenomenon is consistent with the observed global temperature increase and is rated as one of the most reliable natural indicators of ongoing warming trends (Haeberli 2005;IPCC 2007). Variations in glacier mass and, eventually, changes in glacier area are determined by the balance between incoming (mass gain) and outgoing (mass loss) terms. The underlying processes are strongly controlled by atmospheric factors. But the link between glaciers and climate is not a simple onedirectional relationship. Via a number of feedback mechanisms, changes in ice cover can potentially influence the ambient air and thereby modify the local climate. The unique surface characteristics of ice, namely the comparatively high albedo and a surface temperature not exceeding 0°C, influence the components of the surface energy balance. This can result in a cooling of the overlying air mass and a modification of the small-scale circulation (e.g., generation of katabatic winds). In case of the large polar ice sheets a further feedback loop exists via ice thickness, since the simple geometry and shape of an ice sheet can alter the atmospheric circulation (Herterich 2001). In addition to direct physical interactions between surface ice cover and the overlying air mass a further feedback mechanism consists in glacial melt water entering the ocean. This potentially influences oceanic circulations and deep water formation and contributes to sea level rise (IPCC 2007). Also seasonal runoff patterns can be strongly influenced by the presence of glaciers (e.g., Singh and Singh 2001;Beniston 2003;Jansson et al. 2003) and glacier mass variations will therefore induce changes in the hydrological regime. This is likely to have major consequences for water resources management (e.g., Barry 2003;Barnett et al. 2005).
Although glaciers constitute an important component of the terrestrial water cycle and are linked to climate processes via the described feedback loops, they are poorly represented in today's climate models. Both state-of-the-art general circulation models (GCMs) and regional climate models (RCMs) make use of binary static glacier masks which define a land grid point as being either totally glacierized or totally ice free for the entire simulated period. They do neither allow for the simulation of changes in ice area nor for a reliable estimate of changes in the total ice volume. Possible feedbacks via changes in ice extent triggered by changing climatic conditions are neglected. This simplistic approach originally stems from the necessity to represent the large and slowly reacting polar ice sheets as a lower boundary condition. In recent years, the climatic effect of these large ice bodies has been fully acknowledged and the associated feedback mechanisms were investigated by bi-directionally coupling dynamical ice sheet models to GCMs (e.g., Ridley et al. 2005;Winguth et al. 2005;Vizcaino 2006) and by integrating them into Earth System Models of Intermediate Complexity (EMICs; e.g., Wang and Mysak 2002;Kageyama et al. 2004;Calov et al. 2005). Glaciers and ice caps (i.e., surface land ice bodies outside the two polar ice sheets), on the other hand, are usually much smaller than the resolved scale of a climate model and are not accounted for explicitly. To our knowledge no effort has been undertaken so far to fully include these smaller ice bodies in climate models. Such an upgrade of existing models would be highly desirable because the contribution of glaciers and ice caps to global sea level rise is considered as non-negligible (Kaser et al. 2006;Meier et al. 2007) and regional climatic changes in glacierized areas are generally of high interest (for example with respect to effects on river discharge). Therefore, an interactive glacier energy and mass balance scheme including the dynamic adjustment of glacierized areas was developed within the present study. The new scheme was implemented into the RCM REMO where it complements the static glacier mask used so far.
The main characteristics of the new parameterization scheme and the experimental design for investigating its performance are described in Sect. 2. Section 3 then presents modeling results for the second half of the twentieth century with a focus on glacier mass balance and glacier area changes. Also the sensitivity of the results with respect to individual model parameters is investigated. We end with a discussion (Sect. 4) and with concluding remarks on advantages and problems of the new scheme and on possibilities for further model improvement (Sect. 5).
2 Models and methods 2.1 The regional climate model REMO The basic modeling tool for the development and the assessment of the new glacier parameterization is the RCM REMO (REgional MOdel), a three dimensional hydrostatic atmospheric circulation model (Jacob 2001;Jacob et al. 2001). The dynamical core of REMO is based on the former numerical weather prediction model of the German Weather Service (Europa-Modell; Majewski 1991) while the physical parameterizations were adopted from the GCM ECHAM4 (Roeckner et al. 1996).
At the lateral boundaries of the regional model domain, the prognostic atmospheric variables of REMO (horizontal wind components, surface pressure, temperature, specific humidity, cloud liquid water) are relaxed towards a largescale forcing provided by a GCM or by re-analysis/analysis products. At the lower boundary, REMO is forced by the respective land and sea surface characteristics. For this purpose the so-called tile approach is applied: the total surface area of an individual model grid box can consist of a land, water and a sea ice fraction on a subgrid scale (tiles, expressed in percent). This technique allows for a smooth transition from land to sea along the coastlines as well as from ice-covered to ice-free ocean surfaces. The surface fluxes are calculated separately for each tile and are subsequently averaged within the lowest atmospheric level using the respective area contributions as weights. Mountain glaciers are not explicitly accounted for in REMO's land surface scheme. Even their contribution to the mean grid box background albedo (albedo of the snow-free land surface) is neglected since snow-and ice-covered pixels from the underlying high-resolution, satellite-derived albedo dataset are generally excluded. Only the polar ice sheets are represented by a static glacier mask which is allocated at the start of the simulation and which does not change in time.

Overview
In order to account for the energy and mass balance of mountain glaciers, the most recent RCM version REMO 5.3 was extended and modified. These changes primarily concern the physical parameterizations of the model and are represented by three additional subroutines. The new model version (REMO 5.3 plus the glacier extension) is called REMO glacier , hereafter. In order to investigate the influence of glaciers on catchment hydrology, a river routing scheme is additionally coupled to REMO glacier in an offline mode (the HD Model; Hagemann and Dümenil 1998). The present paper, however, puts a focus on the performance of the glacier parameterization itself. The analysis of river discharge will be presented in a subsequent study.
Due to the existence of high-quality observational datasets on glacier characteristics and climatic parameters, the European Alps were chosen as initial test region. Compared to other parts of the world, this region is only slightly glacierized with a total ice surface area in the mid 1970s of about 3,000 km 2 and a total ice volume of about 100 to 130 km 3 (Paul et al. 2004). Currently available RCM resolutions (10 to 50 km) are too coarse to resolve individual mountain glaciers in this region. Hence, the overall behavior of glaciers contained within an RCM grid box has to be described on a subgrid scale in terms of the atmospheric fields as resolved by the RCM. Given the apparent scale mismatch, a spatially resolved simulation of glacier surface energy and mass balance, as it is done using high-resolution distributed models for individual glaciers and the subsequent consideration of ice flow dynamics are not possible. Furthermore, the large amount of data required by such models, for instance information on the specific topographic setting of a glacier or the bedrock topography, can hardly be provided for all existing glaciers contained within a typical RCM domain. An explicit treatment of each individual glacier would also result in a considerable increase in computing time, which is still a limiting factor in regional and global climate simulations. Due to these reasons, the present work has to be seen as an attempt to simulate the overall behavior of glaciers contained in a climate model grid box by strongly simplified and idealized concepts.
In the following section the basic characteristics of the new parameterization scheme will be presented. A more detailed and comprehensive description can be found in Kotlarski (2007).

Fractional surface coverage and cuboid concept
In order to account for mountain glaciers on a subgrid scale, a fourth surface tile ''glacier'' was introduced into REMO in addition to non-glacierized land, water, and sea ice (Fig. 1). The new fraction represents the total area covered by glaciers in the respective grid cell. Its extent is allowed to grow or to shrink depending on the mass balance, but is restricted to the total land surface area of a grid box. Glacier ice is not allowed to expand over water surfaces.
A major simplification of the new scheme concerns the treatment of the glacier population within a climate model grid box. In REMO glacier all individual glaciers located in a specific grid box are represented by one single ice body in form of a two-layer cuboid with a surface area A and a total thickness h (see Fig. 2). The corresponding volume V equals the sum of the estimated volumes of all individual glaciers. Similarly, the surface area A of the cuboid is the sum of all individual glacier areas and determines the glacier-covered fraction of the total grid box area. The idea to group individual glaciers in pools, defined by the extent of a climate model grid box, is supported by the observation that nearby glaciers often show a similar response to a given climatic forcing in terms of their mass balance evolution (e.g., Cogley and Adams 1998;Haeberli et al. 1989;Letréguilly and Reynaud 1990). The glacier cuboids of adjacent grid boxes are independent of each other and do not interact directly. Large-scale ice flow across grid box boundaries is neglected. Given that the size of glaciers in the European Alps (and in many other regions as well) is typically smaller than an RCM grid box, this simplification seems to be justifiable in a first approach.
The altitude of the glacier cuboid surface is assumed to be identical to the mean grid box altitude. It is fixed throughout the entire model simulation and the effect of a possible surface elevation-mass balance feedback is not taken into account. A constant ice density of 880 kg/m 3 , a typical value for glacier ice after Singh and Singh (2001), is assumed for the entire cuboid. In terms of the surface heat budget, glacier ice is taken into account up from a critical ice depth of 0.1 m w.e. In case of large ice depths, the cuboid is split into an upper layer of a thickness of 0.2 m ice representing the interface to the atmosphere or the overlying snow cover, and a lower layer representing the remaining ice volume. Temperature diffusion between the snow, ice and soil layers is calculated based on characteristic values of heat conductivity.

Surface energy balance and snow cover
In each time step the full surface energy balance of the glacierized fraction is computed according to: where Q N is the net radiation (shortwave and longwave) (W/m 2 ), Q H the turbulent flux of sensible heat (W/m 2 ), Q L the turbulent flux of latent heat (W/m 2 ), Q G the ground heat flux (W/m 2 ), Q M the energy consumed by melt of snow or ice (W/m 2 ), dQ snow/ice is the change of the heat content of the upper snow or ice layer (W/m 2 ) The calculation of the individual flux terms follows the parameterizations already applied in standard REMO (DKRZ 1994;Roeckner et al. 1996). Similarly to the non-glacierized land fraction, the glacierized part of a grid box can be covered by a snow layer (Fig. 2). In this case, the residual of the surface energy fluxes is used to warm or to cool an upper snow layer with a maximum thickness of 0.1 m. For calculation of Q N the snow albedo a snow is set to a value between 0.6 and 0.8, linearly depending on the snow temperature (see DKRZ 1994). An ageing effect, i.e., the decrease of snow albedo with snow age, is not considered. It is generally assumed that the snow cover thickness is constant over the entire glacier surface and that the glacier cuboid is either completely covered by snow or completely snow-free. This is a strong simplification regarding the elevation range of real glaciers and the corresponding variability of snow depth.
Once the snow has completely melted, the upper ice layer with an albedo a ice of 0.4 immediately becomes exposed. This value is the mean albedo of clean ice (Paterson 1994) and has to be considered as an upper limit of the surface albedo of mountain glaciers. As in the case of a snow , a ice is assumed to be valid for the entire glacier cuboid (uniform albedo of the entire glacier surface).

Surface mass balance and runoff generation
Mass balance processes are assumed to occur uniformly over the entire glacier cuboid. In a defined balance period, the entire glacier cuboid experiences either net accumulation or net ablation (or, in an unlikely case, a zero balance). The modeled mass balance of the cuboid is assumed to represent the mean specific mass balance b n (m w.e./year) (hereafter simply referred to as mass balance), averaged over all individual glaciers contained in the respective RCM grid box.
The main accumulation process is snowfall onto the glacier cuboid surface which directly adds mass to the system glacier cuboid ? snow cover. The snowfall rate is determined by the atmospheric model component plus possible modifications due to the subgrid redistribution of snow accumulation (see Sect. 2.2.5 and Electronic Supplementary Material). A further accumulation process is rime formation which is directly added to the upper ice layer. The refreezing of rain on the glacier surface is neglected. Snow metamorphism is not calculated Fig. 1 The extended tile approach. Example for a grid box covered by non-glacierized land (45%), glacier (15%), water (30%), and sea ice (10%) Fig. 2 Profile through the land fraction of a partly glacierized grid box in REMO glacier explicitly. Instead, the transformation of snow into ice is parameterized via a snow age threshold. Once the snow layer on any of the two fractions (glacierized, non-glacierized) has reached a critical age of 730 days, it is assumed that parts of the snow volume on this fraction have passed the metamorphism process from snow to ice. This critical snow age is a compromise between the time it takes for snow to be transformed into firn (about 1 year) and the total time necessary for the transformation of snow into ice in temperate regions (5-10 years at minimum; e.g., Hambrey and Alean 2004).
The temperatures of the snow and ice layers are not allowed to exceed the melting point of T melt = 0°C. If this threshold is surpassed as a result of a positive surface energy balance or of heat diffusion, the temperature is adjusted back to 0°C. The released energy is used for ice/ snow melt in the respective layer. Ablation is furthermore possible by surface sublimation (Q L \ 0 in Eq. 1). Internal accumulation, i.e., the refreezing of meltwater within the snow or ice, is neglected as it does not change the total snow and ice mass and its accelerating effect on the elimination of a winter's cold wave in the glacier (e.g., Paterson 1994) could hardly be resolved using a two-layer scheme. Also calving processes and the removal of snow by wind are not accounted for. Meltwater originating from a snow layer on top of the ice cuboid and from the ice cuboid itself either directly contributes to grid box surface runoff (snow and upper ice layer) or to infiltration (lower ice layer).

Atmospheric subgrid variability
On an RCM subgrid level, fundamental differences in the climatic forcing can be expected between glacierized and glacier-free parts. For instance, glaciers generally tend to form in the coldest and wettest parts of terrain. Glacier surfaces can also be subject to pronounced shading effects (location in valley floors), causing a reduced global radiation flux. Furthermore, large parts of a glacier's surface are often situated on inclined slopes with a northern orientation which also results in a reduced exposure to global radiation (e.g., Evans 2006). In order to include some of these effects, the subgrid variability of snowfall and of global radiation within a climate model grid box are explicitly accounted for in REMO glacier by applying simple scaling concepts. In each time step, both quantities are redistributed within an RCM grid cell between the glacierized and the glacierfree surface fraction based on current glacier outlines in the study region, a high-resolution observational precipitation dataset and the application of an offline radiation model which accounts for sloping and shading effects. In the case of snow accumulation, a fixed redistribution factor of j = 2.0 is applied for all grid cells (the glacier cuboid receives the double amount of snowfall compared to the grid box mean). For further details the reader is referred to the Electronic Supplementary Material and to Kotlarski (2007).

Glacier area adjustment
A key element of the new scheme is the calculation of the areal extent of the glacier cuboids (see Fig. 2). In REMO glacier the dynamic response of glacier geometry to a given mass balance forcing is not computed explicitly. Instead, simplifying assumptions are used which relate the calculated change in thickness (the mass balance in m w.e.) to a change in glacier area. Previous studies have shown that mean glacier thickness increases with glacier area (e.g., Müller et al. 1976;Maisch 1992) and therefore glacier volume (mean thickness multiplied by glacier area) could be obtained from glacier area information. Apart from the strong autocorrelation in the statistical relation between area and volume (e.g., Haeberli et al. 2008), a study by Bahr et al. (1997) provides a physical background for such a relationship. Assuming steady-state conditions the glacier volume V (m 3 ) can be related to the glacier area A (m 2 ) via a power law: with a dimensionless scaling coefficient c = 1.375 and a calibrated constant c. The corresponding relation between the glacier area and the mean ice thickness " h (m) is obtained by dividing both sides in Eq. 2 by the area A: The power law relation has its main application in estimating the volume for very large samples ([1,000) of glaciers (e.g., van de Wal and Wild 2001; Raper and Braithwaite 2005). For such global applications, Chen and Ohmura (1990) derived c = 1.357 and c = 0.206 [if A and V are given in (m 2 ) and (m 3 ), or c = 28.5 if A and V in (10 6 m 2 ) and (10 6 m 3 ), respectively] based on data from 63 glaciers in different regions of the world, including 15 glaciers of the Alps. In REMO glacier these values are used to relate the glacier volume V and the mean ice thickness " h to the glacier area A via Eqs. 2 and 3. Note that the parameter c cannot be expected to remain constant under changing climatic conditions (Bahr et al. 1997) as basic changes in glacier geometry (e.g., loss of the tongue) will result in different glacier shapes.
The thickness-area relation from scaling theory (Eq. 3) is similar to the empirical relation derived by Maisch (1992) for a sample of 16 existing glaciers in the Alps and 63 late-glacial ice bodies: Accordingly, the glacier volume V (m 3 ) can be calculated as While the relations (4) and (5) are more representative for glaciers in the Alps, (2) and (3) have the clear advantage of providing a strong monotonic behaviour with a zero thickness for a glacier size of zero. Moreover, the sample of individual glaciers used in this study is quite large ([5,000 in the beginning). Both points justify the choice of Eqs. 2 and 3 as the default method for area adjustment. Methods (4) and (5) are used as an alternative in REMO glacier .
Altogether, three optional area adjustment methods were implemented in the parameterization scheme, controlled by a dedicated switch in the RCM job script. The three methods differ with respect to the volume-area relation applied and to whether the area adjustment is carried out for the entire grid box ice cuboid or for each individual glacier contained in a climate model grid cell (see Kotlarski (2007)

for details):
Method 1 (power law for glacier cuboid): The grid box glacier cuboid is treated as one single glacier and the power law (Eqs. 2, 3) is applied to the entire ice cuboid in each time step.
Method 2 (power law for individual glaciers): After each year of model simulation, the power law (Eqs. 2, 3) is applied to each individual glacier contained within a climate model grid cell. During that year, all individual glaciers are assumed to experience the same mass balance b n as the entire glacier cuboid.
Method 3 (relation after Maisch (1992) for individual glaciers): After each year of model simulation, the relation of Maisch (1992;Eqs. 4, 5) is applied to each individual glacier contained within a climate model grid cell. As in method 2, all individual glaciers are assumed to experience the same mass balance b n .

Experimental design
Four baseline experiments were carried out for the period 1958 to 2003. These include three REMO glacier simulations, differing with respect to the area adjustment method applied (GlacA: Method 1, GlacB: Method 2, GlacC: Method 3) and one control simulation (CTRL) using the standard model version REMO 5.3 without glacier parameterization. Figure 3 displays the common model domain for all experiments. It completely covers the European Alps as well as further parts of central, western and southern Europe. A horizontal resolution of 1/6°on a rotated lat/lon grid, corresponding to a grid box size of approximately 18 km 9 18 km, and 20 vertical levels were used in all simulations. The lateral boundary forcing was directly provided by the ERA40 re-analysis (01/1958to 07/2002Uppala et al. 2005) and the operational analysis of the ECMWF (08/2002 to 12/2003). The internal model time step was 100 s.
For initializing the glacierized grid box fraction and the mean ice depth in each grid cell, glacier area information extracted from the World Glacier Inventory (WGI) was used (WGI 1989;NSIDC 1999). In a first step, glacier areas from this database, which are assumed to represent the conditions in 1973, were scaled in order to obtain glacier areas for the year 1958. For this purpose, the relative area change from 1957 to 1973 for 17 glaciers in the European Alps of different size was related to their area in 1973 by logarithmic regression. The resulting relation between glacier area and relative area change from 1957 to 1973 was then used to derive areas for January 1958 for all glaciers contained in the WGI (see Kotlarski (2007) for details). Subsequently, individual glacier volumes in 1958 were calculated based on the power law (GlacA, GlacB) and the relation of Maisch (1992;GlacC) and were summed up for each REMO grid cell. In most cases, the glacier-covered surface fraction in 1958 computed by this method is smaller than 5%. Still, maximum glacier coverages of more than 45% are obtained in individual high-elevation grid cells.

Glacier mass balance
3.1.1 Entire Alps Figure 4 shows the evolution of the mass balance b n for the entire Alpine region for the baseline experiments as well as for the recent estimate of Dyurgerov and Meier (2005;referred to as D&M 2005, hereafter). The latter was computed by averaging all available glacier mass balance values in the Alps using the respective glacier surface areas as weights. Table 1 (first column) summarizes the results in terms of the mean mass balance " b n , the cumulative mass balance P b n , the standard deviation of b n , and the mean model error.
b n only slightly differs between the experiments. The mean value " b n for GlacA/GlacB/GlacC amounts to -0.08/-0.09/-0.08 m w.e./year, respectively. The general temporal evolution of the simulated mass balance is in good accordance with the observation-based dataset. In both cases positive values appear in the late 1970s and the early 1980s. Before and after that period, mass balances are mainly negative, both in the simulations and in the observation-based dataset. The strong glacier mass loss towards the end of the twentieth century is considerably underestimated by the model. This shortcoming leads to a clear underestimation of the cumulative glacier mass loss over the entire investigated period (Table 1). The rootmean squared error of b n , which can be interpreted as a typical error magnitude, is 0.59 m w.e./year for all experiments. This value is almost as large as the amplitude of interannual mass balance variations in the observationbased time series (0.69 m w.e./year, SD b n in Table 1), indicating a considerable error magnitude. Still, the simulated mass balance in the three experiments lies within ±0.3 m w.e. of the observation-based value in 21/ 22/22 out of 43 years, i.e., during about 50% of the time (number of circles in the gray-shaded area in Fig. 4). The year-to-year variability of b n is basically reproduced by the model, except for the period 1986 to 1995 (correlation coefficient r = 0.66/0.67/0.67 between GlacA/GlacB/ GlacC and observations for the period 1961-2003). However, the amplitude of interannual mass balance variations as expressed by the standard deviation of the detrended time series of b n is underestimated (SD b n in Table 1). This can partly be explained by the small number of individual  Figure 5 presents the mass balance analysis on a catchment scale (see also Table 1). Again, observational values were derived from D&M 2005 by weighting all available individual mass balance values in the four main Alpine catchments by the respective surface areas. A good agreement between simulated and observed mass balances is obtained for the Danube and the Rhone catchment. The mean b n model error ranges from 0.19 to 0.22 and from 0.11 to 0.12 m w.e./year, respectively (Table 1). For the Danube (Rhone) the simulated mass balance lies within ±0.3 m w.e. of the observation-based value during 50% (40%) of the time (gray-shaded area in Fig. 5). Still, the observed mass loss in the last two decades of the twentieth century is clearly underestimated in both regions. A large positive mass balance bias occurs in the Rhine catchment in the period 1980 to 1995. Here, the strong glacier mass loss is not captured by the model which leads to a pronounced mean error of b n and to a positive cumulative mass balance by the year 2003 instead of the observed mass loss (Table 1). The assessment of model performance in the Po catchment is difficult as observation-based mass balance information is only available from 1987 onwards and the validation period is comparatively short. Furthermore, the observational dataset only relies on a very small number of individual glaciers. If the observation-derived mass balance values are trusted in, the model performs relatively well in the last 17 years of the simulation and roughly reproduces the observed glacier mass loss (Table 1). Regarding the amplitude of interannual mass balance variations, REMO glacier underestimates interannual fluctuations in all catchments (SD b n in Table 1). Again, this fact can partly be attributed to the small number of individual glaciers in the observational time series (see above).

Mass balance sensitivities
An important characteristic of the climate-glacier relationship is the mass balance sensitivity of a glacier with respect to variations in a certain climatic parameter. In all catchments and for all experiments, the simulated mass balance shows a clear negative correlation with mean summer (April-September) air temperature (Fig. 6). The correlation coefficients r (ranging from -0.50 to -0.75) agree well with those found by Braithwaite and Zhang (2000) for four Swiss glaciers based on observational data (r ranging from -0.58 to -0.70). The slope of the linear regression curve, expressed by the regression coefficient b, ranges from -0.27 to -0.47 m w.e./°C. It can, with some care, be interpreted as the climate sensitivity of the simulated glacier mass balance with respect to mean summer air temperature. The simulated range of b roughly corresponds to the sensitivities derived by Oerlemans and Reichert (2000) and Reichert et al. (2001) for the two Alpine glaciers Hintereisferner (about -0.35 m w.e./°C) and Rhonegletscher (about -0.34 m w.e./°C). 1 Applying multiple linear regression, Braithwaite and Zhang (2000) derived larger mass balance sensitivities with respect to Table 1 Mean glacier mass balance " b n (m w.e./year), cumulative mass balance P b n (m w.e.) and standard deviation SD bn of the detrended b n time series (m w.e./year) as simulated by the three baseline experiments (GlacA, GlacB, GlacC) and as derived from D&M 2005D&M Experiment 1961D&M -2003D&M 1959D&M -2003D&M 1987D&M -2003 Alps Rhine Danube Rhone Po The mentioned studies are based on modeled mass balance sensitivities with respect to monthly perturbations in temperature. The sensitivity of b n to mean summer temperature, as presented here, was derived by summing up the monthly sensitivities for the months April to September. mean JJA (June-August) air temperature of -0.44 to -0.64 m w.e./°C for four Swiss glaciers. The sensitivities obtained by the application of a degree-day model for the same sample of glaciers are even higher (-0.69 to -0.89 m w.e./°C). Hence, there is an indication that the mass balance sensitivity with respect to summer air temperature is underestimated by REMO glacier .
A significant positive relation exists between the simulated mass balance and mean annual precipitation in all catchments and for all experiments (Fig. 7). The correlation coefficients r range from 0.60 to 0.85. They are larger than those found by Braithwaite and Zhang (2000) for four Swiss glaciers analyzing observational data (0.39 to 0.61). They are also larger than those derived for summer air temperature (see above). The corresponding regression coefficients b, which can be interpreted as mass balance sensitivities with respect to changes in annual precipitation, are largest for the Rhine and the Rhone catchment (0.0324 to 0.0326 and 0.0211 to 0.0220 m w.e./%, respectively) and lowest for the Danube and the Po (0.0176 to 0.0178 and 0.0155 to  Klok and Oerlemans (2004), about 0.025 m w.e./% for Hintereisferner by Oerlemans and Reichert (2000), about 0.026 m w.e./% for Rhonegletscher by Reichert et al. (2001); the last two values have again been calculated by summing up monthly sensitivities]. Hence, in contrast to the temperature sensitivity of glacier mass balance, the dependency of b n with respect to precipitation is overestimated by the model. An increase in annual precipitation leads to a comparatively large increase of glacier mass balance, and the influence of interannual variations of precipitation is exaggerated compared to the influence of variations in air temperature.

Glacier area
After 46 years of continuous model simulation the glacier area in the Alps has considerably decreased with respect to the initial state in all three baseline simulations (Fig. 4, lower panel). The total loss of glacier area from 1958 to 2003 ranges from -17.1% (GlacC) to -23.6% (GlacA). In all experiments a major loss of ice area occurs during the first part of the simulation from 1958 until the mid 1970s as a result of negative mass balances. In December 1972, the total ice area in the three experiments amounts to 2,754, 2,804 and 2,786 km 2 , respectively. These values can be compared to the total Alpine glacier area as derived from the WGI which is assumed to be representative for the year 1973 (see Sect. 2.3). This area is 3,040 km 2 , it is considerably larger than the simulated values. If the scaling of the WGI glacier areas (which has been carried out in order to obtain glacier areas for the year 1958, see Sect. 2.3) is assumed to be valid, the simulated loss of ice area from 1958 to 1972 is too large (-14.5 to -16.0% in the experiments compared to -7.3% obtained from area scaling).
In the following period of positive mass balance years, the simulated glacier area increases and reaches a secondary maximum in 1986. This phase approximately corresponds to the documented stabilization and advance period of Alpine glaciers in the 1970s and 1980s (e.g., Zemp 2006). Towards the end of the simulation negative mass balances again cause a decrease of total ice cover. While the temporal evolution of glacier area is similar in GlacB and GlacC (i.e., in those simulations which rely on an area adjustment method based on individual glaciers), GlacA shows a stronger overall reduction of glacier area from about 1970 onward. The simulated loss of ice area from 1973 to 2000 can be compared to the estimates of Paul et al. (2004) and Zemp et al. (2008). Both studies are based on the WGI and the satellite-derived Swiss glacier inventory of the years 1998/1999 (Paul 2007) and estimate a total loss of Alpine glacier area from 1973/1975 to 1998/ 2000 of about -22%. The simulated area loss in this period ranges from -0.6% (GlacC) to -5.7% (GlacA), i.e., the strong and accelerating retreat of Alpine glaciers in the last two decades of the twentieth century is clearly underestimated. This fact is primarily due to the evolution of glacier Fig. 7 Relation between annual (October-September) precipitation anomaly (%) and mass balance b n (m w.e./year) for the three baseline simulations and the four main Alpine catchments. Values correspond to catchment mean values, weighted by the respective glacierized grid box area. r correlation coefficient (asterisk indicates a significant linear correlation, 99% confidence level), b regression coefficient (slope of line of best fit) area in the Rhine and the Rhone catchment. Here, a net increase of ice area is simulated from the mid 1970s onward (not shown) which is associated with a clear positive bias of the simulated mass balance compared to the observation-based estimates (Fig. 5). Also the total loss of ice area from 1958 to 2000 is underestimated. For the three baseline simulations it amounts to -20.8, -16.5 and -15.5%, respectively, compared to -27.7% as derived from the combination of the area scaling (1958)(1959)(1960)(1961)(1962)(1963)(1964)(1965)(1966)(1967)(1968)(1969)(1970)(1971)(1972) and the estimates of Paul et al. (2004) and Zemp et al. (2008Zemp et al. ( ) (1973Zemp et al. ( -2000. Figure 8 shows the spatial distribution of the glacier area change on a grid box scale. The majority of glacierized grid boxes experiences a pronounced decrease of ice area. In contrast, a number of grid cells in the north-western part of the Alps as well as in the central part show a growth of more than 100% with respect to glacier areas in 1958. By combining the horizontal distribution of glacier area changes with the regional topography (Fig. 3), it becomes evident that those grid cells experiencing a net gain in ice area (and ice volume) are exclusively located along the northern Alpine ridge as defined by the model topography. All grid boxes south of the main ridge show a net decrease of glacier area. The spatial distribution of ice volume changes is very similar to the changes of ice area, i.e., grid cells with a net increase of ice area also exhibit a gain of ice volume and vice versa (not shown).

Simulated climate
As the new parameterization scheme is fully implemented into REMO, the atmospheric forcing for the calculation of glacier mass and energy balance is directly provided by the RCM's atmospheric component. Even if it is assumed that the glacier parameterization scheme itself functions perfectly (general applicability of the cuboid concept, realistic description of subgrid variability of atmospheric parameters, appropriate description of glacier-atmosphere interactions), biases in the atmospheric forcing would still lead to errors in the simulated mass balance. On the other hand, errors in the imposed atmospheric parameters (e.g., precipitation) could compensate for shortcomings of the glacier parameterization scheme itself (e.g., subgrid snow redistribution). The same is true for biases of more than one atmospheric parameter which, in principle, could compensate each other (e.g., too high temperatures vs. too much precipitation). In order to assess the performance of the new glacier parameterization scheme it is therefore necessary to validate the simulated climatic parameters involved in glacier mass and energy balance.

Precipitation
While REMO glacier overestimates precipitation in regions below 2,000 m a.s.l., high-altitude grid cells receive too little precipitation when compared against the two observational datasets CRU TS 1.2 (Mitchell et al. 2004) and ALP-IMP (Efthymiadis et al. 2006). Both datasets are provided on a regular 1/6°grid and were conservatively interpolated onto the rotated 1/6°REMO grid. They do not account for the systematic bias of rain-gauge measurements. The high-altitude precipitation bias is primarily due to an underestimation of precipitation of up to 35% from June to January. In the winter season, this negative bias can reach up to -65% in individual grid cells. Generally, the differences in the simulated amount of precipitation between the three baseline simulations GlacA, GlacB and GlacC are negligible (see Kotlarski (2007) for details). Figure 9 shows the relation between the seasonal precipitation bias of GlacA and the glacierized grid box fraction in the four main Alpine catchments. Most grid cells with a mean glacier cover of more than 10% experience an underestimation of precipitation in summer and autumn (lower panels). In winter and spring (upper panels) the relation is more complex. While the strongly glacierized grid boxes in the Rhone catchment exhibit a large negative precipitation bias of up to 65%, precipitation is too high for most grid cells in the Rhine catchment. These boxes are mainly found along the northern Alpine slopes (see discussion below). A similar overestimation of precipitation in strongly glacierized grid cells is found in the Danube catchment in springtime. Interestingly, the spatial  Fig. 3). The gray color indicates nonglacierized grid cells pattern of the winter precipitation bias (Fig. 10) closely resembles the pattern of glacier area change (Fig. 8). A strong overestimation of precipitation is associated with a growing glacier area and vice versa. The positive bias of glacier mass balance in the Rhine catchment (Sect. 3.1.2) is obviously caused by too high precipitation sums in the model. Figures 9 and 10 also reveal important differences between the two observation-based reference datasets (CRU and ALP-IMP). Still, the general pattern and the magnitude of the precipitation bias are similar for both cases, which indicates a robust model evaluation.

Temperature
The simulated near-surface air temperature in all baseline experiments was validated against two observational datasets: (a) the gridded CRU TS 1.2 dataset, provided on a regular 1/6°grid, and (b) data of individual high-altitude stations of the Swiss Federal Office for Meteorology and Climatology's network (see Kotlarski (2007) for details). In order to account for the elevation difference between the REMO grid cells and the CRU cells/measurement sites, a temperature correction assuming a uniform lapse rate of -0.65°C/100 m was applied prior to the analysis. While the performance of REMO glacier is very good in low-elevation regions problems arise in high-altitude grid cells. Here, the simulated air temperature shows a clear negative bias with respect to observations throughout the entire year. For the experiment GlacA, Fig. 11 shows the relation between the seasonal temperature bias with respect to the CRU dataset and the mean glacier-covered grid box fraction. Biases for GlacB and GlacC are not shown because the temperature differences between the three baseline experiments are negligible (see Kotlarski 2007). Throughout most parts of the year temperature is underestimated with largest negative biases of up to -10°C in autumn and winter (upper left and lower right panels). The negative temperature bias generally grows with increasing fraction of glacier cover (gray lines). This bias can be expected to have a positive effect on glacier mass balance via an increase in the fraction of solid precipitation (snowfall) and a reduced flux of sensible heat towards the glacier surface. In summertime a considerable fraction of glacierized grid cells shows a positive temperature bias. These boxes are primarily located south of the main Alpine ridge in the Po catchment (green color in Fig. 11).

Global radiation
As observation-based gridded datasets were not available, the validation of global radiation relies on quality-checked observational data for individual stations. The mean annual cycle of global radiation was computed for 24 sites in the Alpine region based on two data sources: the Global Energy Balance Archive (GEBA; Gilgen et al. 1998;Gilgen and Ohmura 1999) and the Swiss Federal Office for Fig. 9 Relation between the mean glacier-covered grid box fraction and the seasonal precipitation bias  (%) of GlacA with respect to ALP-IMP (crosses) and CRU (circles). Each marker represents one glacierized REMO grid cell and is color-coded according to the respective catchment. The gray lines denote the mean Alpine-wide precipitation bias for each 5% interval of mean glacierized fraction (average of the individual grid cell biases) with respect to ALP-IMP (solid) and CRU (dashed) Meteorology and Climatology. A detailed analysis can be found in Kotlarski (2007). In summary, REMO glacier shows a good reproduction of the annual cycle of global radiation at all sites. A systematic negative bias seems to prevail in winter and especially in springtime at high-elevation sites. Here, the observed radiation flux is underestimated by up to 80 W/m 2 in the model simulations which is possibly connected to inaccuracies in the simulated cloud cover. At a number of stations, the maximum solar radiation flux in summer is overestimated by several tens of W/m 2 .
The described biases can be expected to directly affect the surface energy balance and, hence, the simulated glacier mass balance. Biases in the summer radiation flux can be thought to have the largest influence since the absorbed fraction of total radiation is larger due to a lower surface albedo (disappearance of snow and exposure of the darker ice surfaces). The overestimation of summertime global radiation will probably lead to a net negative bias of glacier mass balance (enhanced ablation). In contrast, the underestimation of wintertime and springtime global radiation in REMO glacier will be less effective and the associated positive mass balance bias (reduced ablation) will not compensate for the enhanced summer melting.

Climatic effect
In order to investigate the effect of the new glacier parameterization scheme on the simulated climate, the CTRL experiment (using the standard version REMO 5.3) can be compared to the REMO glacier simulations. Such a comparison does not only reveal the (model-specific) effect of including the glacier subgrid parameterization scheme Fig. 10 Mean bias of winter (DJF) precipitation  (%) in GlacA with respect to CRU (left panel) and to ALP-IMP (right panel). Only the central part of the model domain is shown (compare to Fig. 3). The gray color indicates nonglacierized grid cells Fig. 11 Relation between the mean glacier-covered grid box fraction and the seasonal temperature bias  (%) of GlacA with respect to CRU. Each marker represents one glacierized REMO grid cell and is color-coded according to the respective catchment. The gray line denotes the mean Alpine-wide temperature bias for each 5% interval of mean glacierized fraction (average of the individual grid cell biases) into REMO, but also shows the general influence of mountain glaciers on the regional climate. In REMO glacier this effect is primarily due to the special surface characteristics of the glacier-covered grid box fraction. A detailed analysis of these aspects is given in Kotlarski (2007).

Temperature
The inclusion of the glacier parameterization scheme leads to a decrease of near-surface air temperature with respect to CTRL. This effect is restricted to glacierized grid cells and shows a clear dependency on the glacier-covered grid box fraction (largest differences in strongly glacierized grid cells). The strongest influence is exerted in summer and autumn with temperature deviations of up to -7°C. In autumn, more than 50% of the glacierized grid cells show a significant negative temperature difference. This cooling effect is due to the combination of a comparatively high albedo of snow-free glacier surfaces, a glacier surface temperature never exceeding 0°C and the longer persistence of snow on glacier surfaces due to the subgrid snow redistribution. The cooling effect is restricted to the lower troposphere. In the 500 hPa level, no significant temperature deviation is found.

Precipitation
Precipitation is most strongly affected in the summer season and the effect is mainly restricted to the Alpine region (Fig. 12). More than 66% of the glacierized grid cells show a significant negative deviation with a maximum reduction of more than 50% in individual grid cells. No significant differences are found in winter. Compared to air temperature the dependence of the summer precipitation difference on the glacier-covered surface fraction is weak. Furthermore, the spatial pattern of the precipitation signal is generally more diffuse and grid cells in the vicinity of glacierized areas show distinct signals as well. The reduced precipitation in REMO glacier is mainly due to a decrease of summertime convective precipitation (i.e., precipitation generated by the model's convective parameterization scheme) which is reduced by up to 90% in grid cells with a high glacier cover. The influence on the resolved largescale precipitation component is much weaker. Presumably, one reason for the reduction of summer precipitation is the lower surface temperature in glacierized grid cells which reduces the near-surface temperature gradient and hinders the initiation of convection.

Surface evaporation and specific humidity
The supply of humidity by surface evaporation is strongly reduced, with more than 70% of the glacierized grid boxes showing a significant reduction in the glacier experiments (not shown). This leads to a significantly lower near-surface specific humidity in more than 60% of the glacierized grid cells (not shown). The influence on humidity is restricted to the lower troposphere. The reduction of surface evaporation is very likely caused by the existence of non-vegetated snow and ice surfaces in REMO glacier compared to vegetated/bare soil surfaces in CTRL. As a consequence, the local and regional recycling of humidity in summer by precipitation (see above) and evaporation is reduced.

Large-scale flow
The large-scale flow conditions, expressed by the spatial patterns of mean sea level pressure, surface pressure and geopotential height, are not significantly affected by the new glacier parameterization scheme. Obviously, the glacier surfaces in the Alps are too small in order to influence these large-scale parameters. It has to be kept in mind that the RCM domain (Fig. 3) is comparatively small and that the large-scale flow conditions are thus strongly determined by the lateral boundary conditions (which are identical in all experiments).

Model sensitivities
Several additional experiments were carried out in order to assess the main sensitivities of REMO glacier with respect to individual model parameters and regarding the representation of specific processes. The setup of these experiments is identical to the baseline simulation GlacA (area adjustment method 1). The model sensitivities can therefore be derived from a comparison between the sensitivity experiments and GlacA. Formally, also the baseline experiments GlacB and GlacC can be considered as sensitivity studies with respect to the choice of a specific area adjustment method (methods 2 and 3 instead of method 1 as used in GlacA). However, as shown in Sects. 3.1 and 3.2, the mean Alpine glacier mass balance is only slightly affected by the choice of a specific area adjustment method, whereas important differences can occur with respect to glacier area. The sensitivity analysis reveals that glacier mass balance and the corresponding area changes in REMO glacier strongly depend on (a) the chosen value for a number of model parameters and (b) the inclusion of subgrid processes which are so far not part of the standard parameterization scheme. Table 2 gives an overview of all sensitivity experiments and their main results. We here present the most important findings. For a complete description of all sensitivity experiments, the reader is referred to Kotlarski (2007).  The mass balance difference (4th column) refers to Alpine mean values An important control on the simulated glacier mass balance is exerted by the ice albedo (SensA). Replacing the fixed value of a ice = 0.4 in the default setup by a temperature dependent approach (a ice depending on the ice surface temperature and varying linearly between 0.3 for 0°C and 0.4 for -10°C) leads to a lower albedo during melt periods and to increased ablation. This causes a considerable decrease of the mean mass balance " b n . The ice albedo influence is strongest during periods of negative mass balance, which are often associated with below-thanaverage snow accumulation. In these years, the glacier surfaces get exposed earlier and changes of the ice surface albedo are more effective.
An even stronger influence on glacier mass balance is exerted by the implemented subgrid redistribution concepts for solar radiation and snow accumulation. If subgrid radiation scaling is omitted, i.e., if the glacierized grid box fractions receive the mean solar radiation flux as computed by REMO for a non-shaded horizontal plane (see Electronic Supplementary Material), the simulated glacier mass balance strongly decreases (SensB). Also neglecting/ reducing the subgrid redistribution of snow accumulation leads to a strong decrease of " b n via reduced snow accumulation on top of the glacier cuboids (SensC 1 , SensC 2 ). Obviously, the implemented concept of subgrid snow accumulation is a major factor in determining the simulated glacier mass balance.
Modifying the critical snow age for the transformation of snow into ice has only a minor influence on " b n (SensD). On the contrary, a considerable increase of mean glacier mass balance can be obtained by introducing subgrid variability concepts for air temperature which are, so far, not part of the standard parameterization scheme. If the effect of a lower air temperature at the mean glacier altitude of a grid cell on the phase of precipitation is considered, the simulated glacier mass balance increases due to increased snow accumulation (SensE 1 , SensE 2 ). If the computation of the surface sensible heat flux at the glacier surface is based on the air temperature at the mean glacier altitude, the sensible heat flux towards the glacier surface is reduced, which also results in a considerable increase of glacier mass balance (SensF).
A further important parameter is the aerodynamic roughness length of glacier ice. Increasing its fixed value from 0.001 to 0.01 m leads to a larger sensible heat flux towards the glacier surface and to a decrease of " b n via stronger ablation (SensG).

Discussion
The analysis presented above reveals that the broad temporal evolution, the year-to-year variability and the general magnitude of glacier mass balance in the Alps in the second half of the twentieth century can be approximately reproduced by REMO glacier . This is true for both the entire Alpine region and for individual catchments. On an Alpine-wide scale, the simulated mass balance b n lies within a ±0.3 m w.e. range of the observations during 50% of the time with a root-mean squared error of b n of 0.59 m w.e./year. The pronounced mass loss from the 1980s onward is not captured by the model, resulting in a clear positive bias of the cumulative mass balance in most regions by the year 2003. There is some indication that these shortcomings are related to inaccuracies in snow accumulation. The latter is determined by the amount of solid precipitation as simulated by the atmospheric component of the RCM and the subsequent redistribution of snowfall on a subgrid scale. It is the main factor determining the simulated mass balance b n . The sensitivity of b n with respect to precipitation is generally overestimated while the influence of air temperature is comparatively weak. Presumably, this fact is in parts responsible for the underestimation of the observed glacier mass and area loss from about 1980 onward. In the Rhine catchment, two simulations even show a net increase of glacier area during the 46 years of continuous model simulation. This is not consistent with observations and has to be considered as a major shortcoming of the parameterization scheme. The glacier growth is restricted to the north-facing slopes of the main Alpine ridge. Again, the clear dependency of the simulated area changes on the location of the grid cell with respect to the main ridge points to shortcomings in the simulated amount of snowfall. The amount of (solid) precipitation can be assumed to strongly differ between north-west and south-east facing slopes depending on the atmospheric flow conditions, triggering distinct topographic signals of glacier mass balance. A further point which probably contributes to the underestimation of glacier mass loss is the assumption of a uniformly distributed snow layer on top of the glacier cuboids and the neglect of the equilibrium line concept (i.e., the coexistence of accumulation and ablation areas on the same glacier). As a consequence of comparatively large amounts of snow accumulation in the 1970s, the fraction of glacier area which is permanently covered by snow is close to 100% in the Rhine and the Danube catchment in the second half of the experiments (not shown). The glacier cuboids are almost entirely covered by a snow layer and the ice itself does not become exposed during summer and autumn. Hence, ablation is mainly due to snow melt. The enhancement of ablation due to the exposure of the darker ice surfaces after the winter snow layer has melted (snow and ice albedo feedback) does not come into play, resulting in an underestimation of ablation rates and a positive mass balance bias.
The choice of a specific area adjustment method in REMO glacier has a minor influence on the simulated mass balance. Pronounced differences only occur with respect to glacier area which, in turn, determines the intensity of possible feedback mechanisms via changes in glacier extent. The net area decrease is most pronounced in GlacA. Given the general underestimation of glacier area loss, GlacA performs best in this respect which indicates that the consideration of individual glacier areas and volumes (GlacB and GlacC) is not mandatory in the current model setup.
Because the vast majority of Alpine glacier surfaces are located at high elevations, shortcomings in the simulated climatic forcing in high-altitude cells will inevitably influence the model results. One of the most important parameters is winter precipitation which is responsible for the majority of snow accumulation on Alpine glaciers. The validation of precipitation in this season indicates a strong positive bias in the Rhine catchment, especially during the 1970s and 1980s (not shown). This overestimation of winter precipitation very likely causes the positive bias of glacier mass balance in the Rhine catchment for the respective period. In contrast, strongly glacierized grid cells in the other catchments mostly show a negative precipitation bias which, in general, would lead to a negative mass balance bias. However, it is possible that for glacier surfaces a negative bias of wintertime precipitation is partly compensated by the choice of a too large snow redistribution factor j. A further compensating point is the negative temperature bias in most high-altitude grid cells for most parts of the year. This bias generally leads to an overestimation of the fraction of solid precipitation (snowfall) and to an increase of snow accumulation with a positive effect on glacier mass balance. Additionally, the general underestimation of temperature will reduce melt rates. The same is true for the negative bias of springtime global radiation. If the latter is related to an overestimation of cloud cover in REMO glacier the negative bias of global radiation might in parts be compensated by an increase of incoming longwave radiation due to a higher cloudiness. This assumption cannot be verified since the simulated incoming longwave radiation flux has not been validated within the present study due to the lack of observational data on an Alpine-wide scale.
The sensitivity studies carried out reveal that a crucial point in the model formulation is the description of the subgrid variability of global radiation, snow accumulation and air temperature. While in the first two cases simple concepts were incorporated into the scheme, the subgrid variability of air temperature is not yet accounted for explicitly. Given the clear dependency of glacier mass balance on air temperature variations via direct and indirect effects (sensible heat flux and phase of precipitation, respectively) the subgrid variability of elevation and therefore of air temperature should definitely be introduced in a next step. The sensitivity experiments SensE 1 , SensE 2 and SensF propose one possibility of doing so.
Since the explicit consideration of glacier elevations on a subgrid scale with its positive effect on glacier mass balance is not yet incorporated, the general positive mass balance bias of the current scheme indicates shortcomings in the formulation of other processes. One obvious candidate is the subgrid redistribution of snow accumulation. Large differences in the simulated mass balance can be achieved by variations in the controlling parameter j (see SensE 1 and SensE 2 ). A too strong redistribution of snow could easily (over-) compensate for the neglect of subgrid temperature variations (see above). At the same time, a redistribution factor of j = 2.0 in the standard parameterization is not more than a rough estimate and could well be considerably lower or higher. By definition, j accounts for both the subgrid variation of the amount and phase of precipitation and for the redistribution of fallen snow by wind and avalanches. A more physically based concept which discriminates between both features and which considers spatial and temporal variations of j would be desirable. For this purpose a detailed validation of the simulated snow depths both on the resolved and the subgrid scale is necessary, which has so far not been carried out.

Conclusions and outlook
For the first time, the new model system REMO glacier allows for a complete simulation of glacier-climate interactions, including the effect of glaciers on catchment hydrology. The results indicate that observed regional mass balances and glacier area changes can be approximately reproduced based on idealized concepts of glacier-climate interaction. It is shown that realistic results can only be achieved if the subgrid variability of atmospheric parameters within a climate model grid box is explicitly accounted for. On one side, the use of simplified concepts ensures a moderate increase of computing time (the new scheme requires 1% additional computer power) and allows for the application of the model on the scale of entire mountain ranges. The low demand of glacier-specific input data also facilitates the application in regions outside the European Alps. On the other hand, the high degree of simplification can obviously lead to important inaccuracies of the model, e.g., to an underestimation of the strong glacier mass loss towards the end of the twentieth century.
The effect of the new parameterization scheme on atmospheric parameters is restricted to the lower troposphere and to glacierized grid boxes themselves and their direct vicinity. The large-scale flow conditions are not significantly affected. Given this weak influence of the new scheme on the simulated large-scale climate and obvious problems in reproducing observed glacier area and volume changes, the general question arises if the interactive coupling of a glacier scheme to an RCM is meaningful at all. For the European Alps, a coupling is definitely beneficial if an emphasis is placed on the simulation and the analysis of regional water cycles at catchment scales. This is because both precipitation and evapotranspiration are significantly affected by the explicit and interactive treatment of glacier surfaces. The same is true for river discharge, which will be analyzed in more detail in a subsequent study.
In order to obtain more reliable results, a number of deficiencies of the new scheme should be addressed in future studies, especially related to the area-altitude distribution of glaciers on a subgrid scale. The current version should be seen as a first attempt to simulate glaciers as an interacting component of the climate system, which can serve as a basis for future efforts. One of the next steps will be to transfer the approach to regions outside the European Alps where the glacier coverage is more prominent and glaciers can be expected to exert a stronger influence on the regional climate and on discharge conditions and might even be an important factor for global sea level rise. From a general point of view, the spatial transferability of the presented approach is given because the most important processes are represented on a physical basis. Also the amount of glacier-specific input data required by the model is comparatively small. Still, basic information on glacier location, extent and preferably glacier volume is necessary for model initialization and for the subsequent evaluation of model performance. In this respect, an encouraging effort is the establishment of a satellite-derived global inventory of land ice within the GLIMS initiative (Global Land Ice Measurements from Space; http://glims.org). Also the availability of observation-based mass balance time series is a crucial element for model validation. An open question is whether parts of the new parameterization scheme will have to be modified or extended in order to be applicable on a global scale. Internal accumulation (the refreezing of surface meltwater in deeper layers of the glacier) might have to be accounted for if the model is applied in sub-polar regions. For very large glaciers and ice fields a parameterization of ice flow between the RCM grid cells and an explicit consideration of changes in the glacier surface altitude has to be introduced. At the same time, subgrid effects become less significant. Future modifications of the subgrid parameterization will put an emphasis on the incorporation of the subgrid variability of orography and on the representation of the area-elevation distribution of glaciers within a climate model grid cell. This includes the split-up of the glacier surface into different tiles for which the surface mass balance and area changes are computed separately. Only by doing so the mass balance gradient with elevation, which strongly determines the climatic sensitivity of a glacier, and changes in the equilibrium line altitude can be explicitly described. Also the spatial variability of the glacier albedo could be incorporated by such an approach (as opposed to the assumption of a uniform glacier albedo in the current setup).
Despite obvious shortcomings of the current parameterization scheme, some degree of confidence has been gained in its ability to reproduce observed glacier changes. This gives the motivation to investigate the development of Alpine glaciers in a climate change scenario using the current model setup and to challenge the validity of the newly developed parameterizations under climate change conditions. This is especially true for the description of the subgrid variability of global radiation and snow accumulation. Both exert a strong influence on the simulated glacier mass balance but are described by simple scaling techniques based on today's glacier extents. A recent study by Zemp et al. (2006) suggests that by 2100 the major part of Alpine ice cover will have disappeared. The question is if such results can be reproduced by REMO glacier .