Uncertainties in the modelled CO 2 threshold for Antarctic glaciation

. A frequently cited atmospheric CO 2 threshold for the onset of Antarctic glaciation of ∼ 780 ppmv is based on the study of DeConto and Pollard (2003) using an ice sheet model and the GENESIS climate model. Proxy records suggest that atmospheric CO 2 concentrations passed through this threshold across the Eocene–Oligocene transition ∼ 34 Ma. However, atmospheric CO 2 concentrations may have been close to this threshold earlier than this transition, which is used by some to suggest the possibility of Antarctic ice sheets during the Eocene. Here we investigate the climate model dependency of the


Introduction
The first continental-scale Antarctic ice sheet formed during the Eocene-Oligocene transition (EOT) ∼ 34 Ma (Zachos et al., 2001). The extent of Antarctic glaciation prior to this event is disputed (e.g. Miller et al., 2005;Barker et al., 2007b;Gasson et al., 2012). Although various explanations for the cause of Antarctic glaciation have been suggested, such as the formation of the Antarctic Circumpolar Current due to the opening of ocean gateways (e.g. Kennett, 1977;Barker et al., 2007a), arguably the leading hypothesis at present is that Antarctic glaciation was caused by decreasing atmospheric CO 2 concentrations coupled with a favourable astronomical configuration (DeConto and Pollard, 2003). This hypothesis is supported by both ice sheet modelling and climate modelling studies (DeConto and Pollard, 2003;Huber et al., 2004), and proxy records for atmospheric CO 2 (Pagani et al., 2005(Pagani et al., , 2011Pearson et al., 2009). A commonly cited threshold for Antarctic glaciation of 2.8× pre-industrial CO 2 concentration (PIC) (∼ 780 ppmv) is based on the modelling study of DeConto and Pollard (2003), who used an ice sheet model asynchronously coupled to the GENESIS climate model. Proxy records of atmospheric CO 2 suggest that this threshold of 2.8 × PIC may have been crossed at times earlier than the EOT (Beerling and Royer, 2011), raising the possibility of glaciation earlier than this event, during the Eocene . Although other modelling studies have also simulated Antarctic glaciation (e.g. Huybrechts, 452 E. Gasson et al.: Modelled CO 2 threshold for Antarctic glaciation 1993; Langebroek et al., 2009), with the study of Langebroek et al. (2009) suggesting a threshold of ∼ 2.2 × PIC, there has been limited work investigating to what extent the glacial CO 2 threshold is dependent on the climate model used. Here we perform offline ice sheet model (ISM) simulations using the climatology from a variety of GCMs (general circulation models), including the GENESIS GCM used by DeConto and Pollard (2003), to investigate the model dependence of the atmospheric CO 2 threshold for Antarctic glaciation.
The basis for this inter-model comparison is the EoMIP (Eocene Modelling Intercomparison Project) , which collated a number of pre-existing Eocene GCM simulations (Heinemann et al., 2009;Roberts et al., 2009;Lunt et al., 2010;Winguth et al., 2010;Huber and Caballero, 2011). This was an informal inter-model comparison because it was based on a number of independent studies, as a result the GCMs were not set up with identical boundary conditions (such as the astronomical configuration and palaeogeography). Although this precludes a direct assessment of model dependency, it is arguably more faithful to the true uncertainties associated with modelling this period, which has poorly constrained boundary conditions . In addition to the EoMIP simulations, we use Eocene simulations from GENESIS, CESM1.0 (Goldner et al., 2013) and FAMOUS (Sagoo et al., 2013). The aims of this paper are to perform ISM simulations using the climate output from a variety of climate models (HadCM3L, CCSM3, CESM1.0, GENESIS, FAMOUS, ECHAM5 and GISS_ER), compare these results with existing modelling studies, and to diagnose potential differences between the climate simulations used and sensitivity of Antarctic ice sheet growth to the background mean climate states.

Ice sheet model description
We use the Glimmer ISM in this paper. The mechanics of this model are documented in Rutt et al. (2009). Glimmer follows the conventions of a number of previous ISMs (e.g. Huybrechts, 1993;Abe-Ouchi and Blatter, 1993;Ritz et al., 1996;DeConto and Pollard, 2003). It makes use of the shallow ice approximation (SIA), a simplification of the ice sheet physics that significantly reduces computational expense (Hutter, 1983). Although higher-order and full Stokes ice sheet models exist (e.g. Morlighem et al., 2010;Seddik et al., 2012), their computational expense currently prohibits their use for the very long duration (10 4 -10 5 year) ice sheet equilibrium simulations conducted here. For example, Seddik et al. (2012) limited their simulations of the Greenland ice sheet using a full Stokes model to 100 years due to the computational expense of the model. The use of the SIA approximation prohibits the accurate simulation of ice streams or the transfer of mass across the grounding line from ter-restrial ice to floating ice shelves. In this paper we focus on the slow response of the large and predominantly terrestrial East Antarctic ice sheet (EAIS) on long timescales. Because of the lack of necessary dynamics in the ice sheet model used we make no attempt to simulate a marine-based West Antarctic ice sheet (WAIS). The ISM is set up with default settings, which has basal sliding turned off. The ISM has a spatial resolution of 20 × 20 km, and all the simulations are initiated from ice-free conditions.
An offline forcing methodology is used, whereby the climatology from the climate model (surface air temperature and precipitation) is used to force the ice sheet model with no subsequent feedbacks, other than height-mass balance feedback, on the climate system (e.g. Huybrechts and de Wolde, 1999;Lunt et al., 2008;Stone et al., 2010;Dolan et al., 2012). A lapse rate adjustment is made to the temperatures due to the spatial and vertical discrepancy between the GCM and ISM topographies (e.g. Pollard, 2010). All of the GCM simulations prescribe ice-free boundary conditions over the Antarctic . Previous modelling studies suggest that Antarctic glaciation generates a number of feedbacks on the climate system, such as changes in surface albedo, sea-ice and cloud cover (e.g. DeConto et al., 2007;Goldner et al., 2013). Although the lack of these feedbacks will not affect the threshold for the initial accumulation of ice from ice-free conditions, it may affect the rate at which full-scale glaciation occurs. We acknowledge the limitations of our methodology in representing these feedbacks. This methodology differs from that used by DeConto and Pollard (2003), who asynchronously coupled an ice sheet model to a climate model, allowing an approximation of feedbacks from the growth of an ice sheet on the climate system. Because we have included the GENESIS GCM in our intermodel comparison, we can compare our forcing methodology with the more sophisticated asynchronous coupling. The mass balance scheme adopted is the widely used positivedegree day (PDD) method (Reeh, 1991). Alternatives to the PDD method exist, such as physically based energy balance models (e.g. Bougamont et al., 2005), however these are not presently included in the Glimmer ISM. All ISM simulations are set up identically, with only the input climate and GCM topographies differing.

Bedrock topography
The Antarctic bedrock topography used within the ISM needs to be representative of the ice-free conditions prior to the onset of glaciation. There are four bedrock topographies which we use for these simulations (see Fig. 1), our motivation for using multiple bedrock topographies is to explore more fully the uncertainties associated with modelling this period. The first topography used is the presentday Bedmap1 topography (Lythe and Vaughan, 2001) with the ice sheet removed and accounting for isostatic adjustment (e.g. DeConto and Pollard, 2003), which is our default  Lythe and Vaughan (2001), rotated into early Eocene position (TOPO1). (b) A reduced-resolution version of the proprietary topography used by Lunt et al. (2010); we use a higher-resolution version for our ISM simulations than that shown here (TOPO2). (c) Minimum (TOPO3) and (d) maximum extent reconstructed Eocene/Oligocene topography of Wilson et al. (2012). Note the increase in land surface area above present-day sea level, in particular for the West Antarctic. topography (we denote as TOPO1). In addition we use the proprietary topography used by Lunt et al. (2010) for their GCM boundary conditions (here TOPO2) and the two reconstructed topographies of Wilson et al. (2012).
The EOT topographies generated by Wilson et al. (2012) attempt to take into account the erosion, thermal subsidence and plate movements which have occurred since the Eocene (see Luyendyk, 2009 andWilson et al., 2012 for a detailed description of the method). The reconstructions make use of models for sediment erosion and thermal subsidence, constrained by observed sediment volumes deposited around the Antarctic continent. Wilson et al. (2012) generated minimum-extent (we denote as TOPO3) and maximumextent (TOPO4) reconstructions based on different target sediment volumes, due to uncertainties in offshore sediment volumes. Wilson et al. (2012) do not claim that these are accurate reconstructions of the Eocene/Oligocene topography, but argue that they are two plausible end-members. Based on these reconstructions, the accommodation space of the Antarctic continent would have been greater at the EOT than present. The total area above present-day sea level is 12.4 × 10 6 km 2 and 13.1 × 10 6 km 2 for the minimum and maximum reconstructions (Wilson et al., 2012), compared to 10.7 × 10 6 km 2 and 11.1 × 10 6 km 2 for the TOPO1 and TOPO2 reconstructions, respectively.
The majority of the increase in continental area for TOPO3 and TOPO4 is for the West Antarctic. Importantly, Wilson et al. (2012) suggested that during the EOT the West Antarctic continent could have supported a largely continentalbased ice sheet, rather than a marine-based ice sheet as is present today. All of the Eocene GCM simulations available to us have a deglaciated Antarctic and largely submerged West Antarctic. As such, it is possible that the climate would differ if the reconstructions of Wilson et al. (2012) were used for the GCM boundary conditions. Although we will use the Wilson et al. (2012) topographies for sensitivity tests, it is with the caveat that the climate forcing provided to the West Antarctic is from GCM simulations which may have ocean cells over regions which are land in the reconstruction of Wilson et al. (2012). To test the significance of the Wilson et al. (2012) topographies to the formation of the ice sheets at the EOT more accurately, it would be necessary to repeat the GCM simulations using a palaeo-geography which incorporates the Wilson et al. (2012) Antarctic topography.

GCM simulations
The GCM simulations used here are based on a number of previously published independent studies, and as such the GCM boundary conditions are not identical . Although the EoMIP GCM simulations have slightly different boundary conditions, they are broadly similar in that they use an early Eocene palaeo-geography and have prescribed ice-free conditions over Antarctica. Note that EoMIP originally focused on coupled oceanatmosphere GCM simulations only and therefore did not include the GENESIS atmosphere-slab ocean GCM simulations which we have included here. Two separate studies used the CCSM3 model with a slightly different configuration, we denote these as CCSM3_H (Huber and Caballero, 2011) and CCSM3_W (Winguth et al., 2010). We add simulations from two recently published studies using CESM1.0 (Goldner et al., 2013) and FAMOUS (Sagoo et al., 2013); the latter is a reduced complexity version of HadCM3L.
The GCMs used here have been evaluated previously against modern-day observations (e.g. ). It should be noted that modern-day performance may not be relevant to performance under Eocene boundary conditions. Connolley and Bracegirdle (2007) evaluated 4 of the GCMs used here (excluding GENESIS, CESM1.0 and FAMOUS) against 15 other GCMs (used in the IPCC AR4) for their performance compared with Antarctic re-analysis output. They assigned skill scores based on five variables (mean sea level pressure, height and temperature at 500 hPa, sea surface temperature, surface mass balance), giving a skill score between 0 (low skill) and 1 (high skill). Over the Antarctic region (defined as latitudes greater than 45 • S), ECHAM5 had the highest skill score (0.45) of the 15 GCMs based on the 5 chosen variables,   Lunt et al. (2012) for a full description of the simulations. Astronomical parameters: eccentricity (ecc.), obliquity (obl.) and longitude of precession (pre.), with insolation (ins.) for January at 70 • S also shown (Wm −2 ). CS is the modernday equilibrium climate sensitivity for the GCMs, excluding vegetation and chemical feedbacks. with HadCM3 (0.36) and CCSM3 (0.28) 4th and 7th, respectively, and GISS_ER (0.11) 14th. For Antarctic sea surface temperatures the skill of all of the models was low, in part due to the method used to measure skill, however ECHAM5, GISS_ER and HadCM3 were in the top half of the 15 GCMs. HadCM3 had the joint best skill score for surface mass balance over the Antarctic, with CCSM3 and ECHAM5 also scoring highly (> 0.9), however GISS_ER had a low skill score (0.07) (Connolley and Bracegirdle, 2007). The astronomical configuration has been shown to be important for Antarctic glaciation (DeConto and Pollard, 2003;Langebroek et al., 2009); the astronomical configurations vary between the GCM simulations used here, although they are broadly similar (see Table 1 and Fig. S1 in the Supplement). The simulations for HadCM3L, CCSM3_H, CESM1.0 and FAMOUS use the modern astronomical configuration, whereas the ECHAM5 and GISS_ER simulations have greater eccentricity and the GENESIS and CCSM3_W simulations have zero eccentricity. The astronomical configuration used for the GENESIS, GISS_ER and CCSM3_W simulations are likely to be the most favourable for Antarctic glaciation, whereas ECHAM5 has the the least favourable astronomical configuration, based on peak insolation during the austral summer. There are additional simulations for HadCM3L and GENESIS at 2× and 4× PIC, which use astronomical parameters resulting in extremes of summer insolation.
There are additional differences in the GCM boundary conditions. The vegetation prescribed varies, with the GISS_ER and CCSM3_H simulations adopting the vegetation maps of Sewall et al. (2000), CCSM3_W using the vegetation of Shellito and Sloan (2006), the HadCM3L simulation using homogeneous shrubland, and the ECHAM5 simulation prescribing homogeneous vegetation resembling a present-day savanna. All of the simulations have present-day aerosol loading, with the exception of the CCSM3_H simulation which has a reduced aerosol loading. The adoption of this reduced aerosol loading is justified by possible reduced ocean productivity leading to reduced dimethyl sulphide (DMS) production (Huber and Caballero, 2011;Kump and Pollard, 2008). Because of the reduced aerosol load in the CCSM3 simulation of Huber and Caballero (2011), surface temperatures are increased. The global mean surface air temperature of the CCSM3_W 4× PIC simulation is approximately equivalent to the CCSM3_H 2× PIC simulation, largely due to the different approach to aerosol loading . These differences in boundary conditions are important but are representative of plausible boundary conditions; this gives insight into the decisions required when modelling relatively data poor periods, such as the Eocene.
The FAMOUS simulation differs from the other simulations as it was selected from a 100-member parameter ensemble (varying 10 parameters, see Sagoo et al., 2013) based on closest agreement with early Eocene proxy data (Sagoo et al., 2013). The main aim of their paper was to simulate a reduced meridional temperature gradient, which is suggested by proxy data to have occurred in the warmth of the early Eocene (Sagoo et al., 2013). The simulations within EoMIP were also evaluated against proxy data . The EoMIP simulations had closest agreement with proxy data at higher atmospheric CO 2 concentrations. The simulation with the closest agreement with the proxy records was the CCSM3_H simulation at 16× PIC. However, not all of the GCMs were run at the same atmospheric CO 2 concentrations, precluding a direct evaluation of model performance (see Lunt et al., 2012 for a detailed discussion of model performance). Additionally, atmospheric CO 2 is poorly constrained by proxy records in the Eocene (Beerling and Royer, 2011), making assessment of model skill in the Eocene problematic .
The GCM simulations were performed at atmospheric CO 2 concentrations ranging from 1× to 16× PIC (see Table 1). We first perform equilibrium simulations using the climate output at fixed atmospheric CO 2 concentrations. Additionally, for GCMs where simulations were performed at multiple CO 2 concentrations (HadCM3L, CCSM3, CESM1.0 and GENESIS) we perform transient CO 2 experiments by scaling between the simulations following a logarithmic relationship between atmospheric CO 2 and climate (C).

Equilibrium simulations
We first describe results from the equilibrium (50 kyr) ice sheet model simulations using the climate output from the GCM simulations. For the GCMs with simulations performed with multiple astronomical configurations, we select the configuration closest to modern. As can be seen from Fig. 2, the offline simulations using the climate output from CCSM3_H, CESM1.0 and ECHAM5 produce large ice sheets over much of East Antarctica at 2× PIC (10.3-14.6 × 10 6 km 3 ) and GENESIS produces a full continentalsized EAIS at 2× PIC (28.6 × 10 6 km 3 ). However, there is minimal ice in the equivalent 2× PIC simulation using HadCM3L (0.3 × 10 6 km 3 ). Even when using a 1× PIC HadCM3L simulation (not shown) minimal ice forms (1.1 × 10 6 km 3 ). The FAMOUS simulation is completely ice-free at 2× PIC. For CCSM3_H, CESM1.0 and ECHAM5, ice nucleates over Queen Maud Land and the Gamburtsev Mountains. These two smaller ice sheets combine to generate an intermediate-sized ice sheet in the 2× PIC simulations.
The 4× PIC simulations are shown in Fig. 3. There is a relatively large ice sheet for the 4× PIC simulation using CCSM3_H (9.4 × 10 6 km 3 ). The ice sheet in the 4× PIC simulation using CCSM3_H is only ∼ 35 % smaller than for the 2× PIC simulation. This is plausibly a result of the relatively low CO 2 sensitivity of CCSM3 (Huber and Caballero, 2011). This is in contrast to the CESM1.0 simulation, which is mostly ice-free at 4× PIC, likely a result of the higher CO 2 sensitivity of CESM1.0 compared to CCSM3 (although note that GENESIS also has a relatively low CO 2 sensitivity). We also performed offline simulations using the output from CCSM3_H at 8× and 16× PIC (not shown). The simulation with CCSM3_H at 8× PIC generated minimal ice, with a total volume of 0.2 × 10 6 km 3 , and the simulation at 16× PIC was ice-free. This suggests that the glacial threshold for these CCSM3_H simulations is between 8× and 4× PIC. The differences in GCM boundary conditions result in different-sized ice sheets between CCSM_H and CCSM_W at 4× PIC.
Between 4× and 2× PIC a full continental-sized ice sheet forms in the offline simulations using the GENESIS model. This is the same GCM used by DeConto and Pollard (2003) and produces a similar result to their glacial CO 2 threshold. The simulation using the GISS_ER model is for 4× PIC and 7× CH 4 compared to pre-industrial concentrations. Roberts et al. (2009) estimate that the GISS_ER simulation is equivalent to a 4.3× PIC simulation. When we use the climate output from the GISS_ER simulation to force the ISM, it generates a small ice cap over Queen Maud Land, this is a slightly higher volume than the 4× PIC HadCM3L simulation.

Transient simulations
In addition to the equilibrium simulations we next present transient atmospheric CO 2 simulations, in order to better define the CO 2 thresholds. The climate is created by linearly scaling between the GCM simulations at different atmospheric CO 2 concentrations over 1.5 Myr (a rate of CO 2 decrease of ∼ 1 ppm kyr −1 , which is comparable to proxy records for atmospheric CO 2 across the EOT (Pagani et al., 2011)). This is only possible for the GCMs where simulations are available at more than one atmospheric CO 2 concentration, these being HadCM3L, CCSM3_H, CESM1.0 and GENESIS. This scaling is based on the equation for climate sensitivity (e.g. Solgaard and Langen, 2012): where C 2× and C 4× is the climate (temperature and precipitation) for the 2× and 4× PIC GCM simulation, respectively, and CO 2 is the atmospheric CO 2 concentration at the current time step. We are therefore calculating an Earth system sensitivity based on these 2 GCM simulations, this may differ from the climate sensitivities for the GCMs under modern boundary conditions, which are included in Table 1 for reference. The model checks for potential negative values for precipitation resulting from this scaling and resets these to zero. We calculate the CO 2 threshold for the formation of an intermediate (which we define here as 25 m Eocene sea level equivalent (SLE)) and large (40 m Eocene SLE) ice sheet. Ice volumes are converted to Eocene sea levels by accounting for the change in state from ice to seawater and dividing by the total Eocene ocean surface area (372.9 × 10 6 km 2 ; DeConto et al., 2008). In the simulations of Pollard and De-Conto (2005) using an earlier version of the GENESIS GCM with a constant astronomical forcing, the glacial threshold was 2.1× PIC for an intermediate ice sheet and 1.6× PIC for a large ice sheet (shown in Fig. 4). For the equivalent simulations including astronomical forcing, the CO 2 thresholds were higher, at ∼ 3.0× PIC and ∼ 2.8× PIC (Pollard and De-Conto, 2005). Similar results were also found by Langebroek et al. (2009)   forcing, and 1.6× PIC for the constant astronomical forcing experiment (Langebroek et al., 2009). In these transient CO 2 experiments, we scale between 6× and 0.5× PIC over 1.5 Myr using the climate data from HadCM3L, CCSM3_H, CESM1.0 and GENESIS. We interpolate between the simulations at 4× and 2× PIC and then extrapolate for CO 2 values outside this range (for 6-4× and 2-0.5× PIC). Note that by extrapolating we are introducing error, however this error is relatively small compared with the inter-model disagreement (see supplementary information for a comparison of extrapolated climatologies with GCM control climatologies). Because simulations are only available at one atmospheric CO 2 concentration for the other GCMs, we cannot estimate the CO 2 thresholds for these models. However, based on the results of the offline simulations, for the 2× simulation using ECHAM5 an intermediate ice sheet has formed (∼25 m Eocene SLE), suggesting the threshold for a large ice sheet (∼40 m Eocene SLE) is below 2× PIC. For GISS_ER and CCSM3_W, the threshold for glaciation is below 4.3× PIC and 4× PIC, respectively.
In addition to the simulations using a constant astronomical configuration, we perform an experiment including astronomical variability, based on the solutions of Laskar et al. (2004). For this experiment we use the climate output from GENESIS and scale between the simulations with different astronomical configurations using where I is the insolation at 70 • S averaged over the 6 months with peak insolation, and I c , I m and I w are the insolation values at 70 • S from the GCM simulations and C c , C m and C w are the respective climate outputs (Gasson, 2013). The transient CO 2 experiments are shown in Fig. 4. An intermediate ice sheet (25 m Eocene SLE) forms at 3.3× PIC in the experiment using CCSM3_H, 2.9× PIC in the experiment using GENESIS and 2× PIC when using CESM1.0. Again the lack of ice in the experiment using HadCM3L is clearly evident, with a small increase in ice volumes below 1× PIC. A large ice sheet (> 40 m Eocene SLE) forms at 2.8× PIC in the experiment using GENESIS and 1.8× PIC for CESM1.0. Recall that none of these experiments include albedo feedbacks nor feedbacks on precipitation, which may affect the glacial CO 2 thresholds.
The pattern of ice growth varies between GCMs. The CCSM3_H experiment has three distinct steps in ice growth; CESM1.0 has multiple smaller steps, whereas for GENESIS there is one major threshold. The study of DeConto and Pollard (2003), using an earlier version of the GENESIS GCM and an asynchronous coupling method, showed the growth of ice in a series of steps as ice first formed as isolated ice caps in the mountain regions. It therefore appears unusual that our experiment, using a later version of GENESIS, does not show this pattern. However, more recent simulations based on a modified method of that used by DeConto and Pollard (2003) and the same version of GENESIS we use here, also lack the stepped pattern to ice growth (Pollard, personal communication, 2012). Also note the greater ice volume of our GENE-SIS simulations compared with that of Pollard and DeConto (2005) at equivalent atmospheric CO 2 concentrations, this is likely due to the lack of basal sliding in the simulations presented here.
For the GENESIS simulation including a representation of astronomical variability the glacial threshold is at a slightly higher atmospheric CO 2 concentration, the threshold for the growth of a large ice sheet being 3.2× PIC compared with 2.8× PIC for the constant astronomical configuration simulation. This is consistent with the results of DeConto and Pollard (2003) and Langebroek et al. (2009).

Sensitivity to lapse rate and topography
We next present sensitivity tests in order to determine how changing certain poorly constrained parameters affects the glacial CO 2 thresholds. Firstly, we highlight the impact of changing the lapse rate (we use a default value of −7 K km −1 ). The lapse rate has two effects, firstly it allows for the cooling of the ice sheet surface as it rises vertically through the atmosphere. Secondly, the lapse rate is used to scale from the coarse GCM surface topography onto the finer topography used within the ISM. Values for the lapse rate parameter can vary spatially and temporally, largely due to changes in the moisture content of the atmosphere. In a GCM study, Krinner and Genthon (1999) found values for the lapse rate as low as −10 K km −1 for the dry continental interior above continental-sized ice sheets, such as the EAIS. For the moister coastal regions, values as high as −5 K km −1 were found, these values are comparable to empirical results (Magand et al., 2004). As our ISM domain covers the entire Antarctic continent, the default lapse rate chosen is an approximation between these two environments. CCSM3_H and GENESIS, using lapse rates of −6, −7 and −8 K km −1 (see Fig. 5); the default value used in the previous experiments was −7 K km −1 , chosen for consistency with DeConto and Pollard (2003). As can be seen from Fig. 5, the simulations using CCSM3_H are highly sensitive to the value chosen for the lapse rate parameter. With the threshold for the growth of an intermediate ice sheet varying between 1.2× and 5.9× PIC for lapse rates between −6 and −8 K km −1 . With the lower value for the lapse rate, the threshold for the growth of a large ice sheet is crossed at 2.4× PIC. The simulations using GENESIS are less sensitive to the value for the lapse rate, with the threshold for the growth of an intermediate ice sheet varying between 2.6× and 3.3× PIC for the three values for the lapse rate. Similar simulations were also performed using HadCM3L, however these had little impact on the low ice volumes seen in the previous HadCM3L transient CO 2 experiments and are therefore not shown here.
The reason for the strong sensitivity of the CCSM3_H experiment to the lapse rate parameter is due to the Antarctic topography within the GCM. For the simulations using CCSM3_H, the Antarctic topography within the GCM (from the Sewall et al. (2000) palaeo-topography) is significantly lower than the ISM topography. This is evident in the maps shown in Fig. 6. The discrepancy between the GCM and ISM topography for the CCSM3_H simulations exceeds 1 km in certain regions. The Antarctic GCM topography within CCSM3_H (and also ECHAM5, CCSM_W, CESM1.0 and GISS_ER) resembles the present-day Antarctic bedrock topography without isostatic adjustment. Because of this, there is a large lapse rate correction to the surface temperatures as they are scaled from the GCM topography to the ISM topography. This results in the high sensitivity to the value for the lapse rate parameter. This could also explain the GCM results of Huber and Nof (2006), which did not find snow accumulation over the Antarctic in an experiment with an earlier version of CCSM_H. They used the same GCM bound- For the GCM simulations using GENESIS, the GCM topography is much closer to the ISM topography, therefore the ISM simulations are less sensitive to the lapse rate parameter. The Antarctic topography in the simulations using CCSM3, CESM1.0, ECHAM5 and GISS_ER are all significantly less mountainous than the ISM topography and are therefore all likely to be sensitive to the value chosen for the lapse rate parameter. The Gamburtsev mountain range in the centre of the East Antarctic continent is much lower in elevation for these GCM simulations. Although there is uncertainty as to the past uplift history of the Antarctic, the Gamburtsev Mountains are thought to have formed earlier than the Eocene (Cox et al., 2010). This difference in GCM topography over the Antarctic may also affect precipitation patterns, in addition to surface temperatures. Therefore the disagreement between the ISM simulations in Fig. 4 may be due to differences in the GCM boundary conditions, in addition to differences between the GCMs. The differences are therefore a combination of inter-model disagreement and experimental design.
We present further sensitivity tests using four different Antarctic ISM topographies. The ISM topographies we use are TOPO1, the default topography used in the previous experiments; TOPO2, the proprietary topography used by Lunt et al. (2010); and TOPO3 and TOPO4, the minimum and maximum reconstructed topographies of Wilson et al. (2012), respectively (see Fig. 1). Note that all of these topographies are more mountainous than the GCM topography used in the CCSM3, CESM1.0, ECHAM5 and GISS_ER simulations. The glacial CO 2 threshold is sensitive to the choice of Antarctic bedrock topography (Fig. 7). When using the Wilson et al. (2012) topographies (TOPO3 and TOPO4), the onset of glaciation is at a slightly higher atmospheric CO 2 concentration than the default topography (TOPO1). This is especially evident for the experiments using CCSM3_H. This is due to a slightly higher elevation of the mountains in Queen Maud Land and the Gamburtsev Mountains, the regions where ice first nucleates. The difference in mountain elevation is likely a result of the different isostasy models used for our default topography and that used by Wilson et al. (2012). Similar to the previous experiments, a large ice sheet does not form in the CCSM3_H experiments (the lapse rate is −7 K km −1 ). For the GENESIS experiments, the maximum size of the ice sheet varies due to differences in the total Antarctic surface area between the different topographies. For the maximum reconstruction of Wilson et al. (2012) (TOPO4), an ice sheet of 32.5×10 6 km 3 (78 m Eocene SLE) has formed at 2× PIC. This increased ice volume compared to the default topography experiment is largely due to the growth of a continental-based WAIS.

Diagnosing differences between simulations
It is not immediately clear why there is such variation in ice volumes caused by the different climate forcing, in particular why the ISM simulations using the HadCM3L early Eocene simulations of Lunt et al. (2010) and the FAMOUS simulations of Sagoo et al. (2013) should generate such low ice volumes. Although Lunt et al. (2012) noted certain differences between the GCM simulations within EoMIP, their analysis did not identify a disagreement which could explain our ISM results. The variables which are passed to the ISM from the GCM output data are the annual mean air temperature (T a ), annual air temperature half range ( T a ), which is the difference between the warmest month and the annual mean, and the total precipitation (P ). Much of the analysis by Lunt et al. (2012) focused on the annual means from the GCMs. Interestingly, their analysis suggested that when looking at Table 2. Climate variables passed to the ISM from GCM simulations, shown as averages over the East Antarctic continent, with averages at elevations above 1500 m in parentheses. T a is the annual mean air temperature, T a is the annual air temperature half range (difference between the warm month and the annual mean temperature) and P is total annual precipitation. For 2× PIC (upper rows) and 4× PIC (lower rows) simulations the annual mean air temperatures, HadCM3L is cooler than CCSM3_H and ECHAM5 for the 2× PIC simulations. These relatively cool annual mean air temperatures are especially pronounced for the Southern Hemisphere high latitudes. The three variables which are passed to the ISM from the GCM (following lapse rate correction) are summarized in Table 2 as averages over the East Antarctic continent and also as averages over the mountainous regions (> 1500 m), the regions where ice tends to first nucleate. As can be seen from Table 2, HadCM3L has the lowest annual mean air temperature of all of GCM simulations over the East Antarctic continent for the 2× and 4× PIC simulations. FAMOUS is by far the warmest of the simulations at 2× PIC, which explains the lack of ice growth for this simulation.
To investigate the impact of these three climate variables on the ice sheet model results, we use the PDD mass balance scheme to calculate the potential snowmelt (â s ) for various annual mean air temperatures and annual air temperature half ranges. If the total annual precipitation exceeds the potential snowmelt then snow will accumulate. If the total annual precipitation is less than the potential snowmelt than there is no year-to-year snow accumulation and an ice sheet cannot grow. The potential snowmelt is calculated from the PDD sum and the PDD factor for snow (Reeh, 1991): where α s is the PDD factor for snow (3 mm d −1 • C −1 ) and D p is the PDD sum. We use the mass balance scheme to calculate D p using The inner integral in practice is evaluated between 0 and 50 • C, σ T is the standard deviation of temperature fluctuations with a value of 5 • C used, A is the period of the year and T a is the daily surface air temperature calculated using We numerically compute the potential snowmelt (contours in Fig. 8) according to Eqs. (3)-(5) for a range of values for T a and T a . Also shown in Fig. 8 are the values for T a and T a from the GCM simulations, as averages over the Antarctic mountain regions. As can be seen from Fig. 8, the high annual mean air temperatures of the FAMOUS simulation generate a very high potential snowmelt, explaining the lack of ice in this simulation. For HadCM3L, despite the low annual mean air temperatures over the mountainous regions of Antarctica, the potential snowmelt is still relatively high at 2× PIC. This is due to the large annual air temperature half range in the HadCM3L simulations. The potential snowmelt in the HadCM3L 2× PIC simulation is comparable to the CCSM3_H 4× PIC simulation. This CCSM3_H 4× PIC simulation generated a large ice sheet, whereas the HadCM3L simulation did not. The total annual precipitation for the CCSM3_H 4× PIC simulation is approximately double that of the HadCM3L 2× PIC simulation over the East Antarctic. This would suggest that the low precipitation in the HadCM3L simulations is also a significant factor. Based on the Clausius-Claperon relation, the low precipitation is itself likely to be a result of the low annual mean air temperatures. In an idealized simulation where we arbitrarily double the HadCM3L precipitation, a large ice sheet (18.1 × 10 6 km 3 ) forms for the 2× PIC simulation. This ice sheet differs from the other simulations and nucleates on Victoria Land and Wilkes Land, instead of Queen Maud Land and the Gamburtsev Mountains. This would suggest that even with precipitation arbitrarily doubled, the region around Queen Maud Land and the Gamburtsev Mountains remains precipitation-limited for the HadCM3L simulation.
The total annual precipitation and potential snowmelt values we have shown in Fig. 8 are averages over the mountainous regions. As these values vary spatially, ice can grow for simulations where the mean annual precipitation is lower than the potential snowmelt. For example, the potential snowmelt for the 4× PIC CCSM3_H simulation is above the mean annual precipitation for the mountainous regions yet still produced a large ice sheet. Additionally, this data is for ice-free conditions at the first time step, and therefore does not include height-mass balance feedback or ice flow from regions of initial ice nucleation.

GCM seasonality
Given the importance of seasonality to our ice sheet model results, in particular for HadCM3L, we next discuss the different seasonalities of the GCMs. As previously noted, the astronomical configurations are not identical between GCM simulations but are similar; the HadCM3L, FAMOUS, CCSM3_H and CESM1.0 simulations all have a modern astronomical configuration. Maps of annual temperature range (from the warmest month minus the coldest month) are shown in Fig. 9, and we show global maps of seasonality in order to show any potential inter-hemispheric biases.
It is interesting to note that the GENESIS simulations have a relatively high annual temperature range over the Northern Hemisphere, but not the Southern Hemisphere. This pattern is unique to GENESIS amongst the 2× PIC simulations, although the GISS_ER 4× PIC simulation also shows a similar pattern. GENESIS is the GCM used by DeConto et al. (2008) in their study investigating the thresholds for Northern Hemisphere glaciation. Their study suggested that the threshold for Northern Hemisphere glaciation is ∼ 280 ppmv, providing evidence against the early Northern Hemisphere glaciation hypothesis. This hypothesis was based on evidence from ice-rafted debris in the Eocene and Oligocene (Tripati et al., 2005;Eldrett et al., 2007), and discrepancies between benthic δ 18 O and Mg / Ca records across the EOT (Lear, 2000), although this second issue has now largely been resolved (DeConto et al., 2008;Liu et al., 2009;Wilson and Luyendyk, 2009). Given the strong seasonality seen in the GENESIS simulations in the Northern Hemisphere, it would perhaps be interesting to repeat the experiment of DeConto The strong seasonality in the HadCM3L simulations is not just a result of very warm summers, but also cool winters. As can be seen from Fig. 9, for the early Eocene 2× PIC simulations using HadCM3L there is a very large annual temperature range over Antarctica. For HadCM3L, the annual range in surface air temperature over Antarctica exceeds 60 • C in certain regions. The HadCM3L 4× PIC simulation has a slightly lower seasonality than the 2× PIC simulation, but the seasonality is still greater than for any of the other GCMs at 4× PIC. This very large annual temperature range for HadCM3L is also apparent in the high latitude Northern Hemisphere. None of the other GCMs exhibit such a large annual temperature range in both hemispheres. Sensitivity tests using HadCM3L simulations with different astronomical forcing, including a simulation favourable to Southern Hemisphere glaciation (Lunt et al., 2011), did not generate any significant increase in ice volumes (not shown).
To investigate whether the strong HadCM3L seasonality is a result of the early Eocene boundary conditions, or a model bias, we have plotted seasonality maps for modern control simulations from the GCMs in Fig. 10. The seasonality of the ERA-40 data set is also shown. Although the modern control HadCM3L simulation has a relatively high seasonality compared with the other GCMs, especially over northern Asia, it is comparable to the ERA-40 data set. Over Antarctica, which has a large ice sheet in these control simulations, all of the GCMs have a similar seasonality, although HadCM3L is slightly higher than the other models and FAMOUS has a high seasonality over West Antarctica. This suggests that the strong HadCM3L seasonality is mainly caused by the change to early Eocene boundary conditions, although it is interesting that a similar change does not affect the other GCMs.
It is not yet clear why HadCM3L has such a strong seasonality at high latitudes under Eocene boundary conditions; other attempts at understanding why HadCM3L generates such a strong seasonality have included additional HadCM3L simulations using a dynamic vegetation model (TRIFFID) as opposed to the homogenous shrubland used by Lunt et al. (2010) (Loptson, personal communication, 2012; the study of Thorn and DeConto (2006) showed high sensitivity of the Antarctic climate to the polar vegetation cover. In addition, GENESIS simulations have been completed using the proprietary palaeo-geography used by Lunt et al. (2010) (Pollard, personal communication, 2012). These additional GENESIS simulations were also performed with a variety of vegetation types. The HadCM3L simulations with a dynamic vegetation model had an equally strong seasonality, whereas the GENESIS simulations were similar to the standard Eocene/Oligocene simulations. Further diagnostic work is needed to understand why HadCM3L has a strong seasonality under early Eocene boundary conditions, this could include experiments on the East Antarctic ice sheet (similar to the experiments of Goldner et al., 2013) and changes to ocean gateways.

Ice in the Eocene?
Based on previous modelling studies (DeConto and Pollard, 2003;Langebroek et al., 2009), and proxy records of atmospheric CO 2 concentrations (Pagani et al., 2005(Pagani et al., , 2011Pearson et al., 2009;Beerling and Royer, 2011), it is plausible that Antarctica could have been partially glaciated at times during the Eocene. This would support the argument of Miller et al. (2008) that Antarctica experienced ephemeral glaciation earlier than the EOT, based on evidence from the sea level records of Kominz et al. (2008) which show significant fluctuations in the Eocene. The offline simulations undertaken in this paper suggest that the modelled CO 2 threshold for Antarctic glaciation is highly climate-modeldependent. The composite of proxy atmospheric CO 2 records from Beerling and Royer (2011) is reproduced in Fig. 11 for data from 40-0 Ma for comparison with our model results; this includes data from a number of different proxy methods (note that the uncertainty for each of these proxies varies).
The ISM simulations using the climate from HadCM3L  and FAMOUS (Sagoo et al., 2013) do not support the early Antarctic glaciation hypothesis, however, due to the strong seasonality and low precipitation over Antarctica using HadCM3L, there is also no significant glaciation at atmospheric CO 2 concentrations lower than PIC. Given that Antarctica is glaciated today, this result seems unlikely and is also anomalous when compared with previous modelling studies (Huybrechts, 1993;De-Conto and Pollard, 2003;Langebroek et al., 2009) (2005) Beerling and Royer (2011) suggests that atmospheric CO 2 was likely between 4× and 2× PIC throughout much of the mid-to late Eocene (see Fig. 11), although there is significant uncertainty for much of the early Eocene (Beerling and Royer, 2011). With the exception of the experiments using HadCM3L and FA-MOUS, none of the experiments support totally ice-free conditions during the mid-to late Eocene based on current atmospheric CO 2 reconstructions. Although our modelling, combined with the proxy records of atmospheric CO 2 , suggests that isolated ice caps could have existed, we urge caution in assuming that this is correct. This caution is warranted given the significant inter-model disagreement. It seems plausible that a mountainous continent located over the pole would support ice caps. However, there are a number of additional factors which we have not yet fully addressed.
The opening of ocean gateways, in particular the Drake Passage, was proposed as a mechanism for the onset of Antarctic glaciation (Kennett, 1977). The modelling studies of DeConto and Pollard (2003) and Huber et al. (2004), coupled with the synchronous decrease in atmospheric CO 2 at the EOT (Pagani et al., 2011), suggest decreasing atmospheric CO 2 rather than the opening of ocean gateways as the primary mechanism for continental Antarctic glaciation. However, DeConto and Pollard (2003) suggest that the opening of ocean gateways could have lowered the CO 2 glacial threshold. This is because prior to the opening of the Drake Passage and the development of the Antarctic Circumpolar Current (ACC) there was greater oceanic meridional heat transport towards the Southern Hemisphere high latitudes. All of the early Eocene GCM simulations we have used have an open but shallow Drake Passage, resulting in partial development of the ACC. The CCSM3 and CESM1.0 experiments have a closed Tasman Gateway (Huber and Caballero, 2011;Sewall et al., 2000). It is plausible that if the experiments were repeated with a closed Drake Passage then the glacial CO 2 threshold would be lower, potentially below that suggested by the proxy records for the Eocene. In an idealized experiment where the oceanic meridional heat transport was increased by 20 % to represent a closed Drake Passage, DeConto and Pollard (2003) noted a slight lowering of the glacial CO 2 threshold to 2.3× PIC, compared with 2.8× PIC for an open Drake Passage experiment. The GCM simulations used here have a partially opened Drake Passage, so it is possible that the increase in the glacial CO 2 threshold would be less than for the DeConto and Pollard (2003) open/closed experiment, if additional GCM simulations with a closed Drake Passage were undertaken.
Proxy sea surface temperature records suggest that during past warm periods, such as the early Eocene, there was a reduced meridional temperature gradient. During the early Eocene, the high latitudes may have been significantly warmer than present-day (Bijl et al., 2009(Bijl et al., , 2010Hollis et al., 2009;Liu et al., 2009) and the low latitudes only slightly warmer than present-day (Sexton et al., 2006;Lear et al., 2008;Keating-Bitonti et al., 2011). Climate models, including those used here, have had limited success in reproducing this reduced meridional temperature gradient (Roberts et al., 2009;Winguth et al., 2010). For HadCM3L and CCSM3, the best model-data agreement requires high atmospheric CO 2 concentrations, in the range of ∼ 9-18× PIC . These atmospheric CO 2 concentrations appear high when compared with the proxy estimates. However, Huber and Caballero (2011) suggest that this increased radiative forcing is not necessarily just due to atmospheric CO 2 , but could include feedbacks from other greenhouse gases, cloud feedbacks or other unknown factors. This increased radiative forcing could be sufficient to prevent snow accumulation, for example our CCSM3_H simulation at 16× PIC is ice-free. Alternatively, the CO 2 sensitivity could be higher than that suggested by the GCMs, which is particularly low for CCSM3 and GENESIS (Huber and Caballero, 2011). Indeed, simulations using ECHAM5 require only moderate atmospheric CO 2 concentrations (2× PIC) to show reasonable agreement with the sea surface temperature data, a result that is at least in part due to the higher CO 2 sensitivity of ECHAM5 (Heinemann et al., 2009). It is interesting therefore that our ISM simulations using the climate from this ECHAM5 simulation produced a large (10.3 × 10 6 km 3 ) ice sheet. This is perhaps dependent on the large lapse rate temperature correction required from the relatively low Antarctic topography used in the ECHAM5 simulation to the ISM topography we use and the lack of elevation correction for precipitation, which could lead to artificially high precipitation rates. The FAMOUS simulation included here was part of a parameter ensemble of simulations. The ensemble member included here had the best agreement with the proxy data and showed a reduced meridional temperature gradient, although high latitude sea surface temperatures were still lower than suggested by some proxy records (Sagoo et al., 2013).
Our simulations also have relevance to other areas of debate regarding the onset of Antarctic glaciation. There is a large (∼ 1.5 ‰; Coxall et al., 2005) increase in the benthic δ 18 O record at the EOT, caused by deep-sea cooling (Liu et al., 2009;Lear et al., 2010;Pusz et al., 2011) and/or the growth of a continental-sized Antarctic ice sheet (Zachos et al., 2001;Houben et al., 2012). Recent independent estimates suggest that part of this shift was due to ∼ 1.5-5 • C of deep-sea cooling (Liu et al., 2009;Lear et al., 2010). Based on the modelling work of Langebroek et al. (2010) who suggested that the mean isotopic composition of Antarctic ice varies for a small ice sheet (−30 ‰), compared with a large ice sheet (−40 ‰), this would imply that the remainder was due to the growth of an ice sheet with a volume of ∼ 12-44 × 10 6 km 3 . Based on our simulations, the lower ice volume estimate could easily be accommodated on Antarctica, even if the continent was partially glaciated before the event. Our largest ice volume estimate is 32.5 × 10 6 km 3 using the GENESIS simulation at 2× PIC and the upper estimate of Wilson et al. (2012) for the bedrock topography. Therefore if the EOT δ 18 O shift was caused by the growth of an ice sheet of 44 × 10 6 km 3 (i.e. deep-sea cooling was 1.5 • C), it would require ice-free conditions prior to the event and potentially the additional growth of Northern Hemisphere ice sheets.

Conclusions
The inter-model comparison performed in this paper highlights that the modelled Antarctic CO 2 threshold is highly model-and model-configuration-dependent. The threshold for the growth of an intermediate ice sheet (25 m Eocene SLE) varies between 2× and 3.3× PIC (∼ 560-920 ppmv) when using the climate output from GENESIS, CCSM3_H, CESM1.0 and ECHAM5 Eocene simulations, but is not crossed when using the climate output from HadCM3L.
A large part of this disagreement is due to differences in the GCM boundary conditions, in particular the topography over the Antarctic. Some of the pre-existing Eocene GCM simulations we have used here have relatively low topography over the Antarctic. The higher-resolution ISM topographies we use are significantly more mountainous, requiring a large lapse rate correction. Because the lapse rate is a poorly constrained parameter and likely to vary spatially, the lapse rate correction is a large potential source of error. In sensitivity tests, changing the lapse rate between −6 and −8 K km −1 led to glacial threshold varying between 1.2× and 5.9× PIC (∼ 340-1650 ppmv) for CCSM3_H. Future work could involve a repeat of the GCM simulations with identical boundary conditions, which are closer to the ISM topography. We have not investigated ISM dependance in this paper and have used one surface mass balance scheme. It is possible that the CO 2 threshold could also vary if a different ISM or surface mass balance scheme were used. The offline forcing method we have adopted does not take into account feedbacks on the climate system from the growth of an ice sheet, which could affect the glacial CO 2 threshold. However, our results with GENESIS are comparable to the earlier results of DeConto and Pollard (2003) using an asynchronous coupling method.
The simulations using the HadCM3L simulations of Lunt et al. (2010) have relatively low precipitation and a very high seasonality, which results in little snow accumulation, even at low atmospheric CO 2 concentrations. This result is anomalous when compared to the results of the other GCM simulations. When using a FAMOUS simulation which had been tuned to early Eocene proxy data, no ice formed at 2× PIC (560 ppmv). The ISM simulations using the climate output from CCSM3, CESM1.0, GENESIS and ECHAM5, suggests that grounded ice could have existed earlier than the EOT, if current estimates of atmospheric CO 2 are correct. This could support evidence from sea level records (Miller et al., 2005;Kominz et al., 2008). If the Antarctic was ice-free in the Eocene it may suggest that some other mechanism prevented glaciation. For example, it is possible that stronger net radiative forcing (not necessarily due to atmospheric CO 2 ) resulted in warmer high latitudes than shown in the GCM simulations used here. Alternatively, the impact of the opening of ocean gateways and changes in ocean circulation could be greater than suggested by previous studies (DeConto and Pollard, 2003;Huber et al., 2004).