Model description

We use a simple mass-balance model as described by Reference OerlemansOerlemans (2001) that is based on the calculation of the energy balance. Our version of the model is enhanced by an accurate calculation of mean daily potential global radiation using the DEM25 level 2 (a digital elevation model with 25 m spatial resolution from Swisstopo, representing glacier surfaces in the mid-1990s). The model runs at daily steps, and the cumulative mass balance B _cum on day t + 1 is calculated for every time-step and every individual gridcell as:

where t is the discrete time variable, Δ t is the time-step, L is the latent heat of fusion of ice (334 kJ kg⁻¹) and P _solid is solid precipitation in metres water equivalent (mw.e.). Any meltwater is immediately removed from the glacier, and the energy available for melt (E) is calculated as (see Oerlemans, Reference Oerlemans2001):

where d is a reduction factor for incoming shortwave radiation which accounts for cloudiness or haze, a is the albedo of the surface, Q _e is the clear-sky shortwave radiation, and C0 + C T _a is the sum of the longwave radiation balance and the turbulent exchange linearized around the melting point (T _a is in ˚C). As recommended by Reference OerlemansOerlemans (2001), C ₁ was set to 10 Wm⁻² and we used C0 as a tuning factor (giving best fit at −25 Wm⁻² ). For every day, d is calculated by dividing measured and calculated daily means of global radiation. The program operates with three different albedo values, for ice (ai = 0.3), for snow (a_s = 0.7) and for firn (∝ = 0.4). At the start of the calculation the albedo is set to a _i for the entire test site, and for any snow accumulation a _s is used. Accumulated snow is treated with the firn albedo (a = ∝f) when its age exceeds 365 days and after 2 years its albedo is set to the value of ice ( a = ∝i). The mean daily potential global radiation Qe is calculated beforehand for all 365 Julian days from the SRAD code. This code takes into account the effects of slope, aspect and topographic shading as well as a standard atmospheric composition with a clear-sky transmissivity dependent on elevation for every individual gridcell of the respective DEM (Reference Wilson and GallantWilson and Gallant, 2000).

In order to improve the modelling of snow accumulation, snow redistribution by avalanches is included in the model. The parameterization scheme is based upon the multiple flow-direction algorithm developed by Reference Quinn, Beven, Chevallier and PlanchonQuinn and others (1991) and has been adapted for avalanches by S. Gruber (personal communication, 2005). In the model, the percentage of removed snow is zero for cells with a slope of less than 30° and increases linearly to 100% for cells steeper than 60°. First, the transport fraction from each cell to its neighbours is determined. Then, an initial array of transported snow mass is iteratively propagated along its topographically determined flow path. Finally, for cells with a slope of less than 30° a deposition function describes the amount of snow deposited and removed from the propagated mass. When all snow is deposited, the model iteration is stopped. Avalanches are computed for every day when snowfall at 3000 m a.s.l exceeds 0.02 mw.e.; otherwise snowfall is distributed evenly.

Input data

Claridenfirn index measurements are used to verify a 20 year transient model run starting on 29 September 1981 and ending on 17 October 2002. The test site is 18 km by 14 km in size, with Claridenfirn almost in its centre, and is represented by a rectangular section of the DEM25 (Figs 1 and 2). The highest correlation between mean measured winter accumulation on Claridenfirn and average winter precipitation at climate stations nearby exists for Elm at 960m a.s.l. (Reference Müller and KappenbergerMüller and Kappenberger, 1991). However, we calculate daily precipitation for Claridenfirn from the measurements at Elm (25km to the west) and Engelberg (25 km to the east at 1020 m a.s.l.) to have a spatially more homogeneous coverage. Precipitation is distributed to the terrain by means of an elevational gradient that is calculated based on Reference Müller and KappenbergerMüller and Kappenberger (1991) by matching the mean winter precipitation at Elm (1959–83) with the mean winter accumulation measured at the upper stake on the glacier over the same time period. A linear regression gives a correlation coefficient of 0.77 between measured and modelled winter accumulation at the upper stake, with an underestimation of 0.17 mw.e. by the model. However, matching the precipitation measurements with glacier accumulation is not satisfying. The accumulation on the glacier does not depend only on precipitation but is the result of a complex set of variables, including snowdrift or melt events after the autumn stake readings.

Fig. 2.

Modelled mean distributed specific mass balance for all glaciers within the perimeter of test site 1 , calculated with the 20 year transient model run. Stake locations on Claridenfirn are indicated with circles (L is lower stake, U is upper stake). In the background a shaded relief representation of the DEM25 is shown. DEM25 reproduced by permission of Swisstopo (BA057384).

