Introduction

Because there is a strong relation between climate and glaciers, modelling glaciers using historical records of glacier advance and retreat should elucidate past climate change. Many such model studies to simulate glacier behaviour have been carried out, but few field observations on glaciers allow a thorough check of the results of these models. Numerical ice-flow models based on the continuity equation have been developed to simulate geometry changes in response to mass-balance scenarios on Nigardsbreen, Norway (Reference OerlemansOerlemans, 1986), Storglaciären, Sweden (Reference StroevenStroeven, 1996), Unterer Grindelwaldgletscher, Switzerland (Reference Schmeits and OerlemansSchmeits and Oerlemans, 1997), and Hintereisferner, Austria (Reference SchlosserSchlosser, 1997). Unfortunately, the results of these calculations cannot be compared with observations, because (1) observations do not exist or are scarce, or (2) the available observations have been used to calibrate the parameters of the model. In this paper, owing to the wealth of observations, we are able to test the modelling results with the observations. In this way, we underscore the limits of this kind of numerical modelling, particularly for the calibration of flow parameters.

It is rare that mass-balance data, elevation changes, velocities and terminus positions are available for several decades on the same glacier. Glacier de Saint Sorlin is a small glacier (3 km²) located in the Grandes Rousses area, France (45°10′ N, 6°10′ E) (Fig. 1). The mass balance of this glacier has been the subject of several studies (Reference LliboutryLliboutry, 1974; Reference Lliboutry and EchevinLliboutry and Echevin, 1975; Reference Vallon and LeivaVallon and Leiva, 1982; Reference VincentVincent, 1995), and glaciological mass balances and velocities have been measured directly since 1957. In addition, glacier thicknesses are known from boreholes, seismic soundings and gravimetric measurements (Reference BelinBelin, 1962; unpublished information from M. Echevin and from M. Vallon), and altitude changes were determined six times between 1952 and 1997 using photogrammetric or geodetic means.

Fig. 1. Glacier de Saint Sorlin in 1971, showing cross-sections measured by photogrammetry (solid lines) and flowline used in the model (dashed line). The large dots are positions where thicknesses and velocity were calculated in the model (600, 1400, 1800 and 2200 m from col des Quirlies).

In this paper, a numerical ice-flow model is adapted to reconstruct the fluctuations of glacier de Saint Sorlin and to study the large changes in its dynamics. The model is based on the continuity equation integrated over depth and width. Owing to the wealth of observations, it is not necessary to reconstruct mass balance from meteorological data between 1957 and 1997 or to adjust the bedrock profile. Therefore, the measurements of surface velocities and height changes can be fully used to check the results of the model and to gain insight into sliding processes. In the following section, mass-balance observations are thoroughly analyzed in order to obtain a reliable mass-balance series over the entire surface of the glacier. These data are then used to simulate the dynamic behaviour.

Mass-Balance Data

Mass-balance data have been collected using stake measurements scattered over the glacier. As the measurements were not made at the same places each year, we use a statistical method of Reference LliboutryLliboutry (1974), herein called the linear model, to interpolate and extrapolate the data. The linear model assumes that the mass balance can be decomposed into a spatially independent temporal variation, β_t , and a spatial pattern independent of time, α_j , thus:

(1)

where b _j,t Sis the specific balance at site j during year t, α_j is the average balance at that site over the period of record, β is a constant for any given year and ∊ _j,t is an error term. Instead of treating each stake as a site of observation, the surface of the glacier is divided into 0.2 km × 0.2 km squares, within each of which the mass balance is assumed to be constant for a given year (if there is more than one observation within one of these squares, the mass-balance value used is the average of the available measurements). A matrix is then constructed in which each cell represents one of these squares and a particular year. As data are available for only 24% of the cells on the glacier, the Lliboutry algorithm is used to estimate missing values. The standard deviation between reconstructed and observed annual mass balances is 0.28 m w.e. This method can be used:

Fig. 2. Mean specific net mass balance of glacier de Saint Sorlin (m w.e.) over the last 40 years (1957–97).

