Interaction of ice sheets and climate during the past 800 000 years

Abstract. During the Cenozoic, land ice and climate interacted on many different timescales. On long timescales, the effect of land ice on global climate and sea level is mainly set by large ice sheets in North America, Eurasia, Greenland and Antarctica. The climatic forcing of these ice sheets is largely determined by the meridional temperature profile resulting from radiation and greenhouse gas (GHG) forcing. As a response, the ice sheets cause an increase in albedo and surface elevation, which operates as a feedback in the climate system. To quantify the importance of these climate–land ice processes, a zonally averaged energy balance climate model is coupled to five one-dimensional ice sheet models, representing the major ice sheets. In this study, we focus on the transient simulation of the past 800 000 years, where a high-confidence CO2 record from ice core samples is used as input in combination with Milankovitch radiation changes. We obtain simulations of atmospheric temperature, ice volume and sea level that are in good agreement with recent proxy-data reconstructions. We examine long-term climate–ice-sheet interactions by a comparison of simulations with uncoupled and coupled ice sheets. We show that these interactions amplify global temperature anomalies by up to a factor of 2.6, and that they increase polar amplification by 94%. We demonstrate that, on these long timescales, the ice-albedo feedback has a larger and more global influence on the meridional atmospheric temperature profile than the surface-height-temperature feedback. Furthermore, we assess the influence of CO2 and insolation by performing runs with one or both of these variables held constant. We find that atmospheric temperature is controlled by a complex interaction of CO2 and insolation, and both variables serve as thresholds for northern hemispheric glaciation.


Introduction
Earth's climate is characterised by glaciations and deglaciations of the Northern Hemisphere during the past 800 000 years. For this period, the Dome C ice core provides a record of CO 2 , CH 4 , N 2 O and deuterium-based Antarctic temperature (EPICA community members, 2004;Loulergue et al., 2008;Spahni et al., 2005;Jouzel et al., 2007). The Antarctic atmospheric temperature is tightly connected to a reconstruction of sea level over the past five glacial cycles, based on Red Sea bulk sediment analysis (Rohling et al., 2009). Ice coring performed on the Greenland Ice Sheet has provided records of several water isotopes (GRIP members, 1993;NGRIP members, 2004;NEEM community members, 2013). These indicate local temperature variability, although their interpretation is ambiguous (see e.g. different glacial-interglacial reconstructed temperature variability in Kindler et al., 2014, andSimonsen et al., 2011). Greenland and Antarctic climate have been shown to correlate closely to each other (EPICA community members, 2006). Additionally, important information arises from deep-sea sediment records of benthic δ 18 O (Zachos et al., 2008;Lisiecki and Raymo, 2005). These serve as an independent proxy for ice volume only after correction for the also contained deep-sea temperature signal. Köhler et al. (2010) combined data on the radiative forcing of solar insolation, greenhouse gases (GHGs) and modeldeduced ice volume, sea ice, vegetation and dust, to identify the importance of these components in governing temperature variability. This hybrid model-data approach has given a comprehensive overview of the importance of the different factors that played a role in the past glacial-interglacial cycles. However, their inferred global cooling between present day and the Last Glacial Maximum of 6.4-9.6 K was considerably larger than recent proxy-data-reconstructed ones of 5.8 K (Schneider von Deimling et al., 2006), 3.0 K (Schmittner et al., 2011) and 4.0 K (Annan and Hargreaves, 2013). Model studies can help in understanding this mismatch, which is possibly explained by background climate-state dependency of water vapour, lapse rate and cloud feedback strengths (Köhler et al., 2010). In addition, climate models offer the possibility of assessing the importance of its components and the interaction between processes by comparing results with those obtained from runs with one or more variables kept constant. This can also provide insight into the processes governing polar amplification, the ratio of polar to global temperature change.
Modelling ice ages is a classic issue in climate dynamics. Milankovitch recognised that Earth's climate is influenced by changes in solar radiation through three orbital parameters: eccentricity, precession and obliquity (Milankovitch, 1930), all varying with distinct frequencies. In 1941, he was the first to present calculations on the latitude-dependent seasonal solar insolation at the top of the atmosphere based on variations in these parameters (Milankovitch, 1941). Since then, these calculations have been improved (e.g. Berger, 1978;Laskar et al., 2004). Milankovitch proposed that, as a result of changing insolation, summer temperature on the Northern Hemisphere would vary and lead to glacial cycles. Indeed, orbital frequencies were discovered in climate records from marine sediments (e.g. Hays et al., 1976;Shackleton, 2000) and ice cores (Jouzel et al., 1989). The influence of solar insolation changes on the waxing and waning of ice sheets was discussed in many studies of varying complexity (Calder, 1974;Imbrie and Imbrie, 1980;Paillard, 1998) as well as in studies with ice-dynamical models (Pollard, 1978;Oerlemans, 1980Oerlemans, , 1982. However, simulating ice ages solely with solar variability remains difficult because of the lack of spectral power in the eccentricity (100 kyr) band (Imbrie et al., 1992). More recently, research has focused on insolation forcing and internal climate feedbacks between climate, ice sheets, and terrestrial dust (Ganopolski and Calov, 2011) or the lithosphere-asthenosphere system (Abe-Ouchi et al., 2013). Several other studies stress the importance of CO 2 forcing as an important driver of climatic changes (e.g. Paillard and Parrenin, 2004). In addition, Ritz et al. (2011) showed the importance of ocean circulation using a coupled ocean-atmosphere model with ice sheet forcing. Still under debate is the precise role of long-term ice-sheet-climate feedbacks in the climatic response to CO 2 and insolation forcing, which is the focus of this study.
Obviously, land ice and climate are intimately linked. When local temperature is cold enough, winter snow accumulation can exceed summer melt, and as a result an ice sheet is formed. This ice sheet alters the climate in several ways. Two effects are most important on long timescales. Firstly, the ice sheet reflects more shortwave radiation than the underlying ground, leading to lower temperatures. This is known as the ice-albedo feedback. Secondly, the ice increases the height of the surface. The surface is then positioned higher in the atmosphere, where temperatures are lower. This is called the surface-height-temperature feedback. The latter effect is partly counteracted by lowering of the ice sheet bed by the induced ice loading. The cooling, effectuated by the feedbacks, consequently leads to more ice L. B. Stap et al.: Interaction of ice sheets and climate during the past 800 000 years 2137 formation. The coupled nature of land ice and climate calls for combined modelling.
Earlier modelling studies addressing long-term global ice volume variations used stand-alone ice-dynamical models Bintanja and Van de Wal, 2008;De Boer et al., 2010. In these studies, the only information on climate consisted of globally uniform temperature perturbations with respect to pre-industrial (PI) climate. They obtained reconstructions of northern hemispheric temperature, ice volume and sea level from benthic δ 18 O data using an inverse routine. From that, a consistent CO 2 simulation was constructed based on a simple relation between the northern hemispheric temperature and CO 2 data from several proxy records (Van de Wal et al., 2011). These studies, however, lacked information on the meridional temperature distribution. Moreover, the representation of deep-ocean temperatures was very simplistic and CO 2 was not incorporated into the model. Therefore, it was not possible to study the individual contributions of insolation, CO 2 and land ice on temperature and sea level.
In this study, we no longer use the inverse benthic δ 18 O routine. Instead, we aim to quantify the influence of the aforementioned ice-sheet-climate feedbacks by coupling a climate model to a dynamical ice sheet model. Several recent studies have used a similar approach, mostly using three-dimensional ice sheet/shelf models of a single ice sheet, asynchronously coupled to Earth system models of intermediate complexity (EMICs) or global climate models (GCMs). In other cases, the ice sheet/shelf model is forced with look-up tables of climate, calculated beforehand with GCMs. These studies focused predominantly on relatively short timescales, e.g. the last glacial cycle (Helsen et al., 2013;Ganopolski et al., 2010) and the termination of the last interglacial (Herrington and Poulsen, 2011). Alternatively, they addressed a single climatic event, such as the Eocene-Oligocene boundary (DeConto and Pollard, 2003;Wilson et al., 2013), the middle Miocene transition (Langebroek et al., 2009) or the possible future retreat of the Greenland Ice Sheet (Ridley et al., 2005). Transient runs at timescales spanning ∼ 100 kyr (Abe-Ouchi et al., 2013;Ganopolski and Calov, 2011) to ∼1-10 Myr  have also been performed, but to our knowledge no study simulated all the major ice sheets simultaneously in a fully coupled fashion. Here, we combine a zonally averaged energy balance climate model (ZEBCM; Bintanja, 1997) and a onedimensional ice sheet model (ISM;De Boer et al., 2010). The relative simplicity and consequent short computation time of these models enable us to perform a large suite of model tests and 800 kyr long transient simulations of global climate for all five major ice sheets. Furthermore, this approach allows us to study different feedback mechanisms. Our set-up bears resemblance to the LLN 2-D climate model, which also consists of a zonally averaged climate model coupled to a onedimensional ice sheet model. The LLN 2-D model was used to study ice ages (Gallée et al., 1992;Berger et al., 1998;Loutre and Berger, 2000;Pépin et al., 2001) and the past 3 Myr (Berger et al., 1999). However, our model includes the Southern Hemisphere and Antarctica, and we force CO 2 using a longer ice core record. We will compare our results with these studies in Sects. 3.5 and 3.6.
From the results of the coupled model, we obtain a reconstruction of ice volume, sea level and meridional temperature, driven by CO 2 and insolation. Firstly, we will present this reconstruction, along with sensitivity tests to various aspects of the coupled model and a comparison with existing data from proxies and models. Secondly, we will compare simulations with uncoupled and coupled ice sheets to assess the influence of climate-ice-sheet feedbacks. By only passing ice extent or ice height to the climate model, we will differentiate between the effects of the ice-albedo feedback and the surface-height-temperature feedback. Moreover, we will address polar amplification and climate sensitivity over these 800 kyr. Finally, the relative contributions of insolation and CO 2 to climate variability will be studied by performing model runs with one of these variables held constant.

