2.2.1 CESM (CESM1.2_CAM5)

CESM model description

The Community Earth System Model version 1.2 (CESM) is used, which consists of the Community Atmosphere Model 5.3 (CAM), the Community Land Model 4.0 (CLM), the Parallel Ocean Program 2 (POP), the Los Alamos sea ice model 4 (CICE), the River Transport Model (RTM), and a coupler connecting them (Hurrell et al., 2013). In comparison to previous versions of the CESM models that have been used for Eocene simulation, such as CCSM3 (Huber and Caballero, 2011; Winguth et al., 2010; Kiehl and Shields, 2013) and CESM1(CAM4) (Cramwinckel et al., 2018), CESM1.2(CAM5) represents a nearly complete overhaul of the physical parameterizations in the atmosphere model, including new schemes for radiation, boundary layer, shallow convection, cloud microphysics and macrophysics, and aerosols (Hurrell et al., 2013). The new two-moment microphysical scheme predicts both the cloud water mixing ratio and particle number concentration. The new aerosol scheme predicts the aerosol mass and number, and it is coupled with the cloud microphysics, allowing for the inclusion of aerosol indirect effects. The new boundary layer and shallow convection schemes improve the simulation of shallow clouds in the marine boundary layer. These new parameterizations in CAM5 produce a cloud simulation that agrees much better with satellite observations (Kay et al., 2012) and a larger present-day equilibrium climate sensitivity (∼ 4 ^∘C) than previous versions (∼ 3 ^∘C) (Gettelman et al., 2012). CESM1.2(CAM5) reproduces key features of the state and variability of past climates, including the mid‐Piacenzian warm period (Feng et al., 2019), the Last Glacial Maximum (Zhu et al., 2017a), Heinrich events (Zhu et al., 2017b), and the last millennium (Otto-Bliesner et al., 2015; Thibodeau et al., 2018). To make the model suitable for a paleoclimate simulation with a high CO₂ level, the model code has been slightly modified to incorporate an upgrade to the radiation code that corrects the missing diffusivity angle specifications for certain long-wave bands. As a result of the code modification, CAM5 has been re-tuned with a different relative humidity threshold for low clouds (rhminl = 0.8975, versus the default value of 0.8875). These code and parameter changes are not found to alter the present-day climate sensitivity in CESM (Zhu et al., 2019).

CESM model simulations

The CESM Eocene simulations are run at × 1, × 3, × 6, and × 9 CO₂ concentrations (Table 1). The atmosphere and land have a horizontal resolution of 1.9^∘ × 2.5^∘ (latitude × longitude) with 30 hybrid sigma-pressure levels in the atmosphere. The ocean and sea ice are on a nominal 1^∘ displaced pole Greenland grid with 60 vertical levels in the ocean. CAM5 runs with a prognostic aerosol scheme with prescribed preindustrial natural emissions that have been redistributed according to the Eocene paleogeography following the method in Heavens et al. (2012). The vegetation type from Herold et al. (2014) is prescribed in the land model with active carbon and nitrogen cycling. A modified marginal sea balancing scheme was applied for the Arctic Ocean, which removes any gain or deficit of freshwater over the Arctic Ocean and redistributes the mass evenly over the global ocean surface excluding the Arctic. This implementation conserves ocean salinity and is necessary to prevent the occurrence of negative salinity that results from high precipitation and river runoff under warm conditions. A similar balancing scheme has been included for marginal seas in all of the previously published CESM simulations (Smith et al., 2010). The ocean temperature and salinity were initialized from a previous PETM simulation using CCSM3 (Kiehl and Shields, 2013). The sea ice model was initialized from a sea-ice-free condition. All simulations have been integrated for 2000 model years, with the exception of × 1 which was run for 2600 model years.

2.2.2 COSMOS (COSMOS-landveg_r2413)

COSMOS model description

The atmosphere is represented by means of the ECHAM5 (European Centre Hamburg Model) atmosphere general circulation model (Roeckner et al., 2003). ECHAM5 is based on a spectral dynamical core and includes 19 vertical hybrid sigma-pressure levels. The series of spectral harmonics is curtailed via triangular truncation at wave number 31 (approx. $\mathrm{3.75}{}^{\circ }\times \mathrm{3.75}{}^{\circ }$). Ocean circulation and sea ice dynamics are computed by the Max-Planck-Institute for Meteorology Ocean Model (MPIOM) ocean general circulation model (Marsland et al., 2003) that is employed at 40 unequally spaced levels on a bipolar curvilinear model grid with formal resolution of $\mathrm{3.0}{}^{\circ }\times \mathrm{1.8}{}^{\circ }$ (longitude ×  latitude). The coupled ECHAM5–MPIOM model is described by Jungclaus et al. (2006). A concise description of the application of the Community Earth System Models (COSMOS) for paleoclimate studies is given by Stepanek and Lohmann (2012). The COSMOS version used here has proven to be a suitable tool for the study of the Earth's past climate, from the Holocene (Wei and Lohmann, 2012; Wei et al., 2012; Lohmann et al., 2013) and previous interglacials (Pfeiffer and Lohmann, 2016; Gierz et al., 2017) to glacial (Gong et al., 2013; Zhang et al., 2013, 2014; Abelmann et al., 2015; Zhang et al., 2017) and tectonic timescales (Knorr et al., 2011; Knorr and Lohmann, 2014; Walliser et al., 2016; Huang et al., 2017; Niezgodzki et al., 2017; Stärz et al., 2017; Walliser et al., 2017; Vahlenkamp et al., 2018; Niezgodzki et al., 2019). The standard model code of the COSMOS version COSMOS-landveg r2413 (2009) is available upon request from the Max Planck Institute for Meteorology in Hamburg (https://www.mpimet.mpg.de, last access: 10 January 2021).

