Climate model and proxy data constraints on ocean warming across the Paleocene – Eocene Thermal Maximum

Constraining the greenhouse gas forcing, climatic warming and estimates of climate sensitivity across ancient large transient warming events is a major challenge to the palaeoclimate research community. Here we provide a new compilation and synthesis of the available marine proxy temperature data across the largest of these hyperthermals, the Paleocene – Eocene Thermal Maximum (PETM). This includes the application of consistent temperature calibrations to all data, including the most recent set of calibrations for archaeal lipid-derived palaeothermometry. This compilation provides the basis for an informed discussion of the likely range of PETM warming, the biases present in the existing record and an initial assessment of the geographical pattern of PETM ocean warming. To aid interpretation of the geographic variability of the proxy-derived estimates of PETMwarming, we present a comparison ofthis data with the patterns of warmingproduced byhigh p CO 2 sim- ulations of Eocene climates using the Hadley Centre atmosphere-ocean general circulation model (AOGCM) HadCM3L. On the basis of this comparison and taking into account the patterns of intermediate-water warming we estimate that the global mean surface temperature anomaly for the PETM is within the range of 4 to 5 °C.


Introduction
The Paleocene-Eocene Thermal Maximum (PETM; ca 55 Ma) is a relatively short-lived (~200 kyr) transient climatic event characterised by a large negative excursion in stable carbon isotope ratios, a prominent shoaling of the global ocean carbonate compensation depth (CCD) and a widespread positive temperature anomaly Sluijs et al., 2007a). The combination of global warming and the release of large amounts of 13 C-depleted carbon (circa 2000-6000 PgC) (Dickens, 2003;Panchuk et al., 2008aPanchuk et al., , 2008bZeebe et al., 2009;Cui et al., 2011) to the ocean-atmosphere system during the PETM has encouraged analogies to be drawn with modern anthropogenic climate change. As a result the PETM has been the focus of considerable research effort since it was first recognised over twenty years ago (Thomas, 1989;Kennett and Stott, 1991), with the characteristic carbon isotope excursion and associated lithologic and biotic changes being documented at numerous localities across the globe and in disparate environments (see reviews of Sluijs et al., 2007a;McInerney and Wing, 2011).
The PETM is clearly a planetary-scale perturbation of the Earth's carbon cycle, biosphere and climate. There remains, however, considerable uncertainty over the size (Panchuk et al., 2008a(Panchuk et al., , 2008bZeebe et al., 2009;Cui et al., 2011Cui et al., , 2012 and source (Dickens, 2011;Lunt et al., 2011;DeConto et al., 2012) of the PETM carbon release. Consequently the exact magnitude of greenhouse gas forcing arising from this carbon release is poorly constrained. Whereas the extent of the PETM carbon cycle perturbation has received considerable attention (e.g. Dickens, 2003;Zachos et al., 2005), there has been, to date, no compilation and thorough reassessment of the available PETM temperature proxy data. When quoted the PETM sea surface temperature anomaly is usually stated as a range, typically "4-8°C" (Thomas et al., 2002), or more recently "5-9°C" (Zeebe et al., 2009). The robustness, interpretation, quantity and geographic spread of the underlying data are, however, rarely commented upon. Even recent "PETM" climate modelling studies have used a range of published PETM palaeotemperature proxy data as if they were relatively uncomplicated, unbiased and without inherent uncertainty Winguth et al., 2010).
The purpose of this study is to present a full synthesis of all the available PETM marine geochemical temperature proxy data, with the subsequent derivation of site-specific estimates of the PETM temperature anomaly. Given that the palaeotemperature proxy data are both spatially limited and subject to varying degrees of uncertainty and potential biases, we use the results of an early Eocene fully-coupled Atmosphere-Ocean General Circulation Model (AOGCM) to aid in the interpretation of these data and provide a means of estimating the likely global mean PETM temperature anomaly.
The compilation of PETM proxy data is derived from near-surface dwelling planktic foraminifera, sea surface-derived organic lipid biomarkers (TEX 86 ) and benthic foraminifera, mostly from deep-sea drilling cores recovered in intermediate to abyssal water depths (~1000-3000 m). We treat these data, collected over the last 20 years, with a consistent methodology for the identification of pre-PETM and PETM intervals based on carbon isotope stratigraphies, the application of uniform palaeotemperature calibrations and a more robust and transparent assessment of uncertainties associated with temporal averaging and calibrations. We also discuss major uncertainties in these palaeotemperature estimates that cannot yet be reliably quantified, including the preservation-state of foraminiferal carbonate and the application of TEX 86 calibrations from the modern oceans to ancient warmclimates states.
One of the immediate conclusions of this study is that a reliable estimate of PETM near sea surface temperature (SST) change has been established at very few locations. We identify only eight sites with reasonable PETM records of sea surface temperatures: the New Jersey margin sites Bass River and Wilson Lake, DSDP Site 527, ODP Sites 690, 865, 1172, 1209 and the Arctic IODP Site M0004Aabbreviated to ACEX for Arctic Coring Expedition from here on (Fig. 1). Although these sites are geographically disparate, they provide no data from the Indian or Tethys Oceans and most of the tropics. However, within the limitations of the available data some constraint on the plausible range of the global average PETM temperature anomaly is possible. By identifying consistencies in the spatial response of GCM climate simulations and the bestavailable estimates of PETM warming the model can be used as a 'grand interpolator' to produce an estimate of the global mean temperature change. We also hope that a clear presentation of the proxy data and its limitations will justify and encourage future efforts to gather reliable PETM palaeotemperature proxy data across both marine and terrestrial environments and from geographically disparate localities.

Identification of pre-and peak-PETM intervals
Original published marine palaeotemperature proxy data were collated from locations where palaeotemperature proxy data were of sufficient number and quality to resolve at least part of the PETM warming and where these data could be coupled to a carbon isotope stratigraphy that resolved an unambiguous PETM CIE. For each location, we have established two stratigraphic intervals, largely based on carbon isotope stratigraphy, representing two relatively stable climatic states oñ 10 kyr timescales: one immediately prior to the PETM ("pre-PETM") and one at the "peak" of the PETM perturbation ("peak-PETM") . The shape of the PETM carbon isotope excursion (CIE) can vary significantly between sections and the substrate measured (McInerney and Wing, 2011). It is dependent on differential carbonate dissolution, sediment mixing, physiological change of the substrate producer and regional offsets in δ 13 C of the carbon source (Sluijs and Dickens, 2012). Given the great magnitude of the isotopic perturbation, however, in most instances the CIE involves a large and rapid negative excursion in δ 13 C, followed by a short,~80 kyrs, interval that records the most extreme negative δ 13 C values before a long period of recovery in δ 13 C back towards pre-excursion values (Sluijs et al., 2007a).
The use of δ 13 C stratigraphy to identify pre-and peak-PETM intervals, for the purpose of this temperature proxy synthesis, is further complicated by potential diachroneity between the onset of the PETM CIE and PETM-associated sea surface warming, as observed in some highresolution records (Thomas and Shackleton, 1996;Thomas et al., 2002;Sluijs et al., 2007b). Apparent diachroneity on the 1-10 kyr scale has also been observed between proxy records of surface and deep-ocean CIEs and temperature perturbations at two sites on Maud Rise (Thomas et al., 2002;Zachos et al., 2007). To overcome these issues we have independently identified the pre-PETM and peak-PETM conditions for surface ocean proxy data at each locality using the following considerations: 1. The pre-PETM interval should be taken below the onset of the PETM anomaly in all of the surface ocean proxy records at a given locality. Hence, the identification of the onset of the PETM is not based solely on δ 13 C records but also includes a consideration of the positive temperature anomalies, which is not always synchronous with the onset of the CIE. The pre-PETM interval is always taken below these pre-CIE temperature perturbations, as it is being used to characterise pre-PETM "background" conditions. 2. The pre-PETM interval is taken to extend stratigraphically down from the "background" (non-excursion) level nearest the onset of the PETM, over an interval that includes several pre-PETM data points. This interval should exclude any obvious pre-PETM long-term trends and represent conditions prior to the PETM across a stratigraphic interval comparable to that of the peak-PETM. In many cases the pre-PETM interval is constrained at its base by the start of a studied section  Müller et al. (2008). Abbreviated DSDP/ODP site numbers refer to following locations: DSDP Sites 213, 401, 525, 527, 549;ODP Sites 690, 738, 865, 1051, 1172, 1209, 1220, 1258, 1263. Other abbreviations are IODP M0004A (ACEX), Alamedilla (Alam.), Wilson Lake (WL), Base River (BL). Location markers outlined with black circles are those with SST proxy data used in this compilation.
Primary geochemical records across the PETM interval at the New Jersey continental shelf sites Wilson Lake (top) and Bass River (bottom). Wilson Lake planktic foraminifera stable isotope data from Zachos et al. (2006); GDGT data from Zachos et al. (2006) and Sluijs et al. (2007b). Pre-PETM stable isotope data is limited to single sample below 109.95 m where both A. soldadoensis and Subbotina spp. have δ 13 C values~3‰ higher than at the sampling interval above. The PETM onset is more clearly constrained by GDGT data, which shows a rapid warming between 109.93 and 109.8 m. The pre-excursion interval is taken down to the base of high-resolution GDGT data at 112.78 m. The base of the peak-PETM interval is placed at the first planktic foraminiferal δ 13 C excursion value and the beginning of peak warmth indicated by GDGTs at 109.23 m. The top of the peak-PETM is placed at 106.52 m towards the end of the negative plateau in foraminiferal δ 13 C and peak GDGT-based temperatures. Single specimen planktic foraminifera stable isotope data for the Bass River section are from John et al. (2008); GDGT data are from Sluijs et al. (2007b). The onset of the PETM is delineated as the rapid~3‰ negative shift in δ 13 C values of both Acarinina spp. and Subbotina spp. between 357.49 and 357.21 m. The pre-PETM interval is thus defined from 357.58 m down to the limit of high-resolution TEX86 data at 360.0 m. The base of the peak-PETM interval is taken at 357.03 m, above the negative δ 13 C shift in both Acarinina spp. and Subbotina spp. and immediately below the peak positive GDGT-based temperature anomaly. We place the top of the peak-PETM interval at 354.99 m just below the gradual decline in GDGT-based temperatures.
or the start of high-resolution sampling. Establishing strictly comparable pre-PETM stratigraphic intervals between sites is not feasible at present, and would unnecessarily limit the number of data points. 3. The peak-PETM interval is identified as extending stratigraphically upwards from the point at which PETM excursion values have been established in all surface ocean proxy records at a given locality. Again, this overcomes some of the problems of diachrony between palaeotemperature and δ 13 C records during the onset of the PETM perturbation. The length of peak-PETM conditions is harder to constrain but in this study the top of the peak-PETM interval is normally taken at the level where δ 13 C begins to consistently return to more positive values in any phasei.e. at the first sign of δ 13 C recovery as measured in any phase/taxa.
The compilation of the bottom-water temperature records of the PETM uses the existing Nunes and Norris (2006) pre-PETM and "CIE plateau" intervals to represent pre-and peak-PETM conditions. It is also based on their existing compilation of Nuttallides truempyi δ 18 O values, updated with isotope data from other species and the New Jersey shelf sites as well as available benthic foraminifera Mg/Ca data.