The cumulative specific net balance calculated from the linear model has been negative (Fig. 3). Since 1957 the glacier has lost an average of 12.4 m w.e., or 0.31 m a⁻¹.

Fig. 3. Cumulative specific net balance (m w.e.) on glacier de Sarennes and glacier de Saint Sorlin. The dashed line shows the cumulative specific net balance calculated from stakes located in a 0.25 km² area in the lowest part of the glacier. Triangles show changes, relative to 1997, based on geodetic measurements (see text).

In the linear model, we assume that the change in mass balance between any two successive years is the same everywhere on the glacier (i.e. β_t is constant, independent of position). In other words, the balance gradient is independent of time. To test this assumption, some annual mass-balance data are plotted against distance from the terminus to col des Quirlies (Fig. 4). It can be seen that the gradient remains unchanged (within uncertainty of the data) on this relatively small glacier.

Fig. 4. Annual mass balance observed on the flowline between col des Quirlies and the terminus for certain years, and the average over 40 years.

This assumption can be questioned, however, when applied to the full altitude range of larger glaciers (Reference KuhnKuhn, 1984; Reference OerlemansOerlemans, 1993; Reference Vallon, Vincent and ReynaudVallon and others, 1998). Furthermore, measurements in the accumulation zone of glacier de Saint Sorlin have been made for only 11 years. Thus, in order to check these results, mass-balance values have been calculated using a direct glaciological method based on altitude bands. For this purpose, a map with contours of elevation and of net balance is used to determine the average net balance for an average elevation interval of 100 m (Reference PatersonPaterson, 1994, p. 32). Observations were made over the entire glacier surface in only 5 years (about 30 observations over 3 km²), so only these years are available for such a comparison. For these 5 years, it can be seen that the results are in reasonable agreement with the linear model (Table 1). Differences do not exceed 0.5 × 10⁶ m³ w.e., or 0.17 m w.e. averaged over the glacier. Such differences are within the range of uncertainty of the observations.

Table 1. Mass balance (in 10⁶ m³ w.e.) from altitude-band method and from Lliboutry linear model

The raw data used for these two calculations are, however, the same, so the results are not independent. In order to obtain an independent estimate of volume loss, topographic measurements and aerial photographs were used to calculate volume changes. We call this the geodetic method. The dates and origins of the available documents are shown in Table 2.

Table 2. Date and origin of geodetic measurements (photogrammetry and topography)

Instead of mapping contour lines from these photographs, only transverse profiles were constructed (Fig. 1). Thickness variations were then extrapolated over the entire surface. This seemed to be more accurate for the following reasons:

Nevertheless, in this method we do assume that temporal altitude variations on the chosen cross-sections are applicable to a large area up- and down-glacier from the cross-section. The resulting thickness variations are shown in Figure 5; it can be seen that the glacier thickness increased from 1952 to 1984 at the highest elevations, and decreased at lower elevations. From 1984 to 1997, the thickness decreased everywhere.

Fig. 5. Thickness change relative to 1952.

A direct comparison of these results with the mass balance determined from the linear model is not possible because no map is available from 1957, the year of the first mass-balance measurements. To make this comparison, we had to estimate the volume change between 1952 and 1957, using data from glacier de Sarennes (Reference Valla and PiedalluValla and Piedallu, 1997), which is located 3 km away. There is a good correlation between the mass-balance changes of these two glaciers (Reference VincentVincent, 1995). Figure 3 shows the cumulative specific net balance of glacier de Saint Sorlin, together with the changes inferred from volume variations, relative to 1997, obtained by geodetic methods. It is difficult to know whether the slight discrepancies between the geodetic method and the linear model result from: (1) the quality of the first aerial photographs in 1952, for which camera distortions are large; (2) uncertainty in mass-balance measurements (the standard deviation is about 0.20 m w.e. for 1 year, or 1.26 m w.e. for 40 years if the errors are random); (3) the assumptions introduced in the linear model; or (4) dates of the photographs which do not always match those of the mass-balance measurements.