COSMOS model simulations

The COSMOS simulations are carried out at × 1, × 3, and × 4 preindustrial CO₂ concentrations of 280 ppm. The ocean temperatures in the 3 × CO₂ concentration simulation were initialized with uniformly horizontal and vertical temperatures of 10 ^∘C. The initial ocean salinity was set to 34.7 psu. The simulations with 1 × and 4×  CO₂ concentrations were restarted from 3 × CO₂ after 1000 years. All simulations were run with transient orbital configurations until the model year 8000. Subsequently, they were run for 1500 years (to the model year 9500), with fixed, preindustrial orbital parameters. All simulations employ the hydrological discharge model of Hagemann and Dümenil (1998) instead of the river routing provided by Herold et al. (2014).

2.2.3 GFDL (GFDL_CM2.1)

GFDL model description

These simulations use a modified version of the Geophysical Fluid Dynamics Laboratory (GFDL) CM2.1 model (Delworth et al., 2006), similar to the late Eocene configuration in Hutchinson et al. (2018, 2019). The ocean component uses the modular ocean model (MOM) version 5.1.0, while the other components of the model are the same as in CM2.1: Atmosphere Model 2, Land Model 2 and the Sea Ice Simulator 1. The ocean and sea ice components use a horizontal resolution of 1^∘ × 1.5^∘ (latitude × longitude). A tripolar grid is used as in Hutchinson et al. (2018), with a regular latitude–longitude grid south of 65^∘ N, a transition to a bipolar Arctic grid north of 65^∘ N, and with poles over North America and Eurasia. There is no refinement of the latitudinal grid spacing in the tropics. The ocean uses 50 vertical levels with the same vertical spacing as CM2.1. The atmospheric horizontal grid resolution is 3^∘ × 3.75^∘, with 24 vertical levels, as in CM2Mc (Galbraith et al., 2010). This configuration enables relatively high-resolution ocean and coastlines, with the advantage of a faster-running atmosphere. The topography (both land and ocean) uses the 55 Ma reconstruction of Herold et al. (2014), re-gridded to the ocean and atmosphere components. Manual adjustments are made to ensure that no isolated lakes or seas exist and that any narrow ocean straits are at least two grid cells wide to ensure non-zero velocity fields. The minimum depth of ocean grid cells is 25 m; any shallow ocean grid cells are deepened to this minimum depth. In the atmosphere, the topography is smoothed using a three-point mean filter to ensure a smoother interaction with the wind field. This was introduced to remove numerical noise over the Antarctic continent, due to the convergence of meridians on the topography grid. Vegetation types are based on Herold et al. (2014), adapted to the corresponding vegetation type in CM2.1. Aerosol forcing is also adapted from Herold et al. (2014) to the model, and it is a fixed boundary condition. Ocean vertical mixing is identical to that in Hutchinson et al. (2018) – i.e. a uniform bottom-roughness-enhanced mixing with a background diffusivity of 1.0 $\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹.

GFDL model simulations

The model was initiated from idealized conditions, similar to those outlined in Lunt et al. (2017) with reduced initial temperatures: $T\left({}^{\circ }\mathrm{C}\right)=(\mathrm{5000}-z)/\mathrm{5000}\times \mathrm{25}\cos \left(\mathit{\varphi }\right)+\mathrm{10}$ if z≤5000 m, and $T\left({}^{\circ }\mathrm{C}\right)=\mathrm{10}$ if z > 5000 m; here, ϕ is latitude, and z is the depth of the ocean (positive downwards). The initial salinity was a constant of 34.7 psu. The above-mentioned initial conditions were used for the × 1, × 2, × 3, and × 4 CO₂ experiments. These simulations were initially run for 1500 years, after which the ocean temperatures were adjusted in order to accelerate the approach to equilibrium. This adjustment consisted of calculating the average temperature trend for the last 100 years at each model level below 500 m, taking a level-by-level global average of this trend, and applying a 1000-year extrapolation uniformly across the ocean at that level. This choice was based on the observation that all model levels below the mixed layer were consistently cooling at a slow rate, and the rate of temperature adjustment was consistent over a long timescale. After a further 500 years, a second adjustment using the same method was performed. After the second adjustment, all simulations were continuously integrated with no further adjustments for a further 4000 years. Thus, the simulations were run for a total of 6000 years. For the × 6 CO₂ experiment, the initial conditions described above led to transient instabilities due to overheating the surface. Thus, the × 6 experiment was instead initialized using a globally uniform temperature of 19.32 ^∘C. This represents the same global average temperature as in the other experiments and, hence, the same total ocean heat content. For the × 6 CO₂, no stepwise adjustments were made; the model was run continuously for 6000 years.

2.2.4 HadCM3 (HadCM3B_M2.1aN)

HadCM3 model description