Sea surface temperature estimation
Currently available geochemical proxy estimates of near-surface ocean temperatures during the PETM are based on three methodologies: the oxygen isotope composition (δ 18 O) and magnesium/calcium (Mg/Ca) ratios of the calcite tests of planktic foraminifera and the composition of organic tetraether lipids produced by marine Thaumarchaeota (GDGTbased proxies). Each of these methodologies has its own particular advantages and difficulties in determining accurate temperature estimations across the PETM. All of them have been subject to either a recalibration (Kim et al., 2010), a reassessment of the assumed values of seawater chemistry used to calculate temperature (for example Coggon et al., 2010;Tindall et al., 2010;Roberts et al., 2011) or a new understanding of the preservational biases on measured geochemistry (Pearson et al., 2001;Sexton et al., 2006a;Pearson, 2012) since the original publication of much of the data. In particular, the relatively new GDGT-based proxies, dating back to the original TEX 86 proxy (Schouten et al., 2002), have been subject to three significant recalibrations in the past three years (Kim et al., 2010), with the most appropriate version for reconstruction of Paleogene palaeoclimates as the subject of ongoing discussion (Hollis et al., 2012;Bijl et al., 2013). Alongside these recalibrations, new GDGT-based palaeotemperature proxies have been proposed that vary fundamentally from the original TEX 86 proxy in the ratios of tetraether lipids (GDGTs) used Kim et al., 2010).
Given the evolution in the methodological understanding of these approaches, quoting the original temperature estimates and PETM temperature anomalies from the primary literature and making comparisons between studies spread over a decade or more, are somewhat unsatisfactory. Here we present a consistent approach to each proxy methodology, which was guided by three principles: first, to focus on the PETM temperature anomaly rather than the absolute temperatures  Zachos et al. (2003). The pre-PETM interval is taken from the base of the studied section (211.39 mcd) up to 211.27 mcdthe highest level with consistent preexcursion values of δ 13 C from single specimens of M. velascoensis. The peak-PETM interval is taken to extend from the first level with exclusively excursion δ 13 C values recorded in Morozovella velascoensis at 211.24 mcd up to 211.10 mcd. Arguably this may be slightly longer than the period of peak negative δ 13 C values but is necessary to capture the full pattern of increase in Mg/Ca ratios in both M. velascoensis and A. soldadoensis, which slightly lags the peak negative CIE. DSDP Site 527 planktic foraminifera stable isotope data are from Thomas et al. (1999) and Mg/Ca ratios are from Tripati and Elderfield (2004). The pre-excursion interval is defined on the basis of the A. soldadoensis δ 13 C data and is taken to extend from the base of the studied section up to 201.38 mbsf, above which δ 13 C values begin to show a marked negative trend into the main CIE. The peak-PETM interval is taken from 200.88 mbsf through the interval with distinctively low δ 13 C values up to 200.62 mbsf. and climate state (see Hollis et al., 2012;Lunt et al., 2012 for discussion of the background climate state); second, a consideration of the different available proxy methodologies/calibrations where there is active ongoing discussion; and third, to provide a synthesis of the underlying geochemical data, such that it can be easily reanalyzed as and when calibrations develop. Below, we detail our approach to each of the proxy methodologies used.

Oxygen isotopes
For planktic foraminifera-derived δ 18 O temperature estimates we used the calibration of Kim and ONeil (1997) for inorganic calcite, in the form presented in Bemis et al. (1998).  Kennett and Stott (1991); high-resolution single specimen data (top/middle right and bottom panel) from Thomas et al. (2002). Pre-PETM/peak-PETM stratigraphic intervals of Nunes and Norris (2006) are shown as vertical bars to the right of plots; ages in kyrs relative to the onset of the PETM are from Röhl et al. (2007). Pre-PETM interval (172.49 to 170.69 mbsf) of Nunes and Norris (2006) is shortened by the exclusion of the complex transition interval above 170.80 mbsf (Thomas et al., 2002). The base of the peak-PETM interval in this study is taken at 170.54 mbsf, which is the start of consistent excursion δ 13 C values in the planktic isotope data of Kennett and Stott (1991). Specimens of Acarinina in the pre-PETM, and Subbotina in the peak-PETM, with anomalous δ 13 C are circled with a dotted line and have been excluded from the analysis. The top of the peak-PETM interval is taken at 169.84 mbsf based on the end of the distinct plateau interval observed in Subbotina spp. δ 13 C (Kennett and Stott, 1991). At this site the peak-PETM interval is divided into two, a "base PETM" (gold) and "top PETM" (yellow) to allow for a more detailed analysis of the single specimen data of Thomas et al. (2002) see text for further discussion. Single specimen data of Thomas et al. (2002) is plotted for pre-PETM, base-PETM and top-PETM intervals in δ 18 O/δ 13 C space (bottom). Double-headed arrows indicate isotope off-sets between populations of Acarinina and Subbotina through each stratigraphic interval. This equation falls between the high-light and low-light calibrations for the symbiotic planktic foraminifera Orbulina universa (Bemis et al., 1998), and has a very similar gradient to both these equations and other commonly used calibrations (Bemis et al., 1998). For all locations we use the latitudinal correction of Zachos et al. (1994) to estimate the oxygen isotopic composition of local seawater (δ 18 O sw ) and apply a VSMOW to VPDB correction of −0.27‰. As an estimate of calibration error for planktic foraminifera, we use ±0.7°C derived from the low light calibration of Bemis et al. (1998).
Using the δ 18 O composition of planktic foraminifera to estimate the PETM temperature anomaly relies on the unlikely assumption that there is little change in the local oxygen isotopic composition of seawater across this event. An intensified hydrological cycle under warming conditions is likely to increase latitudinal oxygen isotope gradients and continental margin sites may be influenced by increased continental freshwater runoff (Zachos et al., 2006). The inference of significant changes in open ocean salinity across the PETM based on coupled δ 18 O and Mg/Ca thermometry (Zachos et al., 2003;Tripati and Elderfield, 2004), may however, be overestimated due to the damping of the PETM δ 18 O anomaly in poorly preserved foraminifera. In the initial presentation of temperature proxy data, we do not factor in any estimate of the size or direction of PETM salinity perturbationsbased on either modelled δ 18 O sw /salinity or geological evidenceand their potential effects on δ 18 O temperature estimates. Instead this is one of several factors, discussed later, which has the potential to either damp or enhance the proxy estimates of the PETM anomaly.

