Introduction: Ice-Sheet/Climate Relationships

We expect the long-term evolution of the landward margins of the Pleistocene mid-latitude ice sheets over Europe and North America to have been a reflection of climate change, whilst the marine margins may have been strongly influenced by internal dynamic oscillations, at least on short time-scales. Given the capacity of large, slowly varying ice sheets to buffer high-frequency climate oscillations, we might expect fluctuations of the landward margins of ice sheets to be a smoothed, long-term proxy record of atmospheric climate.

It is difficult to use the patterns of ice-sheet fluctuation as a direct proxy for climate that can be readily compared with other palaeoclimate proxy data, because of long ice-sheet response limes. However, numerical models of ice-sheet dynamics are capable of simulating the response of ice sheets to climatic change, permitting such models to be used to translate the geological evidence of ice-sheet fluctuation into more useful indices of palaeoclimate change. Conversely, it should be possible to use models to investigate whether proxy climate sequences derived from other palaeoclimate indicators are compatible with evidence of contemporary ice-sheet fluctuation. We therefore distinguish between a forward approach, in which a climate-change function is prescribed and the ice-sheet response is computed, and an inverse approach, in which ice-sheet behaviour is prescribed and the model is used to infer characteristics of the climate which could have driven the ice sheet.

The Ice-Sheet Model

We use the model described by (Reference Boulton and Payne.Boulton and Payne 1992, Reference Boulton and Payne.1993, Reference Boulton, Payne. and Duplessy1994) which predicts the form, internal velocity field and internal temperature held of an ice sheet in response to the following changes at the surfaces of the ice mass and properties of the bed:

The modelled ice sheet is allowed to evolve through time as a thermomechanically coupled system by solving a vertically integrated ice-mass continuity equation. The velocity and temperature fields are resolved into ten layers, whose resolution increases towards the bed. Parameters controlling the stress/strain response of ice, its dependency on temperature, the amount of strain heating, the response of the lithosphere and the nature of the marine boundary are held constant. Where the basal ice temperature is at the melting point, we have chosen to permit sliding in some runs and to prevent it in others. The sliding velocity is prescribed as a linear function of shear stress. As we are concerned primarily with the landward margin, a simple calving scheme has been used which does not permit the ice sheet to extend into water depths greater than 500 m.

A simple diffusion equation has been used to model the effect of asthenospheric flow on bedrock elevation, using an asthenosphere viscosity, which controls rates of vertical displacement of the bed, that is slightly more viscous than the values obtained In Reference Lambeck, Johnston and Nakada.Lambeck and others (1990). Comparison with sea-level recovery rates in Europe suggests that this compensates for the lack of an elastic lithosphere in the model.

We have modelled along a transect extending from the western continental shelf of Norway to northern Poland through the area occupied by the European ice sheet during; the last (Weichselian) glacial cycle (Fig. 1). This was an approximate flowline through much of the last glacial cycle, apart from during short episodes when ice streams in the southern Baltic crossed the transect from the east.

Fig. 1. The extent of the European ice sheet at the Last Glacial Minimum, its pattern of final decay and the location of the transect to which ice-sheet flowline models have been applied.

The climate drive for ice-sheet fluctuation is provided by independent variation of temperature and mass balance. A prescribed variation in mean annual sea-level air temperature (SLAT) is derived b\ assuming a linear correlation between this temperature and selected European and northeastern Atlantic palaeotemperature records for the glacial cycle (Reference Boulton and Payne.Boulton and Payne, 1992, Reference Boulton and Payne.1993). The vertical variation in air temperature is then derived from an assumed lapse rate of 10°C km⁻¹ (Reference OrvigOrvig, 1970), typical of rates over modern ice sheets. Sensitivity tests (Reference Boulton and Payne.Boulton and Payne, 1993) showed that ice-sheet fluctuations were relatively insensitive to the absolute values of SLAT within reasonable limits suggested by palaeotemperature reconstructions, but were sensitive to major departures of the lapse rate from modem values.

