1. Introduction

The retreat of mountain glaciers around the world is a striking signature of climate change, and it is likely that this retreat is the consequence of a worldwide increase in temperature. However, large regional differences in the behaviour of glaciers are seen. These differences originate partly from different climatic conditions, but there is no doubt that glacier geometry has a strong effect on the climate sensitivity and response time of individual glaciers (Reference Oerlemans, Lisse and BalkemaOerlemans, 2001).

Traditionally, glaciers have been considered as complex systems from which it is difficult to obtain quantitative information about past climates. Recent advances in the modelling of glaciers have changed the picture. A suite of models now exists that allows one to determine climate sensitivity and response time for a glacier of given geometry, ranging from simple analytical considerations (e.g. Reference Oerlemans, Lisse and BalkemaOerlemans, 2001) to full-system glacier flow models (Reference Leysinger Vieli and GudmundssonLeysinger Vieli and Gudmundsson, 2004).

An early attempt to estimate the response time of a glacier to climate forcing was made by Reference NyeNye (1965), based on kinematic wave theory applied to a linearized equation for glacier flow. Subsequently, numerical models based on the shallow-ice approximation (SIA) have been used to study the transient response of glaciers to climate change (e.g. Reference Budd and JenssenBudd and Jensen, 1975; Reference KrussKruss, 1983; Reference OerlemansOerlemans, 1997a). It appeared to be difficult to get a match between the results of a numerical SIA model and Nye’s theory (Reference Van de Wal and OerlemansVan de Wal and Oerlemans, 1995).

A practical method to derive a response time for glacier length was developed by Reference Jóhannesson, Raymond and WaddingtonJóhannesson and others (1989). Their method is based on an estimate of the time it takes to go from one equilibrium glacier geometry to another, given the typical mass-balance conditions. The corresponding response is therefore called a volume timescale. The response time is obtained by dividing a characteristic glacier thickness by the ablation rate at the glacier terminus. This method seems to work for many glaciers, but tends to underestimate the response time when the slope of the glacier is small and the altitude–mass-balance feedback becomes important (Reference OerlemansOerlemans, 1997a, Reference Oerlemans, Lisse and Balkema2001; Reference Harrison, Raymond, Echelmeyer and KrimmelHarrison and others, 2003). However, there is no consistent definition of characteristic ice thickness. The agreement between the Jo´hannesson timescale and other estimates depends on this ice thickness.

An alternative method, more in line with common practice in fluid mechanics, is to estimate a response time by dividing glacier length by a characteristic ice velocity. The characteristic velocity can be taken as the balance velocity at the equilibrium line (Reference Oerlemans, Lisse and BalkemaOerlemans, 2001). Glaciers with a large mass turnover will thus be faster. This is still in line with the Jóhannesson timescale, because such glaciers normally have large balance gradients and therefore high ablation rates at the snout.

Finally, calculations with a higher-order model, in which the full stress field is considered, reveal that models based on the SIA do a good job, but in some cases tend to overestimate the response time (Reference Leysinger Vieli and GudmundssonLeysinger Vieli and Gudmundsson, 2004).

In these and other studies, results are usually compared with results from other models, not with reality. This is because ‘clean’ data describing the response of a glacier to a precisely known change in climate forcing over a period of 100 years or more do not exist. However, in some cases it may be possible to estimate a response time from observed changes in glacier length even without knowing the forcing function a priori. In this work we consider the length records of two pairs of glaciers which have very different geometries but at the same time must be subject to virtually identical climatic forcing.

The first pair of glaciers is formed by Vadret da Palü and Vadret da Morteratsch (Fig. 1a and b; Table 1), Switzerland, which descend from the same mountain. Vadret da Palü flows southward and then eastward from Piz Palü (3900 ma.s.l.), and during the past 100 years the lower part of the glacier has been in steep terrain. Vadret da Morteratsch flows northward from Piz Palü and Piz Bernina (4049m a.s.l.) into a gentle valley. It can be expected that the response time of Vadret da Palü is considerably smaller than the response time of Vadret da Morteratsch.

Fig. 1. Photographs of the glaciers studied. (a) Vadret da Morteratsch looking south; photograph taken in June 2003 by C. Rothenbühler. The maximum stand of the glacier, around 1860, is clearly visible in the landscape. (b) An old photograph of Vadret da Palü (around 1905), from the collection of the Ökologische Gesellschaft München. In 2007 the glacier front is almost on the plateau. Note the large side moraine on the righthand side. (c) Briksdalsbreen, the steep glacier on the left descending from the plateau of Jostedalsbreen; photograph taken by P. Jantunen in June 2001. (d) Nigardsbreen descending from the plateau of Jostedalsbreen in 1948 (oblique aerial photograph by Widerøe).