Mg/Ca ratios
The temperature dependent fractionation of Mg into foraminiferal shell carbonate during growth has become a standard palaeothermometer, especially in Quaternary palaeoceanography, where changes in seawater Mg/Ca ratios (Mg/Ca sw ) diverge little from modern values (Elderfield et al., 2012). In deeper-time applications, the use of foraminiferal test Mg/Ca (Mg/Ca test ) to determine absolute palaeotemperature estimates is problematic due to strongly divergent reconstructions of early Cenozoic Mg/Ca sw (see discussions in Broecker and Yu, 2011;Coggon et al., 2011;Evans and Müller, 2012), the potential for significant variability in the magnesium partition coefficient (D Mg ) with  Tripati and Elderfield (2004); single specimen stable isotope data (top/middle right and bottom left/right) from Kelly et al. (1998). Planktic foraminifera stable isotope analyses across the onset of the PETM are allocated to either pre-or peak-PETM conditions based on their δ 13 C values (N2.5‰ pre-PETM; b2.5‰ peak-PETM). Yellow diamonds on δ 13 C/δ 18 O cross-plots (bottom) represent mean values and the vertical/horizontal bars represent two standard deviations of the data. Pre-and peak-PETM intervals for the analysis of Mg/Ca data are based on the observed PETM Mg/Ca anomaly across all three planktic foraminifera species analysed. We use the intervals from 103.50 to 103.02 mbsf to represent pre-PETM and 102.83 to 102.5 mbsf to represent peak-PETM conditions. These provide reasonable sample populations (n = 6-9), bracket samples with relatively stable pre-or peak-PETM Mg/Ca ratios and avoid the zone of mixing at the onset of the PETM.
changing Mg/Ca sw (Cramer et al., 2011;Evans and Müller, 2012) and the effect of changing oceanic carbonate ion saturation state on Mg/Ca test (see Yu and Elderfield, 2008 and references therein). Temporal changes in D Mg may in fact, be responsible for the divergence in estimates of early Cenozoic Mg/Ca sw derived from coupled δ 18 O-Mg/Ca palaeothermometry assuming invariant D Mg with Mg/Ca sw (Lear et al., 2002;Sexton et al., 2006a;Creech et al., 2010) and nearly all other independent proxy and modelling estimations of the Cenozoic evolution of Mg/Ca sw (Stanley and Hardie, 1998;Lowenstein et al., 2001;Horita et al., 2002;Berner, 2004;Dickson, 2004;Demicco et al., 2005;Coggon et al., 2010;Ries, 2010;Higgins and Schrag, 2012). As Evans and Müller (2012) note, the potentially incorrect assumption of invariant D Mg with Mg/Ca sw can be counteracted by the assumption of unrealistically high Paleogene Mg/Ca sw .
Even given the double uncertainty in both the value of D Mg for planktic/benthic foraminiferal calcite and Mg/Ca sw during the early Paleogene, the form of the test calcite Mg/Catemperature equation is such that temperature anomalies over short timescales (b1 Ma), where Mg/Ca sw is assumed to be constant, remains independent of both D Mg and Mg/Ca sw . Here we use the multi-species sediment trap calibration of Anand et al. (2003), modified to include variable D Mg with changing Mg/Ca sw following the discussion in Evans and Müller (2012), such that ancient calcification temperature can be calculated: where (Mg/Ca) f is the Mg/Ca ratio recorded in ancient foraminiferal calcite; (Mg/Ca) sw and (Mg/Ca) sw′ are the Mg/Ca ratios of modern and ancient seawater respectively; A and B are the standard Mg/Catemperature calibration constants, for the Anand et al. (2003) calibration A = 0.09, B = 0.38. This assumes that there is a power law dependence of test Mg/Ca with changing Mg/Ca sw with power component, H (Evans and Müller, 2012). Although there are, as yet, limited constraints on the appropriate value of H to be used in deep-time palaeothermometry, we use an indicative value of H = 0.42, following the reanalysis of modern culture studies of Globigerinoides sacculifer (Hasiuk and Lohmann, 2010). We use a modern seawater value of 5.15 mol/mol and a latest Paleocene value of 2 mol mol −1 , which is broadly in line with the range of estimates of early Paleogene Mg/Ca sw that are not derived from paired δ 18 O-Mg/Ca palaeothermometry. Where plotted we use the calibration error of ±1.13°C (Anand et al., 2003).
The most recent calibration study proposed two GDGT-based SST proxies; one calibrated to a global core top data set (TEX 86 L ) and the other to a data set that excluded the polar oceans (TEX 86 H ) (Kim et al., 2010 (Bijl et al., 2010). We refer to an extensive review paper for further discussion of GDGT-based proxies .
For ODP Site 1172 and the New Jersey sites, we present GDGTbased SST data based on all three of the most recently proposed and calibrated proxies, TEX 86 H (GDGT index-2), TEX 86 L (GDGT index-1) and the recalibrated form of 1/TEX 86 , as detailed in Kim et al. (2010). For the Arctic Site IODP M0004A we only show TEX 86 ′ (Fig. 5), a slightly modified ratio of GDGT isomers derived to deal with anomalously high concentrations of the GDGT-3 isomer found in the uppermost Paleocene at this site Sluijs et al., 2009  In both cases pre-and peak-PETM intervals are based on δ 13 C stratigraphy. ODP Site 1172 data from Sluijs et al. (2011); ACEX data from Sluijs et al. (2006). Pre-PETM interval in ODP Site 1172 taken between 612.61 and 612.20 mcd where both δ 13 C and GDGT-based temperature estimates are stable and showing pre-PETM values. Peak-PETM interval is taken across the interval with most negative and roughly stable δ 13 C between 611.82 and 611.64 mcd. At the ACEX site the pre-PETM interval is taken between 388.61 and 390.71 mcd at the base of the studied section; this excludes three δ 13 C and GDGT analyses at the top of core 32X, which are within a zone of drilling disturbance and include an anomalously low δ 13 C value. The peak-PETM interval is taken between 386.10 and 383.22 mcd, across the poorly recovered core 31X and into the base of core 30X. This is the interval with relatively stable peak excursion values of both δ 13 C and GDGT-based temperature estimates.
(1) the consistency of behaviour between the various GDGT-proxies and the identification of problematic samples or environments for this approach; (2) a comparison of GDGT-based proxies with other proxy methodologies in the estimation of absolute temperatures; and, (3) the variation between GDGT-proxies in estimating the magnitude of climate anomalies such as the PETM warming. We use calibration errors of ±2.5, 4.0 and 5.4 for TEX 86 H , TEX 86 L and 1/TEX 86 respectively (Kim et al., 2010) and ±2.0 for TEX 86 ′ ).

Bottom water temperature estimation
Geochemical analyses carried out on the test calcite of calcareous benthic foraminifera across the PETM were collated from a total of sixteen sites (DSDP Sites 213, 401, 525, 527, 549;ODP Sites 690, 738, 865, 1051, 1209, 1220, 1258, 1263, the New Jersey drill cores from Wilson Lake and Bass River and the Alamedilla section in southern Spain). Temperature estimates are based upon both oxygen isotope and Mg/Ca analyses. Data is predominantly from the Nunes and Norris (2006) compilation of stable isotope data measured on the cosmopolitan foraminifera Nuttallides truempyi. We have added stable isotope data from the genera Oridorsalis (ODP Sites 1051 and 1263) (Katz et al., 2003;McCarren et al., 2008), Bulimina (ODP Site 865) (Thomas et al., 2000) and Cibicidoides (Wilson Lake and Bass River) (Zachos et al., 2006;John et al., 2008). It is very likely, however, that the "Cibicidoides" records from Wilson Lake, are actually a mixed assemblage dominated by analyses of Anomalinoides (Stassen et al., 2012). We have also included benthic foraminifera Mg/Ca from DSDP Site 527 and ODP Sites 865 and 1209 (Tripati and Elderfield, 2005).
Bottom water temperature estimates from foraminiferal δ 18 O measurements were calculated using the same Kim and ONeil (1997) calibration for inorganic calcite used for planktic foraminifera above. Following Sexton et al. (2006b) and Katz et al. (2003), we correct both Cibicidoides spp. and Nuttallides truempyi to Oridorsalis spp., which is assumed to calcify in equilibrium. The correction factor of +0.28‰ is applied to Cibicidoides spp. δ 18 O values (Katz et al., 2003) and +0.35‰ for N. truempyi (Sexton et al., 2006b); no correction is applied to Bulimina δ 18 O. For all locations we use a uniform oxygen isotopic composition of seawater (δ 18 O sw ) of −1‰ for an ice free world and apply a VSMOW to VPDB correction of −0.27‰.
The Oridorsalis umbonatus calibration of Lear et al. (2002) was used to estimate temperatures from mixed species benthic foraminifera Mg/Ca data. Modified for a variable magnesium partition coefficient with Mg/Ca sw this takes the same form as Eq. (2) above, but with values of A and B of 0.114 and 1.008 respectively. The original Mg/Ca data of Tripati and Elderfield (2005) was already normalised to O. umbonatus by Tripati and Elderfield (2005) using Lear et al. (2002). To determine an estimate of the power law component, H, for O. umbonatus, in the early Paleogene oceans, we follow the method of Evans and Müller (2012) using the coupled δ 18 O-Mg/Ca thermometry at 49 Ma of Lear et al. (2002). We modify their solution slightly by using the original δ 18 O-palaeotemperature estimate of 12.4°C (Lear et al., 2002) and assume a Mg/Ca sw of 2 mol mol −1 rather than 1.6 mol mol −1 (Evans and Müller, 2012). This yields an H value of 0.5, which is only slightly lower than the range 0.52-0.54 derived by Evans and Müller (2012). Again it should be noted that PETM temperature anomalies calculated on the basis of benthic foraminiferal Mg/Ca test are independent of the assumed value of both H and Mg/Ca sw .
The Mg/Ca-based temperature estimates presented here are no more well constrained than previous estimates, and this uncertainty will remain until estimates of early Paleogene Mg/Ca sw begin to converge. They do demonstrate, however, that realistic surface and intermediate water temperatures can be achieved with much lower estimates of Mg/Ca sw than previously used in deep-time palaeoceanography, as long as the effect of a variable magnesium partition coefficient is factored into these calculations.