The Hadley Centre Climate Model (HadCM3) simulations are carried out with the HadCM3B-M2.1aN version of the model, as described in detail in Valdes et al. (2017). Equations are solved on a Cartesian grid with horizontal resolutions of 3.75^∘ × 2.5^∘ in the atmosphere and 1.25^∘ × 1.25^∘ in the ocean with 19 and 20 vertical levels respectively. A few changes are made to the version described in Valdes et al. (2017) to make it suitable for deep-time paleoclimate modelling: (a) a salinity flux correction is applied to the global ocean (at all model depths) in order to conserve salinity; (b) the various modern-specific parameterizations in the ocean model are turned off, such as those associated with the Mediterranean and Hudson Bay outflow and the North Atlantic mixing; and (c) a prognostic 1D ozone scheme is used instead of a fixed vertical profile of ozone. The standard configuration uses a prescribed ozone climatology which is a function of latitude, height, and month of the year that does not change with climate and can become numerically unstable at high CO₂ levels. The prognostic ozone scheme uses the diagnosed model tropopause height to assign three distinct ozone concentrations for the troposphere, tropopause, and stratosphere (2.0 $\times {\mathrm{10}}^{-\mathrm{8}}$, 2.0 $\times {\mathrm{10}}^{-\mathrm{7}}$, and 5.5 $\times {\mathrm{10}}^{-\mathrm{6}}$, in mass mixing ratio, respectively). This allows for a dynamic update of the 1D ozone field in response to the thermally driven vertical expansion of the troposphere. Absolute values for the three levels are chosen to minimize the effects on the global mean and overall tropospheric temperature changes compared with the standard 2D climatology. Concentrations at the uppermost model level are fixed to the higher stratospheric value to constrain the lower bound of total stratospheric ozone. Significant differences to the standard configuration are limited to the stratospheric meridional temperature gradient and zonal winds and are related to the missing latitudinal variations in the 1D field. Although HadCM3 has been used previously to simulate the Pliocene (e.g. Lunt et al., 2008, 2010a), the presented simulations represent the first published application of HadCM3 to pre-Pliocene boundary conditions. However, the lower-resolution HadCM3L model has been previously used to simulate a range of pre-Quaternary climates (e.g. Lunt et al., 2016; Farnsworth et al., 2019a, b).

HadCM3 model simulations

The HadCM3 simulations are carried out at × 1, ×2, and × 3 CO₂ concentrations. Several ocean gateways were artificially widened to allow unrestricted throughflow, and maximum water depths in parts of the Arctic Ocean were reduced. The ocean temperatures were initialized from the final state of Eocene model simulations using HadCM3L. The HadCM3L simulations were set up identically to the corresponding HadCM3 simulations, but with a lower ocean resolution (3.75^∘ × 2.5^∘ as opposed to 1.25^∘ × 1.25^∘). The HadCM3L simulations were initialized from a similar idealized temperature and salinity state as described in Lunt et al. (2017) but with a function that scales with cos ²(lat) rather than cos (lat) and overall reduced initial temperatures to ensure numerical stability in tropical regions. Ocean temperatures below 600 m were set to constant values of 4, 8, and 10 ^∘C (at × 1, ×2, and × 3 CO₂ respectively) based on results from previous Ypresian simulations. The HadCM3 simulations were branched off from the respective HadCM3L integrations after 4400 to 4900 years of spin up and run for a further 2950 years. The initial 50 years of all HadCM3 runs used the simplified vertical diffusion scheme from HadCM3L (Valdes et al., 2017) to reduce numerical problems caused by the changed horizontal ocean resolution. The remaining years of the runs use the standard HadCM3 diffusion scheme (Valdes et al., 2017).

2.2.5 INMCM (INM-CM4-8)

INMCM model description

The INMCM simulations are carried out with the INM-CM48 (INM-CM4-8) version of the model, as described in Volodin et al. (2018). The INM-CM4-8 climate model has a horizontal resolution of 2^∘ × 1.5^∘ in the atmosphere; a total of 17 vertical sigma levels up to a value of 0.01 (about 30 km) are used for the Eocene experiment, and 21 levels are used for the preindustrial experiment. The equations of the atmosphere dynamics are solved by finite difference methods. The parameterizations of physical processes correspond to the INM-CM5 model (Volodin et al., 2017). Parameterization of condensation and cloud formation follows Tiedtke (1993). Cloud water is a prognostic variable. Parameterization of the cloud fraction follows Smagorinsky (1963); cloud fraction is a diagnostic variable. The surface, soil, and vegetation scheme follows Volodin and Lykossov (1998). The evolution of the equations for temperature, soil water, and soil ice are solved at 23 levels from the surface to 10 m depth. The fractional area of 13 types of potential vegetation is specified. Actual vegetation and the leaf area index (LAI) are calculated according to the soil water content in the root zone and soil temperature. This model also contains a carbon cycle and an aerosol scheme (Volodin and Kostrykin, 2016), taking the direct impact of aerosols on radiation into account, as well as the first indirect effect (the influence of aerosols on the condensation rate). The concentration of 10 types of aerosol and their radiative properties are calculated interactively. In the ocean component, the resolution of the INM-CM4-8 model is 1.0^∘ × 0.5^∘ (longitude × latitude) and has 40 sigma levels vertically. Finite difference equations are solved on a generalized spherical C-grid with the North Pole shifted to Siberia; the South Pole is in the same place as the geographical pole.

