A coupled climate model simulation of Marine Isotope Stage 3 stadial climate

. We present a coupled global climate model (CGCM) simulation, integrated for 1500 yr to quasi-equilibrium, of a stadial (cold period) within Marine Isotope Stage 3 (MIS 3). The simulated Greenland stadial 12 (GS12; ∼ 44 ka BP) annual global mean surface temperature ( T s ) is 5.5 ◦ C lower than in the simulated recent past (RP) climate and 1.3 ◦ C higher than in the simulated Last Glacial Maximum (LGM; 21 ka BP) climate. The simulated GS12 is eval-uated against proxy data and previous modelling studies of MIS3 stadial climate. We show that the simulated MIS 3 climate, and hence conclusions drawn regarding the dynamics of this climate, is highly model-dependent. The main ﬁndings are: (i) Proxy sea surface temperatures (SSTs) are higher than simulated SSTs in the central North Atlantic, in contrast to earlier simulations of MIS 3 stadial climate in which proxy SSTs were found to be lower than simulated SST. (ii) The Atlantic Meridional Overturning Circulation (AMOC) slows down by 50 % in the GS12 climate as compared to the RP climate. This slowdown is attained without freshwater forcing in the North Atlantic region, a method used in other studies to force an AMOC shutdown. (iii) El-Ni˜no-Southern Oscillation (ENSO) teleconnections in mean sea level pressure (MSLP) are signiﬁcantly modiﬁed by GS12 and LGM forcing and boundary conditions. (iv) Both the mean state and variability of the simulated GS12 is dependent on the equilibration. The annual global mean T s only changes by 0.10 ◦ C from model years 500–599 to the last century of the simulation, indicating that the climate system may be close to equilibrium already after 500 yr of integration. However, signiﬁcant regional differences between the last century of the simulation and model years 500–599 exist. Further, the difference between simulated and proxy SST is reduced from model years 500–599 to the last century of the simulation. The results of the ENSO variability analysis is also shown to depend on the equilibration.