Daily means of global radiation (to obtain the daily correction factor δ) and temperature (using a constant lapse rate of 0.00625 K m⁻¹) are acquired from the Säntis climate station (2501m a.s.l, 55 km to the northeast; Fig. 1). Despite its long distance from the glacier, this high-mountain station was chosen to avoid extrapolations over large elevation intervals.

Results

The modelled 20 year mean mass-balance distribution for the entire test site is presented in Figure 2. Large bands of avalanche deposits are visible on many glaciers beneath steep mountain slopes, indicating that their dominant source of nourishment is avalanche snow. In general, the modelled mass balances are positive for most glaciers, in contrast to the observed glacier shrinkage in this time period. In view of the simplifications made in the model (the description of the energy balance is very simple, wind redistribution of snow is not included, precipitation is only a function of altitude and cloudiness is a daily constant), the results are nevertheless promising.

To assess the accuracy of the 20 year transient continuous model run, we compare the measurements on Claridenfirn to the modelled values of the corresponding cells of the DEM25. The cumulative mass-balance curves at the positions of the two stakes are depicted in Figure 3a. The final values for both stakes are summarized in Table 1. A more detailed comparison of modelled and measured annual mass balances is depicted in the scatter plot given in Figure 3b. We obtain correlation coefficients of 0.77 for the lower and 0.86 for the upper stake. The mean yearly mass balance is underestimated by 0.09 mw.e. at the lower and overestimated by 0.08 m w.e. at the upper stake (Table 1). As stake measurements are performed at different dates with a varying time-span from year to year, the corresponding time periods are used for the modelled balances. In summary, we find that the simple mass-balance model performs well in the long-term calculation for the two stakes on Claridenfirn, but on the other hand the model run results in positive (up to +1.4mw.e. on some glaciers) and consequently quite unrealistic mass balances for the neighbouring glaciers. Hence, the parameterization of avalanches is an important improvement to the model. Nevertheless, the spatio-temporal distribution of precipitation, as well as the redistribution of snow by wind, are essential to obtain a more realistic mass-balance distribution over large areas and long time-spans. Measurements have shown that even adjacent glaciers have large differences in mean accumulation. In addition, the accumulation pattern of the measured glaciers has revealed a high spatial variability (Reference Machguth, Eisen, Paul and HoelzleMachguth and others, 2006).

Fig. 3.

(σ) Cumulative mass balance for the upper and lower stakes, calculated in a 20 year transient continuous model run from 1981 to 2002. The measured cumulative mass balance is also indicated. The stakes are located well above the equilibrium line, and the measured mass gain does not represent the overall mean yearly mass loss of about 0.7mw.e. observed on Alpine glaciers during the same period. (b) Scatter plot of annual measured and modelled mass-balance values with linear regression trend lines for the lower (2700 m a.s.l.) and upper (2900 m a.s.l.) stakes.

Table 1.

Comparison of modelled and observed values at the two stake locations (2700 and 2900 m a.s.l) on Claridenfirn. The calculation period runs from 29 September 1981 to 17 October 2002

Model description

Curves of cumulative glacier advance and retreat can be converted into time series of temporally averaged mass balance by applying a simple continuity model originally proposed by Reference NyeNye (1960). This approach considers step changes after full dynamic response and new equilibrium of the glacier. Thus, the mass-balance disturbance (db) leading to a corresponding glacier length change (δL) depends on the original length (Lo) and the annual mass balance (ablation) at the glacier terminus (bt): b=bt/L. The dynamic response time (r_r) is hmax/bt (Reference Jόhannesson, Raymond and WaddingtonJohannesson and others, 1989), where hmax is a characteristic ice thickness, usually taken at the equilibrium line where ice depths are near maximum. Assuming a linear adjustment of the mass balance b to zero during the dynamic response, the average mass balance ( b) is roughly assumed to be δb/2. The value so obtained is given in annual ice-thickness change (m w.e.) averaged over the entire glacier surface, and can be directly compared with values measured in the field or modelled values as described in case study 1. The main limitation is the resolution in time. Applying characteristic values for bt and maximum thicknesses of 30 to several hundred metres, the response time is 10−100years (Reference Haeberli and HoelzleHaeberli and Hoelzle, 1995). The calculated mean specific mass-balance values are therefore valid in the range of several decades to a century (Reference Hoelzle, Haeberli, Dischl and PeschkeHoelzle and others, 2003).

For this case study, the mean specific mass balances are calculated for ten selected glaciers of the Swiss glacier-monitoring network according to their individual characteristic response times using the modelled values of bt and db/ dh (Table 2). They are compared to the modelled values of 〈b〉 as obtained with the same distributed mass-balance model as used in case study 1 (Reference MachguthMachguth, 2003). In contrast to the Claridenfirn model run, climatic means of temperature (1961–90) and global radiation (1980–2000) from Meteoswiss are used as input. Seasonal precipitation means (1971–90) are obtained from Reference Schwarb, Daly, Frei and SchärSchwarb and others (2001). The uncorrected precipitation values are increased by 30%. (Reference Schwarb, Daly, Frei and SchärSchwarb and others (2001) recommend a correction of their high-mountain precipitation values by 15–30%.)