INMCM model simulations

The INM-CM4-8 Eocene simulation is carried out at a × 6 CO₂ concentration. The INM-CM4-8 simulation was initialized from a similar idealized temperature and salinity state as described in Lunt et al. (2017), but the initial formula for the ocean temperature is modified as follows: $T=\left(\right(\mathrm{5000}-z)/\mathrm{5000}\times \mathrm{20}\cos (\mathit{\varphi }\left)\right)+\mathrm{15}$, thereby reducing the initial temperatures to ensure numerical stability in tropical regions. The 27 biomes were converted into the 13 model types of vegetation. The duration for the Eocene simulation is 1150 years. Output data are averaged over the years from 1051 to 1150.

2.2.6 IPSL (IPSLCM5A2)

IPSL model description

The Institut Pierre Simon Laplace (IPSL) simulations are performed with the IPSL-CM5A2 Earth system model (Sepulchre et al., 2020). IPSL-CM5A2 is based on the CMIP5-generation previous IPSL Earth system model IPSL-CM5A (Dufresne et al., 2013) but includes new revisions of each components, a re-tuning of global temperature, and technical improvements to increase computing efficiency. It consists of the LMDZ5 (Laboratoire de Meétéorologie Dynamique Zoom) atmosphere model, the Organising Carbon and Hydrology In Dynamic Ecosystems (ORCHIDEE) land surface and vegetation model and the Nucleus for European Modeling of the Ocean (NEMOv3.6) ocean model, which includes the LIM2 sea ice model and the Pelagic Interactions Scheme for Carbon and Ecosystem Studies (PISCES-v2) biogeochemical model. LMDZ5 and ORCHIDEE run at a horizontal resolution of 1.9^∘ × 2.5^∘ (latitude × longitude) with 39 hybrid sigma-pressure levels in the atmosphere. NEMO runs on a tripolar grid at a nominal resolution of 2^∘, enhanced up to 0.5^∘ at the Equator, with 31 vertical levels in the ocean. The performances and evaluation of IPSL-CM5A2 on preindustrial and historical climates are fully described in Sepulchre et al. (2020). Sepulchre et al. (2020) also provide a description of the technical changes that were implemented in IPSL-CM5A2 to carry out deep-time paleoclimate simulations. In particular, the tripolar mesh grid on which NEMO runs has been modified to ensure that there are no singularity points within the ocean domain. Modern parameterizations of water outflows across specific straits, such as the Gibraltar or Red Sea straits, are also turned off.

IPSL model simulations

The IPSL simulations are run at × 1.5 and × 3 CO₂ concentrations. The bathymetry is obtained from the Herold et al. (2014) dataset, with additional manual corrections in some locations, for instance in the West African region, to maintain sufficiently large oceanic straits. Modern boundary conditions of NEMO include forcings of the dissipation associated with internal wave energy for the M2 and K1 tidal components (de Lavergne et al., 2019). The parameterization follows Simmons et al. (2004) with refinements in the modern Indonesian throughflow (ITF) region according to Koch-Larrouy et al. (2007). To create an early Eocene tidal dissipation forcing, the Herold et al. (2014) M2 tidal field (obtained from the tidal model simulations of Green and Huber, 2013) is directly interpolated onto the NEMO grid using bilinear interpolation. In the absence of any estimation for the early Eocene, the K1 tidal field is prescribed as zero. In addition, the parameterization of Koch-Larrouy et al. (2007) is not used here because the ITF does not exist in the early Eocene. The geothermal heating distribution is created from the 55 Ma global crustal age distribution of Müller et al. (2008), on which the age–heat flow relationship of the Stein and Stein (1992) model is applied: $q\left(t\right)=\mathrm{510}\times t^{-\mathrm{1}/\mathrm{2}}$ if t ≤ 55 Ma, and $q\left(t\right)=\mathrm{48}+\mathrm{96}\exp (-\mathrm{0.0278}\times t)$ if t > 55 Ma. In regions of subducted seafloor where age information is not available, the minimal heat flow value is prescribed, which is derived from a known crustal age. The resulting 1^∘ × 1^∘ field is then bilinearly interpolated onto the NEMO grid. It must be noted that the Stein and Stein (1992) parameterization becomes singular for young crustal ages and yields unrealistically large heat flow values. Following Emile-Geay and Madec (2009), an upper limit of 400 mW m⁻² is set for heat flow values after the interpolation procedure. Salinity is initialized as globally constant to a value of 34.7 psu following Lunt et al. (2017). The initialization of the model with the proposed DeepMIP temperature distribution (Lunt et al., 2017) led to severe instabilities in the model during the spin-up phase. Thus, the initial temperature distribution has been modified as follows: $T\left({}^{\circ }\mathrm{C}\right)=(\mathrm{1000}-z)/\mathrm{1000}\times \mathrm{25}\cos \left(\mathit{\varphi }\right)+\mathrm{10}$ if z ≤ 1000 m, and $T\left({}^{\circ }\mathrm{C}\right)=\mathrm{10}$ if z > 1000 m; here, ϕ is the latitude, and z is the depth of the ocean (metres below surface). This new equation gives an initial globally constant temperature of 10 ^∘C below 1000 m and a zonally symmetric distribution above 1000 m, reaching surface values of 35 ^∘C at the Equator and 10 ^∘C at the poles. This corresponds to a 5 ^∘C surface temperature reduction compared with DeepMIP guidelines (Lunt et al., 2017). No sea ice is prescribed at the beginning of the simulations. In IPSL-CM5A2, the NEMO ocean model is inherently composed of the PISCES biogeochemical model. Biogeochemical cycles and marine biology are directly forced by dynamical variables of the physical ocean and may affect the ocean physics via its influence on chlorophyll production, which modulates light penetration in the ocean. However, because this feedback does not much affect the ocean state (Kageyama et al., 2013) and because the early Eocene mean ocean colour is unknown, a constant chlorophyll value of 0.05 g Chl L⁻¹ is prescribed for the computation of light penetration in the ocean. As a consequence, marine biogeochemical cycles and biology do not alter the dynamics of the ocean; as such, the biogeochemical initial and boundary conditions have been kept at modern values. The topographic field is created from the Herold et al. (2014) topographic dataset; LMDZ includes a sub-grid-scale orographic drag parameterization that requires high-resolution surface orography (Lott and Miller, 1997; Lott, 1999). A similar procedure is applied for the standard deviation of orography provided by Herold et al. (2014). Aerosol distributions are left identical to preindustrial values. The × 3 simulation is initialized from rest and run for 4000 years. The × 1.5 simulation is branched from the model year 1500 of the × 3 simulation and run for 4000 years. The × 1.5 and × 3 simulations are identical to those presented in Zhang et al. (2020).