In order to test the linear model further, it has been applied to a small part of the glacier (0.25 km²) close to the terminus; in this way, one can calculate new values for the interannual mass-balance variations (β_t ). Thus, with this small net we evaluate only the temporal term in Lliboutry’s model. The method then relies on the mean mass balance, α_j (obtained from Lliboutry’s model applied over the entire surface), to scale the result to the whole glacier. It can be seen (dashed line, Fig. 3) that the results are similar: the two volume-variation calculations differ by <3 × 10 m w.e. (1 m w.e. in cumulative specific net balance). This is only slightly more than the expected uncertainty in the volume variation over 20 years (0.90 m w.e.).

One might conclude that some stakes in the small terminus area of a glacier, together with maps at the beginning and the end of the period, would suffice to determine the glacier mass balance over a 40 year span with an accuracy which remains within the range of uncertainty This conclusion, however, is valid only for small glaciers like glacier de Saint Sorlin.

Dynamic Behaviour of Glacier de Saint Sorlin

Since 1957, annual velocities of stakes used for mass-balance observations in the ablation zone of glacier de Saint Sorlin have been measured. In two years, 1972 and 1973, measurements were also made in the accumulation zone. Thickness changes are known from geodetic measurements in 1952, 1971,1984,1986,1989 and 1997 (Fig. 5). In addition, a map from 1905 (scale 1: 10 000) is available, although its accuracy is uncertain. Finally, ice thicknesses were measured by geophysical means (see Introduction). Here, we analyze the thickness and velocity variations along a central flowline between col des Quirlies and the terminus. It is immediately apparent (Fig. 6) that the velocities increased significantly between 1960 and 1983 in the central region of the glacier, although the surface geometry remained almost unchanged (the change in slope is about 2% in this region). What triggers such phenomena?

Fig. 6. Measured surface velocities on the flowline between col des Quirlies and the terminus. Only selected years are shown.

Ice-flow model

An ice-flow model is adapted to reconstruct fluctuctions of glacier de Saint Sorlin. This model is based on principles developed by Reference OerlemansOerlemans (1986) and Reference GreuellGreuell (1992). Ice flow is calculated along the central flowline (x axis) (Fig. 1). The change in cross-sectional area with time can be described by the continuity equation:

(2)

where S is the cross-sectional area, U is the mean velocity over the cross-section, B is the mass balance, W is the width of the glacier surface and x is the downstream distance along the flowline. The total ice velocity, U, is the sum of the deformation velocity, u _d, and the sliding velocity, u _s. The mean deformation velocity is calculated from:

(3)

where H is the ice thickness, f is a flow parameter, and τ _b is the driving stress given by

(4)

where ρ is the ice density, g is the acceleration due to gravity and h is the surface elevation. f is taken to be 1.3 × 10⁻²⁴ Pa⁻³ s⁻¹ (Reference LliboutryLliboutry, 1964).

The form of the sliding law is controversial. A number of modellers have assumed that the sliding velocity is proportional to the third power of the driving stress divided by the effective pressure, the ice overburden pressure minus the water pressure at the bed. In addition, as the water pressure has not been measured over long periods, it is assumed this water pressure is a fraction of the overburden ice pressure (Oerlemans, 1998).

Because the physical processes of glacier sliding are no well understood quantitatively, sliding is not taken into consideration in initial calculations with the present mode Nevertheless, in a second step, we introduce sliding velocities from an empirical relation between sliding, shear stress and overburden ice pressure, in order to test the limit of this modelling.

Input data and calibration procedure

Bedrock topography and mass balance are input data and are not used to calibrate the model. The mass-balance value between 1957 and 1997 have been extracted from the analysis in the previous section. Between 1905 and 1957, the glacier d Sarennes mass balance, reconstructed from meteorologice data (Reference MartinMartin, 1978), is used. Although such reconstruction can be questioned (Reference Vincent and VallonVincent and Vallon, 1997), they are probably the best technique available for estimating the mass-balance series prior to the period of observation.