Table 2.

Summary of input data (columns 2–9) used for the model, and output data (columns 10–12) calculated by the model (see case study 2), for the glaciers given in the first column

Table 2 shows that the calculated response times are in general longer than 30 years, except for Oberer Grindel-waldgletscher. The response times as calculated here use the b _t values determined by the mass-balance model from case study 1. These b _t values are much lower than the values calculated in earlier studies (Reference Haeberli and HoelzleHaeberli and Hoelzle, 1995), mainly because of shading effects at the tongue of the glaciers that are now considered. Therefore, the resulting response times are in general too long compared to observations. The calculation of input parameters was not carried out automatically as in previous studies (Reference Haeberli and HoelzleHaeberli and Hoelzle, 1995; Reference Hoelzle and HaeberliHoelzle and Haeberli, 1995; Reference Hoelzle, Haeberli, Dischl and PeschkeHoelzle and others 2003), but individually for each glacier (see input data in Table 2).

Results

Figure 4 shows a comparison between the modelled mean specific mass balances based on the length-change data and the modelled mean specific mass balances based on the mass-balance model described above (time period 1961– 90). Although the general agreement is quite good, there are two cases where the correlation is low. The first case is Tiatschagletscher, which has a complex topography with a cascade-like icefall near the tongue. This complex structure at the tongue produces a special length-change behaviour, with a delayed reaction of the tongue during an advance period. The second case is Silvrettagletscher, which is difficult to model due to contradictory precipitation input data from Reference Schwarb, Daly, Frei and SchärSchwarb and others (2001) and from the Hydrological Atlas of Austria (BLFUW, 2003). We have to keep in mind that the simple model approach applied here has several limitations as well, and that the data should be interpreted with care. Nevertheless the results show a new possibility for mass-balance model validation in areas without direct measurements.

Fig. 4.

Comparison between modelled specific mass balance based on the mass-balance model and modelled specific mass balance based on measured cumulative length changes for the period 1961– 90 for ten glaciers from the Swiss glacier-monitoring network.

Model description and input data

In case study 3, we apply a different distributed mass-balance model which is of intermediate complexity with respect to the energy-balance formulation, and calculate the mass balance for Griesgletscher and Ghiacciaio del Basòdino (for location see Figs 1 and 5a) for the extreme summer 2003 conditions (Reference SchärSchär and others, 2004). Model results are compared with field measurements which are obtained by the direct glaciological method on both glaciers (Reference Herren, Hoelzle and MaischHerren and others, 1999). This model, too, is based on the general approach that the most important variables for Alpine glacier mass balance be included with the highest possible precision (e.g. potential global radiation, albedo, precipitation pattern) while others are treated more generally or are parameterized roughly (cf. Reference OerlemansOerlemans, 2001). The model is a combination of the formulations given by Reference OerlemansOerlemans (1992, Reference Oerlemans1993) and Reference Klok and OerlemansKlok and Oerlemans (2002) and adapted for a forcing by meteorological input data from a climate station (temperature, precipitation, global radiation). Distribution of the input data to the terrain is performed by means of elevation-dependent gradients (using the DEM25) and gridded data fields that are calculated beforehand (Paul and others, in press). This includes high-resolution fields of climatological (1971–90 mean values) annual precipitation totals from Reference Schwarb, Daly, Frei and SchärSchwarb and others (2001), as well as a map of glacier albedo obtained from a Landsat Thematic Mapper (TM) scene acquired on 13 August 2003, utilizing a method developed by Reference Knap, Reijmer and OerlemansKnap and others (1999) and described in more detail in Reference Paul, Machguth and KääbPaul and others (2005).

Fig. 5.

(σ) Annual precipitation sums for the test site from the Reference Schwarb, Daly, Frei and SchärSchwarb and others (2001) climatology (1971–90 mean values) resampled to 25 m cell size. Satellite-derived glacier outlines from 2003 (thick white lines) and 200 m elevation contours from the DEM25 (thin light grey lines) are superimposed. The filled white circle (right, centre) marks the location of the Robiei climate station at 1900 m a.s.l. DEM25 reproduced by permission of Swisstopo (BA057384). (b) Measured global radiation (daily means) at Robiei for 2003, together with modelled potential global radiation from SRAD at the same site before and after correction.