2.2.7 MIROC (MIROC4m)

MIROC model description

The version of the Model for Interdisciplinary Research on Climate (MIROC) used here is MIROC4m, a mid-resolution model composed of atmosphere, land, river, sea ice, and ocean components. Full documentation of the model can be found in K-1 model developers (2004) and is summarized in Chan et al. (2011). The atmosphere has a horizontal resolution of T42 and 20 vertical sigma levels. Details of the land surface model, Minimal Advanced Treatments of Surface Interaction and Runoff (MATSIRO), can be found in Takata et al. (2003). The ocean component is basically version 3.4 of the CCSR (Center for Climate System Research) Ocean Component Model (COCO); the reader is referred to Hasumi (2000) for details. The horizontal resolution is set to 256 × 196 (longitude × latitude), with a higher resolution in the tropics, and the vertical resolution is set to 44 levels, with the top 8 levels in sigma coordinates. Present-day bathymetry is derived from Earth topography five minute grid (ETOPO5) data. For present-day experiments, areas of water such as Hudson Bay and the Mediterranean Sea are represented as isolated basins. As such, ocean salinity and heat are artificially exchanged with the open ocean through a two-way linear damping. This damping and all isolated basins and lakes are removed in the DeepMIP simulation.

MIROC model simulations

Out of the three standard DeepMIP simulations, MIROC is used with a × 3 CO₂ concentration only and run for 5000 model years. The atmosphere is initialized from a previous experiment without ice sheets and with a × 2 CO₂ concentration. For the initial ocean state, salinity is set to a constant value of 34.7 psu, as recommended in Lunt et al. (2017). However, the ocean temperatures are 15 ^∘C cooler than those recommended – i.e. $T\left({}^{\circ }\mathrm{C}\right)=(\mathrm{5000}-z)/\mathrm{5000}\times \mathrm{25}\cos \left(\mathit{\varphi }\right)$ if z ≤ 5000 m, and $T\left({}^{\circ }\mathrm{C}\right)=\mathrm{0}$ if z > 5000 m. Previous MIROC experiments similar to this × 3 CO₂ DeepMIP simulation show that this initialization should be much closer to the final climate state. Simulations were also carried out at × 1 and × 2, but they are not discussed in this paper.

2.2.8 NorESM (NorESM1_F)

NorESM model description

The Norwegian Earth System Model (NorESM) simulations are carried out with the NorESM1-F version of the model, which is described in detail in Guo et al. (2019). The NorESM version that contributes to CMIP5 is NorESM1-M. It has a ∼ 2^∘ resolution atmosphere and land configuration and a nominal 1^∘ ocean and sea ice configuration. In NorESM1-F, the same atmosphere–land grid is used as in NorESM1-M (CMIP5 version), whereas a tripolar grid is used for the ocean–sea ice components in NorESM1-F, instead of the bipolar grid in NorESM1-M. The tripolar grid is also used in the CMIP6 version of NorESM (NorESM2). NorESM1-F runs about 2.5 times faster than NorESM1-M. For the preindustrial simulation, NorESM-F has a more realistic Atlantic Meridional Overturning Circulation than NorESM1-M.

NorESM model simulations

The NorESM simulations are carried out at × 2 and × 4 CO₂ concentrations. The ocean temperatures were initialized from the × 2 CO₂ Eocene simulations with the lower-resolution NorESM-L model (Zhang et al., 2012). The ocean salinity was initialized with constant values of 25.5 psu in the Arctic and 34.5 psu elsewhere. From the initial conditions, the × 2 CO₂ experiment was run for 2100 years in total. The × 4 CO₂ was branched from the end of the 100th year of the × 2 CO₂ experiment and was run for 2000 years. The results from the last 100 years were used in the study. Note that the NorESM simulations were carried out with the Baatsen et al. (2016) paleogeography (based on a paleomagnetic reference frame), not the Herold et al. (2014) paleogeography (based on a mantle reference frame), in contrast to the other simulations described in this paper.