The second pair of glaciers consists of Briksdalsbreen and Nigardsbreen (Fig. 1c and d; Table 1), located in southern Norway. Both glaciers originate from the Jostedalsbreen plateau, (about 2000 m a.s.l.). Briksdalsbreen flows in a northwesterly direction and is quite steep. Nigardsbreen flows in a southeasterly direction and flows into a valley with a gentle slope. The accumulation regions of these glaciers are flat and only about 10km apart.

The length records of the four glaciers used in the present study are shown in Figure 2. Dots show the annual data. There is a gap in the record for Vadret da Palü, partly due to difficulties in defining the position of the glacier snout because of the rugged terrain (Fig. 1b). Although the shape of Vadret da Palü has changed significantly in the last few years, it had a proper tongue for the entire record that is used here.

Fig. 2. Glacier length record for (a) Vadret da Palü and Vadret da Morteratsch and (b) Briksdalbreen and Nigardsbreen. Smooth curves were obtained by Stineman-filtering (Reference StinemanStineman, 1980). Data from the records of the Swiss Glacier Monitoring Network (http://glaciology.ethz.ch/messnetz/index.html), and from Annual Reports of NVE (NorwegianWater Resources and Energy Directorate, Oslo).

During the period considered, all four glaciers have retreated, although the net retreat of Briksdalsbreen is within 500 m. Since 1900, Nigardsbreen and Vadret da Morteratsch have retreated over a distance of about 2000 m.

On a decadal timescale, there are significant differences within the two pairs of length records. A first inspection of the curves indeed suggests that these differences can be explained, at least partly, by a difference in response time and climate sensitivity. The most obvious hint is the fact that Vadret da Palü starts to advance when the retreat of the Vadret da Morteratsch slows (around 1920 and 1975). Similarities can also be seen between the records of Briksdalsbreen and Nigardsbreen if one assumes that Briksdalsbreen responds faster to climate forcing but has a smaller climate sensitivity. It should be possible to take advantage of the differences in length records by reconstructing the climate histories from the length records and then minimizing the difference between these histories by varying the response times. If the reconstructed climate histories are found to be similar for reasonable response times, the method can be considered successful.

We will reconstruct the climate history in terms of the equilibrium-line altitude (ELA), denoted by E _M(t), E _P(t), E _B(t) and E _N(t), for Vadret da Morteratsch, Vadret da Palü, Briksdalsbreen and Nigardsbreen, respectively. Response times are then found by minimizing the difference between E _M(t) and E _P(t), and between E _B(t) and E _N(t).

2. Definition of Response Time

The concept of an e-folding response time formally relates to a linearized system (because then the approach to a new equilibrium is exponential in time). However, glaciers are not linear systems, implying, for instance, that the response time depends on the magnitude of the climatic change imposed on the glacier. This is undesirable, because the response time should be a physical property of a glacier, and be independent of forcing or glacier history. However, when changes in the forcing are sufficiently small, the concept of an e-folding response time is useful.

A first-order linear response equation can be written as

Here L’ and E’ are glacier length and ELA relative to reference values, c (<0) is the climate sensitivity and τ is the response time. The solution (for constant E’ and initial state L′(t = 0) = 0) reads:

This implies that at t = τ about two-thirds of the response has been accomplished.

In the literature, one sometimes encounters the term ‘reaction time ’, loosely defined as the time it takes before a glacier snout reacts to a relatively sudden change in the climate. It should be noted that this is not a well-defined concept, because the reaction time may depend on the glacier/climate history in a non-transparent way.

3. Method

Equation (1) can be applied to the four glaciers by assuming different values for the climate sensitivity (denoted c _i, where the index, i, refers to the specific glacier, i.e. i = M, P, B or N) and the response time (τ _i). To extract the history of the ELA from the glacier length record, ‘backward modelling’ is done, i.e.

In this approach the ELA at any time is determined by the glacier length and the change of glacier length with time (Reference Oerlemans, Lisse and BalkemaOerlemans, 2001, Reference Oerlemans2005).

As discussed in section 1, it is assumed that the pairs of glaciers are subject to the same climatic forcing. This implies that for ideal models with the appropriate response times and climate sensitivities, the outputs (t) and (t), and also (t) and (t), should be identical. However, in practice the outputs will always differ. A possible measure of the differences is

The number of data points in the records is denoted by K or J. For given values of the climate sensitivities, the quantities τ _M, τ _P, c _M, (t) and (t) can be found by minimizing the cost function, ψ_M,P. Extensive experimentation has shown that up to three control parameters can be used to minimize the cost function. It is not possible to derive values for the two response times and the two climate sensitivities. In that case the system runs away in a mode with continuously decreasing sensitivities and increasing changes in the ELA. In the approach adopted here, therefore, τ _M, τ _P, c _M have been defined as the control parameters (and similarly for Briksdalsbreen and Nigardsbreen).

Because a time derivative appears in the righthand side of Equation (3), irregularities in the glacier length records lead to unrealistically large values of E’. This can be avoided by proper preparation of the length record. As a first step a smoothing and interpolation procedure is needed. Generally speaking, glacier length records are characterized by very irregular spacing of data points. Most smoothing methods have problems with this. In principle, one would like to have control of the frequency response of the filter. This is achieved by a Fourier method or Gaussian filtering. However, such methods require further tricks to avoid spurious effects at the beginning and end of a record. Polynomial-fitting methods tend to generate glacier stands outside a realistic range. Here another property of glacier length records is relevant: although older records may have a very few data points indicating maximum glacier stands, the scarcity of data points should not lead to an even larger glacier length in between. Most interpolation methods cannot deal with the fact that interpolated glacier length should not exceed maximum stands. Altogether, after much experimentation, it appears that the method of Reference StinemanStineman (1980) is the most satisfactory (see also http://cran. r-project.org/doc/packages/stinepack.pdf). (The Stineman interpolation method is, for instance, implemented in the Macintosh application KaleidaGraph.) The result is shown in Figure 2 (solid lines).

The optimization procedure is performed using a simple random-walk method. Values of the response times are perturbed at random, and the new values are retained when the corresponding value of ψ is lower than its value in the previous step. In practice, convergence to a stable solution appears to be very rapid. The procedure is repeated for different initial values of the control parameters to ensure a unique solution exists.

The final step is to assign values to the climate sensitivities that are not taken as control parameters, namely c _P and c _B. From a simple analytical glacier model it follows that the first-order sensitivity is given by Oerlemans (2001, ch. 5):

where s is the mean slope along the flowline of the glacier.

This expression is valid for a glacier of constant width resting on a bed with a constant slope. Clearly, in reality, these conditions are not very well met. If the characteristic width of the accumulation zone is larger than the width of the ablation zone, which is the case for all the glaciers considered in this study, c will be larger. Some basic calculations (Reference Oerlemans, Lisse and BalkemaOerlemans, 2001, p. 77) suggest that an increase in c due to converging flow is typically 50% for a glacier with an accumulation basin that is typically four times as wide as the ablation zone.

Estimating the characteristic slope of Vadret da Palü during the past 100 years yields a value of ~1/3. From Equation (6) it follows that c _P ≈ –6 for a glacier of constant width. Although Vadret da Palü does not have a very regular geometry, one can estimate from a topographic map of the glacier (Swisstopo, sheets 1277 and 1278, scale 1 : 25 000) that the width of the accumulation area is typically four times the width of the tongue as it appeared during the past 100 years. The corresponding value of the climate sensitivity is, then, c _P ≈ –9. This is considered to be the best estimate and is used in the following analysis. It implies that a 100 m change in the ELA corresponds to a change in the equilibrium glacier length of 900 m.

The geometric characteristics of Briksdalsbreen (Statens Kartverk, Norway, sheet M711-1318 II, scale 1:50 000) have similarities with those of Vadret da Palü. The glacier is fairly steep and the accumulation basin is much wider than the tongue, although not as extreme as for Vadret da Palü. The resulting best estimate of the climate sensitivity of Briksdalsbreen is c _B ≈ –8.

4. Results

It seems that the random-walk model quickly produces a unique and meaningful solution for the climate history and associated response times. Results are summarized in Table 1 and Figure 3. Figure 3a shows the result for the Swiss glaciers ((t) and (t)). The first thing to note is the good correspondence of the reconstructed ELA for the two glaciers. In fact, the only significant difference is that after 1910 (t) shows a small maximum, whereas the (t) curve is rather flat. The minimum value of ψ is only 7.34 m. The corresponding optimal response times are τ _M = 33.0years for Vadret da Morteratsch and τ _P = 4.4 years for Vadret da Palu. So apparently Vadret da Palü reacts very quickly to changes in the ELA, in view of the steepness of the glacier this is quite possible. Furthermore, c _M = –17.7, implying that the equilibrium length of Vadret da Morteratsch decreases by 1770 m for a 100 m rise of the equilibrium line. This seems to be a reasonable sensitivity for a glacier in a relatively flat valley.

Fig. 3. ELA reconstructed from the glacier length records for (a) Vadret da Palü and Vadret da Morteratsch and (b) Briksdalbreen and Nigardsbreen.

The model predicts an overall increase in the ELA over the past 100 years of about 120 m, which is in broad agreement with earlier studies (e.g. Reference MaischMaisch, 2000). The most obvious feature on a decadal timescale is the distinct maximum of the ELA around 1960. However, values obtained for the last decade are definitely the highest over the past 100 years.

The result for the pair of Norwegian glaciers is shown in Figure 3b. The findings are similar. The reconstructed histories of the ELA match well for realistic response times (τ _B = 5.02 years and τ _N = 34.82 years), and the minimum value of the cost function, ψ, (6.5 m) is slightly smaller than that for the pair of Swiss glaciers (7.34 m). The response time for Briksdalsbreen is small, but one should realize that this glacier is steep and has a very large mass turnover (ablation rates on the snout are up to ~10m of ice). The climate sensitivity found for Nigardsbreen is –26.1, significantly larger than for Vadret da Morteratsch (–17.7). It is likely that the high sensitivity of Nigardsbreen is a consequence of the very small slope of the valley and of the fact that the accumulation basin is much wider than the glacier tongue.

The result for the Norwegian glaciers points to a rise of the equilibrium line during the past 100 years of about 60 m. However, the most prominent feature is the steep rise in the equilibrium line from 1930 to 1955. During recent years the equilibrium line has gone upwards again, but has not yet reached a level comparable to that of the mid-20th century.

Of the four glaciers studied, only Nigardsbreen has been modelled more explicitly with a numerical flowline model (Reference OerlemansOerlemans, 1997b, Reference Oerlemans, Lisse and Balkema2001). This model, based on the SIA, yields a length response time of 68 years. The Jóhannesson timescale (referred to in section 1) for Nigardsbreen is only 15–20 years, which is a consequence of the high melt rates at the glacier tongue. The value found in the present study (34.8 years) is in between. It is not clear why the numerical model produces a value that is twice as large as that found empirically. There are indications that models based on the SIA initially react too slowly to changes in the mass balance (Reference Leysinger Vieli and GudmundssonLeysinger Vieli and Gudmundsson, 2004).

At this point the sensitivity of the results should be considered. To investigate how well the response times are defined by means of the current method, we can look at the behaviour of the cost function for varying response times. The result of such a calculation is shown in Figure 4. When τ _P is varied, τ _M is fixed, and vice versa. Judging from the shape of the cost function, the values of the response times are well defined. A similar exercise for Briksdalsbreen and Nigardsbreen yields a similar result (not shown).

Fig. 4. The cost function, ψ, (defined in Equation (6)) for varying response times for the pair of Swiss glaciers, showing that the optimal values are well defined.

The optimal values found for the response times also depend on the climate sensitivities, ci, of course. Figure 5 shows the optimal values of τ _n, τ _B and c_N for varying c_B. The sensitivity is modest. It should be noted that, for decreasing values of c_B, the cost function, ψ, increases rapidly. The picture for the Swiss glaciers is very similar (not shown).

Fig. 5. Sensitivity analysis for the pair of Norwegian glaciers. The optimal values for (a) τ _B, τ _N and (b) c _N vary smoothly with the prescribed value of c _B.

5. Discussion

In this paper, a method has been proposed to derive glacier response times and to reconstruct ELA histories from glacier length records. The prerequisite is that at least two glaciers can be found that have been subject to almost identical climate forcing. For the glaciers studied, the method works well. The proximity of Vadret da Morteratsch and Vadret da Palü on the one hand, and of Briksdalsbreen and Nigardsbreen on the other, justifies the assumption of similar forcing.

The approach taken here is based on the simplest possible model that can describe a dynamic response to time-dependent forcing, namely the first-order linear response model. This means that, for a pair of glaciers, there are four model parameters. In view of the very good match between ELA reconstructions, i.e. small values of the associated cost function, it makes no sense to go to a more refined model with additional control parameters.

The response times found are perhaps rather smaller than expected on the basis of numerical flowline modelling. Especially for the steep glaciers, Vadret da Palü and Briksdalsbreen, the values are quite small. However, there is no independent evidence that these values are unrealistic. The response times are the outcome of a straightforward procedure and best explain the relation between forcing and response.

The good match between the ELA reconstructions suggests that the choice of pairs of glaciers has been somewhat fortuitous. It would be very interesting to apply the method proposed here to other pairs of glaciers. Even when the glaciers are not as close to each other as those considered in the present study, it should be possible to extract the common part of the climatic forcing. As long as reconstructed ELA histories are well correlated, the analysis makes sense.