The vertical variation of mass balance above and below the equilibrium-line altitude (ELA) has a distinctive form in modern glaciers (Reference Boulton, Smith and MorlandBoulton and others, 1984). We prescribe the form of the mass-balance/altitude curves (Fig. 2) which, for different model runs, vary between continental (low accumulation and ablation) and maritime (high accumulation and ablation) conditions. The gradient of the ELA line along the transect is similar to the calculated modern gradient.

Fig. 2. Empirically based mass-balance patterns used in simulations, showing continental and maritime extremes.

The model is forced by deviations of the elevation of the ELA line from modern values. These alter the altitude at which local mass-balance curves intersect the ice surface at each point along the transect (Fig. 2). A linear correlation is assumed between ELA deviations and palaeotemperatures inferred from selected palaeoclimatic sequences through the last climatic cycle. Parallel variations in ELA and SLAT are therefore produced.

Palaeoglacial and Atmospheric Palaeoclimate Proxy Records

We have summarised ice-sheet fluctuation in Europe through the last glacial cycle by projecting Reference Mangerud and FrenzelMangerud’s (1991) geological reconstruction on to our line of transect (Fig. 5). We recognise that the extent of decay between the glacial maxima is uncertain, that the errors of chains may be significant and that there is an implicit assumption that the most extensive glacial phases are contemporary with the isotopic peaks in the deep-ocean isotopic record.

Fig. 5. Time/distance diagrams showing the extent of modelled ice sheets driven by the temperatures from core SU 90–39. The heavy line shows the approximate projection of Mangerud’s glacial history on to the transect. (a) No-sliding case. The modelled ice sheet reaches the observed LGM position when f₁ = f_1a. With larger values of f₁, the model ice sheet extends further, and becomes unstable when f₁ = f_1a × 1.17. For smaller values of f₁ = f_1a × 0.83, it fails or reach the LGM position, (b) Sliding cases. Here f₁ needs to be larger (f₁ = f_1a × 1.17) to reach the same position as non-sliding cases. because of the fatter profile produced by sliding. When f₁ = f_1a, the modelled ice sheet remains small. Sliding leads, however, to more rapid retreat because of the lower surface profile.

We have used two sources of information about atmospheric conditions in Europe through the last glacial cycle from the European area. Figure 3a shows palaeotemperatures throughout the last glacial cycle from La Grande Pile in northeast France, inferred from a long pollen record (Reference Guiot, Pons, de Beaulieu and Reille.Guiot and others, 1989; Reference GuiotGuiot, 1990). Figure 3b shows a new sea-surface temperature (SST) reconstruction from core SU90–39 at latitude 52°30’N, longitude 22°0’W in the Atlantic. A transfer function relating modern species distributions to sea surface temperatures estimates August SSTs and is constructed using the modern analogue reference method (Reference PrellPrell, 1985) for planktonic foraminifera (Reference Labeyrie, Duplessy and BlancLabeyrie and others, 1987).

Fig. 3. (a) Temperature record in northeast France far the last glacial cycle derived from La Grande Pile pollen record by Reference Guiot, Pons, de Beaulieu and Reille.Guiot and others (1989). (b) August sea-surface temperatures derived from planktonic foraminifera in North Atlantic core SU 90–39.

Forward Experiments: The Ice-Sheet Response to a Given Climate Series

We now test whether the pattern of ice-sheet fluctuation in Europe through the last glacial cycle suggested by Reference Mangerud and FrenzelMangerud (1991) is compatible with contemporary terrestrial and marine palaeoclimate reconstructions.

It is assumed that the ELA and SLAT for our ice sheet are linearly related to inferred palaeotemperature in each experiment, so that

(2)