3.4.1 Proxy datasets

For point-to-point model–data comparisons, we use the SST and surface air temperature (SAT) datasets for the EECO compiled by (Hollis et al., 2019). Following their recommendation, we exclude δ¹⁸O-derived SST estimates from recrystallized planktonic foraminifera because these estimates are generally significantly cooler than estimates derived from the δ¹⁸O value of well-preserved foraminifera, foraminiferal Mg ∕ Ca ratios, and clumped isotope values from larger benthic foraminifera, due to diagenetic effects.

In terms of global metrics, we make use of the global mean near-surface air temperature (GSAT) estimate for the EECO from Inglis et al. (2020), which is based on the Hollis et al. (2019) temperature dataset and also excludes SST estimates from recrystallized foraminifera. The vertical dimensions of the grey filled boxes in Fig. 1a show the 33 % to 66 % and the 10 % to 90 % confidence intervals of this GSAT estimate. For the global mean SST, we again use the GSAT estimate of Inglis et al. (2020) but convert this to global mean SST using a linear function derived from the mean land–sea temperature contrast in all the model simulations shown in Fig. 1:

This SST estimate forms the horizontal dimensions of the grey filled boxes in Fig. 1b.

We make use of SST gradient estimates from Cramwinckel et al. (2018), Evans et al. (2018), and Zhu et al. (2019). Cramwinckel et al. (2018) define a meridional temperature gradient metric as the difference between the tropical mean SST (derived from TEX₈₆) and deep-ocean temperatures (derived from δ¹⁸O), assuming that deep-ocean temperatures are an approximation of high-latitude SSTs. In this way, they reconstruct a metric at 50 Ma of about 20–22 ^∘C (their Fig. 3b). Evans et al. (2018) employ a similar approach but using tropical SSTs and deep-ocean temperatures from Mg ∕ Ca between 48 and 56 Ma; they find a reduction in their metric from the modern to early Eocene of 22 % to 42 %. Given a modern gradient from HadISST (our Fig. 1b) of about 26 ^∘C, this gives an early Eocene estimate of 17.7 ^∘C (20.3 to 15.1 ^∘C). Using a similar approach to Evans et al. (2018), Zhu et al. (2019) find a reduction in their metric of about 20 % to 45 % (their Fig. 1b), giving an estimate of 21 to 14 ^∘C for the early Eocene. Given the differences in methodologies for deriving these estimates, the relatively wide time window of the “early Eocene” for two of the studies, and the differences between the proxy metrics and our modelled metric (defined here as the average SST equatorwards of ±30^∘ minus the average SST polewards of ±60^∘), we take the outer ranges of the three studies as our overall proxy estimate, giving a range from 14 to 22 ^∘C. This range in the meridional temperature gradient estimate forms the vertical edge of the grey filled boxes in Fig. 1b. However, the use of benthic temperatures to approximate high-latitude annual mean surface temperature may result in biases due to the seasonality of deep-water formation (Evans et al., 2018), and we note that a detailed assessment of the meridional temperature gradient implied by proxies, similar to that carried out for global mean temperature by Inglis et al. (2020), would be beneficial for future model–data comparisons.

For CO₂, Anagnostou et al. (2020) give two estimates of EECO CO₂ based on two different calibrations, resulting in 95 % confidence intervals from 1170 to 1830 ppmv and from 1540 to 2490 ppmv. The uppermost and lowermost bounds of these two estimates are close to the × 4 (1120 ppmv) and × 9 (2520 ppmv) simulations; as such, these form the horizontal edges of the light grey box in Fig. 1. A normal distribution in absolute CO₂ would give a corresponding 68 % confidence interval from 1470 to 2170 ppmv, and this forms the horizontal edges of the light grey box in Fig. 1.

Overall, for the purposes of describing the model–data consistency, we use “consistent” to describe a model that sits within the light grey boxes of Fig. 1a and b, and we use “very consistent” to describe a model that sits within the dark grey boxes of Fig. 1a and b.

3.4.2 Comparison with global metrics

For the DeepMIP models, only those that carried out simulations at × 4 CO₂, × 6 CO₂, or × 9 CO₂ are consistent with the CO₂ proxies (i.e. CESM, COSMOS, GFDL, INMCM, and NorESM), and only those that carried out simulations at × 6 CO₂ are very consistent with the CO₂ proxies (i.e. CESM, GFDL, and INMCM). All of these simulations are also consistent with the GSAT proxies, but only COSMOS and GFDL at × 4 are very consistent with the GSAT proxies (inspection of Fig. 1 indicates that CESM would also be very consistent with the GSAT proxies if there was a simulation at × 4). No simulations are very consistent with both the CO₂ and GSAT proxies. Only CESM at × 6, GFDL at × 4 and × 6, and NorESM at × 4 are consistent with the CO₂, GSAT, and meridional temperature gradient proxies.