Zonally averaged energy balance climate model
The zonally averaged energy balance climate model (ZE-BCM) used in this study is described in detail by Bintanja (1997). Briefly, it consists of a coupled atmosphere model, ocean model and sea-ice model, all zonally averaged. It is based on the model of North (1975). The atmosphere model has a latitudinal resolution of 5 • and one layer in the vertical. It has a time step of 0.5 days and includes a radiative transfer scheme, describing radiative forcing by insolation, CO 2 , and cloud cover and optical depth . Turbulent heat transport and divergence of meridional energy transport are parameterised. No hydrological scheme is included. The land cover is subdivided into a forest, grassland and permanent land ice fraction, and affects only the surface albedo. While the glaciated land cover may vary, the ratio of grassland to forest is specified by present-day conditions and remains unchanged; this means that vegetation dynamics are not included. Snow cover on each surface type is calculated as a function of temperature. The ocean model has the same latitudinal resolution and consists of six layers in the vertical. It describes the zonally averaged thermohaline circulation, with water flowing poleward in the upper layer and equatorward in the lowest layer. Upwelling takes place in the entire ocean basin, except in the two latitudinal boxes closest to the pole, where the water sinks. The sea-ice model, with a resolution of 1.25 • , is purely thermodynamical. Where the sea surface temperature drops below the freezing temperature of 35 PSU saline water (271.22 K), the ocean becomes (partially) covered with sea ice of a specified thickness. Similar to the vegetation, the sea ice is subject to snow cover. Radiative forcing per CO 2 doubling (W m −2 ) 3.5 3.7 Enhancement of CO 2 radiative forcing by other GHGs (%) -30 D 0 (north) Diffusion coefficient for the Northern Hemisphere (m 2 yr −1 ) 4.0 × 10 11 2.0 × 10 10 D 0 (south) Diffusion coefficient for the Southern Hemisphere (m 2 yr −1 ) 1.0 × 10 11 0.5 × 10 11 D z Diffusion coefficient for vertical ocean mixing (m 2 yr −1 ) 2.0 × 10 3 5.0 × 10 3 τ cl Cloud optical thickness parameter 3.11 3.5 α SI Maximum albedo of sea ice 0.75 0.65 In an earlier study, the model was used to study presentday (PD) and Last Glacial Maximum (LGM) climate (Bintanja and Oerlemans, 1996). In order to perform transient simulations, we now use insolation and CO 2 as timedependent input variables. Monthly mean insolation follows from Laskar et al. (2004), averaged to 5 • resolution and updated every 1000 years. The EPICA record (Jouzel et al., 2007), interpolated to yearly values, is used as CO 2 input. Cloud cover and optical depth remain at PD value throughout the run. The climate sensitivity of the original model was 1.9 • C per CO 2 doubling (Bintanja, 1997), on the lower end of the 1.5-4.5 • C range reported in IPCC-AR5 (IPCC, 2013). In this study, a higher sensitivity of 2.2 • C is established by two alterations: the radiative forcing induced by a doubling of CO 2 is increased from 3.5 to 3.7 W m −2 , in agreement with Myhre et al. (1998); an additional water-vapour feedback term of 1 W m −2 K −1 is included in the parameterisation of clear-sky upgoing longwave radiation, based on Ritz et al. (2011). To account for non-CO 2 greenhouse gases, the radiative forcing of CO 2 is multiplied by a factor of 1.3. This is based on the mean enhancement of radiative forcing by CH 4 and N 2 O over the last 800 kyr derived from proxy data, as analysed by Köhler et al. (2010). The sensitivity to this factor will be analysed in Sect. 3.3. As in Bintanja (1995), the strength of the thermohaline circulation (C) is made a function of the surface water density difference ( ρ) between the equatorial and polar (60-80 • N and 50-70 • S) waters: where the subscript 0 denotes PI values and ζ is a strengthcontrolling scaling parameter, set to 6. The density difference consists of a temperature-gradient term and a pseudosalinity term. The temperature gradient is calculated from the equation of state for 35 PSU saline water (Gill, 1982). The pseudo-salinity term represents the effect of polar-water freshening by meltwater, by multiplying the change in polar surface-water temperature by a factor of −0.09. The midpoint of the thermohaline overturning is set to 5 • S, in accordance with Trenberth and Caron (2001). When the density difference between south-polar and north-polar waters becomes less than 0.55 kg m −3 , the midpoint is shifted to 10 • S. If the density difference drops below 0.51 kg m −3 , the midpoint is set to 15 • S. This mimics the effect of relative increase in the Atlantic cross-equatorial flow. It returns to 5 • S (10 • S) when the density difference exceeds 0.55 kg m −3 (0.51 kg m −3 ). This midpoint-shift mechanism is included to increase inter-hemispheric heat transfer towards the Northern Hemisphere, leading to lower glacial temperatures in Antarctica that are more in line with the observations (see Sect. 3.2). In Sect. 3.3 we discuss its effect on northern hemispheric temperature, sea level and northern hemispheric overturning strength. Further adjustments with respect to the original setup of Bintanja (1997) have been made to the diffusion coefficients for the Northern Hemisphere (D 0 (north); set to 0.2 × 10 11 m 2 yr −1 ), Southern Hemisphere (D 0 (south); set to 0.5 × 10 11 m 2 yr −1 ) and vertical ocean mixing (D z , set to 5 × 10 3 m 2 yr −1 ). This ensures a better transient temperature response of the ocean than was obtained by using the original values (see Sect. 3.2). Agreement with PD observations is maintained by setting cloud optical thickness parameter τ cl to 3.11 and the maximum albedo of sea ice to 0.65 (see Sect. 3.1). Table 1 summarises the alterations with respect to Bintanja (1997). Ice fraction and surface height are adjusted every 1000 years, based on the output of the ice sheet model. Land fraction is considered constant in time. If, however, the ice fraction exceeds the land fraction, it is set equal to the ice fraction. This means that all ice is considered to be groundbased, in congruence with the ice sheet model. Once the ice fraction returns to smaller values, the land fraction is set back to its original value.