Model derived temperature estimates
Climate model realisations under early Eocene boundary conditions were undertaken with the UK Met Office GCM, HadCM3L. These are the same suite of simulations as described in Lunt et al. (2010) and Tindall et al. (2010). In summary the model is run with a solar constant 0.4% less than the modern, an early Eocene topography and bathymetry of Markwick (2007), no land ice, generic grassland across the land surface, lakes in appropriate regions and a runoff routing scheme based on the model orography. Results presented here are taken from four early Eocene simulations with atmospheric CO 2 prescribed at ×1, ×2, ×4, and ×6 pre-industrial values (i.e., 280, 560, 1120, and 1680 ppm). All early Eocene simulations were initialized from the end of a run with ×3 pre-industrial CO 2 (840 ppmv), which in turn was initialized with zonally-symmetric ocean temperatures and salinity distribution, zero ocean velocity and a modern atmospheric state. After initialization, each early Eocene simulation was run for more than 3400 years. The surface temperatures are based on the average climatology of the final 100 years of these simulations.

Proxy-model comparisons
A potentially significant source of error in deep-time climate model/ proxy data comparisons can arise from divergent or erroneous palaeogeographies of model and proxy data (see discussion in Huber and Caballero, 2011). Relevant to this study are (1) the use of significantly different palaeogeographic reconstructions for the model configuration and proxy data palaeolocations, for example the location of sites relative to continental margins and ocean gateways; and, (2) absolute errors in the palaeogeographic reconstruction of proxy data locationse.g. an error in the latitude, longitude or the assumed representative ocean depth (surface, mixed layer, thermocline) of a proxy measurement. Given the effort required to configure a high resolution AOGCM to palaeogeographies significantly different from the modern, basic model configurations do not update quickly to revised palaeogeographic reconstructions. The model palaeogeography of HadCM3L used here is based on the early Eocene palaeogeographic reconstruction of Markwick (2007) and is the same configuration as used in Tindall et al. (2010), Lunt et al. (2010) and Dunkley .
For consistency with recent Eocene climate modelling studies (Huber and Caballero, 2011;Lunt et al., 2012) all palaeolocations of marine proxy data points are calculated using the GPLATES software and the plate model of Müller et al. (2008). A comparison between the HadCM3L palaeogeography based on Markwick (2007) and the GPLATES reconstruction of continental outlines as well as the locations of proxy data at 55 Ma is shown in Fig. 1. The two palaeogeographies are comparable in all features of the major ocean basins, including the width of key gateways (Tethys, Tasman Gateway, Drake Passage, Panama) and geometry of the central and South Atlantic. There is some discrepancy in the palaeolatitude of the Indian subcontinent, which is reconstructed in a more northerly position under the plate rotations of Müller et al. (2008). Given there are no proxy data for SST conditions in the Indian Ocean, this has no significant effect on modeldata comparisons.
When undertaking the model-data comparison, we explicitly consider uncertainty associated with the assignment of a model depth to proxy data and the potential for latitudinal/longitudinal mismatches between data and model reference frames. This is achieved by comparing each proxy measurement with an array of model temperatures, centred on the reconstructed palaeolocation and most likely depth level. For GDGT-proxies, assumed to represent true sea surface temperatures (Kim et al., 2010), the array is formed of the top two model layers (5 m and 15 m); for mixed-layer dwelling planktic foraminifera species (Acarinina and Morozovella; Norris, 1996) we used data from model layers centred on 25, 35, and 50 m; for thermocline-dwelling foraminifera (Subbotina; Norris, 1996) we use model layers 140 and 200 m.
Although these depths may be somewhat shallow for both mixed-layer and thermocline-dwelling species and do not account for spatial variability in the depth of the thermocline, a more sophisticated comparison is hampered by the loss of model resolution at depths below 200 m. This approach is, however, an improvement on a simple comparison to a single sea surface temperature, by attempting some consideration of the thermal gradients in the upper ocean thermocline. Around any proxy data palaeolocation, the array of model data also includes temperatures from appropriate depth layers in the adjacent model grid cells to the north, south, east and west of the central cell. Where 2-layers are considered, this results in an array of 10 spatial localities (2 depths at 5 points) and 15, when three depths are considered. 100-year seasonal mean temperature data was then extracted from each of these points for each three-month season (JJA, SON, DJF, MAM; see Section 2.4 below).
In addition to the palaeogeographic and depth uncertainties, we explicitly consider model temperatures from all four seasons (JJA, SON, DJF, MAM) to address potential proxy data biases to a particular season of production. This is most likely to affect mid-to high latitude locations where there is a strong seasonality in surface ocean biological productivity that may determine both the flux and species composition of planktic foraminifera (Žarić et al., 2005;Chapman, 2010) and the transport efficiency of Thaumarchaeota to marine sediments (Kim et al., 2010). For each proxy data location an array of data was extracted from the model formed of the three-month seasonal averaged temperatures (JJA, SON, DJF, MAM) at all of the depths and spatial locations described above. In total this produced data arrays of either 40 (for 2 depth layers) or 60 (3 depth layers) model temperatures. At ODP 1172 the model grid box to the south of the site's palaeolocation was terrestrial resulting in a slightly reduced data array (32 points).

Results
The primary geochemical and palaeotemperature proxy data are shown in Figs. 2-6. The pre-and peak-PETM intervals used in this study are also highlighted. Original depth or composite depth scales have been used to present temperature proxy data alongside the associated carbon isotope stratigraphy. Figures also include the extent of the pre-PETM (blue) and PETM (orange) intervals used here. Palaeotemperatures were then calculated from this geochemical data, as outlined above for the different proxies, and mean temperature estimates calculated for the pre-and peak-PETM intervals. These mean values, together with the standard deviations and number of data points in each interval, are listed in Table 1, and provide some quantification of variability and the time resolution of samples within each stratigraphic interval.
In the following we discuss: (1) the complexities of multi-species foraminifera records from ODP Sites 690 and 865; (2) the application of the various GDGT proxy methods across the PETM; and (3) planktic foraminifera preservation and paired δ 18 O and Mg/Ca records from deep-ocean sediments.