Of the pre-DeepMIP simulations in Fig. 1a and b, CCSM_H and CCSM_W at × 8 are also consistent with all the proxy constraints, and CCSM_K is additionally very consistent with GSAT. However, as discussed in Sect. 1, CCSM_K includes somewhat arbitrary modifications to cloud parameters that are designed to enable the model to better fit the Eocene observations; as such, in contrast to the DeepMIP simulations, CCSM_K cannot be considered entirely independent of the temperature data.

Some quantitative metrics for the simulations are presented in Table S2. In this case, the metrics are given for the set of simulations that were carried out at CO₂ concentrations consistent with the proxies.

3.4.3 Comparison with specific locations

The limited range of CO₂ concentrations explored by some models, coupled with the relatively large uncertainties in EECO CO₂ from proxies, means that a model–data comparison of individual model simulations with the site-by-site proxy data can be misleading. Therefore, we only conduct a detailed model–data comparison for models that provide simulations carried out under more than one CO₂ concentration. For these models (CESM, COSMOS, GFDL, HadCM3, IPSL, and NorESM), we apply a global mean scaling factor to the simulated SST and GSAT such that the modelled global means best fit the global mean proxy data. We then compare the spatial patterns in the scaled model outputs with the spatial patterns in the site-by-site proxies. We provide a quantitative metric for the model–data fit and compare this with some idealized temperature distributions to put these metrics in context.

We scale the models by assuming that the spatial pattern of temperature change scales linearly with global mean temperature change and by interpolating or extrapolating to a global mean equal to the estimate from Inglis et al. (2020), i.e. 27 ^∘C for near-surface air temperature, or an equivalent global mean SST given by Eq. (4). This process gives a scaling factor, s, that can be used to create a spatial field of implied temperature, Tⁱ, which is consistent with this proxy-based global mean temperature, 〈T^p〉:

where ^LT and ^HT are the spatial fields of the two model simulations that have global means closest to 〈T^p〉, and 〈.〉 denotes a global mean quantity. This also allows us to calculate an inferred CO₂ concentration, ${\mathrm{CO}}_{\mathrm{2}}^i$:

where ^HCO₂ and ^LCO₂ are the CO₂ concentrations that correspond to the two simulations in Eq. (5). For surface air temperature, this process is equivalent to interpolating or extrapolating the straight lines in Fig. 1a to identify the CO₂ concentration that corresponds to 〈T^p〉.

For CESM and GFDL, the scaling is found by interpolation (s<1.0) because there are simulations that are warmer than 〈T^p〉. For those models where the scaling extrapolates beyond the model simulations (i.e. s>1.0; COSMOS, HadCM3, IPSL, and NorESM), care must be taken due to the assumption of linearity. For HadCM3, IPSL, and COSMOS, this assumption is probably well justified (s=1.51, 1.37, and 1.05 respectively for the surface air temperature scaling), but for NorESM this assumption is probably poorly justified (s=2.02).

For HadCM3, GFDL, IPSL, CESM, COSMOS, and NorESM, the inferred CO₂ concentrations for the surface air temperature scaling are 1030, 1050, 1080, 1130, 1140, and 2270 ppmv respectively. For CESM, COSMOS, and NorESM, these are consistent with the CO₂ proxy estimates of 1120–2520 ppmv (Sect. 3.4.1); the other models have a slightly lower inferred CO₂ concentrations than the proxies indicate. All of these inferred CO₂ values are below the concentration at which the respective models are known to crash (see Sect. 3.1). When the same method is applied to the Lunt et al. (2012) simulations, the inferred CO₂ are all higher than the proxy estimates (2640 ppmv for HadCM3L, 3300 ppmv for CCSM3_H, and 6210 ppmv for CCSM_W). Figure 1 indicates that these relatively cool Lunt et al. (2012) simulations are related to a relatively low climate sensitivity for CCSM3_H and CCSM3_W and to a relatively low response to non-CO₂ forcing for HadCM3L and CCSM3_W.

Figure 5Modelled SST anomalies for the Eocene, relative to the zonal mean of the associated preindustrial simulation. The variable plotted is ${\mathrm{SST}}_{\mathrm{e}}^{\mathrm{m}}-\stackrel{\mathrm{‾}}{{\mathrm{SST}}_{\mathrm{e}}^{\mathrm{m}}}$ in the notation of Lunt et al. (2012). The Eocene simulations have been scaled using a global tuning factor, as described in the text, so that they best fit the global mean SST data inferred from Inglis et al. (2020) (see Eqs. 5 and 4). As such, only models that carried out simulations with more then one CO₂ concentration are shown: (a) CESM, (b) COSMOS, (c) GFDL, (d) HadCM3, (e) IPSL, and (f) NorESM. Also shown are the proxy SST estimates from Hollis et al. (2019) for the EECO, excluding those sites that they identified as being affected by diagenesis.

The scaled SST anomalies, relative to the zonal mean of the preindustrial modelled SST, for CESM, COSMOS, GFDL, HadCM3, IPSL, and NorESM, along with the proxy SST data from Hollis et al. (2019) (relative to the zonal mean of preindustrial observations), are shown in Fig. 5. In general, the models agree reasonably well with the tropical and mid-latitude SST data, but there is a large model–data inconsistency in the southwest Pacific sites around New Zealand and south of Australia, where the modelled anomalies are colder than proxy estimates by 5–10 ^∘C. The reader is also referred to Fig. S7 for the modelled absolute SSTs and absolute SST proxy data.