where f ₁ and f ₂ are constants, Δλ is the change in ELA, and T is the deviation from modern values of the independent palaeotemperature estimate. We have used a value of SLAT at the glacial maximum of approximately −10°C, as suggested by the work of Reference Frenzel, Pesci and VelichoFrenzel and others (1992) for the area of northwestern Europe adjacent to the ice-sheet margin, although the model is relatively insensitive to the precise value. The proportionality constant (f ₁) in the ELA/palaeotemperature relationship is chosen as one which will drive the ice sheet as far as the maximum observed glacial extent along our transect at the time of the Last Glacial Maximum (LGM).

Forward simulations driven by atlantic SST reconstructions

In the initial simulation we assume that sliding does not occur at the ice/bed interface, even though the bed may be at the melting point. Initially, values of f ₁ = f _1a = 10⁶.0 m°C⁻¹ and f ₂ = 15.5 were used. This produces a maximum ice thickness along the transect at the LGM of about 3.4 km and a surface elevation of over 2.5 km. The profile along the transect of the ice-sheet surface and isostaticallv depressed bed, together with the internal temperature distribution, are shown in Figure 4. The coldest temperatures lie in an enclosed area beneath the divide because the ice sheet continues to grow after the temperature minimum. The ice-sheet thickness is large compared with values regarded by Reference Lambeck, Johnston and Nakada.Lambeck and others (1990) as compatible with the relative sea-level data.

Fig. 4. Ice-sheet profile at the LGM using a North Atlantic SST driving signal and showing internal temperature distribution. The completely enclosed low-temperature ice mass beneath the ice divide is a reflection of the fact that ice-sheet growth continued for some time after the start of climatic warming following the coldest part of the glacial period. The Scandinavian mountain mass is shown beneath the ice sheet to the left. The bed shows isostatic flexuring.

The time-dependent simulation of ice-sheet margin fluctuation is shown in Figure 5a compared with the geological reconstruction. In this initial simulation no sliding was permitted, irrespective of basal temperature, The simulated ice sheet has not responded to the assumed climatic forcing as rapidly as the real ice sheet appears to have done: the simulated glacial maxima show large lags; the modelled ice sheet shows little or no retreat dining periods of interstadial warming; and it fails to decay completely during the Holocene. These features tire relatively insensitive to variations between maritime/continental mass-balance extremes (Fig. 2).

Simulations were then also driven by values of f ₁ 0.83, 1.08 and 1.17 times f _1a, where f ₁ = f _1a = 106.0 m °C^−l. A limited glacial peak is achieved at about 60 000 BP in response to the strong forcing peak at 68 000 BP. There is then a long phase of reduced ice extent before the simulated buildup at about 30 000 BP. It is clear that the position of the simulated margin at 30 000–20 000 BP is very sensitive to small increases in forcing. The simulated LGM peak occurs about 10 ka after the geological evidence suggests that it occurred.

The simulations in Figure 5a For the larger values of f ₁ show a more rapid response to the LGM forcing hut overshoot the glacial maximum position. The prescribed mass-balance altitude feed-back is sufficiently strong with this level of forcing that once the modelled ice sheet reaches a critical size its growth is unchecked. The model halts when the domain size is exceeded. Again, none of the simulations succeeds in causing the ice sheet to disappear completely during the Holocene.

We have then introduced sliding at the ice/bed interface into the simulation. We prescribe a basal sliding velocity as a linear function of the shear stress, and allow sliding to occur wherever the basal temperature is at the melting point. The parameter in the sliding function is maximised such that it is the highest the ice sheet will support without producing unstable flow and collapse. Results are shown in Figure 5b for values of f ₁ = f _1a and f ₁ = f _1a × 1.17. The former value is that used to produce the non-sliding response which achieved the LGM extent in Figure 5a. With sliding it generates a much smaller ice sheet and is less responsive to forcing. This is because the lower profile of an ice sheet with sliding generates a smaller net mass balance compared with a steeper-slope, non-sliding ice sheet for the same ELA. Using a value of f ₁ 1.17 times larger, however, drives the ice sheet to simulate the LGM extent, bin again with a lag of 10 ka. The subsequent retreat is greater than in the non-sliding ease because of the lower surface profile.