Multi-species foraminiferal records from ODP Sites 690 and 865
The PETM carbon isotopic excursion was first recognised in cores from ODP Sites 689 and 690, on Maud Rise in the Weddell Sea of the Southern Ocean (Kennett and Stott, 1990). These same cores have subsequently been used for many of the classic studies of the PETM event (Kennett and Stott, 1991;Thomas and Shackleton, 1996;Thomas et al., 2002). The detailed single specimen analyses of Thomas et al. (2002) (Fig. 4)  A second characteristic of the ODP Site 690 record is a distinct shift in the isotopic composition of Acarinina tests to heavier δ 13 C and lighter δ 18 O values at 170.21 mbsf,~25 kyr after the onset of PETM on the orbital age model of Röhl et al. (2007) (Fig. 4). We use this event to separate the PETM into two intervals, "base-PETM" (~0-25 kyr) and "top-PETM" (~25-50 kyr), which together represent the first~30% of the PETM (Röhl et al., 2007), A cross-plot of the single specimen Acarinina and Subbotina δ 18 O and δ 13 C data is shown in Fig. 4. The isotopic gradient between mixed-layer species Acarinina and the deeper-dwelling Subbotina reflects the strength of the thermocline (δ 18 O) and the balance between surface biological productivity and organic matter remineralization with depth (δ 13 C) (John et al., in press). Prior to the onset of the PETM this gradient is about 2‰ in δ 13 C and slightly less than 1‰ in δ 18 O. Immediately after the onset of the CIE, during the base-PETM, the isotopic gradient in both δ 18 O and δ 13 C between Acarinina and Subbotina reduces dramatically to around 0.5‰. This could be due to a transient break down in stratification and enhanced vertical mixing during the base-PETM (Thomas and Shackleton, 1996;Thomas et al., 2000Thomas et al., , 2002. Such a scenario is difficult to reconcile with records of phytoplankton assemblages indicating warminginduced oligotrophy from the very base of the CIE at ODP 690 . Instead, we suggest that this reduction in gradients is largely due to the depth-migration of surface-dwelling Acarinina, towards the habitat of Subbotina during the initial phase of the PETM. This collapse in isotopic gradients may have been enhanced by the loss of photosymbionts by Acarinina during the PETM. Combined depth-migration and loss of symbionts appears to be an integrated ecological response to transient warming events (Kelly et al., 1996;Kelly et al., 1998;Edgar et al., 2013). Symbiont bleaching and migration to deeper habitats both result in negative δ 13 C and positive δ 18 O shifts in planktic foraminifera test calcite. During the MECO bleaching caused a negative δ 13 C perturbation in mid-size range Acarinina tests from ODP Site 748 of~1‰ (Edgar et al., 2013), similar in magnitude to the difference between high and low light experiments on modern symbiontbearing planktic foraminifera (Spero et al., 1997). The effect of bleaching on δ 18 O is more difficult to assess from the geological record, but the same culture experiments indicate that bleaching could cause a positive shift in δ 18 O of~0.3‰ (Spero et al., 1997) due to the loss of CO 2 uptake by the symbionts and the associated carbonate ion effect on oxygen isotopes. During the "top-PETM" interval, above 170.21 mbsf, the δ 13 C gradient between the two genera is restored to~2‰ and the δ 18 O gradient is increased to nearly 1.5‰, as Acarinina morphotypes reestablish populations in the mixed-layer and reacquire photosymbiotic algae. If this interpretation is correct, then using pre-PETM to base-PETM intervals to characterise the CIE and peak-PETM warming is likely to result in an underestimate of warming and an overestimate of the CIE due to the depth-migration of Acarinina in the early stages of the PETM.
The PETM record from ODP Site 865 Allison Guyot in the western Pacific is not ideal, with suggestions that the peak-PETM is missing or condensed (Nunes and Norris, 2006) and evidence for sediment mixing through the onset of the event (Kelly et al., 1998). Nuttallides truempyi isotope data from this site (Bralower et al., 1995) were excluded from the Nunes and Norris (2006) compilation because of their lowresolution and poor stratigraphic constraints. There are, however, more detailed stable isotope records from both Cibicidoides and Bulimina sp. that resolve a clear PETM CIE (Thomas et al., 2000). We have decided to consider the planktic foraminiferal isotope and trace metal data from this site for two reasons: (1) the analysis of numerous single planktic foraminifera specimens across the onset of the PETM (Kelly et al., 1998) allows the clear separation of specimens with preand PETM δ 13 C values (see Fig. 5); and (2) the trace element data from this site (Tripati and Elderfield, 2004) resolve a PETM warming anomaly that can be analysed to give some estimate of pre-and peak-PETM conditions (Fig. 5). Foraminiferal preservation at this site is poor (Plate 1; Kozdon et al., 2011), and all resultant geochemical proxy data should be treated with caution. Due to sediment mixing through the PETM onset, planktic foraminifera stable isotope data are separated into pre-and PETM analyses on the basis of their δ 13 C signature (pre-PETM δ 13 C N 2.5‰; PETM b 2.5‰) rather than stratigraphic position. Considering that the mean δ 13 C values of these two groups differ bỹ 3‰, it is probable that these analyses represent near to peak-PETM conditions at this location. Due to the absence of excursion taxa prior to the PETM, we only use data from the non-excursion taxa to assess the temperature anomaly across the PETM. Excursion taxa could be included to produce generic level isotope anomalies but it is clear that there is a species-level differentiation in δ 18 O isotope values in both genera during the PETM (Fig. 5), which would bias a generic-level estimate of the PETM isotope anomaly. The records of A. soldadoensis and M. velascoensis across the onset of the PETM may themselves be somewhat biased (enhanced CIE/damped Δδ 18 O) by the partial exclusion, depthmigration and loss of symbionts within these taxa during the peak PETM (Kelly et al., 1998). Isotopic analyses from the deeper-dwelling Subbotina spp. show no clear excursion δ 13 C values (Bralower et al., 1995;Kelly et al., 1998), possibly due to their ecological exclusion during peak-PETM conditions (Kelly et al., 1998). These data are not considered further but are available in the supplementary information.
The trace metal analyses at ODP Site 865, based on multiple (10-20) foraminiferal tests picked from each sample interval (Tripati and Elderfield, 2004), are susceptible to sediment mixing and signal attenuation across the onset of the PETM. Pre-and peak-PETM intervals are based on the observed Mg/Ca anomaly across the planktic foraminifera taxa analysed (A. soldadoensis, M. velascoensis and Subbotina sp.) (Fig. 5). All three species show consistently higher Mg/Ca ratios above 102.83 mbsf, and we take the interval between 102.83 and 102.5 mbsf to be representative of "peak-PETM" Mg/Ca ratios and 103.02 to 103.50 mbsf to represent pre-excursion Mg/Ca values. Benthic foraminiferal Mg/Ca data in Tripati and Elderfield (2005) are published as species-specific analyses but the data provided with this are reported as unidentified benthic foraminiferal analyses (NOAA NCDC: ftp://ftp. ncdc.noaa.gov/pub/data/paleo/contributions_by_author/tripati2005b/ tripati2005b.txt). Here they are treated as a mixed benthic assemblage. The PETM excursion in benthic foraminiferal Mg/Ca is more gradual than the CIE onset recorded within Bullimina spp., with peak-PETM Mg/Ca found~15 cm above the most negative δ 13 C values. To maximise the representation of peak-PETM conditions, this interval is selected independently for isotope (102.88-102.98 mbsf) and Mg/Ca records (102.75-102.9 mbsf). For reference the benthic foraminifera extinction level is between 103.03 and 103.01 mbsf (Kelly et al., 1996;depth revision Ellen Thomas pers. comm., 2013).

GDGT-based palaeothermometry
To evaluate the consistency of behaviour between GDGT-proxies at ODP Site 1172 and the New Jersey sites, we have plotted SST estimates from all three proxies against the values of both TEX 86 and GDGT-1 indices (Fig. 7). For plots against TEX 86 , SST estimates based on TEX 86 H and 1/TEX 86both based on the underlying TEX 86 indexsimply follow their respective calibration lines. SST estimates from TEX 86 L , however, are derived from the index GDGT-1, and although showing a strong correlation between SST and TEX 86 also contains a degree of scatter produced by the differential behaviour of TEX 86 and GDGT-1. For plots against GDGT-1, the reverse is true, with TEX 86 L following its calibration and TEX 86 H and 1/TEX 86 displaying a broad trend with some scatter. Ideally, there would be a good correlation between the behaviour of TEX 86 and GDGT-1 across temperature and palaeoenvironments; i.e. a minimum of scatter in the plots shown in Fig. 7. This would reflect a strong and dominant control by ocean temperature on the lipid structures produced by ancient marine Thaumarchaeota. For Bass River, Wilson Lake and ODP Site 1172 the SST estimates from all three proxies behave in a generally consistent manner. The exception to this are four data points below the onset of the PETM in Bass River (duplicate samples from levels 358.83 m and 362.57 m), where, in comparison with surrounding data points, TEX 86 L produces anomalously high SST estimates bỹ 3°C (points with SST estimates of N25°C in the interval prior to the PETM in Fig. 2). In As noted in the original studies, the relationships between GDGTproxies from the Arctic Ocean ACEX site are more complex ). The TEX 86 index does not show a clear warming trend across the onset of the PETM (Fig. S2 of Sluijs et al., 2006). Sluijs et al. (2006) attribute this to anomalously high concentrations of the GDGT-3 isomer in Paleocene (pre-PETM) samples. For a more realistic estimate of the magnitude of PETM warming, they therefore proposed the new proxy TEX 86 ′ that excluded the GDGT-3 isomer. It is important to note that for samples that do not yield high GDGT-3 isomers, TEX 86 and TEX 86 ′ provide almost identical values Sluijs et al., 2009). Moreover, contrary to most other records from across the globe, ACEX TEX 86 and TEX 86 ′ values are within the modern calibration and therefore yield much smaller calibration errors. Therefore, apart from samples in the ACEX record that are accompanied by high terrestrial lipid input as indicated from the BIT index , there is no reason to treat absolute SSTs based on TEX 86 ′ differently than other records. However, SSTs in this region could be skewed to summer values .
Given the relative consistency of GDGT-proxy behaviour at the New Jersey shelf sites and ODP Site 1172, the question remains open as to which index and calibration best represents both absolute SSTs at these sites and the warming anomaly across the PETM. The early Paleogene proxy inter-comparison study of Hollis et al. (2012) in the South west Pacific suggests that TEX 86 L is most likely to represent SSTs in this region but even these estimates may be biased towards summer rather than mean annual temperatures. In Fig. 7, alongside the GDGT-derived SSTs, we plot the distribution of temperature estimates derived from the δ 18 O composition of well-preserved specimens of the mixed-layer planktic foraminifera Acarinina (shown as mean ± 1 standard deviation and maximum/minimum values). These clearly show that TEX 86 L SST estimates are in closest agreement with both pre-and peak-PETM δ 18 O estimates, whilst TEX 86 H and 1/TEX 86 estimates are~5°C warmer.
This temperature offset cannot be explained by poor preservation of foraminifera at this location, which is generally good to excellent (Plate 1) or errors in δ 18 O sw estimates, which in these shelf settings are more likely lower (less saline) than those used here, rather than being the~1‰ heavier required to account for a~5°C cool-bias in δ 18 O SSTs. Large seasonal variations in temperature is a possibility at these locations and is discussed in Section 4.1. Finally, we assess the variability in the PETM temperature anomalies derived from the range of GDGT-proxies. The pattern of GDGT-derived SST anomalies through the PETM at both the New Jersey sites and ODP Site 1172 is relatively consistent between the different GDGTproxies (Figs. 2 & 6). The difference in the magnitude of PETM warming between these proxies is more clearly shown in Fig. 10. For the New Jersey sites and ODP 1172, the warming derived from TEX 86 H is consistently larger than that from 1/TEX 86 , which is simply a result of their respective calibration equations and becomes more pronounced with increasing temperature as their calibrations diverge. In contrast, the anomalies derived from TEX 86 L vary from site-to-site with respect to the other GDGT-proxies, with larger temperature anomalies at the New Jersey sites and a smaller anomaly at ODP Site 1172 (Fig. 10). This is apparent in subtle changes in the relative gradients between TEX 86 L and the other two proxies in plots of SST against GDGT-1 (Fig. 7). This may be indicative of a more systematic control on GDGTproxy gradients with environment or latitude but demonstrating this would require data from more locations.