Introduction
Greenland ice core oxygen isotope records (Dansgaard et al., 1993;North Greenland Ice Core Project Members, 2004) indicate that climate varied between cold (stadial) and mild (interstadial) conditions during MIS 3, approximately 60 ka BP to 30 ka BP. Greenland ice core records show abrupt transitions from cold, stadial to mild, interstadial conditions (known as Dansgaard-Oeschger (DO) events) followed by a slow cooling lasting several centuries and a more rapid drop to cold stadial conditions (Huber et al., 2006). The mechanisms responsible for these abrupt climate shifts are not clear. The shifts could be explained by, e.g. changes in ocean thermohaline circulation, sea ice feedbacks, and/or tropical processes (Clement and Peterson, 2008), although each of these mechanisms has explanatory strengths and weaknesses.
In this article, results of a CGCM simulation of Greenland Stadial 12 within MIS 3, ∼44 ka BP (Svensson et al., 2008), are presented. The purpose of the GS12 simulation is to test the ability of a state-of-the-art climate model to reproduce the climate during an MIS 3 stadial in agreement with proxy data. Until recently (Merkel et al., 2010), no CGCM simulations had been performed for MIS 3 conditions. The present study, however, differs from Merkel et al. (2010) in two major aspects: (1) higher grid resolution of the different components of the CGCM in the present study; and (2) the simulation was integrated for a total of 1538 yr allowing for equilibration of the simulated climate in the present study, as compared to 260 yr of integration in Merkel et al. (2010).
Previous modelling studies of MIS 3 had been performed with an atmospheric general circulation model forced by reconstructed sea surface temperatures (Barron and Pollard, 2002;Pollard and Barron, 2003) and with coupled models of intermediate complexity (Jin et al., 2007;Rial and Yang, 2007;van Meerbeeck et al., 2009). Barron and Pollard (2002) and Pollard and Barron (2003) find that MIS 3 variations in orbital forcing, Scandinavian ice sheet size, and atmospheric CO 2 concentrations could not explain the differences between a cold and a milder MIS 3 state found in proxy data. A limitation of their study is, however, the use of prescribed SSTs which precludes changes in the state of the ocean in response to changes in forcing and boundary conditions. van Meerbeeck et al. (2009) use an Earth system model of intermediate complexity to simulate MIS 3 stadial and interstadial climate. Their study indicates that climate variability during MIS 3 could not be explained by MIS 3 variations in orbital forcing, atmospheric dust concentration and atmospheric CO 2 , CH 4 and N 2 O concentrations.
An advantage of intermediate complexity climate models is that the experiments can be integrated for several thousands of model years using a reasonable amount of computing-time. The CGCM simulation presented here is integrated for ∼1500 model years, which is quite long as compared to other CGCM studies of past climate, but too short to determine if the simulated climate reaches an equilibrium or if there is internal variability on (multi-)millennial time-scales. Trends in atmospheric temperature and oceanic temperature and salinity decrease towards the end of the simulation presented here and we therefore chose to use the term quasi-equilibrium to designate the simulated climate towards the end of our simulation. Authors van Meerbeeck et al. (2009) find that a better agreement with proxy data for stadial MIS 3 climate is obtained with a strong, additional freshwater flux in the midlatitudes of the North Atlantic Ocean, which slows down the AMOC. They conclude that freshwater forcing is necessary to induce the shift from warm to cold interstadials and suggest a change in perspectives, i.e. that stadial climate represents a perturbed climate state rather than a typical nearequilibrium MIS 3 climate. We investigate if this hypothesis is in agreement with the results of the CGCM simulation presented here.
Data from the earlier part of the GS12 simulation presented here was used in dynamical downscaling with a regional climate model for Europe (Kjellström et al., 2010).
Here we present results from the equilibration of the CGCM simulation and the quasi-equilibrium that was reached towards the end of the simulation. The model, forcing and boundary conditions, and the proxy SST data are described in Sect. 2. In Sect. 3, the equilibration of the simulated climate is described in terms of global annual mean quantities and in terms of regional climate. Specifically, the simulated climate from the earlier part of the GS12 simulation is compared to the quasi-equilibrium that was reached towards the end of the simulation. The agreement between simulated and proxy SST and surface temperature is described in Sect. 4. Merkel et al. (2010) investigate El-Niño-Southern Oscillation teleconnections in simulations of MIS 3 stadial and interstadial climate and in a simulation of the LGM climate. They find that glacial boundary conditions induce major modifications to ENSO teleconnections and that the "blueprint" of modern ENSO teleconnections should only be applied with caution to glacial climate periods. Merkel et al. (2010) base their analysis on relatively short simulations, since their MIS 3 simulations were integrated for 260-360 model years starting from their LGM climate simulation. Based on the results presented in the present study, it is reasonable to believe that the simulated climate analysed by Merkel et al. (2010) has not fully reached an equilibrium. We hypothesise that not only the mean state but also the variability is influenced by the equilibration. To test this hypothesis we determine ENSO teleconnections in MSLP for the simulated GS12 climate, and for comparison from simulations of LGM and RP climate with the same model (Sect. 5). Finally, we discuss how the present findings compare to other modelling studies of MIS 3 climate, with focus on the most recent studies by van Meerbeeck et al. (2009), Merkel et al. (2010 and Kjellström et al. (2010) (Sect. 6).

Model and experiments
The National Centre for Atmospheric Research (NCAR) Community Climate System Model version 3 (CCSM3) is a global, coupled ocean -atmosphere -sea ice -land surface climate model (Collins et al., 2006). For the simulations described here, the atmospheric and land components share a horizontal resolution of T42 (equivalent grid spacing of approximately 2.8 • in latitude and longitude). The ocean and sea ice components share a horizontal grid of 320 • 384 points. There are 26 levels in the atmosphere and 40 levels extending to 5.5-km depth in the ocean.
The GS12 climate simulation is compared to an RP (AD 1990) and an LGM simulation. All three simulations were run at the same horizontal and vertical resolution. The RP simulation is the standard control simulation performed at NCAR (DeWeaver and Bitz, 2006), which was run for a total of 1000 model years. Model years 700-999 are used here to represent RP climate.
The LGM simulation is a 1462 yr continuation of the LGM simulation performed at NCAR (Otto-Bliesner et al., 2006). This simulation follows the protocols established by PMIP2. The last 300 yr are used here to represent LGM climate.   This LGM climate has a global surface cooling of 1.1 • C as compared to the LGM climate described by Otto-Bliesner et al. (2006). This difference is discussed by Brandefelt and Otto-Bliesner (2009) who describe the equilibration of the LGM simulation from the quasi-steady state analysed by Otto-Bliesner et al. (2006) to a new quasi-equilibrium. The GS12 simulation was initiated after 900 yr of the LGM simulation, i.e. when the simulated LGM climate had reached a new quasi-equilibrium. The forcing and boundary conditions were changed to represent GS12 conditions and the GS12 simulation was run for a total of 1538 model years. The last 300 yr are used here to represent GS12 climate.

Forcing and boundary conditions
The CGCM determines the coupled evolution of the atmospheric circulation, oceanic circulation, sea ice extent and thickness and the upper portion of the land surface. Due to computational limitations, not all processes that may be of importance are included in the computations. Some components of the system are set constant, under the assumption that these vary on longer time scales than the time scale of the Table 1. Atmospheric GHG concentrations and orbital forcing for the GS12 simulation as compared to the LGM and RP simulations (Spahni et al., 2005;Lüthi et al., 2008;Berger, 1978 dynamical components of the simulation. The constant forcing and boundary conditions in the present study are briefly described here; for a detailed description see Kjellström et al. (2010). Forcing and boundary conditions were chosen to be representative of the year 44 ka BP for the GS12 simulation.

Radiative forcings
The atmospheric GHG concentrations used in the GS12 simulation were chosen according to ice core measurements (Spahni et al., 2005;Lüthi et al., 2008). Atmospheric GHG concentrations were much lower during GS12 than today, but somewhat higher than during the LGM (Table 1). Atmospheric ozone and aerosol concentrations were set to their pre-industrial concentrations in the GS12 simulation as well as in the LGM simulation, following the PMIP2 protocols for the LGM. Aerosol concentrations are likely underestimated since ice core records indicate higher concentrations of mineral dust at high latitudes for large parts of the glaciated time periods (Lambert et al., 2008). The solar constant was set to its present-day value, 1365 W m −2 , in the GS12 simulation as well as in the LGM and RP simulations. The insolation was calculated following Berger (1978). The amount of insolation in Northern Hemisphere mid and high latitudes is mostly controlled by the obliquity (the tilt of the Earth's axis of rotation) and precession (the precession of the Earth's axis of rotation) signals. A major difference in the forcing during MIS 3 as compared to the LGM is the seasonal and latitudinal distribution of the incoming solar radiation. The largest absolute differences between GS12 and the LGM are found in the high latitudes. The Northern Hemisphere summer insolation at 65 • N was 17 W m −2 higher (10 W m −2 lower) at 44 ka BP (21 ka BP) than today (see Fig. 3 in Kjellström et al., 2010). Another difference between GS12 and the LGM is that GS12 minus RP insolation differences are most prominent in the Northern Hemisphere, whereas LGM minus RP insolation differences have similar magnitudes in both hemispheres. The latitudinal gradient in the GS12 minus RP insolation difference is due to the combined effect of obliquity and precession differences. Obliquity was higher during GS12, with large differences between winter and summer www.clim-past.net/7/649/2011/ insolation as a result. Precession influences the direction of the Earth's axis of rotation with respect to the Sun such that the NH (SH) has a stronger seasonal cycle in insolation for the GS12 (RP).

Ice sheets, topography, and bathymetry
Ice sheets were more extensive during the LGM than during MIS 3, which can be inferred from reconstructions of sea level (Lambeck, 2004;Siddall et al., 2008). These reconstructions give a lowering of the sea level of 120 m for the LGM and a lowering of between 45 and 75 m for MIS 3. The extent and topography of the North American, Fennoscandian, and Antarctic ice sheets are not well constrained for the MIS 3 period .
The ice sheet configuration used for the GS12 simulation (see Fig. 2 in Kjellström et al., 2010) is based on Wohlfarth and Näslund (2010). The resulting ice sheets are taken as a combination of a North American MIS 3 ice sheet from a simulation with CLIMBER-2 (Calov et al., 2005), a simulated Fennoscandian MIS 3 ice sheet from SKB (2006) and an Antarctic ice sheet of 14 ka BP from the ICE-5G data (Peltier, 2004). The latter period was taken as a proxy for the Southern Hemisphere ice sheet at 44 ka BP. The CLIMBER-2 model simulates an ice sheet in Alaska for the period around 44 ka BP, but proxy data indicate that no ice sheet was present in this region (see description in Näslund et al., 2008). This part of the CLIMBER-2 ice sheet was therefore removed.
The GS12 simulation presented here is initialised from the LGM simulation keeping sea level and bathymetry unchanged. The impact of this difference (−120 m instead of between −45 m and −75 m as indicated by the reconstructions of sea level) on the position of the coast line is, however, modest at the coarse scale of the global model. The main difference is that the Barents and Kara Seas are mostly land with a sea-level lowering of 120 m, but are below sea level with a lowering of between 45 and 75 m. Although global sea level was higher during MIS 3 than during the LGM, it is not clear what conditions prevailed in this region during GS12. Geological evidence indicate that the Barents and Kara Seas ice sheet was extensive in the period 60-50 ka BP (Svendsen et al., 2004). To determine the influence of the absence of a Barents and Kara Sea on ocean circulation, and specifically on the deep water formation in the Nordic Seas, sensitivity tests are needed. Such tests are, however, beyond the scope of this article. The simulation presented here uses the same setup with a 120 m lowering of the sea level in both LGM and MIS3 climate simulations as van Meerbeeck et al. (2009) and Merkel et al. (2010).

Proxy data
The model results presented here are compared to SST and atmospheric temperature proxy data of GS12 climate. Qualitative, proxy-based information for Dansgaard-Oeschger cycles is also registered in land records (e.g. Wohlfarth et al., 2008;Wohlfarth, 2010), but the chronological control is often poor and quantitative data is lacking why a detailed comparison is limited.

Oceanic temperature and sea ice extent
The simulated GS12 climate is compared to proxy SST data compiled in an extended version of the Voelker and workshop participants (2002) database. The data is listed along with references and errors given by the authors in the Appendix. The 2 σ uncertainty for the proxy SSTs ranges from 0.4 to 8 • C, with a median value of 1.7 • C. This uncertainty range is determined as two times the error given by the authors for each individual proxy SST value plus the standard deviation in those cases when a mean SST value is calculated from the multiple values available for the GS12 interval. The error given by the authors is generally the  (Johnsen et al., 2001;EPICA Community Members, 2006;Jouzel et al., 2007;Kawamura et al., 2007) and from δ 15 N from a Greenland ice core (Huber et al., 2006). The simulated annual mean T 2 m at the closest grid point is given in parentheses for each location. The EDC and Dome F data are given as a difference between GS12 and the last 1000 yr.

Location
Year "analytical error", for example associated with the transfer function method used to transfer primary data to SST. For some past microfossil faunal compositions modern analogs can be scarce or do not exist (non-analog problem), resulting in larger standard deviations and potentially introducing further uncertainty in the proxy SST data. Moreover, the error given by the authors does not include the uncertainty associated with dating and possible migration of a species within the upper water column or seasonal shifts in the months of highest flux. We compensate for all these potential uncertainties by using two times the error estimated for each site. The uncertainty range for the proxy SSTs is comparable to the errors given for LGM SSTs of the gridded multiproxy approach for the reconstruction of the glacial ocean surface (MARGO Project Members, 2009). The error estimates for the MARGO SST data are in the range of 0.8-9.9 • C, with a median value of 1.95 • C. These error estimates represent 1 standard deviation.
To assess the agreement between simulated and proxy SSTs we define δ: (1) at sites where proxy SSTs are available for GS12. The two data types are considered to agree at sites where the simulated SST falls within the uncertainty range of the proxy SST ( SST proxy ), i.e. when δ ≤ SST proxy .
The proxy SSTs are also used to infer sea ice extent during GS12 winter and summer following Sarnthein et al. (2003), who suggested that North Atlantic sites with an SST of more than 0.4 • C (JFM) and 2.5 • C (JAS) would be sea ice free (<50 % SIC).
Proxy temperature representing GS12 conditions is also available for ∼150 m depth at a site located off the coast of Ecuador at 0 • N, 86.6 • W (Pena et al., 2008), and for 2656 m depth at a site located off the coast of Portugal at 37.6 • N, 10.1 • W (Skinner and Elderfield, 2007).

Atmospheric temperature
The simulated temperature is also compared to reconstructed temperature from ice core δ 18 O records from Greenland (Johnsen et al., 2001) and Antarctica (EPICA Community Members, 2006;Jouzel et al., 2007;Kawamura et al., 2007) and a δ 15 N record from Greenland (Huber et al., 2006) representing GS12. The reconstructed temperatures and the location of the ice cores are given in Table 2. We use data from the Greenland Ice Core Project (GRIP; Johnsen et al., 2001) and from the North Greenland Ice Core Project (NGRIP; Huber et al., 2006). For Antarctica we use data from the European Project for Ice Coring in Antarctica (EPICA) from the interior of Dronning Maud Land (EDML; EPICA Community Members, 2006), EPICA Dome C (EDC; Jouzel et al., 2007), and Dome Fuji (Dome F; Kawamura et al., 2007). The available data represents the absolute temperature (GRIP, NGRIP, EDML) and the difference between GS12 and the last 1000 yr (EDC, Dome F). Huber et al. (2006) determine the uncertainty in the NGRIP reconstructed temperature from δ 15 N from ice core records to ±3 • C, representing two standard deviations. Uncertainty estimates are not available for the data from the other ice cores.

Equilibration
As mentioned earlier, the GS12 simulation was initiated from the new LGM quasi-equilibrium described by Brandefelt and Otto-Bliesner (2009). The simulation was initiated with the appropriate increases in atmospheric GHG concentrations, orbital forcing, and decrease in ice sheet topography, but with surface albedo appropriate for LGM ice sheets. The change in the surface albedo (to account for the decrease in the extent of ice sheets) was introduced after 55 yr of integration. This procedure influences the evolution during the first century, but does not have an effect on the long-term evolution of the simulated climate. The stepwise change in the boundary and forcing conditions is visible in the timeseries displayed in Fig. 3. Simulated SST (grey shading) and the difference δ between simulated and proxy SST (coloured circles/squares) for annual mean (top panel), January-March (middle panel), and July-September (bottom panel). The circle/square is surrounded by a black line at sites where δ ≤ (the uncertainty in the proxy SST). Unit • C. Sea ice covered (sea ice free) proxy sites are plotted with squares (circles). The simulated sea ice edge, defined at 50 % areal sea ice cover, is also displayed for GS12 (red line) and RP (green line). The coastline is defined at 20 % of the land fraction. Fig. 1. The initial change in radiative forcing and topography has only a marginal short-term effect on the annual global mean T s , whereas the change in surface albedo gives an increase of the annual global mean T s of approximately 0.5 • C within a few years.
After the change in surface albedo, the annual global mean T s slowly adjusts to the imposed changes in forcing and boundary conditions. From this slow adjustment one may be tempted to conclude that the simulated climate has reached a new equilibrium after 400-600 yr of integration. Although the annual global mean T s only changes by 0.10 • C from model years 500-599 to the last century of the simulation (model years 1440-1539), statistically significant regional differences (at the 95 % level) exist in boreal winter (January-March; JFM). The JFM temperature at 2 m height (T 2 m ) is more than 2 • C higher in the northern North Atlantic region during the last century of the simulation than in model years 500-599 with a maximum of 8 • C located in Southeastern Greenland (Fig. 2). Further, the precipitation rate is increased by more than 20 % in this region with a maximum of 65 % located in Southeastern Greenland (not shown). These differences are associated with a decrease of 20-30 % in sea ice concentration (SIC; Fig. 2).
The cause of these differences is found in the slow adjustment of the ocean, which can be seen in the time series of surface and abyssal salinity (S surf and S abys ) and potential temperature (T surf and T abys ) in Fig. 1. These results are in agreement with Brandefelt and Otto-Bliesner (2009) who show that the evolution of the abyssal ocean temperature and salinity need to be taken into account when determining if a climate simulation has reached a new equilibrium.
Oceanic northward heat transport at 50 • N and AMOC strength increases for the first 1000 yr of the integration leading to a decreased sea ice cover (Fig. 1). After the first 1000 yr changes in T 2 m and SIC are small (Fig. 2). S surf (and S abys ) increases (decreases) by 0.01 (−0.006) psu per century during the last 300 yr of the simulation, but the major adjustment of the salinity field to the change from LGM to GS12 boundary and forcing conditions occurs during the first ∼1200 model years. These results indicate that the simulated climate is approaching a new equilibrium, and we therefore let the last 300 yr of the simulation represent the simulated GS12 climate. To indicate that equilibration may not be fully completed, we use the term quasi-equilibrium. To confirm or disprove if the climate is approaching an equilibrium, the simulation would have to be extended for a millennium or even longer (which would take one calendar year on the computing facility used for these simulations). The difference is displayed based on model years 500-599 (circles) and for the last 300 yr of the GS12 simulation (asterisks). For clarity, the same colour is used for the two differences corresponding to the same proxy site.
In general, the length of the seasons varies with time due to variations in the Earth's orbit. Following Pollard and Reusch (2002), we find that the length of the seasons is close to the present for 44 ka BP, with differences of 0-1 d in the start and length of the seasons. Therefore, present day seasons are used here.
Also displayed in Fig. 3 is the difference between simulated SST and proxy SST, δ (Eq. 1) at sites where proxy SSTs are available for GS12 (see Appendix). Sites where the simulated SST falls within the uncertainty range of the proxy SST (Eq. 2) are plotted with a surrounding black line. Simulated and proxy SSTs agree for 38-45 % of the sites, depending on season. The largest differences (>4 • C) between simulated and proxy SST are found at some coastal sites, which may partly be explained by strong coastal gradients that are not captured by the coarse resolution of the simulation. Further, large differences (>4 • C) are found in the central North Atlantic in winter, indicating too cold surface waters in the simulation. There is a general pattern of positive differences in the tropics and negative differences in the extra-tropics, indicating that the simulated equator-to-pole SST gradient is too large as compared to that inferred from proxy SST. It should, however, be noted that this conclusion is based on only a few proxy SSTs. Barron and Pollard (2002) and Pollard and Barron (2003) performed simulations of MIS 3 stadial and interstadial climate with an atmospheric circulation model forced by estimated SSTs. They find that the simulated air temperatures appear higher than many proxy indicators suggest, and conclude that errors in the specified North Atlantic SSTs is the most likely cause for the discrepancy. The authors sug-gest that this discrepancy may be solved by using a coupled ocean-atmosphere climate model. Authors van Meerbeeck et al. (2009), however, find that this is not the case for their simulation of stadial MIS 3 climate with an Earth system model of intermediate complexity. North Atlantic SSTs are higher than proxy SSTs in their coupled simulation. This is in contrast to our findings of colder simulated North Atlantic SSTs as compared to proxy SSTs.
The difference of our results as compared to van Meerbeeck et al. (2009) is associated with a difference in the AMOC response to the imposed forcing and boundary conditions. In the present simulation the strength of the AMOC is reduced by 50 % as compared to RP conditions, whereas the AMOC was relatively unchanged in van Meerbeeck et al. (2009). Merkel et al. (2010), who perform simulations of MIS 3 stadial climate with the same CGCM as used here, although at a lower resolution, simulate a reduction of the AMOC strength by 40 % in response to MIS 3 stadial forcing and boundary condition as compared to pre-industrial conditions response.
The difference between simulated and proxy SSTs (δ) is reduced from model years 500-599 towards the end of the simulation. To illustrate this reduction, δ is displayed in Fig. 4 for the last 300 yr and model years 500-599 of the GS12 simulation. δ is on average decreased by 0.3 • C, 0.6 • C, and 0.6 • C for annual, NH winter, and NH summer, respectively. The largest decrease in δ from model years 500-599 to the last 300 yr of the simulation is 0.9 • C, 1.3 • C, and 2.1 • C for annual, NH winter, and NH summer, respectively. Simulated SST is lower than proxy SST for most sites in boreal winter and summer, whereas for the annual mean the distribution is more symmetric around zero. The absolute value of δ is decreased from years 500-599 to the last 300 yr of the simulation primarily for sites in which the simulated SST is colder than the proxy SST. At sites where the simulated SST is higher than the proxy SST the difference between the two periods is generally small, though at a few sites δ is increased.
Proxy data of oceanic temperatures below the surface for GS12 are available for two sites (see Sect. 2.3). The simulated temperature at 150 m depth off the coast of Ecuador is 13.1 • C as compared to the proxy temperature 12.8 ± 1.3 • C. The simulated potential temperature at 2656 m depth off the coast of Portugal is −1.9 • C (which gives an actual temperature of −1.6 • C) as compared to the proxy temperature 0.55 ± 0.8 • C. The comparison suggests that the abyssal ocean temperature could be too low in the simulation, but it should be noted that this is based on one single proxy data site. The coarse resolution of the ocean component of the CGCM can partly explain this difference.

Sea ice extent
We use the proxy SSTs to infer which sites were sea ice covered during GS12 winter and summer as described in Sect. 2.3. The sites inferred to be sea ice covered (>50 % SIC) are indicated in Fig. 3 (squares for sea ice covered sites). All sites are ice-free in JAS and only three sites in the northern North Atlantic are sea ice covered in JFM. In JAS the Nordic Seas are ice free (50 % SIC) in the GS12 simulation, whereas the simulated sea ice edge (50 % SIC) is located south of Iceland in JFM (Fig. 3). For most sites proxy SST based sea ice extent is in agreement with the simulated sea ice cover. Differences are found for a few sites. One of the three sites that are sea ice covered in JFM is located just south of the simulated sea ice edge (50 % SIC) and proxy SSTs imply that the Denmark Strait was ice free in JAS during GS12, whereas this area is sea ice covered in the GS12 simulation.
In JFM, a tongue of sea ice stretches east from New Foundland into the central North Atlantic. The proxy SST sites do not cover the region where the simulated sea ice is located, but there are a few sites east and south of this region with substantially lower simulated than proxy SSTs (simulated SST is 4-6 • C lower than the proxy SSTs), indicating that sea ice may be overestimated in this region.
The NH simulated winter ice extent in the RP simulation agrees reasonably well with observations (Holland et al., 2006). The main NH biases occur in winter with too large ice extent in the Labrador Sea and the Sea of Okhotsk, and too small ice extent in the Barents Sea. This bias towards too extensive sea ice is one candidate for resolving the difference found here between simulated and proxy inferred sea ice extent.

Atmospheric temperature
The difference between GS12 T 2 m and RP climate is most pronounced over NH ice sheets, with a maximum of −42 • C over the northern flank of the Fennoscandian ice sheet in boreal winter (Fig. 5). This difference is associated with the elevation of the ice sheet in this region together with an amplification of the topographic wave forced by the Rocky Mountains and the Laurentide Ice Sheet as compared to the RP (Fig. 5). The spatial pattern of the difference between GS12 and RP T 2 m is similar to the PMIP2 simulations of LGM climate (Braconnot et al., 2007), although with a smaller magnitude. The largest differences in T 2 m between our GS12 and LGM climates are found in regions where the surface is ice covered (ice-sheets or sea ice) during LGM but not during GS12 (Fig. 6). GS12 T 2 m is lower than during the LGM around 60 • N in the Pacific and Atlantic regions. This pattern is associated with a lower amplitude topographic wave in the GS12 than in the LGM simulation due to the reduced height and extent of the Laurentide ice sheet in GS12 as compared to LGM (Fig. 6).
Simulated annual mean T 2 m is compared to reconstructed temperature from δ 18 O and δ 15 N from ice core records for Greenland and Antarctica (see Sect. 2.3) in Table 2. Simulated T GS12 −T LGM and T GS12 − T RP differences are in agree-  ment (within ±3 • C) of the corresponding reconstructed temperature differences in Greenland (GRIP). The data from the EDC and Dome F ice cores in Antarctica give the T GS12 −T 1000 difference between GS12 and the last 1000 yr, which we compare to the T GS12 −T RP difference. The simulated T GS12 −T RP differences at these locations are larger by 1.5 to 4 • C than reconstructed T GS12 −T 1000 differences. The main discrepancy is found for absolute temperature (T GS12 ), with significantly warmer conditions in the simulation than in the reconstructions. Reconstructed temperature for GS12 is −50 • C at the location of both the GRIP (Greenland) and EDML (Antarctica) ice cores, and −51 • C at the location of the NGRIP (Greenland) ice core. Simulated T 2 m at these locations is −35 • C, −36 • C, and −36 • C for GRIP, EDML, and NGRIP, respectively. Due to the relatively coarse horizontal resolution in the model, the height of the ice sheet is reduced by up to 1000 m in regions of high topography in the simulation as compared to the actual height. This corresponds to a 6-7 • C difference in temperature, which can partly explain the discrepancy. Interestingly, the simulated T 2 m is higher than the reconstructed T GS12 over Greenland, in contrast to simulated SSTs which are lower than proxy SSTs over the North Atlantic (Fig. 3). Similarly, the few proxy SSTs located in the southern extratropics also indicate that simulated SSTs are lower than proxy SSTs, whereas the simulated T 2 m is higher than the reconstructed T GS12 over Antarctica. This discrepancy calls for further study, with focus on both proxy data and simulations.

Variability analysis
We hypothesise that not only the mean climate state, but also the variability is affected by the equilibration process. To test the importance of the equilibration for variability analysis, we analyse ENSO teleconnections in the GS12 simulation, and for comparison also in the LGM and RP simulations. ENSO teleconnections are chosen for comparison to Merkel et al. (2010), who analyse ENSO teleconnections in simulations of MIS 3 stadial, interstadial, and Heinrich stadial climate. Following Merkel et al. (2010), eastern tropical Pacific variability is identified with the Niño3 index, which is defined as SST anomalies from the long-term time mean, area-averaged over the Niño3 region (150 • W-90 • W, 5 • S-5 • N) (Trenberth, 1997). El Niño (La Niña) years are defined when the JFM Niño3 index exceeds (falls below minus) one standard deviation. ENSO teleconnections of MSLP are determined as the difference between seasonal averages over El Niño years and La Niña years. ENSO teleconnections that are statistically significant at the 95 % level, based on a two-tailed Student's t-test, are analysed. The analysis is based on the last 300 yr of the RP, LGM, and GS12 simulations. January-March is used for consistency with the results for the mean climate presented in the previous section.
ENSO teleconnections in MSLP are displayed in Fig. 7. In agreement with Merkel et al. (2010), we find that glacial (GS12 and LGM) boundary conditions induce major modifications to El-Niño-Southern Oscillation teleconnections as compared to the RP. The amplitude of sea level pressure anomalies in the North Atlantic region associated with El-Niño-La Niña variability is significantly reduced in glacial climates as compared to the RP climate. The strongest re- Isolines for geopotential height are displayed for every 100 geopotential meter. The coastline is defined at 20 % of the land fraction.
duction is found for LGM in agreement with Merkel et al. (2010). As described earlier, the increased extent and height of the North American topography from RP to GS12 and LGM results in an amplification of the NH stationary wave (Figs. 5 and 6). On the contrary, the amplitude of ENSO teleconnections in the North Atlantic region reduces with increasing extent and height of the North American topography. The explanation for these results should be sought in the non-linear interaction between the Rossby waves that communicate the ENSO signal from the Tropics to the midlatitudes and the mid-latitude dynamics (zonal mean, stationary waves, and transient eddies). Such an analysis is beyond the scope of this article. Also in agreement with Merkel et al. (2010), we find that the amplitude of Aleutian low anomalies associated with ENSO variability is reduced in the LGM simulation and enhanced in the GS12 simulation. Southern Hemisphere MSLP teleconnections are relatively similar in the RP, GS12, and LGM simulations presented here. In contrast, Merkel et al. (2010) find that ENSO teleconnections in Southern Hemisphere MSLP are similar in their pre-industrial and LGM climate, while it is significantly different in their MIS 3 stadial simulation.
The ENSO teleconnections in MSLP do not change significantly if model years 1139-1438 are used instead of the last 300 yr (model years 1239-1538) of the GS12 simulation ( Fig. 8; bottom panel). Thus, we conclude that once a quasi-equilibrium has been reached, the results of the variability analysis is not sensitive to which 300 yr period is chosen for the analysis. To test our hypothesis of the importance of equilibration, ENSO teleconnections in MSLP are also determined for a 300-yr period starting after 539 yr of integration in the GS12 simulation (model years 539-838). The signature of ENSO variability with an intensification of the Aleutian low and a change in the spatial pattern over northern North America, Greenland, and the North Atlantic as compared to the RP climate is common for this 300-yr period and the last 300 yr of the GS12 simulation ( Fig. 8; top panel). However, major differences are found in the amplitude of the North Atlantic pressure gradient from the Azores to Greenland, which is 50 % stronger in the period starting after 539 model years than in the last 300 yr. We note that the evolution of the ENSO MSLP teleconnections is not linear in time in our GS12 simulation. The North Atlantic pressure gradient associated with ENSO variability is weak in the initial state (i.e. the LGM climate), relatively strong in the 300-yr period starting after model year 539, and finally relatively weak again during the last 300 yr of the simulation. From this we conclude that equilibration is of importance for variability analysis. Zhang et al. (2008) analyse ENSO at 6 ka BP and 21 ka BP based on 50-100 yr of data from the PMIP2 simulations. Similarly, Merkel et al. (2010) base their analysis on 100 yr of data. In the present study, the analysis is based on 300 yr of data. To test the importance of the length of the time-period analysed for the results, the analysis of ENSO teleconnections was repeated for three 100-yr time-slices of the GS12 simulation (Fig. 9). The signature of ENSO variability with an intensification of the Aleutian low and a change in the spatial pattern over northern North America, Greenland, and the North Atlantic as compared to the RP climate is common for the last 100 yr ( Fig. 9; bottom panel) and the last 300 yr (Fig. 7; bottom panel) of the GS12 simulation. However, differences between the last 100 yr and the last 300 yr occur in LGM −8 −7 −6 −5 −4 −3 −2 −1 0 1 2 3 4 5 6 7 8 Units are hPa, contour interval 1 hPa (excluding the 0 hPa contour). In shading, statistically significant differences (at the 95 % level) are displayed. The coastline is defined at 20 % of the land fraction. the pressure gradient from the Azores to Greenland (which is 50 % stronger for the analysis based on 100 yr of data) and the strength of the Southern hemisphere mid-to high-latitude pressure gradient (which is 25 % weaker for the analysis based on 100 yr of data). Typically, ∼20 yr with high Niño3 index and ∼ 20 yr with low Niño3 index are found in a 100-yr time-slice. These results indicate that 100-yr time-slices are too short for an analysis of ENSO teleconnections, because too few samples are included to determine teleconnections in remote regions with large variability such as the North Atlantic region and the Southern Hemisphere high latitudes. ENSO teleconnections determined from 100-yr timeslices starting after 161 and 539 model years, respectively (Fig. 9), differ significantly from the last 300 yr of the simulation. The major differences occur over North America, Greenland, and the Nordic Seas. The spatial pattern in this region differs somewhat from the RP simulation ( Fig. 7; top panel), but the strength of the pressure gradient from the Azores to Greenland is similar to the RP in both 100-yr time-slices from the first part of the simulation. Thus, if the analysis was based on the 100-yr time-slice of model years 161-260 or 539-638, the (false) conclusion would be that ENSO teleconnections are essentially unchanged from the RP to GS12.

Discussion
The results in the present study are discussed with a focus on comparison to proxy data and other modelling studies of MIS 3 stadial climate (van Meerbeeck et al., 2009;Merkel et al., 2010).

Global mean T s response
The simulated annual global mean T s is 1.3 • C higher during GS12 than during the LGM, due to higher atmospheric GHG concentrations and less extensive continental ice sheets and sea ice cover. A similar difference between simulated stadial MIS3 and LGM conditions of 1.7 • C is found by van Meerbeeck et al. (2009), while Merkel et al. (2010) find no difference in annual global mean T s between stadial MIS 3 and LGM climate. The difference among the three studies in annual global mean T s response to the change from LGM to GS12 forcing and boundary conditions is associated with a difference in the response of the AMOC. The strength of the AMOC is unchanged from LGM to MIS 3 climate in van Meerbeeck et al. (2009) and in the present study, as compared to a substantial reduction in Merkel et al. (2010). We suggest that this difference in AMOC response offers a possible explanation for the difference in the annual global mean T s response. We speculate that the effect of a reduction in AMOC strength (Merkel et al., 2010), if associated with reduced equator to pole heat transport, partly cancels the global annual mean T s rise expected from higher atmospheric GHG concentrations and less extensive continental ice sheets as compared to the LGM simulation. These differences suggest that ocean circulation is important to determine the response of the climate system to changes in the forcing and boundary conditions.

SST response
Simulated SST is compared to proxy SSTs representative of GS12. Simulated and proxy SSTs are considered to be in agreement if the simulated SST is within the error margins of the proxy SST (see Sect. 2). The simulated GS12 climate is in agreement with proxy SSTs in 38-45 % of the proxy sites depending on season. MARGO Project Members (2009) find large discrepancies with respect to LGM SSTs recorded by different microfossil proxies. They find that SST LGM anomalies from the recent past based on foraminifera and dinocyst assemblages are smaller (i.e. warmer) than SST LGM anomalies based on Mg/Ca ratio and alkenones. In the present comparison for GS12 less proxy data is available and thus no firm conclusions can be made for the different proxy data types. Only for NH summer is there a reasonable number of proxy SSTs from planktonic foraminifera (26 SSTs) and Mg/Ca (14 SSTs) to make a comparison. Simulated SST is within the error margins of the proxy SST for 42 % of the foraminifera-based SSTs and in 50 % of the Mg/Ca-based SSTs. Proxy SSTs are higher than the simulated SST in the central North Atlantic in the present study, in contrast to an earlier study of MIS 3 stadial climate in which proxy SSTs were lower than the simulated SST in this region (van Meerbeeck et al., 2009). A possible explanation for the difference between our results and van Meerbeeck et al. (2009) is the response of the AMOC to glacial (MIS 3 stadial and LGM) forcing and boundary conditions as compared to the RP climate. CCSM3 simulates substantially reduced AMOC strength in GS12/LGM as compared to the RP associated with a colder North Atlantic region, whereas the intermediate complexity climate model used by van Meerbeeck et al. (2009) simulates a relatively unchanged AMOC strength in MIS 3 stadial/LGM as compared to the RP. Again, we note that the response of the oceanic circulation is of great importance to determine the response of the climate system to changes in the forcing and boundary conditions.

AMOC response
The AMOC is reduced by 50 % in the GS12 and LGM simulations presented here as compared to the RP (Kjellström et al., 2010;Brandefelt and Otto-Bliesner, 2009). As described by Kjellström et al. (2010), there are some indications from proxy records of a slowdown of the AMOC during LGM. CGCM simulations of LGM climate give very different glacial AMOC, ranging from ∼20 % weakening to ∼20-40 % strengthening of the circulation . Lynch-Stieglitz et al. (2007) review proxy data to determine the current knowledge of the AMOC during the LGM for comparison to CGCM simulations of LGM climate. They find support for an active circulation in the deep Atlantic, but not for a strong, deep AMOC. Quantitative estimates of the oceanic circulation rate for LGM can, however, not be made at present due to insufficient data (Wunsch, 2003). Even less data is available for GS12.
The AMOC reduction in the GS12 simulations presented here is associated with an annual mean sea ice expansion in the Nordic Seas (Fig. 3). Nevertheless, the western parts of the Nordic Seas are ice free in NH summer in the GS12 simulation. These changes present two competing effects on the deep water production in the Nordic Seas. The increased sea ice extent inhibits heat flux from the ocean to the atmosphere, which would cool and thus increase the density of the surface waters in the Nordic Seas. The formation of sea ice during autumn and winter, on the other hand, leads to brine rejection, which results in increased density. We may conclude that the reduction of the AMOC found in the GS12 simulation is in agreement with an increased sea ice expanse and that this effect dominates over the effect of increased brine rejection in the western Nordic Seas. In the coupled climate system the interaction between the ocean, sea ice, and atmospheric dynamics is non-linear and it is therefore difficult to identify cause and effect.
As described in Sect. 4.2, the RP simulation with CCSM3 is biased towards too extensive ice covers in the Labrador Sea and the Sea of Okhotsk in NH winter. This bias may partly explain the difference between the simulated and the proxy inferred sea ice extent in the North Atlantic region. Since all climate models are biased as compared to observations, the results of climate modelling experiments should always be taken with caution. To improve the understanding of the dynamics of the Earth's climate we therefore need a multimodel approach. Authors van Meerbeeck et al. (2009) find that a better agreement with proxy data for stadial MIS 3 climate is obtained with an additional freshwater forcing, which slows down the AMOC. They conclude that their results suggest that freshwater forcing is necessary to return climate from warm interstadials to cold stadials during MIS 3. They therefore suggest a change in perspective, making the stadial climate a perturbed climate state rather than a typical near-equilibrium MIS 3 climate. In contrast to van Meerbeeck et al. (2009), we find that stadial MIS 3 climate in reasonable agreement with proxy SST data can be simulated without additional freshwater forcing. These results indicate that the quasi-equilibrium climate simulated under MIS 3 boundary and forcing conditions, and hence conclusions drawn regarding the dynamics of the MIS 3 climate, are model-dependent.

Equilibration
Climate model simulations can be divided into two categories: equilibrium and transient simulations. Equilibrium simulations are set up with constant forcing and boundary conditions and should be integrated until the simulated climate has adjusted to the new conditions. When the adjustment is completed the system will be in an equilibrium state, or possibly in a state with multiple equilibria, and the unforced variability of the system can be analysed. The initial condition is of importance for the time required to reach the new equilibrium, but is assumed not to influence the final results. Transient simulations on the other hand are started from an initial condition that represents a particular state of the system (e.g. the recent past climate) and integrated forward with prescribed variations in forcing and boundary conditions to simulate the forced climate evolution. Due to limitations in the available computing time, equilibrium simulations with state-of-the-art CGCMs are seldom integrated until the adjustment to the new conditions is completed. As shown in the present study, the equilibration has a significant effect on both the time mean climate and the variability. We therefore argue that the simulated equilibration process should be discussed in studies including simulations of past climate.
MIS 3 stadials recorded in the δ 15 N record from Greenland ice cores lasted for ∼500-1000 yr (Huber et al., 2006). Since this time is shorter than the time required for the ocean to equilibrate it is reasonable to assume that the climate system was not in complete equilibrium during MIS 3 stadials. Yet, the GS12 simulation presented in this study is an equilibrium simulation, which is set up with constant forcing and boundary conditions and integrated to quasi-equilibrium. It can be argued that the simulated climate should not be in equilibrium, since the real climate was probably not in equilibrium during this period. Provided that appropriate variations in the forcing and boundary conditions were used, a transient simulation of MIS 3 climate would be a better choice. However, this would require an appropriate initial condition and a substantial increase in the available computing resources. The purpose of the GS12 simulation presented in the present study is to test the ability of a state-of-the-art climate model to reproduce the climate during an MIS 3 stadial under constant forcing and boundary conditions.
Our GS12 simulation can also be compared to earlier studies that focused on the dynamics of the DO events recorded in Greenland ice cores. Rial and Yang (2007) show that DOlike variations can be simulated under constant forcing and boundary conditions in an earlier version of the intermediate complexity climate model used by, e.g. van Meerbeeck et al. (2009). In the present GS12 simulation, the major adjustment of the simulated climate occurs during the first ∼1200 model years, indicating that the climate is reaching a new equilibrium. There is no indication of large-amplitude variations in the climate simulated under constant GS12 forcing and boundary conditions. However, without continuing the integration for several model millennia, the possibility of multiple equilibria can not be ruled out.

Conclusions
We present results from a simulation of stadial MIS 3 climate performed with a state-of-the-art CGCM at relatively high resolution. The simulation is set up with constant forcing and boundary conditions and integrated for 1538 yr to quasiequilibrium. The aim of the present study is to investigate if a state-of-the-art CGCM can simulate stadial MIS 3 climate in agreement with proxy records of SST. Further, the results are compared to other simulations of MIS 3 climate.
The main findings are: -The simulated GS12 global annual mean T s is 5.5 • C lower than the simulated RP climate with maximum difference over high latitudes and over the continental ice sheets, and 1.3 • C higher than the LGM simulated climate, due to higher atmospheric GHG concentrations and less extensive continental ice sheets and sea ice cover.
-The AMOC is reduced by 50 % in the GS12 simulation as compared to the RP (Kjellström et al., 2010;Brandefelt and Otto-Bliesner, 2009). This reduction is attained without additional freshwater forcing, a method used in other studies to force an AMOC shutdown.
-The simulated SST is in agreement with (i.e. within the uncertainty range of) proxy SSTs in 38-45 % of the sites depending on season. Proxy SSTs are higher than simulated SSTs in the central North Atlantic, in contrast to earlier simulations of MIS 3 stadial climate in which proxy SSTs were found to be lower than simulated SST. This difference to earlier simulations is associated with a difference in the response of the AMOC to the imposed forcing and boundary conditions.
-Reconstructed temperature differences (T GS12 −T LGM and T GS12 − T RP ) from ice cores in Greenland and Antarctica are within ±4 • C of the corresponding simulated T 2 m differences. Absolute reconstructed temperatures in Greenland and Antarctica, however, differ from the corresponding simulated T 2 m by 15-16 • C. Part of this discrepancy is explained by the coarse horizontal resolution of the simulation, which results in a flatter topography than in reality.
-The simulated GS12 climate is in better agreement with proxy SSTs for MIS 3 stadials than the stadial MIS 3 climate simulated by van Meerbeeck et al. (2009). They find that the climate simulated with an additional freshwater forcing in the North Atlantic, that slows down the AMOC, is in better agreement with proxy data for stadials. Authors van Meerbeeck et al. (2009) therefore suggest that stadial climate should be viewed as a perturbed climate state rather than a typical, near-equilibrium MIS 3 climate. We show that the simulated MIS 3 climate, and hence conclusions drawn regarding the dynamics of the MIS 3 climate, is highly model-dependent.
-Although the annual global mean T s only changes by 0.10 • C from model years 500-599 to the last century of the simulation, differences of 2-8 • C in the seasonal mean T 2 m and 20-65 % in precipitation rate exist over the northern North Atlantic region (Fig. 2). Further, the difference between simulated and proxy SSTs in the North Atlantic region is reduced for the last century of the simulation as compared to the period (model years 500-599) described by Kjellström et al. (2010) (Fig. 4).
The reduction amounts to on average 0.3-0.6 • C, with the largest reduction (∼1-2 • C) found for the sites in the North Atlantic.
-Equilibration is also shown to influence the climate variability. ENSO teleconnections in MSLP determined for a 300-yr time-slice starting after 539 yr of the GS12 simulation are significantly different from the last 300 yr of the simulation. Further, we find that the results are not robust if the analysis of ENSO teleconnections is based on 100-yr time-slices, as in an earlier study of MIS 3 climate (Merkel et al., 2010).
-The simulation presented here is integrated for ∼1500 model years with constant forcing and boundary conditions to quasi-equilibrium. We find that this setup is useful, although MIS 3 climate may not have reached equilibrium. This setup allows us to test the ability of stateof-the-art climate models to reproduce the climate during a MIS 3 stadial under constant forcing and boundary conditions. Further, the results may be compared to earlier studies of the dynamics of DO events. -There is no indication of large-amplitude variations in the simulated climate under constant GS12 forcing and boundary conditions. However, without continuing the integration for several model millennia, the possibility of multiple equilibria can not be ruled out.
This article is limited to the study of stadial MIS 3 climate, the challenge of understanding the large differences between stadial and interstadial MIS 3 climate as inferred from proxy data remains.

Sea surface temperature proxy data
This Appendix presents the SST proxy data used for comparison to the CGCM results. Data for annual average conditions and Northern Hemisphere winter (January-March)  Oppo and Sun (2005) and summer (July-September) average conditions are included. The proxy data was compiled as an extension to the Voelker and workshop participants (2002) database. The data are given below separately for average, NH winter, and NH summer. References are given for each record in Tables 7, A, A, and A. The SST values are given as a central value plus/minus one standard error (1s), which is generally the analytical error. When appropriate the SSTs given were calculated as mean values over different estimates. In these cases the error given is the sum of the analytical error and the standard deviation of the different SST estimates included in the mean. As described in Sect. 2.2.3, we use a 2s error.
Also included are data from 4 records of February and August SSTs (MMD99-2331, MD95-2042, ODP Site 1060, and MD04-2845). These are used here for comparison to the seasonal data. Table A5. Proxy based SST data for annual mean conditions with ±1σ error (given either by the author(s) of the respective data or estimated by the authors of this paper).