Ice sheet model
The ice sheet model (ISM) consists of five hypothetical, axisymmetrical "continents" on which ice sheets can grow. These represent the five major ice sheets of the Cenozoic: the Eurasian (EuIS), North American (NaIS), Greenland (GrIS), West Antarctic (WAIS) and East Antarctic (EAIS) ice sheets. The continents are initially cone-shaped, thus being uniquely described by the initial height of the centre (H cnt ) and the slope (s) of the initial bed. These quantities determine the maximum size of the ice sheet and its sensitivity to temperature (Table 2). For the GrIS, WAIS and EAIS we used linear initial slopes as in De Boer et al. (2010), while for the EuIS and NaIS they are parabolic so as to realise better temperature sensitivity (see Sect. 3.2). Each continent has an equally spaced grid (Table 2). It is forced by the mass balance, which is the annual sum of monthly values of accumulation minus ablation.
Because the ZEBCM does not include a hydrological scheme, accumulation is calculated internally by the ISM. It is determined as a temperature-dependent fraction of precipitation P (Bintanja et al., 2002). The precipitation parameterisation is based on the Clausius-Clapeyron equation, taking into account the fact that precipitation increases at higher temperatures (T ) by 4 % per • C. Furthermore, precipitation decreases when the radius of the ice sheet (R) grows (Oerlemans, 2004): (2) P 0 and R c are a reference precipitation and critical radius, chosen to match PD observations of precipitation and its sensitivity to changes in radius (Table 2). Ablation (M) follows from an energy-balance model approach (Van den Berg et al., 2008), incorporating the effects of monthly mean temperature and local radiation (Q), varying over time: Q is taken at 65 • N for the EuIS and NaIS, 70 • N for the GrIS, and 65 • S for the EAIS and WAIS. C abl determines the threshold for the onset of ablation. It serves as a tuning parameter for ice sheet inception (Table 2). When the mass balance becomes positive, an ice sheet forms. The flow of the ice is then described by the continuity equation: where H is the ice thickness, r is the distance from the centre and B = P − M is the mass balance. Mean horizontal velocity (U ) results from the driving stresses, as obtained from the commonly used shallow ice approximation (SIA) and a Weertman-type sliding velocity at the bed. Bedrock adjustment is taken into account by applying  Table 2). Furthermore, glacial inception in East Antarctica starts when the CO 2 concentration is close to 700 ppm (at PD insolation), in accordance with DeConto et al. (2008); West Antarctica and Greenland obtain their maximum size during the LGM. The North American and Eurasian ice sheets do not have a geographically imposed maximum size. The tuning parameters are chosen to ensure a total sea level drop of approximately 120 m at the LGM, and the disappearance of the North American and Eurasian ice sheet for pre-industrial conditions.

