1. Introduction
For a long time, glacier variations have attracted the attention of people living in mountainous regions. Variations in the frontal position of a number of large valley glaciers have been reconstructed with the aid of descriptions of historical documents, drawings and paintings, dating of end moraines and, in more recent times, measurements of the distance between the glacier fronts and benchmarks. Some long records of front positions of glaciers in Europe have been summarized in Reference Oerlemans,Oerlemans (1988). Of course, the earlier parts of the records are more uncertain but the extreme positions are fairly reliable (due to the presence of trimlines and end moraines). The general impression is that these glaciers have behaved quite similarly. The most striking feature is the large retreat that has taken place since the middle of the last century. In fact, glacier retreat seems to have been a worldwide phenomenon for the past hundred years or so (e.g. Reference Porter,Porter, 1986; Reference Oerlemans,Oerlemans, 1994). This suggests that the cause of the retreat is a climatic change on a global scale.
Several attempts have been made to simulate historical glacier fluctuations over the last three to four centuries using numerical ice-flow models (Nigardsbreen: Reference Oerlemans,Oerlemans, 1986; Glacier d’Argentière: Reference Huybrechts., deNooze, Decleir and Oerlemans,Huybrechts and others, 1989; Rhonegletscher: Reference Stroeven,, van de Wal, Oerlemans and Oerlenmans,Stroeven and others, 1989; Hintereisferner: Reference Greuell,Greuell, 1992). In these studies, climatic series from nearby stations and/or climatic records based on proxy data were used as input. Simulations for the whole period were not very successful. However, a simulation of the variations in the length of Hintereisferner over the last 100 years turned out to be relatively good (see Reference Greuell,Greuell, 1992).
As Greuell demonstrated, failure of the simulations cannot be blamed on uncertainties and assumptions in the flow model but it is more likely to be due to poor quality of the forcing. The reason for this is two-fold: (i) climatic series and/or climatic reconstructions are not very reliable, and (ii) an inadequate description of the relation between climate and mass balance leads to an erroneous mass-balance reconstruction. Recently, progress has been made concerning this last point and we now have mass-balance models available which are based on the energy balance of the ice-snow surface (Reference Greuell, and OerlemansGreuell and Oerlemans, 1986; Reference Oerlemans,Oerlemans, 1992).
The response of the position of a glacier front to climatic change is the outcome of a process that consists of several steps (see, for instance, Reference Meier., Wright, and FreyMeier, 1965; Reference Paterson.Paterson, 1981; Reference Furbish, and Andrews.Furbish and Andrews, 1984). Modelling this response requires an approach in which the two main modules, the mass-balance model and the ice-flow model, are compatible. The mass-balance model should translate changes in meteorological conditions into changes in specific balance. This serves as forcing for the ice-flow model, which brings about changes in glacier geometry in the course of time. The mass-balance model is coupled to the flow model in a simple way: a polynomial fit is made to the calculated mass-balance profiles and included in the flow model (section 5.2).
The purpose of this study is to simulate the historical length variations in Unterer Grindelwaldgletscher, starting from climatic records. The main reason for choosing this glacier is that it has the longest-known historical record. Moreover, a relatively large amount of climate data is available, including meteorological data from the high-altitude station of Jungfraujoch (since 1938) and from the Grindelwald weather station (since 1914). One drawback is the complicated geometry of the glacier. Another is that there are no mass-balance observations and no velocity or ice-thickness measurements for Unterer Grindelwaldgletscher. The only way to overcome the absence of such information is by calibrating carefully.
In the next section, we describe Unterer Grindelwaldgletscher. A short description of the ice-flow model follows in section 3. The input data for the mass-balance model will be given in section 4. Some sensitivity experiments with the ice-flow model and with the mass-balance model will be presented in sections 5 and 6. Then, we come to the ultimate goal of this paper: the numerical simulation of the historical length variations in section 7. We conclude by investigating the sensitivity of the glacier to greenhouse warming (section 8).
2. Unterer Grindelwaldgletscher
Unterer Grindelwaldgletscher (Figs (1) and 2; Table 1) is situated in the Bernese Oberland, Switzerland (46 °35′N, 8 °05′ E). The distance from the glacier’s head and glacier length are taken along the flowlines shown in Figure 2. In the upper reaches, the glacier consists of two branches, called Fieschergletschcr and Ischmeer. Below about 1700 m, the glacier consists of one branch. In this paper, we shall call this branch together with Fieschergletscher the main stream. The surface slope of the main stream of Unterer Grindelwaldgletseher varies considerably; the maximum slope is about 54% near the top and the minimum about 3% at about 7.5 km from the glacier’s head. The length of the Ischmeergletscher branch is about 6.4 km. Its head is at a height of 3293 m and its average exposure is west to northwest. The ice flux from the Ischmeergletscher branch into the main stream is believed to be significant. The average exposure of Fieschergleischer is east to northeast.
Fig. 1. Photograph of Unterer Grindelwaldgletscher in 1858. From Reference Zumbühl,, Messerli and PfisierZumbühl and others (1983).
Fig. 2. Map of Unterer Grindelwaldgletscher showing present (1987) extent, maximum (1600) length and the central flowlines of the main stream and the Ischmeergletscher branch with grid-point spacing of 100 m. The 13 basins which deliver ice to the main stream and to the Ischmeergletscher branch are indicated by circled numbers.
Table.1. Topography of Unterer Grindelwaldgletscher (data 1987)
The glacier is assumed to be temperate throughout, This assumption is not completely correct but reasonable in view of the presently available ice-temperature data from Alpine glaciers (e.g. Reference Hooke,, Gould and BrzozowskiHooke and others, 1983; Reference Haeberli, and HoelzleHaeberli and Hoelzle, 1995).
During the past five centuries, the length of the glacier has varied between 8.2 and 10 km (Fig.3). The Neoglacial maximum was reached in 1600. From then, up to the beginning of the 19th century, the glacier retreated slightly (on average). Between 1814 and 1855, a new advance occurred; after that, the glacier retreated rapidly, although there were two small advances, one between 1916 and 1934 and the other between 1970 and 1981. The last-published quantitative observation was made in 1983.
Fig. 3. Observed length of Unterer Grindelwaldgletscher from 1534 to 1983 (marked by dots). (From Zumbü and others (1983).)
A big disadvantage of basing a study on Unterer Grin-delwaldgletscher is the absence of mass-balance data and of velocity and ice-thickness measurements. The only data available are climatic data from the Jungfraujoch and Grindelwald stations and topography from maps. Therefore, it is not possible to calibrate the mass-balance and ice-flow models independently of each other; only by coupling the two models can a calibration be performed. The calibration is done by adjusting unknown or uncertain parameters within their range of uncertainly (Table 2). Note that time-dependent simulations have been performed in order to calibrate the coupled models. In the tuning procedure, the root-mean-square difference between the calculated surface elevation of the model glacier in 1987 and the observed surface elevation derived from 1:25 000 topographic maps, issued in 1980 and 1987, has been minimized. The advantage of this “dynamic” calibration is that no assumption of steady state has to be made in order to compare simulated and observed surface elevations. In addition, the root-mean-square difference between the simulated and observed front positions has been minimized by adjusting the annual precipitation before 1865, which was assumed to be constant, and by perturbing the mass balance with a constant value for the whole period.
Table.2. Tuning parameters in the coupled model system
A second disadvantage of choosing Unterer Grindelwaldgletscher is that its geometry does not resemble that of an ideal “model glacier”: it has a large variety of surface slopes and many basins that deliver ice to the main stream and to the Ischmeergletscher branch. Thus, it is questionable whether a one-dimensional ice-flow model can adequately describe its flow.
However, since the ice-flow model and mass-balance model have been tested thoroughly for other glaciers, for which more measurements are available (see, e.g., Reference Oerlemans,Oerlemans, 1992), it is worth trying to simulate the historical front positions of Unterer Grindelwaldglelscher by coupling the two types of model.
3. Description of the Ice-Flow Model
In the ice-flow model the dynamics of a glacier are calculated along a central flowline. The model is similar to the flow models used by Reference Oerlemans,Oerlemans (1986) and Reference Greuell,Greuell (1992).
The model involves: (1) parameterization of the three-dimensional geometry; (2) the continuity equation; (3) calculation of ice velocities; and (4) calculation of fluxes from the Ischmeergletscher branch into the main stream and from the smaller basins into the main stream and the Ischmeergletscher branch.
3.1. Parameterization of the Three-Dimensional Geometry
The model is basically one-dimensional (flowline along the x axis) but the three-dimensional geometry is implicitly taken into account by the parameterization of the cross- sectional geometry at each grid point. The cross-sectional profile is assumed to he trapezoidal. It is determined by four parameters (Fig.4); the bedrock height b, the width of the valley bottom W, and the complements of the slope angles on the right (αR) and left (αL) sides of the valley. The parameters αR, αL and the valley width at the surface Ws were taken from topographic maps in such a way that the area elevation distribution is not distorted too much.
Fig. 4. Parameterization of the cross-sectional geometry. (Adapted from Reference Greuell,Greuell (1992).)
The shape of the glacier itself is then determined by the bed geometry and the ice thickness H. Because there are no measured ice-thickness profiles for Unterer Grindelwaldgletscher, a preliminary ice thickness was estimated assuming a constant basal stress τb:
Here ρ is the ice density.(900 kg m−3), g is the acceleration of gravity and h is the surface elevation. The latter was also taken from the topographic maps. The bedrock height b then can be computed by subtracting the estimated ice thickness from the observed glacier-surface height at each grid point.
The assumption used in Equation (1) is that the product of ice thickness and surface slope is constant; this would be true if ice deformed as a perfectly plastic material. The value of 1.5 bar was found after a tuning procedure, with the model. In the timing procedure the root-mean-square value σ was minimized:
Here, h c, is the calculated height of the model glacier surface, h t is the observed (1987) surface height (read from the topographic maps), i is the index of the grid points and N is the total number of grid points.
However, the constant basal shear-stress assumption is problematic. Reference Haeberli, and SchweizerHaeberli and Schweizer (1988) showed with numerical calculations, using historical data from Rhone-gletseher, that stresses must generally be higher in Steep zones than in flat areas in order to enable balanced flow. Therefore, it is justified to adjust the bedrock in the following way. The bedrock was smoothed by using a five-point filler and the smoothed bedrock was then adjusted further manually until the model produced a glacier surface close to the 1987 profile. This reduced σ further.
Expressions for the width of the valley bottom W and the cross-sectional area of the glacier S are:
3.2. The Continuity Equation
The dynamic behaviour of the glacier is described in terms of changes in ice thickness. These have to be constant, so a conservation equation for ice density is taken to be constant, so a conservation equation for ice volume can be used:
In Equation (5), U is the mean ice velocity in the cross-section and B is the annual mass gain. Thus, the first term on the righthand side represents the divergence of the volume flux and the last term stands for net accumulation at the surface.
Substitution of Equation (4) into Equation (5) yields:
where μ = 0.5(tan(αR) + tan(αL)).
3.3 Calculation of Ice Velocities
The depth-averaged ice velocity (U) is determined by internal deformation (U d) and sliding (U s). For the calculation of (U), the following equations are used (Reference Budd,, Keage and BlundyBudd and others, 1979; Reference Paterson.Paterson, 1981):
Here, τd is the driving stress and f 1 and f 2 are flow parameters. Note that the contribution that the basal water pressure makes to the effective basal normal stress is neglected in this study, although it is recognized that this pressure can affect the sliding velocities, particularly in summer (e.g. Reference HodgeHodge, 1976; Reference Budd,, Keage and BlundyBudd and others, 1979). However, in the light of other approximations made, a more refined treatment was not attempted.
The flow parameters f 1, and f 2 are not known accurately. They depend on bed conditions, debris content and crystal structure of the basal ice layers. Therefore f 1 and f 2 have been used as tuning parameters. The values adopted here are:
These values for the flow parameters are within the range of the values used in similar studies (e.g. Reference Budd,, Keage and BlundyBudd and others, 1979; discussion in Reference Greuell,Greuell, 1989).
3.4. Calculation of fluxes from the Ischmeergletscher branch and the smaller basins into the main stream
On the valley walls of Unterer Grindelwaldgletscher (above the ice streams) there are 13 ice basins, which deliver ice to the main stream and to the Ischmeergletscher branch (Fig.2). Their contribution has to be taken into account, as these basins gradually add part of their volume to the main stream. To calculate the magnitude of the contribution, the following algorithm was developed.
First, the ice volume of the basin V bas(t) is calculated according to:
Here A(hk ) and B(hk ) are respectively the area and the mass balance of the elevation interval (100 m in this study) centred around hk ; k is the number of the elevation interval. The first term on the righthand side thus represents the “volume gain” of the basin (per year). Furthermore, a characteristic time-scale for the ice Ilow (t*) has been introduced. The reason for this is that the mass flux due to ice flow from the slopes may lake several years. The time-scale for the ice flow is defined as follows:
Here, L is a characteristic length scale which is chosen to be half the length of the tributary basin and u* is a characteristic ice velocity which depends on the mean ice thickness and the mean surface slope of the basin (analogous to Equation (7)). Therefore, the value for t* differs from basin to basin. Here, it varies from 2 years for a short basin with a large mean surface slope to 32 years for an oblong basin with a small mean surface slope (Table 3). Note that area, length, mean ice thickness and mean surface slope of the tributary basins are taken from topographic maps ( perfect plasticity is assumed in order to derive ice thickness; Equation (1)). The second term of Equation (9) therefore represents the volume loss of the basin (per year due to the mass flow from the basin into the main stream.
Table.3. Parameters the calculation of the ice volume in the basins. The elevation interval is 100 m. The basin numbers refer to the numbers in Figure 2
Equation (9) can easily be put into an incremental form following a forward differention scheme:
Here, n is the number of the time step. From Equation (11), it is clear that the volume loss of the basin, ΔV bas, due to the mass flow from the basin into the main stream is:
The volume added to the main stream ΔV main stream has to be converted into added ice thickness ΔV at a certain grid point according to (Fig.4):
Here, ΔS is the change in cross-sectional area at a certain grid point, j is the number of grid points minus 1, over which the added volume has to be divided, and Δl is the grid point interval. It may be obvious that the determination of the value of j is subjective.
Now, we introduce the ice flux Φ, which is defined as:
The ice flux from the Ischmeergletscher branch into the main stream is treated as follows. The ice flux between the two last grid points of the lschmeerglctscher branch is calculated according to Equation (14). We assume this is a good approximation of the ice flux from the Ischmeergletscher branch into the main stream. This flux is multiplied by Δt to yield the ice volume that is added at every time step. If the elevation of the glacier surface is higher in the Ischmeergletscher branch than in the main stream, the volume will be added to the main stream.
Again, the volume added to the main stream has to be converted into added ice thickness ΔH at a certain grid point according to Equation (13). Analogous arguments hold for the (equal) ice volume that is removed from the Ischmeergletscher branch.
3.5 Numerical Details and Stability of the Scheme
In the flow model, a staggered grid is used to evaluate the ice flux divergence dΦ/dx along the flowlines. The ice-flow model uses an explicit scheme for time integration, implying that the time step has to be rather small. We have determined experimentally the maximum time step for which the scheme is stable and the solution is accurate (i.e. independent of the time step). The maximum time step is of the order of 10−2 year. In this study, a time step of 0.02 year is used.
4. Input Data for the Mass-Balance Model
The ice-flow model is forced by a numerical mass-balance model, which is based on the energy balance of the ice/snow surface and on the accumulation rate. As the mass-balance model employed here is that used by Reference Oerlemans,Oerlemans (1992), we refer to that paper for a description of it.
The following meteorological input data are needed to run the mass-balance model (Table 4): temperature and humidity at standard measuring height, cloudiness and precipitation. Air temperature T is generated by assuming sinusoidal shapes for seasonal and daily variations:
Here, M T is the annual mean temperature reduced to sea level, γM is the lapse rate, Nday is the Julian day number and t is local time in h. The precipitation rate was kept constant through the year but was allowed to depend on altitude.
Table.4. Input parameters for the mass-balance calculations. When the values differ between the main stream and the Ischmeergletscher branch, a* is used tn indicate the main stream and a† to indicate the Ischmeergletscher branch. The meteorological input data (temperature, humidity and cloudiness) are from the high-altitude meteorological station of Jungfraujoch. They can be found in Schweizerische Meteorologische Anstalt (1865-1992)
The largest uncertainty generally concerns the altitudinal gradients of the input parameters, especially precipitation. To determine the altitudinal gradient of precipitation, we used precipitation data from the weather station of Grindelwald (1040 m a.s.l.), as well as measured precipitation (from totalizing rain gauges) and accumulation data for Jungfraufirn (from Reference Braun, and AellenBraun and Aellen, 1990). Because of the large uncertainly, we used this gradient as a tuning parameter in the mass-balance model by coupling the latter to the ice-flow model. It turned out that different altitudinal gradients of precipitation had to be taken for the main stream and for the Ischmeergletscher branch (Table 4) in order to yield reasonable results. In other words, reasonable results are obtained only if the equilibrium-line altitudes are different for the main stream and for the Ischmeergletscher branch (see section 5.1). For the lapse rate, we used a value of-0.007 K m−1.
The other meteorological input data in Table 4 are from the high-altitude meteorological station of Jungfraujoch (3580 m a.s.l.), which is situated a few kilometres southwest of the glacier’s head. The data from the Grindelwald and Jungfraujoch stations can be found in Schweizerische Meteorologische Anstalt (1865-1992).
The slope and exposure at each grid point are read from the topographic maps. The average exposure for the main stream is northeast and for the Ischmeergletscher branch is west-northwest. The mean slope of the main stream is 20% and the mean slope of the Ischmeergletscher branch is 15%.
5. Sensitivity Experiments with the Mass-Balance Model
5.1 Reference Calculations
We now wish to explore how sensitive the mass-balance model is to variations in model parameters. Using the input parameters from TabIe 4, we calculated the reference mass-balance profiles for both the main stream and the Ischmeergletscher branch (Fig.5). The differences between the two curves result from the fact that the altitudinal gradients of precipitation had to be unequal, as mentioned above, in order to obtain agreement between model and observation, and from the fact that surface slope and exposure are different for the main stream and for the Ischmeergletscher branch. The main stream has an equilibrium-line altitude of 2750 m. whereas the equilibrium-line altitude of the Ischmeergletscher branch is 2797 m. As no mass-balance measurements have been made at the Unterer Grindelwald-gletscher, we cannot compare these results with observations. However, we have some confidence in the result, because the meteorological input from the Jungfraujoch station is very detailed. On the other hand, the uncertainty in the altitudinal gradients is large, especially with regard to precipitation.
Fig. 5. Calculated reference mass-balance profiles for the main stream and the Ischmeergletscher branch. The profiles are unequal due to differences in exposure, surface slope and the altitudinal gradient of precipitation. Parabolic fits to the curves an indicated by dashed lines.
5.2 Sensitivity of Mass-Balance Profiles to Climatic Change
There are many model constants and parameters that can be varied but we will restrict our sensitivity analysis here to changes in (seasonal, i.e. spring/simmer/autumn/winter) temperature and precipitation. In this context, the word sensitivity is used to denote the change in mass balance or equilibrium-line altitude for a certain change in climatic conditions.
The results of a series of experiments for the main stream, in which these quantities were systematically varied, are summarized in Figure 6. The results for the Ischmeergletscher branch are similar and are therefore not reproduced.
Fig. 6. Sensitivity of the mass balance of the main stream to changes in seasonal temperature and precipitation, as indicated by the labels. Parabolic fits to the curves are indicated by dashed lines.
It is difficult to explain all the details, because of the many feedback mechanisms involved. Hence, we will give a short description of the most important features. The general impression we get from Figure 6 is that the response to a change in a climatic parameter is more or less asymmetric around the reference state.
Concerning temperature, it is obvious from Figure 6 that the mass balance is most sensitive to a change in summer temperature, as expected. However, sensitivities to changes in spring and autumn temperature are also relatively high, and therefore cannot be neglected. The same applies to changes in winter temperature for the lower elevations; (under 2000 m), where the sensitivity is comparable to the sensitivity to changes in spring and autumn temperature. This is because some melting takes place during the winter at these elevations. At the higher elevations, the mass balance does not change because temperatures in winter are lower and the fraction of precipitation falling in solid form remains the same regardless of the prevailing winter temperature. For temperature changes in spring, summer and autumn, the change in mass balance is small at high elevations but increases significantly (especially for changes in summer temperature) near the equilibrium-line altitude (due to the albedo feedback); at lower elevations, the change in mass balance tends to decrease slightly.
With regard to precipitation, the mass balance is most sensitive to a change in winter precipitation, as expected. The sensitivity to a change in summer precipitation is the lowest and is only comparable to changes in precipitation in the other seasons at elevations above 3000 m. The overall precipitation picture can be described as follows. At the highest elevations, virtually all precipitation falls as snow and the change in mass balance simply equals the change in precipitation. At somewhat lower altitudes, the change in mass balance tends to decrease until the equilibrium line is reached. There, a (sharp) increase occurs due to the albedo feedback (especially for changes in spring and winter precipitation). At lower elevations the change in mass balance decreases rapidly to (almost) zero. At elevations of 1300-1500 m, only a small fraction of the annual precipitation falls as snow. At these elevations, therefore, changes in precipitation hardly affect the mass balance.
The sensitivity of the Ischmeergletscher branch is compared with that of the main stream in Table 5. Changes in the equilibrium-line altitude E, due to perturbations in annual mean temperature, mean summer temperature and annual precipitation are shown. The values are averages of the positive and negative perturbation experiments. The results for the two branches do not differ very much, although Ischmeergletscher seems to be the more sensitive. The mean value dE/dT s is 50 m K−1 which is more than half the mean value dE/dM T (92 m K−1), implying that a change in mean summer temperature has a relatively large effect on the change in equilibrium-line altitude compared to a change in annual mean air temperature. The mean value of dE/dP is 4.5 m %−1. Thus, a 20% change in annual precipitation has roughly the same effect on the equilibrium-line altitude as a 1 K change in annual mean air temperature.
Table.5. Sensitivity of equilibrium-line altitude, E, with respect to changes in annual mean temperature MT, mean summer temperature Ts and annual precipitation P
Like Reference Oerlemans, and HoogendoornOerlemans and Hoogendoorn (1989), we conclude that the response of a mass-balance profile to climatic change cannot be expected to be independent of altitude. The results in Figure 6 suggest, further, that the sensitivity to changes in seasonal temperature and precipitation is so large that all seasonal variables have to be included in an expression for the mass balance. Moreover, the sensitivity of the mass balance to annual mean temperature is found to be highest at lower elevations. This means that the glacier front will be very sensitive to climate change (Reference Oerlemans,Oerlemans, 1992). We shall investigate this point further in the next section.
The mass-balance model is coupled to the flow model in the following way. First, we fit the calculated mass-balance altitude profiles with a parabolic relation, B = a 1 + a 2 h + a 3 h 2, and then fit the coefficients a 1, a 2 and a 3 with parabolic functions of a perturbation in a climatic parameter dp c, thus: a 1 = b l + b 2dp c + b 3dp c 2. The contributions from the different climatic parameters can simply be added for a certain coefficient. Finally, we included the resulting expression in the ice-flow model. Thus, the feedback from changing surface elevation to the mass balance is taken into account (Reference Oerlemans,Oerlemans, 1992). The parabolic fits to the calculated mass-balance profiles are shown in Figure 6.
The advantage of this approach is that the mass-balance model need not be re-run when the glacier geometry and the climate have changed. Climatic records are needed as input for the ice-now model in the case of climate forcing (section 7). A disadvantage is that the feedback from changing glacier geometry to the absorbed short-wave radiation is not accounted for. This effect may be important but, in view of other approximations made and the larger computation time that would be required, this effect is neglected.
6. Basic Sensitivity Experiments with the Ice-Flow Model
6.1. Steady-State Length as a Function of Seasonal Temperature and Precipitation Anomalies
We now wish to explore how sensitive the ice-flow model is to variations in input parameters. In Figure 7a, the steady-state length of the glacier is depicted as a function of an anomaly in the seasonal air temperature. It will be seen that, for example, a 1 K increase in summer temperature would lead to a reduction in steady-state length from 8.75 km to 7.15 km. Note that the sensitivity of glacier length to a change in summer temperature is highest and the sensitivity of glacier length to a change in winter temperature is almost zero. Anomalies in spring and autumn temperatures affect the steady-state length in a similar way, although the sensitivity to a change in spring temperature is somewhat higher. This is consistent with Figure 6.
Fig. 7. Steady-state length of the main stream vs anomalies in: (a) seasonal mean air temperature and (b) seasonal precipitation, as indicated by the insets.
The greatest sensitivity shown in Figure 7a is to decreases in summer temperature in excess of 0.6 K. This is a result of the elevation mass-balance feedback. Since the slope of the bed decreases significantly when the glacier is longer than approximately 9.75 km and thus reaches the valley floor (Fig.8a), the very effective mechanism of the elevation-mass-balance feedback causes ice thickness to grow exponentially (see also Equation (1)). Therefore, the model glacier length shows a large increase as well.
Fig. 8. Simulated (1987) ice-thickness profiles for: (a) the main stream and (b) the Ischmeergletsrher branch of the model glacier. They are compared with the observed (1987) surface profile of Unterer Grindelwaldgletscher. Note the difference in scales.
If the extreme glacier stands in 1600 and 1970 (Fig.3) had been steady states, it follows from Figure 7a that they would have corresponded to a difference of 0.9 K in the mean summer temperature.
If the steady-state length were controlled entirely by seasonal precipitation, the change in length with change in precipitation would be as shown in Figure 7b. The length is most sensitive to changes in winter precipitation and least sensitive to changes in summer precipitation, as expected. The sensitivity to anomalies in spring precipitation is higher than to changes in autumn precipitation. This is due to the albedo feedback: high amounts of snow in spring cause the albedo to increase, as a result of which the ablation season starts later.
We can conclude that the sensitivity of glacier length to a change in climatic conditions is higher when the terminus is on a slightly inclined part of the bed (Figs (7) and 8a; Reference Oerlemans, and Oerlemans,Oerlemans, 1989). One should bear in mind that, in spite of the high sensitivity of glacier length to a change in climatic conditions, the relative changes in ice volume are much smaller since the bulk of the glacier is above 2500 m.
In view of the results reported in this section, it appears that the historical variations in the length of Unterer Grin-delwaldgletscher, which are of the order of 1.8 km (Fig.3), can be explained by relatively small climatic changes.
6.2 Response Times
The response time, τ, is defined as the time that elapses between a perturbation and the time when the glacier reaches a length of:
Here, L 0 is the initial steady-state length and ΔL is the difference between the initial and the final steady-state lengths.
If one starts with zero ice volume, our modelling suggests that it takes 150-450 years for a steady state to be reached. This (theoretical) growth time is significantly larger than the response time. This result demonstrates that a simulation of the historical record must start at a time well before the beginning of the observed record.
We next investigated the response time by perturbing the reference steady state with various (abrupt) changes in the mass balance. Figure 9 shows an example of the response (in terms of glacier length) after various perturbations in the mean summer temperature are imposed on the flow model as step functions. The response times range from 34 to 45 years. The larger values correspond to negative changes in the mass balance (i.e. positive perturbations of temperature).
Fig. 9. Reaction of the glacier-front position to a stepwise change in mean summer temperature of given magnitude: the response time τ is given in years.
7. Simulation of the Historical Length Variations
7.1. Records used for the Reconstruction of Climatic Data
A continuous precipitation record is available from the weather station of Grindelwald for the period from 1914 onwards. The temperature record from the high-altitude meteorological station of Jungfraujoch is available only from 1938 onwards. Therefore, the records for these locations for the period before these dates (Fig.10) have been reconstructed using records from other stations and climatic data for Switzerland based on proxy data (Reference Pfister,Pfister, 1984). The climatic records from the Beatenberg, Grindelwald, Jungfraujoch and Säntis stations (Schweizerische Meteoro-logische Anstalt, 1865-1992) and from the Basel station have been used.
Fig. 10. Records used for the reconstruction of: (a) seasonal mean temperatures for Jungfraujoch and (b) seasonal sums of precipitation for Grindelwald.
7.2. Reconstruction of Climatic Data
7.2.1 Reconstruction of the Seasonal Mean Air Temperatures
We reconstructed the seasonal air temperatures for Jungfraujoch for the period 1755-1882 using the seasonal air temperatures for Basel. We did this by calculating four linear regression equations, one for each season, between the seasonal air temperatures of both stations from the 49 years of common measurements (1938-86). The correlation coefficients between the Basel record and the Jungfraujoch record are in the range 0.66-0.88. The Basel record was used for the reconstruction because it was the longest one.
We used the data from the high-altitude meteorological station of Säntis (2490 m a.s.l.) to reconstruct the seasonal air temperatures for Jungfraujoch for the period 1883-1937. From the 55 years of common measurements (1938-92), we computed regression equations between the seasonal air temperatures for the stations Jungfraujoch and Säntis. Our main reason for choosing the Säntis record for the reconstruction was that the correlation coefficients turned out to be high (0.83-0.97).
For the period 1530-1754, we reconstructed 10 year mean values of seasonal air temperature for Jungfrattjoch (extended with the reconstruction for 1755-1937) using 10 year mean values of seasonal air temperatures for Switzerland (Pfisler, 1984). Again, four linear regression equations between the 10 year mean values of seasonal air temperatures for Switzerland and Jungfraujoch were computed. The correlation coefficients are in the range 0.79-0.87.
Because we need temperature anomalies as a forcing for the flow model, the mean values of the seasonal air temperatures for the period 1951-80 were subtracted from the (partly reconstructed) seasonal temperature records for Jungfraujoch.
7.2.2. Reconstruction of Seasonal Precipitation
It is more difficult to reconstruct precipitation records, for a number of reasons (Reference Greuell,Greuell, 1992): precipitation is correlated over much shorter distances than temperature; records are more liable to contain inhomogeneities due to changes in surroundings and equipment, and long and well-investigated records are scarce.
We began by collecting precipitation records from stations situated in the vicinity of Grindelwald. The precipitation record of Beatenberg (1150 m a.s.l.) turned out to have the highest correlation with the Grindelwald record. Therefore, this record was used to reconstruct the seasonal precipitation record of Grindelwald for the period 1865-1913. The correlation coefficients are in the range 0.83-0.91.
No measured precipitation data were available for the period 1530-1864. Therefore, an attempt was made to reconstruct 10 year mean values of precipitation for Grindelwald using 10 year mean values of precipitation for Switzerland (Reference Pfister,Pfister, 1984). However, the correlation coefficient for annual precipitation obtained in ihis way was too low (0.49). Thus, the annual precipitation during the period 1530-1864 was assumed to be constant and was adjusted until the root-mean-square difference between the simulated and observed front positions reached a minimum. The figure, obtained in this way, appeared to be equal to the mean annual precipitation for Grindelwald during the period 1865-1978.
Because precipitation anomalies are needed as a forcing for the flow model, we calculated them using the mean values of the seasonal precipitation for the period 1951-80. The anomalies are expressed in per cent relative to the respective mean values.
Figure 11 depicts the seasonal temperature and precipitation forcing functions (together with smoothed values) that were obtained. The most striking feature of the temperature curves in Figure 11a is the pronounced warming that has occurred in the present century. As far as the precipitation curves are concerned, it appears from Figure 11b that seasonal values reached a maximum in the 1970s and 1980s.
Fig. 11. Variation in seasonal mean temperature at Jungfraujoch (reconstructed from 1530-1937). b. Variation in seasonal precipitation in Grindelwald (reconstructed from 1865-1913); a constant annual precipitation value has been used for the period 1530-1864 (see text). The plots also show smoothed values.
7.3. Result of the Simulation
As noted earlier, because of the long response time of the glacier, the integration should start well before the beginning of the observed record of front variations. Thus, we started the model in the year AD 1000 with zero ice thickness.
For the period 1000-1992, the model was forced with (dTs = 0.4 K, dT sp = dT au = dTw = 0 K and dP = 0% (all relative to the 1951-80 mean). For the period 1530-1992, the seasonal mean temperature record Jungfraujoch and seasonal precipitation record for Grindehvald (Fig.11) were used as additional forcings. This yielded a steady state with a length of 9.55 km for the period 1000-1529 (which is slightly longer than the median of the length observations for the period 1534-1983).
As an illustration of the type of temporal variation in the mass balance that resulted from this climatic forcing, curves of mass-balance anomalies for two different surface elevations (1500 and 3000 m ) are shown in Figure 12a.
Fig. 12. a. Mass-balance anomalies for two different surface elevations (1500 and 3000 m), resulting from the climatic forcing; the anomalies are relative to the corresponding steady-state mass-balance values for the period 1000-1529. Smooth-curve fits are also indicated, b. Simulation of historical variations in the front of Unterer Grindelwaldgletscher (main stream) compared with observations. Forcing functions are seasonal temperature and precipitation anomalies.
Figure 12b shows the result of a model simulation for the period from 1530 onwards, with the seasonal mean temperature and precipitation anomalies as forcings. From the smooth curve fits in Figure 12a, we can deduce that the mass-balance anomalies reach maximum values around 1840 and minimum values around 1950. For comparison, the model glacier reached its maximum length around 1870 and its minimum length around 1980 (Fig.12b), suggesting a 30 year lag. The simulation is generally consistent with the observed retreat during the late 1800s and early 1900s, with the observed interruption of this retreat by an advance between 1920 and 1940 and with the observed advance from 1970 onward. However, there is a 1-2 decade lag between the observed behaviour and the modelled one. The root-mean-square difference between the simulated and observed front positions is 0.28 km. The observed total retreat of the main stream between 1534 and 1983 was about 1.0 km, whereas the model calculates a retreat of 1.3 km.
On the whole, the response time of the model glacier is somewhat larger than the response time of Unterer Grindelwaldgletscher (with a difference of about 10 years). An example is the simulated maximum around 1860 which lasts longer than the observed maximum around that time. This is probably due to the fact that there are no records for precipitation prior to 1865 and therefore precipitation cannot be reconstructed. Because the model glacier has a long response time, it does not react rapidly to the large negative precipitation anomalies which lead to a retreat.
The magnitudes of the simulated length variations generally correspond to those of the observed front variations. One exception is the magnitude of the simulated advance between 1570 and 1600 which is half the magnitude of the observed advance which led to maximum lengths at the beginning of the 17th century. Another exception is the magnitude of the simulated advance between 1920 and 1940, which is twice as large as that of the observed advance between those years. A possible explanation for these discrepancies is that the mass-balance reconstruction is incorrect for these periods. This could well be the case for the first period, for which there is a lack of detailed climatic data. The mass-balance anomalies around 1570 are indeed smaller than the anomalies around 1840 (Fig.12a). Another possible explanation for the poor quality of the simulation in the second period is the occurrence of ice avalanches from an isolated tributary glacier (i.e. point 2 in Fig.2). These ice avalanches fill parts of the gorge and cause an advance of the tongue. Thus, they have nothing to do with the dynamics of the main glacier. The re-advances around 1920 and 1980 as calculated for the model may, there-fore, overestimate the dynamic response of the main glacier (personal communication from W. Haeherli and H. Gudmundsson, 1996).
In Figure 8, the simulated surface elevation of the model glacier in 1987 is compared with the observed surface elevation of Unterer Grindelwaldgletscher (1987 glacier stand). We are not unduly surprised to find that the simulated ice-thickness distributions correspond well to the observed shape, because the model was tuned with respect to the observed surface elevation. The root-mean-square error (defined by Equation (2)) is 13 m for the main stream and 14 m for the Ischmeergletscher branch. The length of the model glacier in 1987 is 8.25 km, which is 0.2 km shorter than observed. Despite these slight anomalies, it can be coucluded that the “dynamic” calibration procedure worked well.
The simulation indicates that, on the whole, the coupled model system, as well as the reconstructions of seasonal temperature and precipitation, are reasonably accurate and realistic.
8. Prediction of Future Front Variations
In this section, we will investigate how large the retreat of Unterer Grindelwaldgletscher is likely to be in the next century if two different greenhouse-gas emissions scenarios are applied. These scenarios are taken from the report of the UN Intergovernmental Panel on Climatic Change (IPCC; Reference Bruhl,, Houghton,, Jenkins and EphraumsBruhl and others, 1990). One category of emission scenarios, which is referred to as policy scenarios, represents a broad range of possible controls to limit the emissions of greenhouse gases. Estimates of the change in global mean surface air temperature due to man-made forcing were made using a box-diffusion-upwelling model (Reference Bruhl,, Houghton,, Jenkins and EphraumsBruhl and others, 1990). In this experiment, estimates of the temperature rise for the two extreme policy scenarios BaU Business-as-Usual) and D (early controls on greenhouse-gas emissions, are used.
From 1993 on, the ice-flow model is forced by the predicted temperature scenarios (values relative to the 1951-80 mean), whereas the seasonal precipitation values are set at the mean values of 1951-80. For the period prior to 1993, the forcing functions used are the same as those used in section 7.3. It is assumed that the predicted temperature scenarios refer to changes in annual mean air temperature.
As an illustration of the type of temporal variation in the mass balance that resulted from the climatic forcing, curves of mass-balance anomalies for two different surface elevations (1500 and 3000 m) and for the two scenarios BaU and D (from 1993 on) are shown in Figure 13a. From 1993 onwards, the mass-balance variations are due entirety to temperature variations, so they are larger at 1500 m than at 3000 m (in accordance with Figure 6).
Fig. 13. a. As Figure 12a, but now extended with mass-balance anomalies, calculated from estimates of temperature rise, based on the tun (extreme) IPCC “policy” emissions scenarios BaU and D (Reference Bruhl,, Houghton,, Jenkins and EphraumsBruhl and others, 1990). b. Prediction of future length (solid lines) and volume (dashed lines) of Unterer Grindelwaldgletscher based on the two different scenarios, BaU and D.
The model predicts that the glacier length will continue to increase until a maximum of 9.25 km is reached in 2017 (Fig.13b). This initial increase is a result of the high amounts of precipitation in the 1970s and 1980s (Figs 11b and 13a). Note that the length response time is larger than the volume response time. It is predicted that after 2017 the glacier will show a continuous retreat. By 2100, it will have a length of 4.95 km in the case of the BaU scenario and a length of 6.65 km in the case of scenario D. Thus, the total magnitude of the retreat (1983-2100) will range from 1.6 to 3.3 km. The total volume decrease in this period will range from 0.57 to 1.25 km3.
Note that several assumptions have been made. First, it is assumed that the enhanced greenhouse effect will not alter the total amount of precipitation. It should be emphasised that most climate models show an increase in global mean precipitation with increasing global mean temperature due to the fact that a warmer atmosphere will contain more water vapour. However, because climate models are still unable to predict the regional distribution of precipitation and because we want to get an idea of the degree of glacier retreat due to a temperature rise alone, we have assumed here that total precipitation will not change. If the amount of annual precipitation were to increase, this would lead to more accumulation in the higher parts of the accumulation zone (Fig.6). This accumulation would then counteract the effect of increased melting in the ablation zone. Total retreat would therefore be less.
Secondly, it is assumed that local air temperatures will increase by the same amount as global mean air temperatures. Because the regional temperature distribution is also unknown in the current state of climate modelling, this is the best assumption one is able to make.
9. Conclusions and Discussion
In this study, the variations in length of Unterer Grindelwaldgletscher which have occurred in the past have been simulated by coupling a mass-balance model to an ice-flow model. As climate forcing, we used (partly reconstructed) seasonal temperature and precipitation records. The root-mean-square difference between simulated and observed front positions turned out to be 0.28 km. Total retreat was computed fairly accurately (error of 0.3 km). Note that the mass-balance variations themselves were not tuned to the resulting length variations. They were determined only by the mass-balance model, the polynomial fits to the calculated mass-balance profiles and the climatic data. Further-more, the simulated glacier geometry for 1987 fits the observed geometry for that year reasonably well. Therefore, we can conclude that the “dynamic” calibration worked well.
Thus, the simulation of the historical variations in the length of Unterer Grindelwaldgletscher is successful and in fact better than the simulations for Nigardsbreen (Reference Oerlemans,Oerlemans, 1986), Glacier d’Argentière (Huybreehts and others, 1989) and Rhonegletscher (Reference Stroeven,, van de Wal, Oerlemans and Oerlenmans,Stroeven and others, 1989). This is probably due to the good quality of the climatic series and reconstructions and the fact that the relation between climate and mass balance was handled using a mass-balance model based on the energy balance of the ice-snow surface.
Furthermore, we conclude that glaciers with a complicated geometry can also be modelled by calculating explicitly the fluxes from smaller basins into the main stream, using a method that includes the residence time of the ice in these basins.
Since the simulation of the historical front variations was a success, future length and volume variations can be predicted. The model experiments confirm that the glacier is very sensitive to greenhouse warming. For a BaU scenario (from Reference Bruhl,, Houghton,, Jenkins and EphraumsBruhl and others, 1990) only 29% of the 1990 ice volume would remain in 2100.
Acknowledgements
We are grateful to W. Greuell for providing the climatic record from Basel and for his useful comments on tin earlier version of this paper. Furthermore, we wish to thahk W. Haeberli, H. Gudmundsson, J. L. Fastook and R. LeB. Hooke for comments that helped us improve the manuscript.