Influence of planktic foraminifera preservation on palaeotemperature estimates
A comparison of planktic foraminiferal δ 18 O and Mg/Ca-derived temperatures estimates clearly shows the now widely recognised cool-bias in δ 18 O-derived SSTs from deep-ocean sediments (Pearson et al., 2001;Sexton et al., 2006a). This is particularly clear in the low and mid-latitudes with δ 18 O derived temperature~6-10°C cooler than Mg/Ca temperatures (DSDP Site 527, ODP Sites 865 and 1209; Fig. 8). Although there is significant uncertainty in both the estimate of early Paleogene seawater Mg/Ca and the appropriate value of the magnesium partition coefficient, the values used here produce palaeotemperature reconstructions from well-preserved Eocene planktic foraminifera that are not inconsistent with δ 18 O reconstructions (cf. Sexton et al., 2006a; reconstructed palaeotemperature using these parameters and shell calcite Mg/Ca ratio of 4.5 mmol mol −1 is 31.9°C). Given the extensive work on the effects of diagenetic recrystallisation on δ 18 O composition of planktic foraminifera (Pearson et al., 2001;Sexton et al., 2006a) and a preservation mode typical of other mid-to low-latitude deep-ocean locations (Plate 1), we consider δ 18 O derived temperature estimates from all of these sites to be substantially cool-biased. Trace metal concentrations, and associated Mg/Ca-derived palaeotemperatures, appear to be less impacted by diagenetic recrystallisation but there is little understanding of the processes behind this observation. The impact of diagenesis on δ 18 O records of transient temperature perturbations, such as across the PETM, is more difficult to quantify.
One argument against a diagenetic bias on estimates of warming, proposed by Zachos et al. (2003), is that the addition of deep-ocean diagenetic calcite to foraminifera tests is likely to operate uniformly across the whole PETM sequence, producing an offset in absolute temperature estimates but not biasing the magnitude of short-term perturbations. Diagenetic processes are now known to operate via pervasive micronscale recrystallisation throughout the foraminiferal test, rather than the discrete addition of diagenetic calcite (Sexton et al., 2006a). The resultant δ 18 O composition of the recrystallised test is thus on a mixing line between its primary value and in situ deep-ocean diagenetic calcite, depending on the extent of recrystallisation (Pearson et al., 2001). In tests that have undergone significant recrystallization, this will act to damp a primary anomaly in δ  Fig. 8. a and b. Site-specific response of HadCM3L model to increased GHG forcing. Model responses are plotted as the mean (line and crosses) and maximum to minimum range (shaded area) of a temperature array that includes both spatial (one grid cell laterally and 2 or 3 depth layerssee text for details) and seasonal variability. Proxy temperature estimates are plotted for comparison, including both 95% confidence intervals derived from data variability within the pre-and peak-PETM intervals (black error bars) and indicative calibration uncertainty for each proxy methodology (blue/red error bars). In all cases the calibration uncertainty is significantly larger than the variability of temperature estimates within the pre-and peak-PETM intervals.
δ 18 O, it is notable that the sites with the worst planktic foraminiferal preservation (ODP Sites 865 and 1209) show the smallest perturbations in δ 18 O, which are inconsistent with Mg/Ca derived temperature anomalies from the same locations. As a conservative interpretation, we consider these offsets between δ 18 O and Mg/Ca derived temperature anomalies to be primarily due to poor test preservation rather than any major changes in surface ocean δ 18 O sw and salinity (Zachos et al., 2003;Tripati and Elderfield, 2004).

Discussion
In assessing the PETM temperature anomaly, we take an approach that uses the characteristic response of an early Paleogene climate model to CO 2 -driven warming to inform the interpretation of PETM climate proxy data and act as an interpolator of sparse data. The use of model simulations to interpret the palaeoclimate proxy data aims to better quantify the global mean surface temperature response on the basis of spatially limited proxy data. In this section we first describe the model's site-specific response to CO 2 forcing and use this to derive an indicative "PETM" temperature anomaly before discussing the proxy records of PETM warming informed by the model responses.

Site-specific responses of HadCM3L to CO 2 forcing
The site-specific SST responses of HadCM3L, at CO 2 forcings of 1×, 2×, 4× and 6× pre-industrial CO 2 (280 ppm), are shown in Fig. 8. These temperature responses are taken from the data arrays outlined in Section 2.4, which include spatial, depth and seasonal variability. The mean and maximum to minimum range of these arrays, are shown alongside proxy temperature estimates for pre-and peak-PETM intervals. Notable features of these plots are: (1) the sitespecific nature of warming trends with increasing CO 2 ; (2) the range and variability of seasonality between sites; and, (3) the first-order comparison between model and proxy SST estimates.
Comparisons should be made with the warming trend observed between the highest CO 2 simulations (4× and 6× CO 2 ), where absolute temperatures are closest to those of the pre-PETM proxy palaeotemperature estimates. The site-specific patterns of warming between these simulations reflect the more general response of the model δ 18 O from planktic foraminifera; mixed layer data from Acarinina spp. only; sub-thermocline data is from Subbotina spp; ODP 690 data is from Kennet & Stott (1991) and Thomas et al. (2002).
Mg/Ca from planktic foraminifera; for the mixed layer data plotted in the order Acarinina and Morozovella spp. (left to right); sub-thermocline data is from Subbotina spp. climatology, which depending on the accuracy of the model response and fidelity of the proxy records, may be present in palaeoclimate proxy data. Briefly, these are (1) a reduction in the magnitude of SST warming of tropical to mid-latitude sites in the open ocean, specifically DSDP Site 527 and ODP Sites 865 and 1209, in response to the same CO 2driven radiative forcing at higher CO 2 concentrations; (2) high SST seasonality at the New Jersey sites on the north-west Atlantic continental margin; (3) relatively high temperature sensitivity to CO 2 -forcing at the Southern Ocean ODP Site 690; and, (4) seasonally asymmetric response of the Arctic Ocean site to increasing CO 2 with more pronounced warming in the summer season. Given the caveats on the proxy data discussed in various sections above, a comparison between proxy and model temperature estimates (Fig. 8) clearly shows, (1) the cool-bias in δ 18 O SST estimates at DSDP Site 527 and ODP Sites 865 and 1209; (2) due to the large model seasonality of SSTs at the New Jersey sites, all GDGT-and δ 18 O-based based proxy estimates for the pre-PETM fall within the temperature envelopes of the 6×CO 2 model simulation; (3) GDGT-proxies are on the order of 10°C warmer than even the warm-season model temperatures at the high-latitude sites ODP 1172 and ACEX. Point 2) is important to consider if using proxy inter-comparisons to identify the most appropriate version of GDGT palaeothermometry for deep-time applications. For 6×CO 2 the modelled seasonal range at the New Jersey margin sites is~10°C, which is at the low end of the~15°C or greater seasonal SST range over much of the New Jersey Shelf today (Yoder et al., 2002). If these model results are genuinely indicative of the seasonal SST range of this region during the early Paleogene then any of the GDGT-based proxy estimates could be reconciled with δ 18 O-based SST estimates by seasonal offsets in the productivity of planktic foraminifera and marine archaeota.
Explaining the large mismatch in GDGT-proxy and model SSTs at ODP Site 1172 in the high southern latitudes and the ACEX site in the Arctic, is a critical challenge for understanding early Paleogene warm climates. Although this mismatch between warm proxy data and cooler model SSTs is a persistent feature across early Eocene climate models (Lunt et al., 2012), the scale of this offset is likely exacerbated by particular features of the HadCM3L configuration and climatology. ODP Site 1172, in the western Tasman Sea, is the coldest modelled region at this latitude in the Southern Ocean (Fig. 9), which is also a feature of other AOGCM simulations (Hollis et al., 2012) but may be highly 6x -4x normalised to CO 2 doubling; proxy ΔT 6xCO 2 (1680ppm) Fig. 9. Early Eocene HadCM3L GCM simulations of sea surface temperature compared to pre-and peak-PETM proxy SST estimates. 4× pCO 2 simulations compared to pre-PETM (top left) and peak-PETM (top right) proxy SST estimates. 6× pCO 2 simulations compared to pre-PETM (bottom left) and peak-PETM (bottom right) proxy SST estimates. SST difference 6x-4x simulations normalised to a doubling of pCO 2 (centre) compared against proxy estimates of the PETM temperature anomaly (ΔT). Proxy estimates plotted are as follows: ACEX -TEX 86 ′; Wilson Lake -1/TEX 86  dependent on regional bathymetry and oceanography (Bijl et al., 2013). Eocene simulations (1500 ppm CO 2 ) with the intermediate-complexity UVic model, however, have both higher Southern Ocean SSTs in general, of the order 10-14°C, and do not show this cooler region on the east Australian margin (Sijp et al., 2011). The authors attribute this to a proto-East Australian Current, even with a closed Tasman Gateway, that warms this region by 2-4°C (Sijp et al., 2011) but this remains inconsistent with biogeographic patterns . This difference in model response between intermediate-complexity and AOGCMs is important for understanding heat transport to the Southern Ocean in warm climate states; the SST fields of intermediate-complexity models are in better agreement with available proxy data but the reasons for this are unclear. Even assuming the~10°C warmer Southern Ocean SST fields of UVic and a summer bias of~5°C, model-proxy agreement for ODP Site 1172 cannot be fully achieved. In the Arctic, a strong seasonality in archaeal lipid production would place ACEX GDGT-proxy estimates close to warm month mean terrestrial temperatures of over 18°C (Eberle et al., 2010). Together with reconstructions of the Arctic zero degree isotherm (Markwick, 2007), these re-emphasize that, even at 6× CO 2 , these Eocene HadCM3L simulations have a persistent, year-round cold-bias in the Arctic on the order of 5-10°C.