Forward simulation driven by palaeotemperature reconstructions from La Grande Pile, France

This record is given as temperature departures from the present; we therefore use a value of f ₂ = 0.0 to equate it with the other signals. We have initially used the same value of f ₁ as in the first simulation (f _1a).

Figure 6a and Figure 6b show the simulated ice-sheet response to ELA Forcing based on the Grande Pile data for the no-sliding condition and the most extreme value of the sliding parameter. Figure 6a shows the response of the ice sheet to Forcing with no basal sliding. The forcing signal produced when f ₁ = f _1a, is the maximum the model will support before it produces an unstable ice sheet. When using the initial value of f ₁ = f _1a × 1.08 the modelled ice sheet becomes too large. Much-reduced ice volumes result if f ₁ is set smaller than f _1a.

Fig. 6. Forward simulation for La Grande Pile palaeotemperature sequence, (a) No-sliding cases. We were unable to produce the peak in ice extent at the LGM using this record. When f₁ = f_1a, the modelled ice sheet approximately matches the geological record but misses the extremes of growth and decay. If f₁ is set larger (f_1a × 1.08) in an effort to generate the peaks, the model becomes too large, (b) Sliding cases. The behaviour is similar to that foe the no-sliding cases other than that more forcing is acquired for similar-sized ice sheets. As before, the modelled ice sheet has a very flat response to the signal.

Figure 6b shows the response of the ice sheet to forcing when maximal basal sliding is included. With a value of f ₁ = f _1a, the ice sheet again remains relatively small and of roughly constant size through the glacial cycle. If f ₁ is increased to f _1a × 1.17 (Fig. 6b) there is a relatively good fit to the average size maintained by the ice sheet during the periods 100000–80000 and 50 000–12 000 BP although the peaks at 60 000 and 18 000 BP are not simulated. Moreover, there is a substantial decay during the early Holocene, although even under conditions of the softest bed compatible with dynamic stability, the ice sheet does not decay completely in the early part of the present Inter-glacial. Moreover, there is not n sufficiently dramatic reduction in ELA at 70000 BP to produce the 70000–60000 BP glacial maximum, nor after 30 000 BP to produce the LGM peak. Overall, although the Grande Pile temperature record between 115 000 and 85 000 BP produces a reasonable lit with the geological record, there is insufficient variability in the forcing after 85 000 BP to reproduce the principal features of the ice-sheet response to climate. If the value of f ₁ is increased, the modelled ice sheet again grows too much.

Our interim conclusion is that the greater variance in the SST records compared with the French pollen sites, and the more differentiated forcing it produces, are a better guide to the ice-sheet surface climate required to generate the geological record of ice-sheet fluctuation. None of these records, however, is able to generate a forcing able to destroy the ice sheets if we assume a linear relation between palaeotemperature records away from the ice sheets and ice-sheet surface climate.

Inverse experiments to establish the climate compatible with a given pattern of ice-sheet fluctuation

In these experiments, we have sought to establish, through the medium of the model, the nature of the ELA forcing functions which would be compatible with the geological evidence of ice-sheet fluctuation. We have not employed a formal inverse technique, but have run multiple forward simulations so as to match the geological record.

We used both a non-sliding and a maximally sliding model in these experiments. The match between model simulations and geological data, together with the ELA forcings required to generate the simulations, are shown in Figure 7a-b. These ELA forcings and the ELA variations derived from palaeotemperature records are compared in Figure 8.

Fig. 7. Inverse model for (a) m sliding, compared with geological reconstruction, and (b) extreme sliding. The high-frequency fluctuation of the ice margin compared with the no-sliding case reflects dynamic oscillation of the ice sheet (Reference Boulton, Smith and MorlandBoulton and Payne, 1994).