The RMS skill score of the scaled absolute simulations, relative to the SST proxies, σ_s (^∘C), is 7.0 for NorESM, 9.6 for GFDL, 9.7 for CESM, 10.5 for HadCM3, 10.7 for IPSL, and 12.0 for COSMOS. Note that the NorESM score is not directly comparable to the others because the NorESM simulation, and the proxy data it is compared with, are on a paleomagnetic reference frame rather than a mantle reference frame (Fig. 5f). For comparison, the GFDL skill score is 7.3 when calculated on the paleomagnetic reference frame. Note that we calculate all RMS scores from a specific point-to-point comparison of models and data, not from zonal means.

To put these numbers in context, we also calculate the same skill score for some idealized temperature distributions (on the mantle reference frame), expressed as anomalies relative to the zonal mean of the preindustrial observations. This approach is similar to that used by Hargreaves et al. (2013) in the context of Quaternary model–data comparisons. Our idealized temperature distributions are (i) a constant value of zero (i.e. no change from the zonal mean of the preindustrial simulation), (ii) a non-zero constant value (C), and (iii) a function $f\left(\mathit{\varphi }\right)=A+B(\mathrm{1}-\cos \mathrm{2}\mathit{\varphi })$. For the constant value, C, we choose a value that is equal to an estimate of global mean SST change from the proxies. This estimate of SST change is scaled from the proxy-based estimate of GSAT, 〈T^p〉=27 ^∘C, using the scaling in Eq. (4), minus the preindustrial global mean SST. For the function f(ϕ), we choose A and B such that the global mean is equal to C, and the polar amplification metric, defined as the average SST equatorwards of ±30^∘ minus the average SST polewards of ±60^∘, is equal to our central estimate, i.e. 18 ^∘C minus the preindustrial value (see Sect. 3.4.1).

Figure 6(a–d) Zonal mean SST (solid lines) and near-surface air temperature (dashed line) anomalies, relative to the zonal mean of the associated preindustrial simulation, for the scaled version of the (a) CESM, (b) COSMOS, (c) GFDL, (d) HadCM3, (e) IPSL, and (f) NorESM models. Also shown are the EECO SSTs and error bars from Hollis et al. (2019), expressed as a difference relative to the zonal mean of the preindustrial observations. Panels (g)–(i) are the same as panels (a)–(f) but instead of a model we show idealized temperature distributions of (g) 0, (h) C, and (i) $A+B(\mathrm{1}-\cos \mathrm{2}\mathit{\varphi })$. All plots also show the proxy SST estimates from Hollis et al. (2019) for the EECO, excluding those sites that they identified as being affected by diagenesis (black circles with uncertainty bars). The modelled and idealized SSTs at the specific location of the proxies (red squares) are also displayed.

Download

These idealized functions are shown in Fig. 6g–i as zonal means, along with the scaled DeepMIP models (Fig. 6a–f), and are all expressed as anomalies relative to the zonal mean of the preindustrial simulation. The global anomaly of zero relative to the zonal mean of the preindustrial simulation is associated with an RMS skill score σ_s=20.1, the global mean constant temperature anomaly (C) is associated with σ_s=11.6, and the f(ϕ) temperature profile is associated with σ_s=9.0 (σ_s=7.5 on the paleomagnetic reference frame). This means that all of the models apart from COSMOS can be considered as having some skill in capturing the first-order patterns of SST change (because the skill score of those models is better than that of the global constant), but only NorESM has skill in capturing the second-order, more regional temperature patterns (because the skill score of the other models is worse than that of the f(ϕ) distribution when calculated on the appropriate reference frame). However, the performance of the scaled NorESM simulations should be viewed with some caution because of its relatively high scaling factor (s, in Eq. 5).

So far, this analysis has focused on SSTs, but we also carry out a comparison with terrestrial near-surface air temperature data (SAT), even though it is generally less well constrained in age than SSTs and, as such, likely represents a wider range of climate states. The absolute SAT model–data comparison for each DeepMIP simulation is shown in Fig. S8. For those models that carried out more than one CO₂ simulation (CESM, COSMOS, GFDL, HadCM3, IPSL, and NorESM), Figs. S9 and S10 show the SATs from the scaled models in comparison with terrestrial proxy data.

These figures show that both models and SAT terrestrial proxies show a similar amount of polar amplification. In particular, the southwest Pacific site SATs are better simulated in the models than the SSTs; the RMS skill score decreases in the southwest Pacific by 30 % on average for the SATs compared with the SSTs across the ensemble. This implies that there may be inconsistency between marine and terrestrial temperatures in either the proxies or models in this region. This discrepancy could be related to a potential summer bias in mid- and high-latitude SST proxies (Hollis et al., 2012; Davies et al., 2019). An alternative hypothesis is that the discrepancy is related to Red Sea-like features of GDGT (glycerol dialkyl glycerol tetraethers) distributions in high SST samples from the southwest Pacific and Wilkes Land that appear to amplify proxy SSTs where isoGDGT_RS>30 (Inglis et al., 2015), which is an idea supported by recent work in the context of the Cretaceous (Steinig et al., 2020). However, the discrepancy may also be caused by physical processes that are not captured by any of the models.