Proxy-and model-based estimates of PETM warming
To estimate the site-specific CO 2 -driven warming in the model, we calculate the temperature change between the 4× and 6× simulations normalised to a doubling of CO 2 concentrations. This metric is an approximation of the temperature anomaly between the 4× simulation and the predicted response of a model run at 8× CO 2 , if the temperature sensitivity of all sites to GHG-forcing was constant across this range of CO 2 . This is represented by the dashed line extensions to 8× CO 2 in Fig. 8 and the geographical pattern and magnitude of warming is shown in Fig. 9.
Proxy-derived PETM temperature anomalies, estimated as the difference between mean pre-and peak-PETM temperatures, are shown in Fig. 10 and Table 1. These are grouped by location and ranked by the size of anomaly. There is significant variability in temperature anomalies between different proxy methods at almost all locations, whilst the model responses are relatively uniform across depths (surface, mixed-layer, thermocline). The variability in the temperature anomaly estimated using the different GDGT proxies (discussed in Section 3.2), which are derived from the same underlying set of GDGT analyses, must arise from the form of the proxy and particular calibration used. The only site with a reasonably good independent proxy estimate of warming to compare to the GDGT-proxies is Bass River, where there are sufficient pre-PETM Acarinina and Subbotina δ 18 O-derived temperature estimates to define an anomaly. Here, the δ 18 O-derived anomaly (3.0 and 3.7°C for Acarinina and Subbotina respectively) is smaller than any of the GDGT-proxy anomalies (4.3, 5.4 and 5.7°C for 1/TEX 86 , TEX 86 H and TEX 86 L respectively). The difference between δ 18 O and Mg/Ca-derived temperature anomalies at deep-ocean sites ODP 865 and 1209 of~1-2°C, is more easily accounted for by foraminiferal test recrystallization significantly damping the PETM perturbation in δ 18 O. This was originally ascribed to major changes in open ocean salinity across the PETM (Zachos et al., 2003;Tripati and Elderfield, 2004), but we attribute this to poor preservation. The magnitude of the uncertainty of these estimates also approaches the difference between them. The final marked variation between proxy estimates of the temperature anomaly is at ODP Site 690, where, at 9°C, the temperature anomaly between the pre-PETM and top-PETM interval in the single specimen Proxy-model ΔT residual o C ACEX S u b . Fig. 10. Proxy estimates of the PETM temperature anomaly based on peak-PETM minus pre-PETM temperatures for each site and proxy methodology. Colour coded by broad geographical regions. Error bars on proxy estimates represent 95% confidence intervals based on the sampling frequency and variability in temperature estimates during both the pre-and peak-PETM stratigraphic intervals; they do not include any representation of proxy calibration uncertainty but are presented to indicate the reliability of the primary proxy data sets. Site-specific and depth-specific (surface, mixed layer, deepsee text for details) warming between 4× and 6× HadCM3L model simulations normalised to a doubling of CO 2 are shown as horizontal grey bars. Top panel shows the residual between these proxy estimates of PETM warming and the model-derived warming.
Acarinina spp. analyses of Thomas et al. (2002) is over 3°C larger than any other δ 18 O-derived estimate at this site. The possible upwards depth migration, from the thermocline to mixed layer, and reacquisition of photosymbionts is the likely explanation for the coordinated shift in δ 18 O and δ 13 C at the transition into the top-PETM interval. This is, however, unlikely to account for the large anomaly as measured from the pre-PETM state to the top-PETM state because in both intervals it seems likely that Acarinina dwelt in the mixed layer and was photosymbiotic. A related ecological explanation is possible, such as the post-PETM recolonization of the mixed-layer by a different population of Acarinina with a significantly lighter δ 18 O signature or altered seasonality. Although this region does show one of the strongest warming responses in HadCM3L (Fig. 9), a strong vertical stratification is not reproduced, with a differential warming between the mixedand deep-layers of HadCM3L of less than a degree (4.9 and 4.0°C respectively). Given the uncertainty in the ecological controls on the Acarinina isotope record, and the consistent model prediction that mixed-and deep-layer temperature anomalies should be similar within a degree, we consider the deep-dwelling Subbotina record to provide a more accurate estimate of warming at ODP Site 690, in the region of 5.2-5.8°C.
Considering the specifics and limitations of the proxy data for each site, we suggest that the following estimates of the PETM anomalybased on the available datato be the most robust for each site: Mg/Ca at DSDP Site 527 (~2-3°C) and ODP Sites 865 (~1-2°C) and 1209 (~4°C); TEX 86 ′ at ACEX (~4°C); 1/TEX 86 at the New Jersey margin sites (~4 and~6°C); TEX 86 L at ODP Site 1172 (~4°C); and, Subbotina spp. δ 18 O at ODP Site 690 (~6°C). The estimate from ODP Site 865 is certainly compromised by very poor foraminifera preservation. Across all the sites, however, these figures are somewhat lower, than previous estimates of warming at the PETM, for example the "5-9°C" of Zeebe et al. (2009). There is also a spatial variability in the proxy-derived estimates of warming, which could be accounted for in three ways; (1) the larger estimates of warming (e.g. ODP 690 and New Jersey sites) are closer to the global response, whilst other sites are biased downwards by poorpreservation, dissolution or incomplete records; (2) the larger estimates of warming are biased upwards, for example by salinity effects on δ 18 O; or, (3) the differential reflects the genuine spatial pattern of warming at the PETM. Given the consistency of warming estimates between inorganic and organic proxies at New Jersey, we consider the second explanation unlikely. Obtaining higher quality proxy records in regions that, in this analysis, show small magnitudes of warming, would help differentiate between explanations 1 and 3.
In the absence of such new data, here we compare the proxy and model-derived geographical patterns of warming (Fig. 10). The ranking of sites by the mean size of the proxy-derived PETM temperature anomaly, is largely consistent with site-to-site variation in model-derived warming estimates. At the level of accuracy of the data, this is best de- South Atlantic/Arctic sites versus N. Atlantic and Southern Ocean sites in Fig. 10). This spatially heterogeneous response to GHG-forcing is to be expected from our general understanding of the climate system (Christensen et al., 2007) and is similar to other climate models of this period (Winguth et al., 2010;Huber and Caballero, 2011). That the climate models produce a very similar spatial pattern of warming as the proxy records is a critical point. It questions the assumption that the highest recorded temperature anomalies are the most representative of the global response, with other records biased downwards by stratigraphic incompleteness or proxy data quality. Instead, this model-data comparison implies at least some genuine spatial heterogeneity in the patterns of warming. Palaeotemperature proxy data acquisition in regions where the model predicts the smallest degrees of warming would test this hypothesis. In Fig. 11, we further investigate the consistency of the spatial patterns of warming between proxy and model estimates, and the implications for estimates of global mean warming at the PETM warming. Here we continue to assume a constant sensitivity of the HadCM3L model, at all locations analysed, beyond the 6× simulation. With this assumption we ask at what level of forcing, and hence of global mean temperature change in the model, are the sum of square residuals minimised between the proxy-and model-estimates of warming. This is undertaken with three sub-sets of the near-sea surface proxy temperature data presented here, first with all the data, second without planktic foraminifera δ 18 O from sites most likely to be compromised by diagenetic recrystallisation (DSDP Site 527 and ODP Sites 865 and 1209) and finally without any data, δ 18 O or Mg/Ca, from these same sites. The striking feature of this analysis is that the inclusion or exclusion of this datathose with the smallest proxy-derived temperature anomalieshas almost no impact on the best-fit model forcing and associated mean modelled global warming. The best-fit model warming to all of the data in this analysis is 4.3°C, with modification in data selection only having an impact of ±0.05°C. It might be thought that the exclusion of sites with a small proxy temperature anomaly would increase global mean warming in the best-fit climate model analysis. Instead, consideration of the site-by-site proxy and model temperature anomalies (Fig. 11 middle row) and the proxy-model residuals ( Fig. 11 bottom row), shows that there is very similar spread of proxy-derived temperature anomalies around the model-derived anomalies at both the low end of responsiveness (low proxy anomaly and low modelled anomalies) and the high end of responsiveness (high proxy anomaly and high modelled anomaly). The exclusion of sites with low proxy data temperature anomalies thus has very little effect on the overall best fit. Further, on this basis of this analysis and the available data, it is very difficult to argue for a global mean temperature anomaly in excess of 5°C.