Coupling scheme
A two-way coupling of the ZEBCM and ISM is established as follows (see Fig. 1): first the ZEBCM is run for 1000 years, enough time for the atmosphere and ocean surface to adapt to the distribution of surface type and surface elevation. Thereafter, the ISM receives temperature data from the ZEBCM and is run for 1000 years. Next, the ZEBCM is forced with the new CO 2 and insolation data and the updated ice volume and height data obtained from the ISM. The mutual communication between the ZEBCM and ISM continues at 1000year time intervals. The coupling time step serves to obtain the desired 1000-year output resolution, as the physics in our model are not detailed enough to simulate sub-millennial climate dynamics. Although it is possible to use a shorter coupling time step, it slows down the model and does not significantly change the results.  The temperature data going from the ZEBCM to the ISM consist of mean monthly atmospheric temperatures in • C (reduced to sea level), averaged over the areas where the ice sheets grow (Table 3). The mass balances of the northern hemispheric (NH) ice sheets are forced with these monthly values, corrected for height with a lapse rate of 6.5 • C km −1 . The volumes (V ) of the different ice sheets are transferred to the ZEBCM. The relation for the volume-to-area ratio of a perfectly plastic ice sheet V = C vol × A 5/4 (Bahr et al., 1997) is then used to calculate ice extent (A). The ice-sheetdependent gain factor C vol is calculated from the data of De Boer et al. (2013) ( Table 3). The NH ice sheets are centred at the locations specified in Table 2 and grow circularly. The influence of the position of the ice sheets in the ZEBCM on the climate will be studied in Sect. 3.3. The sum of the southern hemispheric (SH) ice sheet areas is considered to be a single circle, centred at the South Pole (90 • S).
The PD topography in the ZEBCM is taken from Becker et al. (2009). This topography is assumed to be consistent with the ice sheet heights calculated from a steady-state run of the coupled model, forced with PD insolation and preindustrial (PI) CO 2 (280 ppm). For the NH ice sheets, height changes with respect to PI are determined from the ISM. To calculate the height changes for Antarctica, the SH ice sheets are considered to be a single perfectly plastic ice sheet. The total area is used to get the ice sheet radius, and the ice thickness profile from the edge to the centre follows from h(x) = 2τ 0 x ρg (Oerlemans, 2011), where h(x) is the thickness at location x, τ 0 the yield stress (calculated by matching PD ice volume and PD ice area), ρ the density of ice and g the gravitational acceleration. Integration of this formula per latitude band yields the mean zonal thickness.