The total annual precipitation shows a local maximum near the Robiei climate station (1900 ma.s.l.) with decreasing values towards higher elevations (Fig. 5a). In order to retain the complex precipitation pattern observed throughout the region, we normalized the precipitation data by the value found at the location of Robiei. The daily precipitation at each cell is then obtained by multiplication of the measured precipitation at Robiei with the normalized grid and finally increasing it by 20%, as proposed by Reference Frei and SchärFrei and Schär (1998) for elevations above 1500m a.s.l. After a correction has been applied to the global radiation from SRAD, the cloud factor is obtained for the entire test site from the ratio of the measured to the potential global radiation at Robiei (Fig. 5b).

An albedo of 0.75 is used for freshly fallen snow (with an exponential aging curve) and TM-derived glacier albedo is used if snow depth in the model becomes zero. The model takes into account the extremely low glacier albedo values (<0.2) observed at the end of summer 2003 (Reference Paul, Machguth and KääbPaul and others, 2005). A constant value of 0.1 is added to all albedo values in order to account for effects of the non-Lambertian reflectance characteristics of ice and snow (Reference Greuell and de Ruyter de WildtGreuell and de Ruyter de Wildt, 1999). Some other variables required in the model (e.g. pressure, humidity) use a fixed climatologic mean value at a certain elevation and are extrapolated from the DEM25 for each pixel according to their gradients.

The model starts on Julian day 271 (1 October 2002) with zero snow depth at each cell, and calculates cumulative mass balance for each cell in daily steps until day 635 (30 September 2003) is reached. Snow redistribution by wind and avalanches is not yet included in this model, but seems to be an important factor for the modelled mass balance of many smaller glaciers (see case study 1). Mass balances and mass-balance profiles at 50 m elevation bins are obtained by Geographical Information System (GIS)-based modelling, i.e. intersecting the recent (2003) satellite-derived glacier perimeter (Reference Paul, Kääb, Maisch, Kellenberger and HaeberliPaul and others, 2002) with the DEM25 and the obtained mass-balance distribution.

Results and validation

The modelled mass-balance distributions for Griesgletscher and Ghiacciaio del Basòdino are depicted in Figure 6a and b, respectively, with highlighted glacier areas and approximate positions of the ablation stakes (numbered). While mass-balance distribution on the glacier tongue is governed by the albedo pattern, the potential global radiation exerts a major influence in all other regions. In particular, regions on steep slopes which are prone to heavy topographic shading exhibit positive balances, as avalanches are not included in this model. In consequence, these regions are actually not covered by glaciers (see Fig. 6a and b). The modelled (measured) mean specific mass balances (in mw.e.) for Griesgletscher and Ghiacciaio del Basòdino in the 2002/03 balance year are: –2.63 (–2.52)/ –2.08 (–2.04). These results are promising with respect to the parameterizations made in the model and in view of the fact that no local tuning is applied.

Fig. 6.

Modelled mass-balance distribution for the 2002/03 balance year with stake positions (numbered), 100 m elevation contours from the DEM25 (black) and highlighted satellite-derived glacier areas for the year 2003: (σ) for the region of Griesgletscher, and (b) for the region around Ghiacciaio del Basòdino. DEM25 reproduced by permission of Swisstopo (BA057384). The legend given in (σ) is also valid for (b).

A more detailed comparison of the modelled and measured values is presented in Figure 7a, showing the mass-balance profiles of both glaciers. Apart from the overall good agreement of the two curves for both glaciers, deviations are visible as well. They can be explained partly by the strong smoothing of the measured curves, which is due to the interpolation techniques applied. While the values for Griesgletscher are determined by a regression function (Reference Funk, Morelli and StahelFunk and others, 1997), the 100 m mean values for Ghiacciaio del Basòdino are obtained by manual interpolation (personal communication from A. Bauder, 2005). It is difficult to say which profile is more realistic. The interpolated curves might smooth out local topographic effects (that are clearly visible in the modelled curves) and even fail under extreme conditions (Reference Paul, Machguth and KääbPaul and others, 2005). Assuming that the measured values are correct, the model overestimates ablation in the lower parts of Griesgletscher and underestimates ablation for the upper parts, and the opposite for Ghiacciaio del Basòdino.

Fig. 7.

(σ) Modelled mass-balance profiles sampled at 50 m elevation bins for Griesgletscher and Ghiacciaio del Basòdino, in comparison to measurements (regression curves) in the 2002/03 balance year. (b) Scatter plot of measured vs modelled mass-balance values at the stake positions (see Fig. 6).

In order to avoid any interpolation effects, we also compared the original stake readings with the modelled values at their locations (Fig. 7b). The scatter plot displays a linear regression trend line which indicates a general overestimation of the ablation in the case of Griesgletscher by the model and an under- (over)estimation for the lower (upper) parts of Ghiacciaio del Basòdino. However, considering the large number of local effects controlling the values at individual stakes (e.g. wind drift) and the general assumptions made in our model (e.g. clouds at Robiei = clouds at Griesgletscher), the model results are satisfactory.