The three-dimensional geometry is taken into account by allowing flowlines to converge linearly in the accumulation area: a convergence parameter, p, allows W (Equation (2) to change along the central flowline, thus:

(5)

where W ₀ is the initial width at col des Quirlies, p is the convergence parameter, x is the distance from col des Quirlies, and L is the distance over which the flowband is reduced (1 km). p is the only adjustable parameter in our approach. It remains constant with time.

In the ablation zone, the flowlines are inferred from velocity data and are known (and thus modelled) to be more-or-less parallel. In the accumulation zone, velocity measurements have been carried out over only 2 years, and p has been roughly adjusted so that the simulated longitudinal profile (from steady-state conditions) fits the observed 1971 profile and simulated velocities on the flowline fit observed velocities. It was found that in order to achieve this, the distance between flowlines had to be decreased 50% over the first kilometer. With regard to calibration of p, we note that:

1. (1) This pattern of W(x) is consistent with the assumption that ice flows normal to contour lines. In addition, further experiments show that temporal variations of calculated velocities are not strongly affected by the value of p.

2. (2) Under steady-state conditions, using this value of p and a (steady-state) mass-balance pattern obtained by adding 0.31 m w.e. to the pattern in Figure 2, the maximum discrepancy in the surface elevation reaches 25 m in the center of the glacier where the ice is about 120 m thick. This difference may be due to the approximate knowledge of the bedrock topography, the non-steady-state conditions of the glacier in 1971 or the principles of the model (in which ice flow is calculated along a central flowline).

For model runs, the initial surface profile is taken to be the bedrock profile. The first aim is to reconstruct the 1905 glacier (the first available map). For this purpose, the present-day annual mass balances shown in Figure 2 are increased by 0. 55 m w.e., and the mass balance is assumed to be constant with time. After about 500 years of simulated time, an equilibrium, close to the 1905 state, is reached. We thus, implicitly, assume that the 1905 glacier was close to a steady state. This assumption is most probably wrong, but without consequence because the period analyzed (1957–97) is relatively far from this date. According to Reference Jóhannesson, Raymond and WaddingtonJóhannesson and others (1989), the response time-scale for adjustment of a glacier to changes in mass balance is given by the maximum thickness divided by the balance rate at the terminus. For glacier de Saint Sorlin, the response time-scale would thus be about 60 years. Therefore, it can be assumed that changes in mass balance before 1905 do not significantly influence either glacier geometry or velocities between 1957 and 1997.

Results and comparison with the observations

Shown in the lefthand column of Figure 7 are modelled changes in surface elevation at four points along the flowline from col des Quirlies, based on the model output. Also shown are observed elevations at those points at six to eight different times. In the righthand column of the figure, velocities calculated at various locations along the flowline are compared with observed velocities at the same locations. The modelled thickness variations are very close to the measured ones except for surface elevations at the highest and lowest locations in 1905. For 1905, it is not surprising to see some discrepancies: either mass-balance data or map accuracy could be responsible for these, as could the assumption of a steady state in 1905.

Fig. 7. Variations in thickness and velocity through time (solid lines), calculated from the model without sliding, for locations at 600 m (a), 1400 m (b), 1800 m (c) and 2200 m (d) down-glacier from col des Quirlies along the flowline. The dots are observations. The dashed lines without dots are results from the model after introducing sliding

The modelled velocity variations between 1957 and 1997 (solid lines on right side of Fig. 7) are not, however, in accordance with observations. Figure 7 shows that:

1. (1) One can reconstruct thickness variations which fit well with observations, although the dynamic processes are misunderstood.

2. (2) Between 1960 and 1983, observed velocities increased by 80%, although the surface elevation remained almost unchanged (plus or minus a few meters) over most of the glacier, and decreased 15 m in the vicinity of the snout (Fig. 7d).

3. (3) Present velocities are higher than the velocities in 1960, although the thickness decreased everywhere.

4. (4) The time series of the differences between observed surface velocities and calculated deformation velocities (Equation (3)) are very similar on the lower part of the glacier (Fig. 8). Higher on the glacier (Fig. 7a), the observations are not as complete so no conclusions can be drawn.

Fig. 8. Difference between measured surface velocities and velocities calculated from internal deformation alone for locations 1400, 1800 and 2200 m from col des Quirlies. Summer ablation on glacier de Sarennes is shown at top.

Additional modelling tests showed that it is not possible to improve the results by changing f in Equation (3) or by adding sliding velocities which depend upon empirical relations between sliding, shear stress and overburden ice pressure (Reference PatersonPaterson, 1994, p. 157). These changes can increase the velocities but do not improve the agreement with observed changes in velocity.

The velocity differences shown in Figure 8 are probably a good estimate of the sliding velocities over the lower half of the glacier. Accordingly, these differences were introduced into the model as sliding velocities in order to study the sensitivity of the thickness to the change in velocity. Of course, in this experiment the resulting surface velocities (dashed lines without dots in Fig. 7) now fit observed velocities. Because u _d decreases rapidly as H decreases, thus compensating for the increase in u _s (Equations (3) and (4)), adding this sliding component to the model results in only small changes in surface elevation; the maximum difference is 7 m.

The inferred sliding velocities were compared with summer ablation on glacier de Sarennes (Fig. 8, top). No relation is evident on time-scales of 1 or 2 years. Nevertheless there is a trend toward increasing ablation from 1970 to 1990, accompanied by a clear increase in sliding. Both then decrease after 1990.

Finally, the measured terminus position (Fig. 9) shows a pattern of change similar to that of the observed thickness changes on the lower part of the glacier: the average retreat rate was about 10 m a⁻¹ between 1952 and 1970, and 3.5 m a⁻¹ between 1970 and 1997.

Fig. 9. Measured terminus positions since 1905.

Discussion and Conclusions

Mass-balance and velocity measurements carried out on glacier de Saint Sorlin over 40 years provide insights into dynamic processes.

First, the mass-balance analysis suggests that reliable annual mass-balance measurements from only a few stakes on a reduced part of the glacier can be used to estimate interannual mass-balance changes over a few decades on a small glacier like glacier de Saint Sorlin. Moreover, the time series of cumulative specific net balance can then be determined using maps made at the beginning and end of the study period. For larger glaciers, this method is not appropriate because the mass-balance gradient may vary with elevation; thus observations must be extended to higher elevations (Reference Haeberli, Haeberli, Hoelzle and SuterHaeberli, 1998; Reference Vallon, Vincent and ReynaudVallon and others, 1998).

Dynamic changes on glacier de Saint Sorlin since 1957 are significant. A simple numerical ice-flow model, based on the continuity equation, simulates dynamic changes during this period. The reconstructed thickness variations agree with observations, although velocities do not.

The present study shows the limits of this kind of modelling, particularly as concerns calibration of flow parameters from geometric changes. In this kind of model, the total ice velocity is expressed as the sum of the deformation velocity and the sliding velocity. The deformation velocity is expressed as a function of the thickness and the driving stress (Equation (3)) with a flow parameter. With regard to the form of the sliding law, a number of modellers use an empirical relation between sliding velocity, basal shear stress and effective pressure (Reference PatersonPaterson, 1994, p. 157) with a second flow parameter. Modelling tests of the present study show that such empirical relations for sliding do not improve the agreement between modelled and observed changes in velocity. Moreover, we found that, over 40 years, the surface geometry change is not very sensitive to changes in velocity (experiments with doubling the first flow parameter or with introducing the large inferred sliding velocity). As a result, it is not useful, and indeed it is misleading, to calibrate flow parameters from terminus fluctuations in this kind of model, for two reasons:

1. (1) The calibration of the flow parameters from the historical glacier-length record, called the dynamic calibration procedure (Reference OerlemansOerlemans, 1997; Reference Wallinga and van de WalWallinga and others, 1998), requires >50 years of observations (glacier-length and mass-balance records), as has been shown in the present study. Thus, because direct mass-balance observations are not available prior to 1945 anywhere in the world, it seems impossible to fix these parameters correctly.

2. (2) The present study shows clearly that sliding conditions can change significantly over 20 years and it is misleading to hold the sliding-law parameter constant in time in an empirical relation between sliding, shear stress and overburden ice pressure.

Computed velocities, determined from the calculated driving stress alone, differ significantly from observed velocities. Similar results were obtained from 22 years of velocity measurements on Nisqually Glacier, Washington, U.S.A. (Reference MeierMeier, 1968). On glacier de Saint Sorlin, comparison between the reconstructed deformation velocities, calculated from the driving stress, and the observed surface velocities provides an estimate of the sliding speed. It is noteworthy that the inferred temporal variations in sliding speed are similar over the bottom half of the glacier for the last 40 years (Fig. 8). Moreover, as in Meier’s study, the inferred sliding velocity cannot be described by classical analysis (Reference WeertmanWeertmann, 1957) or empirical relations in which the sliding velocity depends on the driving stress (Reference PatersonPaterson, 1994, p. 157). Nor can the increase in velocity between 1960 and 1985 be explained by secular change in the local driving stress.

It seems, in addition, that the increase in velocity between 1960 and 1985 is not related to the quantity of water coming from surface ablation. Basal water-pressure fluctuations might explain the observed acceleration, but we do not have such measurements over such a long period of time. Thus, it is not possible to draw conclusions about the influence of the basal water pressure. The thorough studies carried out over some months on Findelengletscher, Switzerland (Reference Iken and BindschadlerIken and Bindschadler, 1986), or Storglaciären (Reference JanssonJansson, 1995) relate the enhanced basal sliding in summer to the high basal water pressure. At the surging-glacier end of the spectrum, the detailed observations of the surge of Variegated Glacier, Alaska, U.S.A., in 1982/83 (Reference KambKamb and others, 1985) show that the flow acceleration is caused by the build-up of high water pressure in the basal passageway system. Glacier de Saint Sorlin shows a different pattern of flow acceleration, with a gradual rise in basal sliding over two decades. This acceleration is also probably linked to enhanced basal water pressure, triggered by changes in the geometry and water-transport characteristics of the basal conduit system.

These conclusions about the glacier de Saint Sorlin velocities are valid as long as longitudinal stress gradients can be neglected. H. G. Gudmundsson (personal communication, 1999) notes that a flow model of this type cannot be used to calculate surface velocities accurately because it does not incorporate such stress gradients. The volume time-scale is predicted correctly, but the diffusion and propagation time-scales are not (Reference Gudmundsson, Raymond and BindschadlerGudmundsson and others, 1998). However, the integral of the velocities over depth and across the width of the glacier can be correct. In fact, the flux must be more-or-less correct because the flux gradient is determined by the mass balance (personal communication from H. G. Gudmundsson, 1999).

With regard to a 100 year record of velocity observations on Hintereisferner, Reference Span, Kuhn and SchneiderSpan and others (1997) noted that an increase in ice thickness alone could not explain an enormous increase in surface velocity in about 1919 (from 30 m a⁻¹ to 125 m a⁻¹ over 5 years). In search of possible causes of this flow acceleration, Span and others mention a positive feedback process in which increased velocities result in increased crevassing which in turn increases water input to the glacier. In higher parts of glacier de Saint Sorlin, the mean longitudinal strain rate did not exceed 0.02 a⁻¹ (Fig. 6), which is less than the value of 0.03 a⁻¹ at which crevasses form (Reference Meier, Kamb, Allen and SharpMeier, 1974). In the lower part, crevasses are not widespread and changes are not visible in aerial photographs. Therefore, the increase in observed velocities on glacier de Saint Sorlin cannot be explained by this mechanism.

Flow processes of glaciers must be understood in order to simulate glacier response to climate change. For example, the response time depends on flow processes; with constant climate similar to present conditions, our model predicts that the glacier would reach an equilibrium state in 2140, and the average recession rate of the terminus would be about 4 m a for the next 50 years (in the model, the equilibrium state is considered to be reached when thickness variations are <0.1 m). If the velocities are doubled (the present study shows that the uncertainty in calculated velocity can reach 100%), the equilibrium state would be reached in about 2170 and the recession rate would be about 1.5 m a⁻¹ for the next 50 years.