Introduction
The waxing and waning of huge ice sheets in the Northern Hemisphere are among the most prominent changes in the Earth’s environment during the Pleistocene. Understanding the mechanisms responsible for changes in global ice volume and their 100 000 year periodicity are important topics in palaeoclimatology.
It is widely believed that the succession of glacials and interglacials is ultimately caused by changes in the solar insolation and the seasonal and lateral distribution of this insolation (e.g. Reference Imbrie, Berger, Imbrie, Hays, Kukla, Saltzman and Reidel Publishing Co.Imbrie and others, 1984; Reference ShackletonShackleton, 2000). The climate system itself responds in a non-linear way to these changes, due to internal feedback mechanisms, the non-linear growth and decay of ice sheets and possibly climatic thresholds. To investigate this non-linear response, independent time series of important climate quantities are required, that can be used in coupled numerical models. Among these quantities are sea level, temperature, the strength of the ocean circulation and the atmospheric concentration of greenhouse gases (e.g. Reference Bintanja, van de Wal and OerlemansBintanja and others, 2005a).
One of these climate quantities, temperature, can be estimated on glacial timescales using the famous relation between stable-isotope (δ18O and δD) content along ice cores and depositional temperature, as first described by Reference DansgaardDansgaard (1964). This approach is commonly used to generate temperature forcings for numerical-model simulations of the ice sheet from which the ice core was obtained (e.g. Reference Ritz, Fabre and Letréguilly.Ritz and others, 1997; Reference Huybrechts and De Wolde.Huybrechts and De Wolde, 1999; Reference Van de WalVan de Wal, 1999). Although popular, the method of using ice-core records for temperature estimates has its pitfalls, as pointed out by Reference CuffeyCuffey (2000). One of Cuffey’s main concerns is that necessary temperature corrections for changes in the ice sheet’s surface elevation can only be accounted for when the evolution of the ice sheet is simulated simultaneously with the temperature estimate. In other words, difficulties with the interpretation of the ice-core records arise because the ice core itself is part of a climate feedback. Another complication with the use of ice-core records is that they reflect mostly local estimates. The applicability of temperature forcings generated from ice cores is questionable when the temperature forcing is used on a continental scale for the simulation of ice sheets other than the ice sheet from which the ice core was obtained, because of large spatial variability.
Against the background of these limitations, Reference Bintanja, van de Wal and OerlemansBintanja and others (2005b) presented a new and independent method to estimate temperature on glacial timescales, based exclusively on sea-level observations. This temperature time series is believed to be representative of a substantial part of the Northern Hemisphere (>40º N). The method is based on two assumptions: (i) that sea-level changes during the last ice age have been caused to a large extent by the growth and decay of large ice-sheet complexes over North America and the Eurasian continent (e.g. Reference ClarkClark and others, 1993; Reference Huybrechts and T’siobbel.Huybrechts and T’siobbel, 1997) and (ii) that the evolution of these ice sheets is primarily influenced by temperature changes. These two assumptions imply a strong causal relation between changes in sea level and temperature variations during the last glacial cycle. Consequently, sea-level observations can be used as a constraint to reconstruct temperature using inverse modelling techniques. This requires sea-level observations, an ice-sheet model, a mass-balance model and a climate description. Figure 1 clarifies the process by showing the relation between these different models.
Overview of the models and variables involved in the inverse problem. Variables (or quantities) are SL (sea level), T(temperature), AT (temperature anomaly), P (precipitation), A (accumulation), M (ablation), h (surface elevation), α (albedo) and Q (insolation), functions of x (position), t (time), m (month, seasonality). Precipitation and temperature are present-day monthly fields (x,m,0). Arrows indicate the input and output of the different models. Time series of AT is the reconstructed temperature time series.
In this paper, we examine the effect of several physical processes that control the ice-sheet evolution and are therefore important when interpreting the temperature reconstructions. In our opinion, these processes are the key to understanding the relation between ice volume (or sea level) and temperature. The complexity of this relation is illustrated in Figure 2. This figure shows the reconstructed temperature by Reference Bintanja, van de Wal and OerlemansBintanja and others (2005b) as a function of sea level. A remarkable feature of this graph is that there are periods in which sea level drops while temperature rises, and also periods in which sea level rises while temperature drops. Another important observation from Figure 2 is that temperatures in the period following the Last Glacial Maximum (LGM) are higher for similar values of sea level than the temperatures in the period prior to the LGM. We show that most of these peculiarities can be explained by the fact that the ice sheets are not in equilibrium. In fact, discrepancies in the timescales for the changes in temperature and the mechanical and thermodynamical response of the ice sheet to these changes underlie this transient behaviour.
Reconstructed temperature vs sea level by Reference Bintanja, van de Wal and OerlemansBintanja and others (2005b). Arrows indicate the time direction in the experiment. Markers and labels indicate 10 kyr periods.
Physical mechanisms behind the complicated sea-level–temperature relation shown in Figure 2 are examined by comparing one dimensional (1-D) and three-dimensional (3-D) model simulations of ice sheets in the Northern Hemisphere over the last glacial cycle (described in the next section). It is shown that these simple models can produce reasonable temperature reconstructions over the last glacial cycle. Furthermore, we show that comparing the results of the 1-D model with those of the 3-D model provides useful insight, due to the simplicity of the 1-D model. Interactions among land surface and ice geometry, ice flow and ice surface mass balance are further explored by driving the 1-D and 3-D models with hypothetical stepwise changing sea-level patterns.
Models and Methods
The different models used in our study and their relations to each other are presented in Figure 1. The 3-D model, the inverse method and the surface mass-balance model are the same as in the earlier study by Reference Bintanja, van de Wal and OerlemansBintanja and others (2005a). Therefore, we only give a brief overview of these models and the results obtained with them. For a detailed description, the reader is referred to the papers by Reference Van de WalVan de Wal (1999) and Reference Bintanja, van de Wal and OerlemansBintanja and others, (2002, Reference Bintanja, van de Wal and Oerlemans2005a). Wherever our 3-D model is compared to other 3-D ice-dynamical models, it is called the ‘Bintanja model’.
3-D ice-sheet model
The dynamics of an ice sheet are described by the fundamental equations for conservation of mass, momentum and energy from classical continuum mechanics, together with an equation of state relating the density of the ice to the local pressure and temperature. Usually the equation of state is covered by the assumption of incompressible ice. The constitutive equation relating the strain rate to the stress field is Glen’s flow law. These formulations lead to a set of prognostic equations for ice thickness and temperature (e.g. Reference HuybrechtsHuybrechts, 1990).
In the Bintanja model the so-called shallow-ice approximation is followed (Reference HutterHutter, 1983). This approximation is suitable for ice masses with small aspect ratios and small surface and bedrock slopes. An implication is that longitudinal stresses can be neglected, which simplifies the prognostic equations describing the ice dynamics. More specifically, the local change in ice thickness can be described by a diffusion equation, with a diffusivity largely dependent on the surface slope. The shallow-ice approximation is widely used in ice-sheet models of the 1990s (e.g. Reference Huybrechts and T’siobbel.Huybrechts and T’siobbel, 1997; Reference Ritz, Fabre and Letréguilly.Ritz and others, 1997; Reference Marshall, Tarasov, Clarke and PeltierMarshall and others, 2000).
The energy-balance equation accounts for conductive and advective heat fluxes as well as for strain heating by deformation of ice. It governs the temperature distribution in the ice sheet, which feeds back to the ice dynamics due to the Arrhenius temperature dependence of the viscosity in Glen’s flow law and the temperature-dependent conditions for sliding at the bed. This is also a very common approach in ice-sheet modelling (see the studies referred to above). In the Bintanja model, sliding is set to occur when the basal temperature is within 1 K of the pressure-melting point.
Adjustment of the bedrock to the time-evolving ice load is an important process that has to be incorporated in an ice-sheet model, because it directly affects the surface elevation of the ice sheet and consequently the ice sheet’s surface mass balance. The Bintanja model uses the elasticlithosphere–relaxing-asthenosphere (ELRA) model of Le Reference Meur and HuybrechtsMeur and Huybrechts (1996), with global uniform values for the flexural rigidity of the lithosphere and the relaxation time of the asthenosphere.
The ice-shelf model is taken from Reference Oerlemans and van der VeenOerlemans and Van der Veen (1984) and adapted by Bintanja and others (2002). The model specifies the thickness of the ice shelves as a function of the distance to the grounding line and the ice thickness at the grounding line, ignoring dynamical aspects. Although simple, this model satisfactorily simulates the occurrence of ice shelves in narrow embayments and allows the ice sheet to expand over shallow seas. In contrast to some of the dynamical ice-shelf models, the frontal position is not explicitly prescribed.
The governing equations are solved using a finite-difference scheme on an equidistant rectangular grid of 20 km resolution. The domain contains the Eurasian and North American continents (excluding Greenland). The time-step is taken to be 1 month.
The Bintanja model has previously been used for the simulation of the behaviour of the Greenland ice sheet (Reference Van de WalVan de Wal, 1999). Furthermore, Reference Bintanja, van de Wal and OerlemansBintanja and others (2002) showed that the model is able to simulate the global evolution of ice sheets over the last glacial cycle in reasonable agreement with other studies (SPECMAP (spectral-mapping project) sea-level curve: Reference Imbrie, Berger, Imbrie, Hays, Kukla, Saltzman and Reidel Publishing Co.Imbrie and others, 1984; Reference Tarasov and PeltierTarasov and Peltier, 1997).
1-D ice-sheet model
The 1-D model simulates the time evolution of a hypothetical axially symmetric ice sheet. The ‘continent’ is initially cone shaped; the bed has a constant negative slope in the radial direction. The initial height at the centre of the bed and the bed slope are free parameters in the model. Although the geometry is simple, we do not assume a fixed ice-sheet surface profile, but calculate changes in the ice-sheet geometry using a prognostic equation for temporal changes in the ice thickness.
The model contains those processes that are considered necessary for the simulation of ice-sheet evolution on a continental scale. These processes are, following Reference Tarasov and PeltierTarasov and Peltier (1997), the feedback between mass balance and surface elevation, the mass-balance–albedo feedback and the adjustment of the bedrock to changes in the ice load.
The 1-D model is applied along a flowline from the centre to the edge of the domain. The domain contains 200 equally spaced (12.5 km) gridpoints. For every gridpoint, the ice thickness H and the bed elevation b are explicitly calculated at every time-step t (where the steps are 1 month apart). The rates of change of these quantities are given by the following prognostic equations:
In these equations, U is the mean horizontal ice velocity, B is the local ice equivalent surface mass balance, r is the distance from the centre, τ b = 3 kyr is the relaxation timescale of the asthenosphere, k is the ratio of the densities of the bedrock material and ice (set equal to 3) and b 0 is the initial local ice-free bedrock elevation.
The prognostic equation for the ice thickness (Equation (1)) is a continuity equation. With the assumption of a constant ice density, this equation is equivalent to conservation of volume. As in the 3-D model, the shallow-ice approximation is used to obtain an expression for the mean horizontal ice velocity due to the deformation of ice. The total horizontal velocity U is the sum of the deformation velocity and a sliding velocity, based on a Weertman-type sliding law:
where n is the exponent in Glen’s flow law, set to the commonly used value 3, and τ d the driving stress, which is proportional to the ice thickness H and the gradient of the surface elevation. The deformation parameter f d = 1.0 × 10¯15 Pa–3 a–1 is related to the proportionality constant between the deviatoric stress and the strain rate in Glen’s flow law (a measure of the viscosity of the ice). The sliding parameter fs = 3.0 × 10–11 Pa–3m2a–1 connects the sliding velocity to the driving stress.
The important adjustment of the bed elevation to time-varying ice loads is described by an expression representing the bed’s tendency to reach local isostatic equilibrium (Equation (2)). This equation is a commonly used description of bedrock adjustment in 1-D ice-sheet models (e.g. Reference Van der Veen and RotterdamVan der Veen, 1999).
Finally, in the 1-D model there is no interaction between the ice sheet and sea level. The cone-shaped ‘continent’ has no marine margins and the ice sheet can simply continue to grow below sea level. Thus, calving does not occur and there is no formation of ice shelves.
Surface mass-balance model
The model for the surface mass balance, defined as the difference between accumulation and ablation, is the connection between the temperature and the changes in the total ice volume and the sea level. Therefore, it plays a crucial role in the temperature reconstructions. We use two different mass-balance formulations. We comment on the suitability of these formulations in the discussion of the experiments in which the results with the different formulations are compared.
In the first formulation, which is only used in combination with the 1-D model, surface mass balance is a function of surface elevation. This formulation is referred to as the mass-balance-height (MBH) formulation. The positive feedback mechanism between surface mass balance and surface height is directly incorporated in this formulation. In the second formulation, used in combination with both the 1-D and 3-D models, accumulation and ablation are calculated independently. This formulation incorporates the local surface albedo, and is therefore referred to as the mass-balance-albedo (MBA) formulation. Obviously, this formulation incorporates the positive feedback mechanism between surface mass balance and albedo. The mass-balance–height feedback is incorporated in the second formulation in an indirect manner, because the ablation is a function of local surface air temperature and hence influenced by the local surface elevation.
The MBH formulation is adopted from Reference OerlemansOerlemans (2004b). Surface mass balance increases linearly with surface elevation, up to the critical height h c. If the surface elevation exceeds this height, mass balance remains constant. The continental mean temperature reduced to sea level, 7sl, is tied to the critical height using an empirical relation based on mass-balance studies in different parts of the world, including Severnaya Zemlya, Greenland, Svalbard, Southern Norway and the Alps (see Reference OerlemansOerlemans, 2004a, and the references therein).
In this equation, T sl is in °C, hc in m. Accumulation (i.e. the mass balance at the critical height) is coupled to the temperature and the radius of the ice sheet, implying drier conditions for lower temperatures and larger ice sheets (Reference Bintanja, van de Wal and OerlemansBintanja and others, 2002). The mass balance at the critical height determines the mass-balance gradient.
Ablation and accumulation are calculated independently in the MBA formulation. The ablation (M) is calculated as a linear combination of surface air temperature (T) and absorbed solar radiation (Q):
where α is the surface albedo, which depends on the depth of the snow layer on the ice sheet. The coefficients c 1, c 2 and c 3 (c 1 = 0.04 mw.e. a–1 °C–1, c 2 ¼ 5.13 × 10–3ma–1 (Wm–2)–1 and c 3 ¼ –0.32 m a–1) are derived from minimizing errors between mass-balance observations and model runs over Greenland (see Reference Bintanja, van de Wal and OerlemansBintanja and others, 2002). Accumulation is calculated from the precipitation rate using a temperature-dependent snow fraction and a correction for temperature changes: precipitation increases by 4% per degree temperature rise. This commonly used parameterization is based on the Clausius–Clapeyron relation between temperature and saturation vapour pressure (e.g. Reference Tarasov and PeltierTarasov and Peltier, 1997).
Climate description
Climate is described differently in the 1-D and 3-D models. These differences cannot be avoided because realistic climate fields are used in the 3-D model and this is impossible in an axially symmetric 1-D model. Climate is embodied in three quantities: (i) local surface air temperature, (ii) uncorrected precipitation rate and (iii) insolation. In the following paragraphs, the differences between the 1-D and the 3-D model with respect to these quantities are discussed.
-
i For both the 1-D and the 3-D model, temperature is the target of our inverse calculations. However, the nature of the temperature reconstructed with the inverse method differs between the two models. The 1-D model reconstructs the continental mean surface air temperature reduced to sea level, T sl. Local and monthly surface air temperatures T sa are calculated using a constant lapse rate γ in order to correct for changes in surface elevation h, and a monthly variation δ m:
(6)
The monthly variations in temperature are represented using an annual sinusoidal variation with its maximum in August, superimposed on the 100 year mean (reconstructed) temperature. The amplitude of the seasonal cycle in temperature is constant in time and over the entire domain. We acknowledge that annual ranges in temperature are usually smaller in maritime than in continental environments. Therefore, a possible improvement in the temperature calculation would be to have a seasonal cycle that decreases with the distance to the centre of the domain. However, the 1-D model is only used for conceptual simulations, and therefore we prefer our more transparent formulation.
Temperature anomalies, not absolute temperatures, are obtained when the inverse method is applied in combination with the 3-D model. This anomaly is constant over the entire domain. The anomaly is added to the present-day monthly mean surface air temperature field (US National Centers for Environmental Prediction re-analysis; Reference KalnayKalnay and others, 1996). The resulting temperatures are corrected for changes in the local surface elevation relative to the present-day surface elevation. The reader is referred to Reference Bintanja, van de Wal and OerlemansBintanja and others (2005b) for a discussion of the uncertainties in the temperature reconstructions associated with the important assumption of the spatially constant temperature anomaly.
-
ii Precipitation is corrected for temperature changes in both the 1-D model and the 3-D model. In the 1-D model, the uncorrected precipitation rate is constant in time. This is in contrast to the 3-D model, which uses spatially and monthly variable precipitation fields. Variations in precipitation over time are generated via the temperature change only.
-
iii The different geometries of the 1-D and 3-D models require different treatments of the variation in insolation. In the 3-D model, local insolation is calculated for every month as a function of latitude, including Milankovitch variations following Reference BergerBerger (1978). In the 1-D model, present-day monthly insolation at 65° N is used for the calculations of the mass balance.
Inverse method
In order to reconstruct temperature using sea-level observations an inverse method is employed. The essence of the inverse method is that a reconstructed temperature time series is built up from the start of the simulation by simultaneously running the ice-sheet model in forward mode and calculating temperature variations anticipating sea-level changes to come (see Reference Bintanja, van de Wal and OerlemansBintanja and others, 2005a). More specifically, temperature T(t) is calculated every 100 modelled years (the resolution of the sea-level interpolation):
where a is a constant tuning parameter;ΔSL100 is the difference between the sea level as calculated from the total ice volume in the model and the observed sea-level value 100 years later, and 〈T〉 N is the mean temperature over the last N periods of 100 modelled years. With this method, we run the ice-sheet model in the normal forward mode for successive periods of 100 years. At the end of a period, the volume of the ice sheet is translated into a sea-level value and Equation (7) is used to calculate the new temperature for the next period. The reconstruction is started with no ice as the initial condition.
The values for the parameters a and N in the inverse method are optimized by minimizing the root-mean-square difference between the observed and modelled sea level. The physical meaning of the parameters is hidden in the complex interaction between the different timescales in the physical system. For example, a too low value for a results in a too slow response of the ice sheet, while a too high value leads to large temperature fluctuations in the reconstruction. The parameter a thus relates the timescale of the physical system to the temporal resolution of the sea-level observations. The self-consistency of the method is best indicated by the close match between observed and modelled sea level. In all our simulations, this difference never exceeded 0.5 m for the optimized values for a and N.
A necessary and very important assumption in the inverse problem is that changes in global sea level are proportional to calculated changes in the modelled ice volume. Sea-level variations during the last ice age were caused to a large extent by variations in ice volume on the continents of North America and Eurasia. Reference Bintanja, van de Wal and OerlemansBintanja and others (2002) calculated that the ice volume on these continents accounted for 86% of the total difference between LGM global ice volume and present-day global ice volume. The remaining 14% of the total changes in ice volume relative to the LGM were caused by variations in ice volume in South America (4%), Tibet (2%), Greenland (4%) and Antarctica (4%). In the experiments with the axially symmetric 1-D model, ice volume equals the integral 
and is assumed to be a constant fraction of the global ice volume. This fraction is equivalent to the average contribution of the ice volume on North America and Eurasia to global ice volume at the LGM being 43%.
Input sea-level observations
We use the SPECMAP sea-level curve (Reference Imbrie, Berger, Imbrie, Hays, Kukla, Saltzman and Reidel Publishing Co.Imbrie and others, 1984) for reconstructions with the 1-D model. This record, based on analyses of oxygen isotope fluctuations in benthic foraminifera in marine sediments collected in the Atlantic and Pacific Ocean, is one of the first detailed records of sea level over the last glacial cycle ever produced. Nowadays, other sea-level records are available as well, based on analyses of coral terraces (e.g. Reference Lambeck and ChappellLambeck and Chappell, 2001) or inferred from isotopes and salinity changes in the Red Sea (Reference SiddallSiddall and others, 2003). Results in this paper do not qualitatively depend on the choice of the sea-level curve and any record with sufficiently fine resolution could be used.
Temperature Reconstructions and Discussion
Several temperature time series are reconstructed, using the MBH and MBA formulations, with different values for the tuning parameters a and N in Equation (7). We define the ‘reference reconstruction’ as the reconstruction obtained with those values for a and N that resulted in the sea-level pattern with the smallest root-mean-square (rms) deviation between the modelled and the observed (input) sea-level pattern. For the reference reconstruction, the rms deviation is smaller than 0.5 m.
The temperature reconstructions presented in Figure 3 show similarities as well as differences between the results obtained with 1-D and 3-D models. The most obvious agreement between the results with the two models is the timing of the temperature maxima and minima. Furthermore, the large variations in temperature correspond to large variations in sea level. Eye-catching differences between the results with the two models are the 5–7ºC shift of the temperature reconstructions of the 1-D model, in comparison with the results of the 3-D model, and the differences in the amplitude of the temperature variations.
Reference experiments. Reconstructed temperature with the 1-D model (solid lines) in time for (a) the MBH and (b) the MBA formulation, in comparison with the reference reconstruction of the 3-D model (dashed lines), and the input (SPECMAP) sea-level variations (dash–dot line, right axis).
The observed similarities in the results reflect the fact that both models respond in a consistent way to the imposed sea-level variations. The differences in the amplitude and offset of the response can be explained by the different design of the two models. The differences in the model design include both the differences in geometry and the differences in climate description.
As an example of the difference in geometry, the mean elevation of the land in the 3-D model is well defined, whereas in the 1-D model it depends on the rather arbitrary initial height at the centre of the domain and the slope of the bed. Changes in the initial elevation of the centre of the domain obviously lead to a shift in the reconstructed temperatures. It can be shown that increasing the slope of the bed leads to a larger amplitude in the temperature reconstructions.
Results with the 1-D model are sensitive to the uncorrected accumulation rate, because this parameter affects the ability of the ice sheet to grow rapidly. A lower uncorrected accumulation rate leads to a larger amplitude in the reconstructed temperatures. The same holds for the amplitude of the annual temperature variation. Increasing the amplitude leads, most importantly, to higher summer temperatures and consequently to lower overall reconstructed temperatures with a larger amplitude.
Another, more specific and important difference between the (reference) reconstructions of the 1-D model and the 3-D model is observed in the period 20–10 kyrBP. Temperatures in the reconstruction with the 3-D model are found to be about 3°C higher than in the experiments with the 1-D model (relative to the temperature differences in the rest of the experiment).
In this section, we show that systematically examining the importance of other differences between the 1-D and 3-D models provides valuable insight into the models as well as the feedback mechanisms between the ice-sheet evolution and its geometry.
Equilibrium states
The ice sheet is obviously not in an equilibrium or quasiequilibrium state at any time during the glacial period, while growing or decaying towards larger or smaller volumes. Delays between temperature and sea-level changes are caused by the discrepancy in the timescales for changes in mass balance (negligible) and ice thickness (thousands of years) and the timescale associated with the adjustment of the bedrock (thousands of years). These transient effects explain most of the complexity of the relation between sea level and temperature, as illustrated in Figure 2.
Processes other than the phase difference between temperature forcing and ice-sheet response may contribute to the complexity of the sea-level–temperature relation. We explore this possibility using sea-level patterns, which allow the ice sheet to evolve towards equilibrium states. These sea-level patterns consist of periods in which sea level remains constant (Fig. 4). The periods are 25 and 50 kyr for the 3-D and the 1-D model, respectively. The only possible solution for a non-changing sea-level value in time is an ice sheet in steady state; the mass balance integrated over the ice sheet is zero and the temperature is constant. This temperature is referred to as the equilibrium temperature for a given value of sea level.
Stepwise changing sea-level pattern as used for experiments with steady-state ice sheets. The time interval between the steps is 25 kyr. In experiments with the 1-D model 50 kyr periods are used.
Figure 5 shows the equilibrium temperature as a function of sea level for the equilibrium experiment conducted with both the 1-D model (MBH formulation) and the 3-D model. In the experiment with the 1-D model, a linear dependence of the equilibrium temperature on sea level is found. The result is also symmetric: equilibrium temperatures in the first half of the experiment, during which sea level is decreasing, are identical to those in the second half of the experiment, during which sea level is increasing. These results are not observed in the experiment with the 3-D model. Instead, a non-linear relation is found, in which a clear hysteresis is observed. These discrepancies may be explained by several differences between the 1-D and 3-D models. In the following subsections the importance of some of these differences – mass-balance formulation, thermodynamics, geometry and spatial climate variations is discussed.
Equilibrium temperature as a function of sea level, using the MBH formulation for (a) the 1-D model and (b) the 3-D model. Arrows indicate in which way sea level is stepwise changed.
The hysteresis in our experiments is not identical to the hysteresis described elsewhere, for example by Reference Letreguilly, Huybrechts and ReehLetréguilly and others (1991). They investigated steady-state properties of the Greenland ice sheet and found that the temperature– ice-volume relation shows a hysteresis, but in their reconstructions the starting point was always either the present-day ice sheet or the ice-free Greenland bedrock. In our experiments, we simulate transitions from one equilibrium state to another. In this case, the occurrence of hysteresis is not an evident result but rather a useful point of comparison in evaluating various aspects of the different models.
The importance of mass-balance formulation
Biases in the two different mass-balance formulations used in combination with the 1-D model can be explored using the step-function sea-level time series. Assuming that the parameters are chosen appropriately, the MBA formulation should be more realistic than the MBH formulation. The latter does not account for seasonal variations in temperature and insolation, or albedo changes due to the covering of the bedrock by snow and ice. We repeat the equilibrium experiment with the 1-D model using the MBA formulation (Fig. 6) instead of the MBH formulation (Fig. 5a). In this experiment, the linearity in the results is clearly not obtained, but neither is the hysteresis that we find in the experiment with the 3-D model (Fig. 5b).
Equilibrium temperature as a function of sea level, for the 1-D model using the MBA formulation. Arrows indicate in which way sea level is stepwise changed. Note: for SL = –50 m in the first half of the experiment (decreasing sea level) and for SL = –70 m in the second half of the experiment (rising sea level) no equilibrium state is reached. Dashed lines are used for transitions from one equilibrium state to another, where the equilibrium state in between is missing.
We conclude (i) that the high level of idealization in the MBH formulation is the most likely cause of the linearity in the result of the equilibrium experiment and (ii) that the hysteresis is the result of geometric differences between the 1-D and 3-D models. In particular, the axial symmetry in the 1-D model excludes multiple solutions for steady-state ice sheets.
The importance of thermodynamics
The importance of including thermodynamics in the ice-sheet model can be tested with the 3-D model. Here, we consider two experiments: the reference equilibrium experiment for the North American continent and an experiment with the same domain in which the dependence of ice dynamics on temperature is switched off. We find that the temperature range is smaller in the second experiment and that both experiments show a hysteresis (Fig. 7).
Equilibrium temperature as a function of sea level, for the North America reference experiment (thermodynamics included; solid line) and the experiment in which thermodynamic effects are excluded (dashed line). Arrows indicate in which way sea level is (stepwise) changed.
The smaller range of temperature in the ice sheet modelled without thermodynamics is explained as follows. For the experiment excluding thermodynamics, we chose the flow parameter value 6.0 × 10¯17 Pa–3a–1 from Reference Van de WalVan de Wal (1999). This value leads to a similar volume for the Greenland ice sheet at the present-day conditions as was calculated from the 3-D thermodynamic ice sheet with the typical values for Glen’s flow parameter = 1.14 × 10¯5 Pa–3a–1, enhancement factor = 3 and creep activation energy = 60 kJ mol–1 for ice above –10°C (see e.g. Reference Marshall, Tarasov, Clarke and PeltierMarshall and others, 2000). The constant flow parameter leads to a constant stiffness of the ice, which results in an underestimate of the stiffness in the central parts and an overestimate in the marginal zone and hence a flatter ice sheet. A flatter ice sheet is more sensitive to climate change than an ice sheet with steeper surface gradients, and hence smaller temperature ranges are required to obtain equally large steps in sea level.
The importance of realistic geometry
There are two types of differences between the realistic geometry used in the 3-D model and the highly idealized cone-shaped continents used in the 1-D model. First, the geometry is different on a local scale. Specifically, there are significant differences in both the local gradients in surface elevation and the hypsometry between the 1-D and 3-D models. Local gradients are crucial for the ice dynamics (e.g. Reference Oerlemans and BalkemaOerlemans, 2001), while the hypsometry is important for the mass balance (e.g. Reference Marshall and ClarkeMarshall and Clarke, 1999). Second, the geometry is different on a continental scale. Different ice-sheet complexes on the Eurasian and North American continents contributed to the sea-level variations during the last glacial cycle, each of them formed by the merging of smaller ice sheets. The number of ice sheets that form on the continents, the locations at which these ice sheets form and how they interact (i.e. merge or separate) may all be important to the temperature reconstruction. The consequence of merging and separation of ice sheets can be a redistribution of the total ice volume over the different ice sheets. The 1-D model is unable to simulate this process. This is the key to understanding the different results from the 1-D and 3-D models.
The importance of a realistic ice distribution between the two different continents is investigated with the 3-D model. Reference Bintanja, van de Wal and OerlemansBintanja and others (2005b) presented an overview of ice distributions at the LGM, simulated with the 3-D model using the inverse method. These ice distributions are not fully in agreement with geomorphological evidence and results from rebound studies (Reference ClarkClark and others, 1993; Reference PeltierPeltier, 1994). This problem can be solved to a large extent by imposing a hypothetical temperature difference between the continents of Eurasia and North America. This adjustment has little effect on the reconstructed temperature time series. Significant deviations of about 2°C from the reference temperature time series emerge only after the LGM. For this reason, we conclude that the distribution of ice over the different continents is not important for the temperature reconstruction.
The importance of a realistic ice distribution on a single continent is examined in this study. The hysteresis in the temperature reconstructions with the 3-D model (Fig. 8) can be understood by considering the corresponding distribution of ice sheets (Fig. 9). In the experiment represented by Figures 8 and 9, sea level is stepwise changed and its variations are assumed to be proportional to ice volume variations on the North American continent. It is evident that the ice distributions do not change symmetrically about the inflection in the sea-level time series. The most remarkable of these changes in the ice distribution is found in the eastern part of the continent. Three small ice sheets form in regions with high accumulation and low temperature; they expand when the temperature decreases and merge, forming one bigger ice sheet. The newly formed highland is preserved under the improved conditions for glaciation – due to elevation and albedo changes – when temperature increases again. The merging of small ice sheets appears to reduce the climate sensitivity of the total ice volume. This explains the hysteresis in the sea-level– temperature relation.
Equilibrium temperature as a function of sea level for the experiment with the 3-D model, in which sea-level changes are assumed to be proportional to ice volume variations over North America. Arrows indicate in which way sea level is stepwise changed.
Ice distribution over North America for a stepwise changing sea-level pattern. Ice distributions are plotted at the end of certain sealevel steps. Sea-level (SL) values are indicated. The top row (a) refers to the first half of the experiment, the bottom row (b) to the second half.
This conclusion is supported by the observed discrepancy between the reconstructed temperatures after the LGM in the experiment with the 3-D model, in which a longitudinal temperature gradient was applied, and the reference experiment with the 3-D model (see second paragraph of this subsection). The merging of two large ice-sheet complexes on the North American continent (see Reference Bintanja, van de Wal and OerlemansBintanja and others, 2005b, fig. 6) causes a reduction of the climate sensitivity and therefore higher temperatures after the LGM for the same sea level. Furthermore, the hysteresis explains the observed temperature difference after the LGM between the 1-D and the 3-D model (see end of introduction to this section, above).
The importance of spatial climate variations
The importance of the connection between ice-sheet geometry and spatial climate variations is further evaluated by imposing spatial homogeneous fields for potential temperature and precipitation in the 3-D model. For this experiment, we use the same stepwise changing sea-level pattern as in the previous experiment. The resulting temperature reconstruction and ice distributions are shown in Figures 10 and 11. In this scenario, glaciation initiates in the western highlands. Merging and separation of ice sheets are observed during glaciation and deglaciation (Fig. 11). Without the interaction between ice-sheet geometry and temperature, the total ice volume is not redistributed in this process. In the end, those areas that have initially the best conditions for glaciation remain glacierized for the longest periods. The hysteresis has almost disappeared in the experiment with the homogeneous climate (Fig. 10).
Equilibrium temperature as a function of sea level for the experiment with the 3-D model (North America), in which spatial climate variations are excluded. Arrows indicate in which way sea level is stepwise changed.
Ice distribution over North America for a stepwise changing sea-level pattern for the experiment in which spatial climate variations are excluded. Ice distributions are plotted at the end of certain sea-level steps. Sea-level (SL) values are indicated. The top row (a) refers to the first half of the experiment, the bottom row (b) to the second half.
Conclusions
A new method to estimate continental mean ice age temperatures was developed and presented by Reference Bintanja, van de Wal and OerlemansBintanja and others (2005b). Here we expand on that work to evaluate the importance of climate feedbacks and mass-balance formulations to the relationship between sea level and temperature.
Using the 1-D model, we are able to show that the oversimplified mass-balance–height formulation leads to a linear relation between sea level and temperature in experiments with steady-state ice sheets. In this formulation, surface mass balance is a function of surface elevation only and it lacks a seasonal cycle and the important feedback between mass balance and surface albedo.
Results with the 3-D model show an interesting hysteresis in the sea-level–temperature relationship, which is not observed in experiments with the 1-D model. We conclude that this hysteresis (asymmetry between the equilibrium temperature in warming and cooling climates) occurs when ice volume is redistributed as a consequence of the processes of merging and separation of ice sheets. The effect of the merging of ice sheets in a certain area leads to a reduction of the climate sensitivity of the total ice mass in that area. Hysteresis does not occur if the topography does not support a relatively flat interior dome, formed by the merging of small ice sheets, that tends to preserve itself under warming conditions.
The results show the importance of correctly modelling the locations at which ice sheets form and how they interact. For example, the hysteresis effect provides a direct explanation for the result shown by Reference Bintanja, van de Wal and OerlemansBintanja and others (2005b), that reconstructed temperatures are about 2°C higher after the LGM if the climate setting causes the large ice sheets on the North American continent to merge. We have examined some of the physical mechanisms behind the geometric effects, i.e. the importance of mass-balance formulation, thermodynamics, realistic geometry and spatial climate variations. The understanding of these mechanisms is the starting point for further improvements of the temperature reconstructions. In particular, incorporating changes in the spatial distribution of temperature and precipitation derived from global circulation models run for glacial conditions (including the dynamic effect rather than expressing climate relative to present day) may lead to a better representation of the ice distribution over the last glacial cycle and therefore an improved temperature reconstruction.
Acknowledgements
Large parts of this work have previously been described in a Master thesis by F.W. at the University of Groningen. We thank H. Meijer for suggestions and improvements to that thesis. We are grateful to H. Oerlemans for useful comments. We were able to improve the manuscript with the constructive suggestions of three anonymous reviewers. During his research period at Utrecht University, F.W. was able to join a summer school on glaciers and ice sheets in Karthaus, Italy, in September 2003 and the European Geosciences Union 1st General Assembly, Nice, France, in April 2004. Both very instructive experiences were financially supported by the Spinoza Award that H. Oerlemans received in 2001 from the Netherlands Organization for Scientific Research (NWO).