Fig. 8. ELA forcings used in various runs. The inverse signals (a and b) are those required to drive the model to achieve a reasonable match with the geological evidence. The others (c and d) are the ELAs used to drive the model derived by applying the temperature/ELA relationships to the palaeotemperature curves. In each case they show the ELAs used in no-sliding cases such that the ice sheet achieves a maximum size but does not become unstable. This means that each forcing signal produces a roughly comparable ice volume at maximum, bill other features differ. The palaeo-signal from core SU 90–39 is much closer to the derived “inverse” signals than the pollen record from La Grande Pile.

Conclusions

The ice-sheet model has been used to explore the compatibility of selected continental and marine palaeotemperature records with a geological reconstruction of the fluctuations of the margins of the European ice sheet through the last glacial cycle by using a numerical model of ice-sheet dynamics. We draw the following conclusions:

• (a) When the ice-sheet model is driven by northeast Atlantic SSTs or continental pollen palaeotemperatures, the ice sheet responds slowly to forcing and does not exhibit the strongly varying behaviour of Mangerud’s geologically based reconstruction. This may reflect:

  errors in the geological reconstruction of ice-sheet

  fluctuation;

  errors in palaeotemperature estimates;

  greater climatic variance over the ice sheet than in extra-glacial areas.

• (b) Early Holocene ELAs must have been significantly higher than modern values in order to produce complete degradation. The conventional view of Holocene palaeotemperatures (e.g. Reference MörnerMørner, 1980) suggests a thermal optimum at about 6000 BР. However, a recent detailed reconstruction by Reference Seret, Guiot, Wansard, de Beaulieu and Reille.Seret and others (1992) suggests that the Holocene thermal optimum may indeed have been reached between 7000 and 8000 BP.

• (c) The poor match in amplitude variation between simulations driven by palaeotemperature/ELA forcings and the geological record implies that the amplitude of temperature change beyond the ice sheet is much smaller than on the ice-sheet surface and that there may be a nonlinear relationship between ice-sheet climate and extra-glacial climate. Such a non-linearity is likely to be derived from feed-back processes between the ice sheet and atmosphere. The ocean SST record is, however, much closer to what is required to explain ice-sheet fluctuation than the continental pollen-based record. However, the depleted continental faunas and floras of cold periods are probably much less effective as palaeoclimate indicators than those of warm periods. It may therefore be that extra-glacial climates were more severe than currently recognised.

• (d) The periods of strong degradation at 85 000–80 000, 15 000 and 10 000 BP indicate anomalies between the glacial and palaeotemperature records. The glacial record requires an abrupt rise in the FLA to above modern values. It is possible that the 85 000–80 000 BP de-glaciation could be explained by shifts in the three-dimensional geometry of the ice sheet which our two-dimensional simulation cannot match, rather than dramatic warming. The known geometry of the ice sheet alter 15 000 BP makes this unlikely for the second and third periods. It is possible that there was very strong warming at these times, which the inertia of oceanic and floral response failed to record. Indeed, palaeotemperature reconstructions by Reference CoopeCoope and others (1977) based on fossil coleoptera suggest this. The other possible explanation for these apparent anomalies is that the ice sheet underwent dynamically driven collapse during these periods.

• (e) Ice-sheet surface climate shifts must lead the glacial maxima by several thousand years. Geological reconstructions in which the precise timing of glacial events is largely timed by correlations with records elsewhere should take this into account.

• (f) Easy sliding does not produce high rates of ice-sheet buildup. Compared with a non-sliding ice sheet, it is associated with a lower ice-sheet profile, which produces a smaller net positive mass balance and therefore a smaller buildup rate. Stronger climate forcing is therefore needed to maintain the same rate of buildup for a sliding ice sheet. The converse occurs during decay, however. A lower profile gives rise to a more negative mass balance, a greater rate of decay and an earlier final degradation.