Present-day validation and reference experiment
First, an equilibrium PI/PD model test was performed by running the coupled model for 100 kyr, starting with no ice, PI insolation and PI CO 2 of 280 ppm as input. After 50 kyr, the input was replaced with PD values (350 ppm CO 2 ), but the ice sheet forcing in the ZEBCM was kept constant at the final PI level, because PD land ice is not in long-term equilibrium with PD climate. Year 50 000 represents the equilibrium climate for PI, and year 100 000 for PD. The resulting equilibrium meridional temperature distribution for both climates is plotted and compared to PD ERA-40 reanalysis data for 1970-2000 (Uppala et al., 2005) in Fig. 2a. The difference between PI and PD temperatures is rather small: 1 K globally averaged, though locally up to 4 K. The model generally captures the observed temperature distribution and is  (Laskar et al., 2004, red), CO 2 (Jouzel et al., 2007, green) input, modelled reference NH temperature (orange), sea level relative to PD (blue), and global temperature (black). globally averaged only slightly too warm: 288.1 K, whereas 287.7 K is observed. At 10-50 • N the model produces moderately higher temperatures than the observations, while at 50-90 • N they are lower. In the Southern Hemisphere, temperatures are locally somewhat overestimated at 52.5 • and 72.5 • S. Local differences are most likely due to missing atmosphere and ocean dynamics. In this respect, the most important deficiency of the model is the assumption of infinite zonal exchange of energy between ocean and land atmosphere. This may lead to underestimation (overestimation) of the seasonal temperature cycle at land (ocean)-dominated lat-itudes. Tests with finite wind strength, however, did not produce better results. The results of our model set-up, where ice volume may vary over time, are in general accordance with PD model tests performed with the old set-up of the ZEBCM forced by observed PD ice (Bintanja, 1997). Overall, the model performs reasonably well in simulating steadystate PI and PD climate.
Thereafter, a reference run of the coupled model was performed. The model was spun up for 386 kyr, with 1072-786 kyr ago insolation data. CO 2 was gradually lowered from 730 to 280 ppm. During the last 100 kyr of the spin up, an artificial interglacial-glacial CO 2 cycle was simulated by a cosine, with CO 2 values ranging from 280 to 180 ppm. After the spin up, the model was run for 786 kyr. CO 2 data of the EPICA record (Jouzel et al., 2007), interpolated to yearly values, were used as input, along with 5 • resolution insolation data of Laskar et al. (2004). The 65 • N insolation and CO 2 input are plotted in Fig. 3, as well as the modelled NH temperature (40-80 • N mean), sea level and global temperature. Figure 2b shows the modelled meridional temperature difference between PI and LGM (19-23 kyr ago mean), compared to the zonally averaged data from a recent LGM reconstruction (Annan and Hargreaves, 2013, AH13). As can be seen, the modelled distribution lies within the uncertainty band of the reconstruction everywhere except at the North Pole. The modelled temperature difference gradient is larger then the reconstructed one, especially in the Northern Hemisphere. The PI-LGM global temperature difference calculated by the model is 4.9 K; this is higher than the 4.0 K found by Annan and Hargreaves (2013) A comparison of the modelled sea level with an earlier modelling study (De Boer et al., 2010) and Red Sea sediment data ( (Rohling et al., 2009, R09) is shown in Fig. 4. In general, the glacial-interglacial cycles are simulated well by the model. There is reasonable agreement between model and data (correlation coefficient of 0.71, RMSE = 31.8 m; Fig. 4a). The last glacial cycle (125 kyr ago to present) is particularly well modelled (correlation coefficient of 0.86, RMSE = 22.5 m; Fig. 4b). Although the general trend is simulated more accurately by De Boer et al. (2010), the strength of the ∼ 40 kyr fluctuations in the sea level data at 90 kyr ago and 60 kyr ago is captured better by our coupled model. Around 330 to 260 kyr ago, modelled glaciation takes place much later than in the reconstructed sea-level record. This mismatch will be discussed further in Sect. 3.6. The model seems to slightly underestimate the sea level drop at 250, 350 and 440 kyr ago. This may at least in part be attributed to an underestimated sensitivity of the Antarctic ice sheets to temperature changes. Recent estimates of the LGM Antarctic ice volume difference show an increase of ∼8-10 m s.l. with respect to PI (e.g. Maris et al., 2014). In our model, a ∼ 4 m s.l. size increase in the West Antarctic Ice Sheet during glacials is almost entirely offset by a decrease in the East Antarctic Ice Sheet due to decreased precipitation. This stability is a consequence of the simplified geometry of the onedimensional ice sheet model and the lack of ice shelf dynamics, as modelling studies using more sophisticated ice sheet/shelf models show larger Antarctic Ice Sheet glacialinterglacial differences (De Boer et al., 2013;Pollard and DeConto, 2009).

Comparison of reconstructed temperatures
A comparison is made between the modelled Greenland temperature anomaly and a 2 kyr smoothed data reconstruction (Johnsen et al., 1995) based on δ 18 O from the GRIP core (GRIP members, 1993) (Fig. 5a). Modelled Greenland temperature is assumed to be the simulated 70 • N temperature at sea level, corrected for the height of the Greenland Ice Sheet centre. The overall agreement between model and data is decent; the model captures the main trend of the data (correlation coefficient of 0.82). The rapid, sub-millennial fluctuations, however, are not present in the model output be-cause of lacking short-timescale forcing variations. Furthermore, the model seems to underestimate the amplitude of the temperature record by approximately 5 K. Some uncertainty resides in the interpretation of the data, as Greenland temperature records inferred from different water isotope data show different glacial-interglacial variability (Kindler et al., 2014;Simonsen et al., 2011). However, most likely this underestimation of temperature variability is predominantly a consequence of the zonal averaging in our model, which distributes temperature anomalies over the entire meridional band. In reality, local temperature differences can be much larger than the zonal average, especially over ice sheet areas, which are prone to large albedo and height changes.
In Fig. 5b, the modelled Antarctic temperature anomaly (60-90 • S surface temperature) is shown, together with the temperature reconstruction based on the deuterium isotope record of the Dome C ice core (Jouzel et al., 2007). The modelled temperature shows rapid fluctuations of larger amplitude than the data. These are associated with shifts of the ocean-overturning midpoint. These shifts have a profound effect on Antarctic climate. By increasing inter-hemispheric heat transfer, they lead to lower Antarctic glacial temperatures that are more in line with observations. However, the shifts have a discrete character as, restricted by the model resolution, the midpoint moves 5 • at once. This leads to overestimation of temperature fluctuations. Nonetheless, a good coherence of model output and data is visible, both in timing (correlation coefficient of 0.73) and in amplitude (RMSE of 2.5 K). During the past four interglacials this agreement is less good; the data shows up to 4 K higher temperatures, while the model only calculates conditions similar to PI. This underestimation was also found by Ganopolski and Calov (2011), who attributed it to not including West Antarctic Ice Sheet changes. Although we do include Antarctic volume changes in our model, the combined Antarctic Ice Sheet shows too little variability to produce interglacial temperatures markedly higher than PI(see also Sect. 3.1). Figure 5c shows the deep-sea temperature anomaly, as calculated by the model and as reconstructed from the Mg / Ca ratio of a sediment core (Elderfield et al., 2012). The second vertical layer of the ocean model is assumed to represent the deep ocean. The sediment core was drilled as part of the Ocean Drilling Program (ODP) at 41 • S (leg 181, site 1123). Two model records are displayed: the mean deep-sea temperature anomaly of the total ocean (red line) and the zonal average at 42.5 • S (blue line). While the model output at 42.5 • S corresponds reasonably well with the data (correlation coefficient of 0.59), the amplitude of the anomaly is smaller than observed. The aforementioned stability of the Antarctic Ice Sheet probably plays a role in this underestimation. Moreover, the zonally averaged ocean is only a crude representation of the real ocean system, which has three main basins, and where the deep-ocean temperature is also determined by the zonally varying circulation of deep water. Also visible in Fig. 5c is that the total ocean anomaly differs from the zonal

Sensitivity analysis
Many assumptions and choices have been made to develop the set-up of the coupled model. Some of these may have a large impact on the model results. To test the sensitivity of the model to the aspects introduced in this study, we performed additional runs with one aspect changed at a time. Many options were tested, but most were found to be of minor importance. Here, we present the tests on the most important aspects: (1) the formulation of ocean overturning, (2) the position of the ice sheets in the ZEBCM, and (3) the radiative forcing strength of non-CO 2 greenhouse gases.
The historical strength of the NH ocean overturning is not well constrained and varies significantly between different models (e.g. Bakker et al., 2013;Weber et al., 2007). In our reference run, the strength is 9-11 Sv in fully glaciated conditions and 14-15 Sv during interglacials. The strength of our SH ocean overturning is ∼ 13 and ∼ 27 Sv in glacials and interglacials, respectively. Both ranges are similar to the results of Ritz et al. (2011), who used a more detailed ocean model. In our model, overturning strength is variable in two ways (see Sect. 2.2): it is dependent on the equator-to-pole waterdensity difference, and on the pole-to-pole water-density difference through the shifting midpoint. A model run is performed in which the overturning circulation is fixed to its PD strength and pattern (Fig. 6, blue lines). At 620 to 500 kyr ago, glaciation is initially stronger, but ceases earlier than in the reference run. Aside from this period, the temperature and sea level reconstructions look very similar to the reference experiment. In a model run where the overturning strength is varied depending on equator-to-pole water density difference but the midpoint is fixed to its PD position (Fig. 6, black lines), the overturning strength is more stable than in the reference run. The 620 to 500 kyr ago period looks similar to the fixed PD-overturning case. This indicates that the delayed major glaciation in the reference run is caused by the shifting midpoint. In the fixed PD-overturning case, NH temperatures fluctuate slightly less. This is reflected in a smaller LGM-PI temperature difference in the NH. Shifting the overturning midpoint only has a very small effect on this difference. For Antarctic temperatures, the opposite holds. Varying ocean strength increases the SH LGM-PI difference only a little, while the midpoint shift more than doubles it. The difference between NH and SH (polar) temperatures in sensitivity to oceanic heat transfer in the ZEBCM was earlier noted by Bintanja and Oerlemans (1996). They used manually set reduced SH overturning strength to obtain a lower Antarctic LGM temperature. In a more physical way, our approach of changing ocean overturning strength and a shifting midpoint, mimicking the effect of relative increase in the Atlantic cross-equatorial flow, also leads to good agreement with the data (see Sect. 3.2). The NH ice sheets grow circularly in the ZEBCM, starting from a fixed centre. In reality, growth is more or less southward, starting at higher latitudes. The present-day situation, with a matured Greenland Ice Sheet, but no Eurasian and North American ice sheets, is accurately modelled in terms of temperature (Fig. 2a). Nonetheless, it is desirable to know the influence of the positioning of the NH ice sheets in the ZEBCM. We conduct a model run with 5 • more northerly, Figure 7. Modelled LGM-PI atmospheric temperature difference, for the reference run (red), a run with the NH ice sheets placed 5 • more to the north (blue), and a run with the NH ice sheets placed 5 • more to the south (black).
as well as a run with 5 • more southerly, NH ice sheet centres. The further the ice sheets are located towards the equator, the larger the effect of increased albedo becomes, since more radiation is reflected. Consequently, temperatures are colder for more southern ice sheet positions. Moreover, the LGM-PI temperature difference minimum shifts by 5 • to 65 • N (Fig. 7). As is shown in Sect. 3.4, this is caused by the surface-height changes. In contrast to the overall colder temperatures, this minimum shift cannot be undone by a different tuning of the ISM. The set-up of the reference run, the Greenland Ice Sheet centred at 65 • N, and the Eurasian and North American ice sheets at 60 • N results in the best correspondence to the LGM-PI difference found by Annan and Hargreaves (2013), while the ice sheets are still located in geographically realistic positions.
Lastly, our assumption of enhancing the radiative forcing of CO 2 by a factor of 1.3, to account for non-CO 2 GHGs, is tested. When we do not enhance CO 2 radiative forcing, temperature deviations are smaller (not shown), as expected. However, the simulated evolutions of the global temperature and sea level are only different from the reference reconstruction from 620 to 500 kyr ago. In this case, there is no glaciation at all during the whole period. In addition, a model run is performed, using the actual enhancement of CO 2 by CH 2 and N 2 O, from records of Loulergue et al. (2008) and Spahni et al. (2005) following Köhler et al. (2010). The radiative forcing in this case is almost the same as in the reference case. Not surprisingly, there is also little difference in the reconstructed NH temperature and sea level. Even in the period 120 to 80 kyr ago, the onset of the last glacial period, where the radiative forcing is considerably stronger than in the reference, the NH ice sheets do not grow more quickly. This indicates that the build-up of ice sheets is limited by precipitationrather than the temperature decrease. As could be expected, the LGM-PI temperature difference is increased by taking non-GHGs into account (not shown). Conversely, using a CO 2 enhancement factor of 1.3 has a similar effect to using real radiative forcing data of the most important non-GHGs.

Climate-ice-sheet feedbacks
To quantify the effect of long-term ice sheet variability on the climate, the fully coupled model runs are compared to (partly) uncoupled runs (Fig. 8). The ISM is still forced with ZEBCM temperatures, but no information (blue line), only surface-height change (orange line), or only ice volume (black line), is transferred back to the ZEBCM. The missing information needed to calculate ice sheet size (albedo) or surface height, or both (see Sect. 2.2), is prescribed as PD values. Figure 8 shows that ice-sheet-climate feedbacks do not alter the evolution of global temperature fluctuations, nor do they delay glaciation. However, the amplitude of the global temperature and sea level cycles is reduced to about half the size. Temperature minima in the NH are reduced even more severely. This becomes most apparent from the LGM-PI difference, which is augmented by the ice-sheet-climate to about 5 times the uncoupled-case difference (Fig. 9a). Icesheet-climate feedbacks are most prominently affecting the NH. This is not surprising, since there ice sheet variability is much larger than on the SH. In Figs. 12 and 13, it is also visible that the ice-albedo feedback has a larger influence than the surface-height-temperature feedback. Uncoupling ice albedo has almost the same effect as completely uncoupling the ice sheets, while climate variability is only slightly reduced by uncoupling surface height changes. Meanwhile, the distinct local minimum in the LGM-PI temperature difference at around 62.5 • N is caused by surfaceheight changes, as it is only present in the runs with varying surface height.
A quantitative measure of the influence of ice sheets on the climate is given by the ratio of global temperature change in the fully coupled case to that in the uncoupled-ice case ( Fig. 9b; red dots). As can be seen, this ratio increases approximately linearly with the total ice amount present on Earth. A linear least-squares fit through all data points with a larger than 0.5 K global temperature decrease in the uncoupled-ice case renders a slope of 0.04. This means that global temperature change is enhanced by 4 % per 10 6 km 3 of ice, which is equivalent to 1 % per 0.63 m sea level. This is mainly established by the ice-albedo feedback (Fig. 9b, black dots), rather than the surface-height-temperature feedback (orange dots). In the coupled-albedo case, the enhancement percentage is 3 % (1 % per 0.84 m s.l.), while in the coupled height it is only 0.7 % per 10 6 km 3 of ice (1 % per 3.6 m s.l.).
The effect of ice sheets on atmospheric temperature is mostly localised to the region where they grow. This is de- Figure 8. Modelled global temperature (top), sea level relative to PD (middle), and NH temperature (bottom) for the fully coupled reference run (red), a run with surface height and albedo in the ZE-BCM prescribed as PD values (ice uncoupled, blue), runs with only surface-height changes (height coupled, orange), and runs with only ice volume (albedo coupled, black) transferred back to the ZEBCM. The missing information on albedo or surface height, or both, is prescribed as PD values. The temperature data come from the ZEBCM, and the sea level data from the ISM. picted by a comparison of NH to global temperature change in Fig. 10a and b. The ratio between the two, known as (northern) polar amplification, is approximately constant in all cases. When ice is uncoupled, a polar amplification factor of 1.43 is deduced from a linear least-squares fit. When ice is coupled, this factor increases by 94 % to 2.78. Coupling only the ice albedo has a larger effect (factor of 2.29, 60 % increase) than coupling just the surface height (factor of 1.83, 28 % increase). The ice-albedo and surface-heighttemperature feedback also strengthen each other. This can be deduced from both the enhancement of global temperature change and polar amplification. In the case of polar amplification, the sum of their individual contributions would lead to an 88 % increase, which is ∼ 7 % smaller than their pooled effect.

Relative influence of CO 2 and insolation
Finally, we tested the relative contributions of CO 2 (Fig. 11) and insolation (Fig. 12) by keeping one of the variables LGM-PI atmospheric temperature difference, for the same runs as in Fig. 8. (b) Ratio of global temperature anomalies from the fully coupled reference run (red), the height-coupled run (orange) and the albedo-coupled run (black) to those from the ice-uncoupled run, as a function of ice volume in the (partially) coupled run. The temperature data come from the ZEBCM, and the sea level data from the ISM. Neither the 280 ppm CO 2 case (Fig. 12, orange line) nor the 125 kyr ago insolation case (Fig. 11, black line) exhibit inception of the Eurasian and North American ice sheets. This means that both variables show a threshold inception value. Likewise, both variables play a role in deglaciation, because at fixed PD insolation and at fixed 180 and 230 ppm CO 2 (Fig. 12, black and blue lines), these ice sheets do not disappear. The importance of transient effects is demonstrated by the low sea level stand of −40 m at PD, when PD insolation is used (Fig. 11, blue line). From the global temperature records of the aforementioned cases with no NH glaciation (Fig. 12, orange line; Fig. 11, black line), it becomes clear that the effect of CO 2 is more global than the effect of insolation. Greenhouse gases affect both hemispheres in a similar manner. In contrast, the insolation signal is dominated by precession, the effects of which are opposite for the two hemispheres. When the NH ice sheets are present, temperature fluctuations are much larger in both the Figure 11. Modelled global temperature (top), sea level relative to PD (middle), and NH temperature (bottom) for the fully coupled reference run (red), a run with insolation fixed at PD value (blue) and a run with insolation fixed at the precession maximum of 125 kyr ago value (black). fixed CO 2 and the fixed-insolation cases because of the icesheet-climate feedbacks. Insolation affects ice volume both indirectly, through temperature in the ZEBCM, and directly, through the mass balance in the ISM (Eq. 3). The indirect effect is more important, as fixing only the ISM input has almost no effect (not shown). The reference sea level is correlated more strongly to CO 2 (coefficient of 0.84) than to insolation (coefficient of 0.35). This indicates that CO 2 levels provide the envelope of sea level. Insolation changes have a damping or amplifying effect, as well as being decisive for the timing of glacial inception. A contrasting conclusion was drawn based on work with the LLN 2-D climate model (Loutre and Berger, 2000). Using this model, glacialinterglacial variability in good coherence with − but slightly smaller than -observations, was simulated with a constant CO 2 level of 210 (Gallée et al., 1992;Berger et al., 1998;Pépin et al., 2001), comparable to our simulation with 230 ppm CO 2 . However, keeping insolation constant at PD level and using the Vostok ice core record as forcing for CO 2 , they hardly found any variability in temperature and sea level (Loutre and Berger, 2000;Pépin et al., 2001). This result is very different from our fixed PD insolation case, where considerable variability is found (Fig. 11, blue line). In addition, Figure 12. Modelled global temperature (top), sea level relative to PD (middle) and NH temperature (bottom) for the fully coupled reference run (red), a run with CO 2 fixed at 180 ppm (full glacial value, blue), a run with CO 2 fixed at 230 ppm (intermediate value, black), and CO 2 fixed at 280 ppm (interglacial (PI) value, orange). keeping CO 2 high at 290 ppm did not prevent glaciation of the NH in the LLN 2-D model (Berger et al., 1998), as is the case in our 280 ppm experiment (Fig. 12, orange line). These differences may be attributed to the LLN 2-D model having a small sensitivity to GHG changes (Gallée et al., 1992), due in part to ignoring non-CO 2 GHGs. In this analysis, insolation and CO 2 are considered as separate forcings, while in reality CO 2 levels are part of the Earth system's response to external forcings (e.g. through the chemical weathering feedback Walker et al., 1981).

Discussion
By using an energy balance climate model coupled to a onedimensional ice sheet model, we reconstructed global mean sea level and a zonally averaged air temperature over the past 800 kyr. The sensitivity tests presented in Sect. 3 show that SH temperature is largely dependent on the parameterisation of the ocean overturning. In contrast, NH temperature and sea level are affected less and are more robust. Furthermore, the influence of non-CO 2 GHGs is shown to be very important for the strength of glaciation. They largely determine the amplitude of sea level fall during the whole reconstructed period. However, the reconstructedtemporalevolution of tem- perature and sea level are only affected in the period 620 to 500 kyr ago. During this period, glaciation relies on sufficiently low NH temperatures, generated by both non-CO 2 forcing and ocean overturning. Multiplying the radiative forcing of CO 2 by the mean amplification over the past 800 kyr (factor of 1.3) according to Köhler et al. (2010) has to a large degree the same impact as explicitly modelling non-CO 2 GHG forcing. Finally, the sensitivity analysis shows only a minor influence of the position of the ice sheets in the ZEBCM with respect to the resulting temperature and sea level reconstruction. The set-up in this study is chosen to ensure that the simulated global sea level is in the best possible accordance with the Red Sea record of Rohling et al. (2009), and the modelled temperature difference between LGM and PD agrees within uncertainty limits with the data reconstruction provided by Annan and Hargreaves (2013). Overall, the modelled sea level is in reasonable agreement with the data, except at around 330 to 260 kyr ago. During this period, NH glaciation is simulated much later than in the data record. This is most likely a consequence of large hysteresis in the ISM induced by the simplified geometry. If Eurasia and North America are deglaciated entirely in the model, glacial inception on the NH only takes place when a drop in CO 2 is accompanied by a large enough drop in NH summer insolation. This is demonstrated by separate hysteresis tests, where CO 2 is decreased in steps of 10 ppm per 30 kyr from 280 to 180 ppm and then increased in the reverse manner. Insolation is kept at constant level of 335 (high NH summer insolation), 326 (intermediate) and 325 (low) kyr ago. In the reference case (REF), it is first decreased in steps of 1 kyr from the high to the low extreme, and then increased back to the initial level, in phase with CO 2 . An index is assigned to the input values, linearly ranging from favourable (180 ppm CO 2 (and low NH summer insolation for REF); index 0) to unfavourable (280 ppm CO 2 (and high NH summer insolation for REF); index 1) for NH glaciation. From comparison of the reference case to the fixed-insolation cases, it is apparent that hysteresis is reduced when NH summer insolation and CO 2 are in phase (Fig. 13). A more realistic topography than used here would ensure that some high parts remain glaciated (or at least glaciate much earlier) and cause a positive feedback when the CO 2 level decreases, thereby reducing the hysteresis. Due to high CO 2 and summer insolation in the NH, the Eurasian and North American ice sheets completely disappear at 334 kyr ago. During the period 334 to 280 kyr ago, NH summer insolation does not drop to low enough levels to realise inception. This only occurs at 280 kyr ago, when the CO 2 level is very low. The NH ice sheets then grow quickly because of the low CO 2 levels and the ice-sheet-climate feedbacks. Although it cannot be excluded that the importance of radiation is overestimated in our model, it only leads to a model-data mismatch from 334 to 280 kyr ago. Removing this period from our data analysis does not significantly alter the resulting strength of the icesheet-climate feedbacks (4 % per 10 6 km 3 of ice) and polar amplification factor (2.78), which we consider our most important findings.
In our model, the influence of ice on the climate is about equally large, as was obtained by Gallée et al. (1992) with the LLN 2-D model. This is judged by considering the LGM-PI temperature difference, averaged over the whole Northern Hemisphere (0-90 • N). By including ice-sheet-climate interactions, this difference is enhanced by a factor of 4.3 in their model, and by an almost equal factor of 4.5 in ours. However, in the LLN 2-D model the influence is less localised to the regions where the ice sheets grow. This may be caused by a lower albedo of ice in our model (0.80 as opposed to their 0.85), but this seems unlikely because land albedo is also lower in our model, offsetting the land-ice albedo difference. More likely, differences in atmospheric and oceanic transport play an important role. A thorough analysis, however, is hampered because the LLN 2-D model does not include a Southern Hemisphere. Furthermore, Gallée et al. (1992) compared runs with only varying insolation, while we compare runs with CO 2 changes also included.
As mentioned before, land ice has an influence on climate through albedo and surface height. However, the complete temperature response to land ice is also determined by changing sea ice and snow cover, which are fast feedbacks. To quantify this influence we did two further runs of the model with sea ice and snow cover kept constant: one with coupled ice sheets, and one with land-ice cover and surface height fixed at PD levels (Table 4). Compared to the completely fixed run, coupling land ice leads to an increase in global temperature anomalies by 0.7 % per 10 6 km 3 of ice (1 % per Figure 13. Modelled global sea level in hysteresis tests, where CO 2 is decreased in steps of 10 ppm per 30 kyr from 280 to 180 ppm and then increased in the reverse manner. Insolation is kept at constant levels of 335 (high NH summer insolation, blue line), 326 (intermediate, black line) and 325 (low, orange line) kyr ago. In the reference case (REF, red line), it is first decreased in steps of 1 kyr from the high to the low extreme, and then increased back to the initial level, in phase with CO 2 . An index is assigned to the input values, linearly ranging from favourable (180 ppm CO 2 (and low NH summer insolation for REF); index 0) to unfavourable (280 ppm CO 2 (and high NH summer insolation for REF); index 1) for NH glaciation.
3.6 m s.l.). Coupling only sea ice and snow (the land-ice fixed case) increases global temperature anomalies by 2.2 % per 10 6 km 3 of ice (1 % per 1.2 m s.l.). Land-ice cover is taken from the ISM output, which is not used to force the ZEBCM in this case. Therefore, the effect of including only fast ice feedbacks is more than three times as important as the effect of including only land-ice changes. The true importance of the land ice feedback only appears when we compare the fully coupled results to those of the completely fixed run. The temperature anomalies are then increased by 5.2 % per 10 6 km 3 of ice (1 % per 0.48 m s.l.), implying that the temperature response is more than twice as large then for the case where only the sea-ice and snow variation are taken into account. This result is much larger than the product of the separate effects, which is caused by positive reinforcement of the effect of land ice, as well as sea ice and snow, on atmospheric temperature.
Our modelled polar amplification of the LGM-PI surface air temperature difference is consistent with a recent LGM data reconstruction (Annan and Hargreaves, 2013) within the uncertainty range provided by that study, and the results of seven models presented in the 2013 IPCC report (IPCC, 2013). Using a GCM, Singarayer and Valdes (2010) performed two tests on polar amplification over the past 120 kyr. In one test, they varied insolation and CO 2 but kept surface elevation and ice extent at PD levels. In the other test, insolation and CO 2 were still varied but they used the ICE 5G ice sheet reconstruction of Peltier (2004) to force their model. Both in the coupled-and uncoupled-ice case, the ratio of Greenland temperature change to global temperature change corresponds to our results. In the uncoupled-ice case, our meridional LGM-PI temperature difference distribution, normalised by global temperature change, is also in broad agreement with the results of an ensemble of GCMs (Holland and Bitz, 2003). Although they based their distribution on the difference induced by a doubling of CO 2 with respect to PD (keeping ice sheets constant at PD size), the results are comparable because of the normalisation. Our model shows an ∼ 2.2-fold global warming temperature difference at the North Pole; this is within the 1.5-to 4.5-fold amplification range in Holland and Bitz (2003). The agreement with these results, obtained with more sophisticated models, provides confidence in our resulting polar amplification factors.

Summary and conclusion
We have coupled a zonally averaged energy balance climate model (ZEBCM) to a one-dimensional ice sheet model (ISM). The coupled model was used to obtain an 800 000year record of meridional atmospheric and ocean temperature, as well as global eustatic sea level. The sea level record was tuned to match a Red Sea sediment data reconstruction (Rohling et al., 2009) as closely as possible, with a correlation coefficient of 0.71 and an RMSE of 31.8 m. The simulated Antarctic temperature, Greenland temperature and deep-sea temperature are in reasonable agreement with proxy-data reconstructions. Our model shows that NH polar temperatures are mainly determined by the waxing and waning of the Eurasian and North American ice sheets, while the Antarctic Ice Sheet is more stable and SH climate is heavily influenced by the strength of ocean overturning. We realise that all our results are model-specific and only large-scale features are significant because of limitations imposed by the zonal averaging and rather coarse resolution (5 • ) of the ZE-BCM as well as the lack of detailed physics. Furthermore, the ISM has a highly simplified geometry and misses ice shelf dynamics. However, this study significantly improves on earlier stand-alone ISM results, forced by an inverse δ 18 O routine (De Boer et al., 2010). By including the ZEBCM, a more comprehensive view of global climate is acquired, as it is no longer an anomaly with respect to the PD climate. Now, we calculate not only the meridional temperature distribution but also more realistic deep-sea temperatures. Deepsea temperature fluctuations are tightly coupled to fluctuations of atmospheric temperatures in the same region, reduced approximately by a factor of 5. Switching CO 2 and L. B. Stap et al.: Interaction of ice sheets and climate during the past 800 000 years insolation forcing on and off enabled us to study the relative influence of these two variables. Atmospheric temperature on the NH is controlled by a complex interaction of CO 2 and insolation. Both variables, as well as transient effects, govern the evolution of ice sheets on the Northern Hemisphere. By comparing the fully coupled results to results from the ZEBCM forced with PD albedo and/or surface height, we have quantified the importance of the ice-albedo feedback and the surface-height-temperature feedback. The combined long-term ice-sheet-climate feedback amplifies global temperature anomalies from PI with 4 % per 10 6 km 3 of ice (1 % per 0.63 m s.l.). In full glacial conditions (∼ 40 × 10 6 km 3 of ice), this leads to an increase in the temperature change by a factor of 2.6. Furthermore, it almost doubles polar amplification, from 1.43 to 2.78. The dominant feedback is the icealbedo feedback, while the surface-height-temperature feedback has only minor influence. Compared to the results of a model run where also the fast ice feedbacks (sea ice and snow) are fixed, the amplification factor of global temperature anomalies increases to 5.2 % per 10 6 km 3 of ice (1 % per 0.48 m s.l.) in the fully coupled run. This is more than twice as large as the effect of only including sea-ice and snow cover variation (2.2 % per 10 6 km 3 of ice or 1 % per 1.2 m s.l.).
In future research, the use of the coupled model will be extended into the Pliocene (5 Myr ago to PD), and the early to middle Miocene (23 to 13 Myr ago). During the former period, the Antarctic Ice Sheet was more variable in size, while the latter period is characterised by glaciations and deglaciations of the Northern Hemisphere. A marine benthic δ 18 Oisotope model will be incorporated in order to allow for an independent comparison with palaeo-proxy-records (Lisiecki and Raymo, 2005;Zachos et al., 2008;Cramer et al., 2009).