Proxy records of intermediate water warming
In addition to the SST estimates presented above, an indication of the magnitude of PETM temperature anomaly can be derived from intermediate to deep-water temperatures based on benthic foraminiferal δ 18 O and Mg/Ca. These records can be compromised by the effects of deepocean carbonate dissolution, incomplete recovery of the PETM and the ecological exclusion of calcareous benthic foraminifera at the peak of the event (Thomas, 1998(Thomas, , 2003. Depth-dependent variations in the severity of dissolution and foraminifera assemblage changes can introduce a systematic depth bias to the magnitude of the recorded PETM temperature anomaly (McCarren et al., 2008). Spatial variability in the extent of carbonate dissolution (Panchuk et al., 2008a(Panchuk et al., , 2008bZeebe et al., 2009) may introduce a similar systematic bias. Estimates of warming based on benthic foraminiferal Mg/Ca ratios may also be  Stott et al. (1996); Stott (1990, 1991); Kelly et al. (2005); Thomas et al. (2002); Nunes and Norris (2006) 11.9 0.69 16 18.1 0.27 2 6.2 All measured on Nuttalides truempyi unless stated otherwise. NB: Used O. umbonatus calibration of Lear et al. (2002) for all benthic Mg/Ca data. Mixed benthic Mg/Ca has already been normalised to O. umbonatus by Tripati and Elderfield (2005).
damped by 1-2°C due to reductions in the carbonate saturation state (Tripati and Elderfield, 2005). The PETM temperature anomaly was estimated from the benthic foraminiferal proxy data using the same methodology as for the SST proxy estimates (Table 2; Fig. 12). The minimum temperature anomaly is recorded at the deep, northwest Pacific site ODP 1209 on Shatsky Rise (~2°C), with relatively small warming also recorded in the eastern equatorial Pacific (ODP 1220). Sedimentation rates are very low at ODP 1209, making it unlikely that peak-PETM warming is fully recorded at this site. Both of these sites were also located at palaeodepths of greater than 2500 m Lyle et al., 2002) and together may reflect a genuine damped response of Pacific intermediate and deep ocean to global warming. In contrast, sites on the New Jersey shelf, in waters depths of b150 m (Gibson et al., 1993;van Sickel et al., 2004), show some of the highest temperature anomalies. At such shallow palaeodepths, benthic foraminifera temperature records will be dominated by the seasonal dynamics of stratification and mixing of shelf-waters, rather than the wider patterns of warming in North Atlantic intermediate waters.
From the North Atlantic to the Southern Ocean the majority of sites, across latitudes, record warming on the order of~5-6°C, including the North Atlantic (ODP 1051), the tropical South Atlantic (ODP 1258) and the Southern Ocean (ODP 690 and 738). This intermediate water warming is slightly higher than the estimate of global mean surface temperature warming presented above. This is, however, likely driven by enhanced warming in the high latitude regions from which these intermediate waters are sourced. In this context, 5-6°C intermediate water warming is actually remarkably consistent with the~6°C warming of the Southern Ocean simulated with a doubling of GHG forcing. In two regions, the northeast and southeast Atlantic, the estimates of warming from multiple sites are less coherent. First, in the northeast Atlantic, the temperature anomalies from DSDP Site 401 and Alamedilla are~2°C less that that recorded at DSDP Site 549 (~5-6°C). All three of these sites are estimated to have been in bathyal depths of 2000-3000 m during the PETM (Shipboard Scientific Party, 1979;Shipboard Scientific Party, 1985;Lu et al., 1996). None of these sites have ideal records: at DSDP 549 there is~1 m of carbonate dissolution during the event (Nunes and Norris, 2006); foraminiferal preservation at Alamedilla is poor, such that planktic foraminifera record almost no δ 18 O excursion across the PETM (Lu et al., 1996); and rotary drilling at DSDP Site 401 likely prevented full recovery of the PETM interval. DSDP Site 549 was excluded from the Nunes and Norris (2006) compilation because of anomalously low δ 13 C values attributed to diagenetic overprinting. We suggest that the best available constraint on PETM warming of northwest Atlantic intermediate waters is thus a minimum estimate based on DSDP Site 401 (ΔT N3°C).
In the southeast Atlantic, there is significant variability in benthic foraminifera δ 18 O anomalies between sites drilled in the vicinity of Walvis Ridge. PETM records from DSDP Leg 74 Sites, 525 (modern water depth 2467 m) and 527 (water depth 4435 m) are relatively low resolution but both show clear negative δ 18 O anomalies, equivalent to~5°C of warming. In contrast, large negative δ 18 O (and δ 13 C) anomalies, Bottom: estimates of pre-and peak-PETM temperatures including 95% confidence intervals based on sampling density and variability within each stratigraphic interval. Middle: PETM temperature anomaly derived from these proxy data plotted in rank order and colour coded by geographical region, plotted with 95% confidence intervals as above. Simple mean of all of these benthic foraminifera-derived temperature anomalies (5.4°C) shown for reference. Top: palaeolatitude of each of these sites. equivalent to~9°C of warming, are recorded in the benthic foraminifera Nuttallides trumpeyi and Oridorsalis umbonatus at the more recently drilled ODP 1263 (water depth 2712 m). These anomalies exceed the δ 18 O response at intermediate water far-field sites to the south (ODP 690) and north (ODP 1258) as well as the Leg 74 sites nearby. McCarren et al. (2008) attribute this strong response to a potential change in intermediate water sources from the southern high latitudes to warm, nutrient-rich, oxygen depleted waters. Such a change is supported by the circulation responses of HadCM3L at high-levels of GHGforcing  and independent redox proxies from the South Atlantic (Chun et al., 2010). McCarren et al. (2008) show that this large isotope anomaly can be damped with increasing palaeodepth and intensity of carbonate dissolution. It is likely that the smaller anomaly recorded at the older, single hole at DSDP Site 525, at a similar water depth to ODP 1263, is due to incomplete recovery of the PETM, given similar problems at ODP Site 1263 that were only recognised with the drilling of multiple holes (Shipboard Scientific Party, 2004).
In the pre-PETM interval, benthic foraminiferal δ 18 O and Mg/Ca data also provide a constraint on the minimum (cold season) surface temperatures in the principal regions of high latitude deep-water formation. All δ 18 O-based pre-PETM intermediate water temperature estimates from the Atlantic sector of the Southern Ocean (ODP 690) and South Atlantic (DSDP 525, 527, ODP 1263) are all 12°C ± 1°C. Given the pattern of austral winter mixed layer depths in the 6 × CO 2 model simulation (Fig. 13), with the major locations of deep-water formation in the Weddell Sea and Atlantic-Indian Ocean sectors of the Southern Ocean, these benthic proxy temperatures (~12°C) should reflect winter-time SSTs in these regions. For the 6 × CO 2 model runs these austral winter SSTs are of the same order or slightly warmer than reconstructed intermediate water temperatures in the Atlantic and Southern Ocean. The agreement of both modelled winter surface water and proxy intermediate water temperatures in the Atlantic and Indian Ocean sectors of the Southern Ocean and between modelled mean annual and proxy (δ 18 O) surface water temperatures from ODP Site 690 is a positive step forward in addressing problems of highlatitude warm climate model-proxy mismatches. Most notable in these model simulations is the N10°C gradient in winter temperatures between the Indo-Atlantic and the (higher latitude) Pacific sectors of the Southern Ocean. The most significant model-data mismatches, including the GDGT-data from ODP Site 1172 discussed herein and the recently published early Eocene terrestrial temperature estimates from the Wilkes Land margin (Pross et al., 2012;Bijl et al., 2013) are both located in this overly cool Pacific-sector in the model representations of the Southern Ocean.

Conclusion
This study presents a detailed reassessment of the existing marine palaeotemperature estimates across the PETM hyperthermal event. We are able to apply a consistent and up-to-date set of calibrations to all published data. This provides the basis for a more rigorous comparison of the PETM temperature anomaly between records published over the past two decades, from widely distributed locations. On the basis of these recalculated PETM palaeotemperature estimates, we find a range of estimated near-sea surface PETM temperature anomalies, from almost no change (δ 18 O-based estimates from ODP Site 865) up to a maximum estimated warming of 9°C (TEX 86 L at Wilson Lake and Acarinina δ 18 O at ODP Site 690). Such face value estimates of PETM temperature change must be treated with the utmost caution. All of these proxybased estimates of the PETM temperature anomaly are subject to multiple uncertainties and potential biases. Three critical uncertainties in proxy temperature estimations are (1) the impact of foraminiferal carbonate recrystallisation on δ 18 O and Mg/Ca temperature estimates, especially in deep-ocean sediments; (2) the unconstrained variability in local seawater δ 18 O across the event, critical for δ 18 O-derived palaeotemperatures; and, (3) ongoing significant discrepancies in deep-time, warm-climate SST estimates derived from different GDGT-based proxies and calibrations.
The fundamental contention of this study is that the variability in proxy estimates of PETM warming between different sites contains some component of genuine spatial variability in warming at the PETM. Here we have tried to deconvolve the relative influence of proxy bias and the genuine patterns of warming by comparing the proxy-derived pattern of PETM warming to the location-specific response of early Eocene AOGCM simulations to greenhouse-gas forced warming. This comparison clearly shows that the first-order spatial pattern of warming is consistent between the proxy-derived and model-derived temperature anomalies. Further, that taking all proxy records into account, it appears unlikely that the global mean surface temperature anomaly across the PETM is in excess of 5°C. Although this estimate of mean global warming at the PETM may appear to be relatively modest in relation to some of the larger PETM temperature anomalies recorded at individual sites, as an estimate of the global surface temperature response it represents a remarkably large climate anomaly. For climate perturbations occurring on the~10 kyr timescale it is, to date, unequalled in the Cenozoic. If this assessment of the global temperature anomaly proves to be relatively robust then, with refined controls on the PETM carbon cycle perturbation and carbon input, a more accurate estimate of early Paleogene climate sensitivity is within reach.