4.1.1 Theoretical background of oxygen isotopes

Oxygen isotope paleothermometry is based on the temperature-dependent fractionation of the stable isotopes of oxygen (¹⁶O and ¹⁸O, expressed as δ¹⁸O) into biogenic calcite from its growth medium, typically seawater. Pearson (2012) provides more details on the development of the δ¹⁸O proxy. The temperature dependence of δ¹⁸O is a combination of thermodynamic and vital effects (e.g., Kim and O'Neil, 1997; Bemis et al., 1998). Most δ¹⁸O paleotemperature reconstructions use the shells of foraminifera, a group of protists that live and calcify either as plankton in the ocean's surface mixed layer, thermocline or sub-thermocline depths, or as benthos on or just below the sea floor. Temperature is calculated by an empirical calibration to quantify the fractionation between the δ¹⁸O of ancient seawater and biogenic calcite (Bemis et al., 1998).

4.1.2 Strengths of oxygen isotopes

Foraminiferal oxygen isotopes have been the primary proxy for reconstructing ocean temperatures spanning the past ∼120 million years, largely due to the low analytical cost, small sample size requirements and relative ease of measurement. Analytical uncertainty on δ¹⁸O measurements is typically small, ±0.1‰ (equivalent to <1 ^∘C). Moreover, the theoretical basis for temperature-dependent fractionation is firmly tied to field- and laboratory-based relationships between foraminifer test δ¹⁸O values and temperature (e.g., Kim and O'Neil, 1997; Bemis et al., 1998; Lynch-Stieglitz et al., 1999). The range of planktic foraminiferal depth habitats also allows for the reconstruction of water column profiles and thermocline structure (e.g., Birch et al., 2013; John et al., 2013), which can be compared to modeled upper ocean structure (e.g., Lunt et al., 2017). The relatively short life span of planktic foraminifera, 2 weeks to 1 month for most modern species, can also provide constraints on paleoseasonality (Pearson, 2012). Further, based on the assumption that deep-water formation is largely focussed at high (subpolar) latitudes throughout the Cenozoic (Cramer et al., 2011), benthic foraminiferal δ¹⁸O values can (arguably) be used as a SST proxy in these areas.

4.1.3 Weaknesses of oxygen isotopes

The dissolution and subsequent replacement of primary biogenic calcite by inorganic calcite (recrystallization/diagenesis) in pelagic, carbonate-rich sediments during early diagenesis is known to shift planktic foraminiferal calcite to higher values (Schrag et al., 1995), with the effect that most low- and mid-latitude δ¹⁸O-derived Paleogene SSTs from deep-ocean carbonate-rich successions may be systematic underestimates (Pearson et al., 2001, 2007; Tripati et al., 2003; Sexton et al., 2006; Pearson and Burgess, 2008; Kozdon et al., 2011; Edgar et al., 2015). The effect of seafloor recrystallization and diagenesis on SST estimates will be proportionally less significant in areas where cooler surface waters are closer to deep ocean temperatures, which may be the case in some Paleogene high-latitude or upwelling regions, and is considered to be insignificant for benthic foraminifera because temperatures during early diagenesis are very close to growth temperatures (Schrag, 1999; Edgar et al., 2013; Voigt et al., 2016). During diagenesis, foraminifer tests may also become overgrown and/or infilled with calcite precipitated from sediment pore waters. Where marine sediments are exposed on land, those secondary precipitates can incorporate oxygen from isotopically light meteoric waters and hence yield artificially warm apparent temperatures. For these reasons, temperature estimates from foraminiferal δ¹⁸O are only considered to be reliable where a good state of preservation has been confirmed by scanning electron microscope (SEM) examination and illustration. Well-preserved (non-recrystallized or glassy) late Paleocene to early Eocene foraminifera have been reported from low-permeability clay-rich facies in shallow marine or hemipelagic settings from Tanzania (Pearson et al., 2004; 2009; Sexton et al., 2006), the New Jersey margin and California (Zachos et al., 2006; John et al., 2008; Makarova et al., 2017), New Zealand (Hollis et al., 2012), and Nigeria (Frieling et al., 2017).

One key assumption for δ¹⁸O-based temperatures is that foraminifera precipitate their tests in isotopic equilibrium with seawater. However, foraminiferal physiology (e.g., respiration, metabolism, biomineralization, photosymbiosis) and ecology (e.g., depth migration during life cycle, seasonality), often termed “vital effects”, commonly lead to isotopic offsets that can bias temperature reconstructions (e.g., Urey et al., 1951; Birch et al., 2013). The biology and ecology of foraminifera may also vary in either time or space. For instance, the depth habitat of a species may change in response to rapid environmental change, or due to evolution within a lineage. Foraminifera that host algal photosymbionts may also be subject to bleaching events, which can substantially alter the test micro-environment within which calcification, and the associated isotopic fractionation, occurs (Wade et al., 2008; Edgar et al., 2013; Luciani et al., 2016; Si and Aubry, 2018). In some sites there is evidence that planktonic foraminifera and other eukaryotes disappeared from the record during peak PETM warming, possibly because environmental conditions became too extreme (Aze et al., 2014; Frieling et al., 2018), in which case peak conditions would go unrecorded by this proxy. Enhanced dissolution in the PETM would have the same effect.

Calculating ancient SST from foraminiferal δ¹⁸O requires an estimation of the oxygen isotopic composition of seawater (δ¹⁸O_sw) at the time of precipitation. This is not straightforward because δ¹⁸O_sw varies spatially in the surface ocean, largely following patterns of salinity (Zachos et al., 1994; Rohling, 2013), and temporally due to changes in the cryosphere (Broecker, 1989; Cramer et al., 2009). Large and spatially variable changes in the intensity of the hydrological cycle are inferred across the PETM (Bowen et al., 2004; Zachos et al., 2006; Pagani et al., 2006; Carmichael et al., 2017); hence it is unlikely that δ¹⁸O_sw at any single location remained constant through the Paleocene–Eocene interval. Continental margin settings, where foraminifera are typically well preserved, may be particularly sensitive to changes in δ¹⁸O_sw related to the hydrological cycle.

Finally, culture studies and field observations demonstrate that the δ¹⁸O value of foraminiferal calcite decreases as the pH of the culture medium increases (Spero et al., 1997; Bijma et al., 1999; Zeebe, 1999, 2001; Russell and Spero, 2000). This “pH effect” can influence ancient SST reconstructions if seawater pH varied rapidly or was significantly different from today (Zeebe, 2001). Studies suggest that early Paleogene surface ocean pH was as much as ∼0.5 units lower than today (Penman et al., 2014; Anagnostou et al., 2016; Gutjahr et al., 2017), which implies that background δ¹⁸O-based SST estimates for this time interval could be ∼3 ^∘C too low, or even more so for the PETM when pH may have declined by a further ∼0.3 units (Uchikawa and Zeebe, 2010; Aze et al., 2014).

4.1.4 Recommended methodologies for oxygen isotopes

Here we outline our recommendations for generating δ¹⁸O data from fossil foraminifera, and for converting δ¹⁸O values into temperature estimates. We have compiled available planktic foraminiferal δ¹⁸O data from 10 DSDP, ODP and IODP sites and 9 onshore sections (Supplement Data File 3). Using the methods outlined below, we have calculated SSTs and compiled a summary of proxy-specific SST estimates for each time slice (LP, PETM, EECO). SST estimates are based on species that are inferred to have inhabited near-surface waters or the mixed layer. Data for deeper-dwelling (thermocline) planktic species are also included in Supplement Data File 3. For compilations of benthic foraminiferal δ¹⁸O, see Zachos et al. (2008), Cramer et al. (2009, 2011) and Westerhold et al. (2017, 2018), although these data are not considered here.

Depth ecologies of all Paleogene species, as inferred mainly from stable isotope evidence, have been compiled by Aze et al. (2011). The principal groups used for sea surface temperature reconstruction are open-ocean mixed-layer species with or without algal symbionts and high-latitude species (ecogroups 1, 2 and 5 of Aze et al., 2011). The main mixed-layer genera for the DeepMIP time slices are Morozovella and Acarinina but other relevant groups are Igorina, Planorotalites and Pseudohastigerina. Different species of Morozovella and Acarinina may exhibit consistent offsets indicating a degree of depth stratification within the mixed layer and upper thermocline, possibly related to sinking at the time of reproduction. Hence SST reconstructions from different “mixed layer” species may vary by up to several degrees Celsius. Because of this, combining various species of the same genus in analyses (e.g., measuring Acarinina spp.) is likely to produce underestimates of SST. In this compilation, we derive average SST estimates from “mixed layer” species for each site in a given time slice. This is to promote consistency and aid intersite comparisons but is considered to be a conservative approach to estimating SSTs. For future work we recommend more detailed investigation of the various important species and identification of those species which most faithfully record the warmest upper ocean mixed layer conditions.

To expand on the δ¹⁸O database, we propose a three-pronged approach. First, analysis of new and classic carbonate-rich sites containing recrystallized foraminiferal tests should involve the novel technique that uses a secondary ion mass spectrometer (SIMS). Pioneering SIMS studies by Kozdon et al. (2011, 2013) suggest that areas furthest from the test exterior are less susceptible to diagenetic overprinting. These areas yield SSTs up to ∼8 ^∘C warmer than conventional analyses from the same sample and are in better agreement with δ¹⁸O-based SSTs from glassy tests (Kozdon et al., 2011). On a cautionary note, a recent study by Wycech et al. (2018) has reported an offset between SIMS and traditional isotope-ratio mass spectrometer (IRMS) analyses for modern foraminifera, with SIMS δ¹⁸O values being ∼0.9‰ lower. Further study is needed to determine if this offset also affects fossil foraminifera. The SIMS technique provides hope for recovering reliable SSTs from recrystallized foraminiferal tests but it is time intensive and is not a practical approach for the analysis of all samples. Thus, our second recommendation is that whole specimen analyses are undertaken together with the SIMS analysis, in order to constrain the magnitude of diagenetic bias on IRMS δ¹⁸O values. For lithologically uniform sediments, quantitative estimates of this diagenetic bias in a few widely spaced samples could provide calibration points for higher resolution data generated by conventional whole shell methods. Whole specimen isotopic analyses should be species-specific and use a prescribed size fraction (e.g., 250–300 µm or 300–355 µm) to minimize variations in vital effects (Birch et al., 2013). Our third recommendation is that new sites that contain glassy foraminifera are sought out to provide the material for whole specimen analysis, as well as SIMS analysis of selected samples to further validate the two methods. All samples included within the database are categorized as either glassy or recrystallized, using the criteria of Sexton et al. (2006) and Pearson and Burgess (2008), as a guide to reconstructed SST reliability.

Empirically derived δ¹⁸O–temperature calibrations can differ by several degrees Celsius, , although offsets decrease with increasing temperature (Bemis et al., 1998; Pearson, 2012). Use of multiple equations may capture a range of plausible temperature values but our preference is the calibration of Kim and O'Neil (1997) for inorganic calcite, which is appropriate for the Paleogene because it is based on inorganic calcite precipitated from water temperatures between 10 and 40 ^∘C. Both epibenthic and asymbiotic planktic foraminifera yield values close to the resulting regression (Bemis et al., 1998; Costa et al., 2006). Field or laboratory studies only include calcite precipitated up to 30 ^∘C (e.g., Bemis et al., 1998; Lynch-Stieglitz et al., 1999), yet Paleogene SSTs likely fall close to or above the upper limit of these studies. The recommended calibration, as modified by Bemis et al. (1998), derives temperature from Eq. (1):

where T is the water temperature in ^∘C and δ¹⁸O_C and δ¹⁸O_sw are the δ¹⁸O of calcite (‰ VPDB) and ambient seawater (‰ VSMOW), respectively. This equation may overestimate SST for symbiont-hosting planktic foraminifera by ∼1.5 ^∘C based on the consistent 0.3 ‰ offset observed between Orbulina universa grown under high-light compared to low-light conditions (Spero and Williams, 1988; Pearson, 2012). This offset is likely caused by algal photosymbionts modifying the pH in the calcifying microenvironment (Zeebe et al., 1999). Orbulina universa is inferred to share a similar ecology to the dominant Eocene genera Morozovella and Acarinina typically analyzed for Paleogene SSTs (e.g., Shackleton et al., 1985; D'Hondt et al., 1994) and is currently the best modern analogue for which calibration data are available. However, we do not recommend applying a symbiont correction to these genera because of uncertainties in photosymbiont activity levels in Paleogene conditions.

It is assumed that these δ¹⁸O calibrations are insensitive to evolving seawater chemistry, unlike some trace element proxies (e.g., Evans and Müller, 2012). These calibrations do, however, require an appropriate estimate of δ¹⁸O_sw at the time of test formation. Changes in global mean δ¹⁸O_sw are largely driven by continental ice volume and isotopic composition. Various estimates have been proposed for the adjustment of the global ocean value for ice-free conditions of the early Paleogene: Shackleton and Kennett's (1975) initial estimate of −1.00 ‰ can be compared with the more recent estimates of −1.11 ‰ of L'Homme et al. (2005) and −0.89 ‰ of Cramer et al. (2011). Differences arise from uncertainty surrounding the mean δ¹⁸O value of the modern ice caps (see Pearson, 2012 for discussion). Pending resolution of these discrepancies we here use the value −1.00 ‰ for early Paleogene δ¹⁸O_sw under ice-free conditions, which is the mean of the L'Homme et al. (2005) and Cramer et al. (2011) estimates and identical to the estimate of Shackleton and Kennett (1975) upon which many historical temperature estimates have relied. A correction for the local effects of salinity and hydrology on δ¹⁸O_sw should also be incorporated into SST estimates where possible. Zachos et al. (1994) developed a correction using present-day δ¹⁸O_sw latitudinal gradients, which is widely used despite acknowledged shortcomings. The correction is applied globally but is only based on Southern Hemisphere data and does not incorporate local freshwater runoff effects on continental margins or the influence of boundary currents. It also ignores expected changes in the relationship between latitude and δ¹⁸O_sw through time. Isotope-enabled climate models have been used to generate predictions of early Paleogene surface ocean δ¹⁸O_sw distributions, following an assumption of mean average δ¹⁸O_sw conditions (e.g., Tindall et al., 2010; Roberts et al., 2011), but the δ¹⁸O_sw fields generated are dependent on the surface climatology of the model concerned (Hollis et al., 2012). The use of model-derived δ¹⁸O_sw to generate proxy estimates of SST also introduces a problematic model-dependency within the proxy dataset. Here, we chose to update the approach of Zachos et al. (1994) using the more highly resolved global gridded (1^∘ × 1^∘) dataset of modern δ¹⁸O_sw produced by LeGrande and Schmidt (2006). We calculated δ¹⁸O_sw for each site by relating the site's paleolocation to 10^∘ latitudinal bins in the modern dataset (see Supplement Data File 3). Median values and 95 % confidence intervals were calculated using paleolocations derived from both paleomagnetic and mantle-based reference frames (see Sect. 7.2). We encourage future studies to improve the empirical fit to the modern spatial variability in δ¹⁸O_sw but also recognize that the modern system only provides a first-order estimate of spatial δ¹⁸O_sw patterns in deep time. With improved data coverage, paired δ¹⁸O–Mg∕Ca and δ¹⁸O–Δ47 analyses hold promise for direct reconstructions of ancient δ¹⁸O_sw variability.

We do not recommend applying a pH correction to δ¹⁸O data because of significant uncertainty in the magnitude of the effect below pH 8 (Uchikawa and Zeebe, 2010), which is inferred to be the upper limit for most of our Paleogene records (Anagnostou et al., 2016; Penman et al., 2014; Gutjahr et al., 2017). Further, the pH–δ¹⁸O_C sensitivity of asymbiotic planktic and benthic foraminifera is not well known (e.g., Anagnostou et al., 2016; McCorkle et al., 2008; Mackensen, 2008).

4.2.1 Theoretical background of Mg∕Ca ratios

The sensitivity of foraminiferal Mg∕Ca ratios to temperature has a basis in thermodynamics (Lea et al., 1999) through the exponential temperature dependence of any reaction for which there is an associated change in enthalpy. However, most species do not conform to theoretical calculations (Rosenthal et al., 1997), being characterized by a temperature sensitivity of Mg incorporation around 2–3 times greater than that for inorganic calcite. The reasons for this remain elusive (e.g., Bentov and Erez, 2006), necessitating empirical calibration of foraminiferal Mg∕Ca to temperature (e.g., Anand et al., 2003).

4.2.2 Strengths of Mg∕Ca

As with δ¹⁸O paleothermometry, Mg∕Ca paleothermometry can be applied to both benthic and planktic foraminifera and be used to constrain the past thermal structure of the water column (Tripati and Elderfield, 2005). The long residence time of Ca and Mg in the ocean means that the Mg∕Ca of seawater (Mg∕Ca_sw) can be treated as constant over short timescales (<10⁶ years), which removes one source of uncertainty in calculating relative changes in temperature (e.g., Zachos et al., 2003; Tripati and Elderfield, 2004). A major advantage of Mg∕Ca is its use in paired measurements with δ¹⁸O on the same substrate, which allows for the deconvolution of δ¹⁸O_sw and temperature effects on measured foraminiferal δ¹⁸O. As Mg∕Ca measurements are increasingly routine in many laboratories, are relatively inexpensive to produce, and have low sample-size requirements, high-resolution time series are readily achievable. Technological developments in laser-ablation inductively coupled plasma mass spectrometry (ICPMS) and other spatially resolved methodologies, such as electron probe microanalysis (EPMA) and SIMS, allow for Mg∕Ca measurements on or within a single individual test, providing new information on foraminiferal ecology and short-term environmental variability (e.g., Eggins et al., 2004; Evans et al., 2013; Spero et al., 2015). Planktic foraminiferal Mg∕Ca may also be more robust to shallow burial diagenetic recrystallization than δ¹⁸O, based on values measured in the same material (Sexton et al., 2006).

4.2.3 Weaknesses of Mg∕Ca

A long-standing challenge for the deep-time application of the Mg∕Ca temperature proxy is that the seawater Mg∕Ca ratio (Mg∕Ca_sw) influences shell Mg∕Ca, but there is still debate over how Mg∕Ca_sw has varied through time (Horita et al., 2002; Coggon et al., 2010; Broecker and Yu, 2011; Evans and Müller, 2012). Nonthermal influences on foraminiferal Mg∕Ca ratios can also be difficult to account for, including pH, bottom water carbonate saturation state and sample contamination (Barker et al., 2003; Regenberg et al., 2014; Evans et al., 2016b). Determining the reliability of Mg∕Ca paleotemperatures requires an understanding of these challenges, and the use of independent paleoenvironmental proxies where available, such as boron isotopes to constrain carbonate system parameters (Anagnostou et al., 2016).

Impact of foraminiferal preservation on the Mg / Ca paleothermometer. Many foraminifer tests from deep ocean sediments are affected by diagenetic alteration (Pearson et al., 2001; Edgar et al., 2015), with the type and extent of alteration controlled by original test morphology, taphonomic processes, the characteristics of the host sediment and burial history. For example, planktic foraminifera tests can undergo immediate post-mortem or post-gametogenic dissolution as they sink into deeper, less carbonate saturated waters (Brown and Elderfield, 1996). In general, MgCO₃ is more soluble than CaCO₃ such that dissolution tends to decrease foraminiferal Mg∕Ca. This may result in artificially low Mg∕Ca-based temperature estimates if unaccounted for (e.g., Rosenthal and Lohmann, 2002; Regenberg et al., 2014; Fehrenbacher and Martin, 2014), although not in all cases (e.g., Sadekov et al., 2010; Fehrenbacher and Martin, 2014). The effects of dissolution can be minimized by selecting sites with relatively shallow paleodepths, above the calcite lysocline (e.g., <2000 m). Foraminiferal shells can also be subject to diagenetic overgrowths of various mineral phases, including oxy-hydroxides and authigenic carbonates, depending on seafloor and sub-seafloor conditions (e.g., Boyle, 1983). Recent work suggests that some textural recrystallization of planktic foraminiferal tests may occur in semiclosed chemical conditions, potentially allowing original geochemical signals to be retrieved using microsampling techniques (e.g., Kozdon et al., 2011, 2013, see Sect. 5.2.5).

Challenges for benthic foraminiferal Mg / Ca paleothermometry. The temperature sensitivity of Mg incorporation into benthic foraminiferal calcite varies between species, necessitating genus- or species-specific temperature calibrations (Lear et al., 2002). Fortunately, some extant species are common throughout the Cenozoic (e.g., Oridorsalis umbonatus – Lear et al., 2000) and offer a means to develop calibrations for coeval extinct species. There is also no consensus as to whether relationships between benthic Mg∕Ca and temperature are best described by linear or exponential fits (Cramer et al., 2011; Evans and Müller, 2012; Lear et al., 2015); we recommend that calibrations are applied with caution where Mg∕Ca ratios are outside the range for which temperatures have been empirically determined.

Present-day benthic foraminifera appear to increase their discrimination against magnesium when calcifying in waters with very low carbonate ion saturation state ($\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$); a relationship that has been empirically quantified for some species (Elderfield et al., 2006; Rosenthal et al., 2006). Measurements of benthic foraminiferal B∕Ca in tandem with Mg∕Ca can be used to identify temporal variations in $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ (Yu and Elderfield, 2007) and may provide a means of correcting for this secondary effect, although it can be difficult to identify the threshold for a $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ influence within downcore records (e.g., Lear et al., 2010). Mg∕Li ratios may have a more consistent empirical relationship with temperature than Mg∕Ca (Bryan and Marchitto, 2008). However, limited understanding of how the Mg∕Li seawater ratio has varied over geological time means that this proxy can only be used as a guide to relative temperature change in deep time studies (Lear et al., 2010). The calcification of endobenthic foraminifera species within buffered porewaters may make them relatively insensitive to variations in bottom water $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ (Zeebe, 2007; Elderfield et al., 2010). However, the saturation state of porewaters is dependent on many factors and likely also varies through time (Weldeab et al., 2016). We recommend that tandem trace metal ratios that are sensitive to carbonate saturation state (e.g., B∕Ca, Li∕Ca) are examined to assess downcore $\mathrm{\Delta }{\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$ variations, even in infaunal records (e.g., Lear et al., 2010, 2015; Mawbey and Lear, 2013).

Evidence from multiple proxies indicates that early Paleogene Mg∕Ca_sw was significantly lower than the modern value (Horita et al., 2002; Coggon et al., 2010; Lear et al., 2015; Evans et al., 2018a). Correcting benthic Mg∕Ca data for secular changes in seawater chemistry is complicated by the fact that the benthic foraminiferal magnesium partition coefficient (D_Mg = Mg∕Ca${}_{\mathrm{CALCITE}}/{\mathrm{Mg}/\mathrm{Ca}}_{\mathrm{sw}})$ decreases with increasing Mg∕Ca_sw according to a power function (Ries, 2004; Hasiuk and Lohmann, 2010; Evans and Müller, 2012; Lear et al., 2015). Moreover, the sensitivity of shell Mg∕Ca to changes in Mg∕Ca_sw, specifically the curvature of the power relationship that relates the two factors, appears to be genus-specific. Lear et al. (2015) and Evans et al. (2016b) argued for a low sensitivity of shell Mg∕Ca to Mg∕Ca_sw for Oridorsalis and the endobenthic genus Uvigerina. In contrast, a higher sensitivity is thought to characterize the epibenthic genus Cibicidoides/Cibicides (Evans et al., 2016b), which is also widely used in paleoclimate studies. Evaluating this aspect of benthic foraminiferal geochemistry is challenging and has yet to be assessed in other widely utilized species. As such, best practice would be to report to what extent Mg∕Ca temperatures would change when considering the uncertainty in the slopes of these seawater–shell Mg∕Ca relationships.

Challenges for planktic foraminifera Mg / Ca. The relationship between planktic Mg∕Ca ratios and temperature is species- or group-specific (e.g., Regenberg et al., 2009), such that species-specific calibrations should be used whenever possible. Nonetheless, many planktic foraminifera conform to a broader Mg∕Ca–temperature relationship (Elderfield and Ganssen, 2000; Anand et al., 2003), which is one method by which modern calibrations can be applied to extinct taxa. As for benthic foraminifera, when working with pre-Pleistocene samples the control exerted by changes in seawater elemental chemistry over geological time must also be considered. Culture experiments in modified seawater demonstrate not only that Mg∕Ca_sw impacts planktic foraminifera shell chemistry (Delaney et al., 1985), but also that the slope of the Mg∕Ca–temperature relationship may be sensitive to Mg∕Ca_sw (Evans et al., 2016b). In addition, several other nonthermal controls on Mg∕Ca should be considered when interpreting data from planktic foraminifera. Culture and core-top studies demonstrate a relatively minor salinity effect (Kısakürek et al., 2008; Hönisch et al., 2013) in which, for example, a 2 PSU salinity increase results in a temperature overestimate of ∼1 ^∘C. In contrast, the carbonate system has been shown to have a large influence on Mg∕Ca in several species (Lea et al., 1999; Russell et al., 2004; Evans et al., 2016a). Lower pH and/or [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] results in higher shell Mg∕Ca; for example, a 0.1 unit pH decrease results in a temperature overestimate of ∼1 ^∘C. The effect of the carbonate system on planktic Mg∕Ca has been identified in sediment-trap as well as culture studies (Evans et al., 2016a; Gray et al., 2018), with Gray et al. (2018) demonstrating that the widely used Mg∕Ca–temperature sensitivity of ∼9 % per degree Celsius in Globigerinoides ruber is an artifact of the covariance of temperature and pH, through the temperature effect on the dissociation constant of water. Specifically, the pH of seawater decreases with increasing temperature, resulting in an increase in the incorporation of Mg into planktic foraminiferal calcite due to both processes (see Evans et al., 2018b). The secondary pH effect accounts for around one-third of the observed increase in Mg∕Ca ratios at higher temperatures, leaving a primary Mg∕Ca “temperature only” sensitivity of 6 % per degree Celsius (Gray et al., 2018), which is significantly lower than that widely utilized. Whilst these factors may be accurately accounted for in the recent geological past (Gray and Evans, 2019), it is challenging to account for them in deep time because high-resolution pH records are scarce (see below for detailed recommendations).

4.2.4 Recommended methodologies for Mg∕Ca

Here we outline our recommendations for generating Mg∕Ca data from fossil foraminifera and for converting these data into temperature estimates. We have compiled available planktic Mg∕Ca data from five DSDP and ODP sites and six onshore sections and, using the methods outlined below, calculated SSTs and associated uncertainties for the late Paleocene and early Eocene (Supplement Data File 4). SST estimates are based on species that are inferred to have inhabited near-surface waters or the mixed layer. Data for deeper-dwelling (thermocline) planktic species are also included in Supplement Data File 4 but are not discussed here. For compilations of benthic foraminiferal Mg∕Ca ratios, see Cramer et al., 2011).

Sample preparation. Foraminifera need to be thoroughly cleaned prior to analysis. Clay and organic contaminants are removed using a short oxidative procedure, whereas removal of metal oxide contaminants requires a longer procedure including a reductive step (Boyle and Keigwin, 1985). These two procedures result in offsets in Mg∕Ca values, which must be corrected when making comparisons to other records (e.g., Barker et al., 2003; Yu et al., 2007). Cleaning efficacy is assessed using Al∕Ca, Mn∕Ca and Fe∕Ca ratios (e.g., Boyle, 1983; Barker et al., 2003). In some cases, a simple threshold value may be used to screen samples, e.g., Al∕Ca > 80 µmol mol⁻¹ (Mawbey and Lear, 2013). In many cases, the threshold depends on contaminant composition such that Mn∕Mg or Fe∕Mg ratios may be a more useful indicator (Barker et al., 2003). Cleaned foraminifera are commonly dissolved in acid and analyzed by ICPMS, taking into account the dependence of measured Mg∕Ca on analyte concentration (the matrix effect – Lear et al., 2002), although other techniques may also be employed (see below). Prior to crushing, several representative specimens should be selected for SEM analysis to record the extent of textural recrystallization on broken chamber walls. Sr∕Ca values routinely collected alongside Mg∕Ca data can provide one means of monitoring the impact of recrystallization on test geochemistry. Inorganic calcite tends to have lower Sr∕Ca and higher Mg∕Ca than foraminiferal calcite (Baker et al., 1982), which leads to inverse relationships in diagenetically altered downcore records.

Recommended steps. The following are recommended steps in converting planktic Mg / Ca ratios to temperatures:

• i.

  The possible presence of dissolution should be determined. In the recent geological past, a dissolution correction can be based on the bottom water saturation state at the sample site (Dekens et al., 2002), but uncertainties in ocean carbonate chemistry currently preclude this approach in deep-time studies. We recommend that planktic foraminiferal test weights are reported along with Mg∕Ca data to assess the potential impact of dissolution (e.g., Rosenthal and Lohmann, 2002). Alternatively, chemically resistant domains within individual tests can be selected for analysis (see above).

• ii.

  A salinity correction should be applied if there is independent evidence that the sample site experienced substantial deviations from normal salinity, such as the large changes in the hydrological cycle inferred for the PETM (Zachos et al., 2003). Normalizing culture data of three modern species (compiled in Hönisch et al., 2013; Allen et al., 2016) to the Mg∕Ca observed at a salinity of 35 psu for each species results in the following multi-species salinity sensitivity:

  $$\begin{array}{}\text{(2)}& \begin{array}{rl}\mathrm{Mg}/{\mathrm{Ca}}_{\mathrm{CORRECTED}}& =(\mathrm{1}-(\mathrm{salinity}-\mathrm{35})\times \mathrm{0.042}\pm \mathrm{0.008})\\ & \times \mathrm{Mg}/{\mathrm{Ca}}_{\mathrm{MEASURED}}\end{array}\end{array}$$

  (see the Supplement for further details). This sensitivity of 4.2±0.8 % per PSU is in good agreement with global sediment trap and plankton tow data for G. ruber (Gray et al., 2018).

• iii.

  A correction for past changes in the carbonate system should be applied. This is complicated, however, because pH and [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] are not only driven by long-term changes in the carbon cycle but also by factors such as the temperature effect on the dissociation constant of water (K_W). Therefore, whilst the results of Gray et al. (2018) indicate that the Mg∕Ca–temperature sensitivity in the modern ocean is 6 % per degree Celsius when the effects of pH and temperature are fully deconvolved, this can only be applied if the control of temperature on K_W (and therefore pH) is accounted for (ideally through δ¹¹B-derived pH reconstructions using the same material). For instance, if the pH reconstruction available for a given interval was determined at a different site, a temperature sensitivity of 6 % per degree Celsius should be applied only if the difference in temperature between sites can be estimated, so that the temperature-driven intersite pH gradient can be accounted for (pH differences between sites may also exist for other reasons). In practice, this requires that the equations for pH and temperature be solved iteratively, given that Mg∕Ca is sensitive to both factors. We recommend differing approaches for Mg∕Ca data treatment depending on whether a δ¹¹B pH record is available for the same site (see step v below).

  Where δ¹¹B-based pH reconstructions are available, Mg∕Ca ratios should be corrected for pH. However, this correction should be considered with caution until the controlling carbonate system parameter on foraminifera Mg∕Ca is identified, particularly given that pH and [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$] may be decoupled over geological time (Tyrrell and Zeebe, 2004). If no pH reconstruction is available for the site of interest, then the temporally and spatially closest data should be used and uncertainties in applying this considered. Based on a linear fit through culture data from three modern species (Evans et al., 2016a), the correction is as follows:

  $$\begin{array}{}\text{(3)}& \begin{array}{rl}\mathrm{Mg}/{\mathrm{Ca}}_{\mathrm{CORRECTED}}& =(\mathrm{1}-(\mathrm{8.05}-\mathrm{pH})\times \mathrm{0.70}\pm \mathrm{0.18})\\ & \times \mathrm{Mg}/{\mathrm{Ca}}_{\mathrm{MEASURED}}.\end{array}\end{array}$$

  Note that Mg∕Ca may relate nonlinearly to pH outside the range 7.7–8.4, and more complex relationships have been suggested (Russell et al., 2004; Evans et al., 2016a). This sensitivity of $-\mathrm{7.0}\pm \mathrm{1.8}$ % per 0.1 pH unit is in agreement with the $-\mathrm{8.3}\pm \mathrm{7.7}$ % derived from a global sediment-trap and plankton tow G. ruber dataset (Gray et al., 2018), although we recommend applying the culture-derived expression in deep time because it is calibrated over a much wider pH range.

• iv.

  A correction for Mg∕Ca_sw is usually applied to the pre-exponential component (B) of the relevant Mg∕Ca–temperature calibration of the form Mg∕Ca = Bexp^AT (Hasiuk and Lohmann, 2010; Evans and Müller, 2012):

  $$\begin{array}{}\text{(4)}& B_{\mathrm{CORRECTED}}=({\mathrm{Mg}/{\mathrm{Ca}}_{\mathrm{sw}}}^H/{\mathrm{5.2}}^H)\times B_{\mathrm{MODERN}},\end{array}$$

  where B_MODERN is the pre-exponential coefficient derived from modern calibrations, H is the nonlinearity of the relationship between shell and Mg∕Ca_sw, which may be estimated from culture experiments under variable seawater chemistry (Evans et al., 2016b; Delaney et al., 1985), Mg∕Ca_sw is that of the time interval of interest, and 5.2 is the modern seawater Mg∕Ca ratio in moles per mole (mol mol⁻¹). However, the observation that the slope of this relationship is sensitive to Mg∕Ca_sw in culture experiments means that equations describing the change in both constants (A and B) have been reported for modern taxa (Evans et al., 2016b), with the implication that the use of modern calibrations may underestimate relative temperature changes during the early Paleogene. For estimates of early Paleogene Mg∕Ca_sw we recommend the use of a relatively high-precision, million-year (average) resolution Mg∕Ca_sw reconstruction from the coupled analysis of Mg∕Ca and clumped isotopes in foraminifera (Evans et al., 2018a).

• v.

  Mg∕Ca is converted to temperature using an exponential calibration equation:

  $$\begin{array}{}\text{(5)}& T=\ln \left(\mathrm{Mg}/{\mathrm{Ca}}_{\mathrm{CORRECTED}}/B_{\mathrm{CORRECTED}}\right)/A,\end{array}$$

  where A and B_CORRECTED are derived from species or group-specific calibrations. The most appropriate calibration for extinct Eocene species should be chosen based on similarities to extant species in terms of shell chemistry or, for example, position within the water column and the presence or absence of symbionts. Best practice would be to report the sensitivity of a reconstruction to the choice of calibration equation. Importantly, the sensitivity factor A should be modified depending on how the carbonate system correction described above is performed. If pH is explicitly accounted for at the site of interest through δ¹¹B then the 6 % per degree Celsius sensitivity of Gray et al. (2018) should be applied. However, if the best available pH reconstruction is from a different site or time interval or derived from a model which represents the global mean (e.g., Tyrrell and Zeebe, 2004), then we recommend applying the apparent sensitivity derived from culture (Kısakürek et al., 2008; Evans et al., 2016a) as this indirectly accounts for the effect of temperature on K_W.

These recommendations have been applied to the Mg∕Ca analyses included in the DeepMIP database (Supplement Data File 4). We have compiled Mg∕Ca data for planktic foraminifera from five DSDP or ODP sites and six onshore sections. SST was derived from Mg∕Ca ratios as follows: 1000 random draws were performed of salinity (33–37 psu), seawater Mg∕Ca (within the 95 % CI given by Evans et al., 2018b), pH (±0.2 units) and the Mg∕Ca–pH sensitivity (anywhere between 0 %–8.8 % per 0.1 unit, i.e., anywhere between not sensitive at all to the upper confidence interval on the modern culture calibrations) for each data point. Calibration uncertainty was assessed by randomly choosing either the laboratory calibrations of Evans et al. (2016b), which define a Mg∕Ca_sw-dependent Mg∕Ca–T sensitivity, or the modern calibration with an “H-factor” applied to the pre-exponential constant (Evans and Müller, 2012), which maintains the modern Mg∕Ca–T sensitivity in deep time. The uncertainties on each data point are then the 97.5 and 2.5 percentiles of these 1000 sets of assumptions. The “best estimate” SSTs are the 50 percentiles of the subset of these 1000 draws that use the calibrations of Evans et al. (2016b), which is preferred because the available evidence suggests that the Mg∕Ca–T sensitivity varies as a function of Mg∕Ca_sw. Note that the data are subject to revision following replication of that study, and that using the 50th percentile of all estimates including both calibrations would result in overall cooler SST. Analytical uncertainty is not considered significant given the magnitude of these uncertainties.

Intra-test Mg / Ca analysis by LA-ICPMS and EPMA. In contrast to solution inductively-coupled-plasma mass spectrometry, which enables high throughput of pooled, dissolved foraminifera, highly spatially resolved techniques such as laser-ablation ICPMS (LA-ICPMS) facilitate targeted in situ analysis of carbonates (Eggins et al., 2003) and allow intra-specimen preservation to be assessed (e.g., Creech et al., 2010; Evans et al., 2015). Small samples such as foraminifera may be analyzed without embedding or sectioning, and high-sample throughput means that the technique is relatively inexpensive. Laser spot sizes are typically 20–80 µm in diameter, with 5 µm possible (Lazartigues et al., 2014), enabling repeat measurements of individual chambers or seasonality reconstruction in long-lived organisms with incremental growth layers (Bougeois et al., 2014; Evans et al., 2013). Because each laser pulse removes less than 100 nm of material on carbonates (Griffiths et al., 2013), depth-profiling through the sample has an effective resolution of <0.5 µm when using a fast wash-out ablation cell (Müller et al., 2009). Therefore, element profiles through foraminifera chamber walls not only facilitate the characterization of intra-specimen preservation, which can be assessed by the simultaneous collection of Mg∕Ca, Sr∕Ca, Mn∕Ca, Al∕Ca and Fe∕Ca ratios, among others, but also allow diagenetically affected areas, such as surface overgrowths, to be excluded from the measurement used to calculate the oceanographic variable of interest (e.g., Hollis et al., 2015; Hines et al., 2017). The disadvantage of LA-ICPMS is that overall data acquisition and reduction is typically more time-consuming compared to solution-based techniques and there is relatively large intra-specimen variability.

Electron probe microanalysis (EPMA) is a microanalytical technique based on the detection of element-characteristic X-rays produced by bombarding the sample with an accelerated and focused electron beam. For quantitative analyses, the intensity of the element-specific X-rays is compared against those of the same elements from chemically well-characterized standards. However, the sensitivity of EPMA is limited, and typically only Mg∕Ca and Sr∕Ca ratios in foraminiferal shells can be quantitatively measured (e.g., Brown and Elderfield, 1996; Anand and Elderfield, 2005). One advantage of this technique is the high spatial resolution; typical beam spot sizes for quantitative analyses in foraminiferal shells is ∼2 to 10 µm (e.g., Hathorne et al., 2003). In addition, semiquantitative elemental distribution maps can be acquired with submicron resolution for a larger suite of elements (e.g., Pena et al., 2008), allowing for the identification of diagenetic phases.

Sample preparation for EPMA is time-consuming, as tests need to be embedded in epoxy and polished until suitable cross sections for analysis are exposed. This method is unique in allowing SEM imaging of chamber walls prior to analysis, including an assessment of contaminant and diagenetic phases. The same epoxy mounts can also be used for in situ δ¹⁸O analysis by SIMS, so that chemical and isotopic information can be deduced from adjacent domains just microns apart (Kozdon et al., 2013; Wycech et al., 2018). Semiquantitative maps of intrashell Mg∕Ca variability with submicron resolution are possible with EPMA (e.g., Eggins et al., 2004), but for quantitative Mg∕Ca measurements, an electron beam with a diameter of several microns is required, which averages between the high- and low-Mg bands or domains observed in many planktic (e.g., Eggins et al., 2004; Spero et al., 2015) and benthic (e.g., de Nooijer et al., 2014) foraminifera. The average of multiple Mg∕Ca measurements by EPMA from an individual shell is typically comparable to solution-phase data that consume the whole shell (e.g., Hathorne et al., 2003). A major advantage of this method is that it is nondestructive.

4.3.1 Theoretical background of clumped isotopes

The carbonate clumped isotope thermometer is based on the temperature-dependent proportion of ¹³C−¹⁸O bonds in carbonate minerals (Ghosh et al., 2006; Eiler, 2007). The proxy has a direct basis in thermodynamics (Schauble et al., 2006; Hill et al., 2014) and has been applied to a wide range of marine and terrestrial samples, from foraminifera to paleosol carbonates (e.g., Tripati et al., 2010; Snell et al., 2013; Douglas et al., 2014). The zero-point energy of atomic bonds decreases with the mass of the atoms involved, favoring bonds between the rare, heavy isotopes. However, this effect decreases with increasing temperature, leading to the theoretical and observed decrease of “clumping” of heavy isotopes with increasing formation temperature of the mineral (Eiler, 2007). Excess abundance of ¹³C−¹⁸O bonds is abbreviated to Δ₄₇ and refers to the over-abundance of CO₂ with the composition ${}^{\mathrm{13}}\mathrm{C}{-}^{\mathrm{18}}\mathrm{O}{-}^{\mathrm{16}}\mathrm{O}$ relative to a stochastic distribution of all isotopes (Eiler and Schauble, 2004). Δ₄₇ is measured on an isotope ratio mass spectrometer after acidification of the carbonate, in a very similar way to classical δ¹⁸O measurements. The only difference is that the abundance of mass 47 is recorded in addition to the traditional masses 44–46; δ¹⁸O and δ¹³C are obtained as by-products of the measurements and are needed to calculate the expected stochastic abundance of ¹³C−¹⁸O bonds in the respective sample, which is then compared to the observed abundance to calculate Δ₄₇.

4.3.2 Strengths of clumped isotopes

There are three key strengths to this carbonate-based paleothermometer: (1) both theory and empirical studies demonstrate that the isotopic composition of water exerts no measurable control on the clumped isotope signature (Ghosh et al., 2006; Schauble et al., 2006); (2) the technique involves the simultaneous measurement of Δ₄₇ and δ¹⁸O, enabling the independent and direct calculation of ancient δ¹⁸O_sw; (3) many, although not all, biogenic carbonates and inorganic precipitates fall on the same calibration line (e.g., Ghosh et al., 2006; Zaarur et al., 2013; Tang et al., 2014; Tripati et al., 2015). Molluscs and brachiopods (Came et al., 2007; Eagle et al., 2013a; Henkes et al., 2013), foraminifera (Tripati et al., 2010; Grauel et al., 2013; Evans et al., 2018a; Peral et al., 2018; Piasecki et al., 2019), paleosols (Passey et al., 2010), land snails (Eagle et al., 2013b), and other forms of carbonate (Eiler, 2007; Kele et al., 2015) all appear to be reliable archives for the measurement of Δ₄₇. Clumped isotope analysis of speleothems and certain coral species (Ghosh et al., 2006; Tripati et al., 2010, 2015; Saenger et al., 2012; Affek and Zaarur, 2014; Loyd et al., 2016; Spooner et al., 2016) as well as other taxa (Davies and John, 2019) are more uncertain and require further study.

4.3.3 Weaknesses of clumped isotopes

Whilst clumped isotope thermometry has been successfully applied to a wide range of sample types, there are several challenges associated with paleoclimate reconstructions. Of these, the most fundamental is the low abundance of doubly substituted (“clumped”) carbonate, which typically makes up only ∼46 ppm of the total CO₂ produced from a sample. Precision is therefore ultimately limited by our ability to cleanly measure mass 47 CO₂. In practice, this implies relatively large sample masses, typically ∼10 mg of material (>500 planktic foraminifera), although recent advances in instrumentation have seen this reduced by a factor of 5–10 in some laboratories (Meckler et al., 2014; Müller et al., 2017). Analytical precision for these sample sizes limits the accuracy of the technique to ±2–3 ^∘C (1 sigma), which can be improved by performing a greater number of replicate measurements (e.g., Huntington et al., 2009; Thiagarajan et al., 2011; Tripati et al., 2014), with an obvious trade-off between sample size and precision. The presence of organic carbon, which can contribute to mass 47, may necessitate stringent sample cleaning procedures. In common with many proxies, the potential for seasonal growth of some archives must be considered. This is particularly the case for molluscs. Careful sample selection and geological context are of critical importance when designing studies and interpreting clumped isotope data.

As for δ¹⁸O and Mg∕Ca, the preservation of foraminifera and other carbonates is a key issue that must be addressed. While the impacts of dissolution are not known, recrystallization at different temperatures or the addition of secondary diagenetic calcite precipitated after deposition can bias clumped isotope measurements of planktic foraminifera (Shenton et al., 2015; Stolper et al., 2018). However, recent work has shown that glassy foraminifera and some nonglassy planktic and benthic foraminifera appear to yield reliable clumped isotope data for paleoceanographic reconstructions (Leutert et al., 2019). Solid-state reordering within the calcite mineral will alter the isotope ordering, although only in samples that have experienced burial temperatures well above 100–150 ^∘C for over ∼10 Ma (Passey and Henkes, 2012; Henkes et al., 2014; Shenton et al., 2015).

Several empirical calibrations of Δ₄₇ to temperature have been developed, with differences partly related to laboratory-specific analytical methodology; these include acid digestion temperature, instrumentation, ion corrections and standards. Some discrepancies observed in earlier studies may be due to the limited temperature range and replication when calibrating biogenic carbonates, whereas calibrations with larger temperature ranges (>90 ^∘C) across different types of carbonates agree well (Bonifacie et al., 2017; Fernandez et al., 2017; Kelson et al., 2017). Material-specific calibrations have also been suggested, for example, for marine molluscs (Eagle et al., 2013b; Henkes et al., 2013). Accurate absolute temperature reconstructions depend on empirical calibrations being developed or checked in each laboratory. Differences between laboratories have been attributed to a range of factors (Dennis et al., 2011; Wacker et al., 2014; Defliese et al., 2015; Daëron et al., 2016; Schauer et al., 2016) but are not yet fully resolved. With the use of rigorous standardization procedures to correct for instrument drift and more consistent methodologies overall, it is hoped that calibrations between instruments and labs will converge (e.g., Bernasconi et al., 2018).

4.3.4 Recommended methodologies for clumped isotopes

Here we outline our recommendations for generating clumped isotope data from fossil shells, and for converting these data into temperature. The current dataset for the targeted age range is limited to five early Eocene onshore sections (Supplement Data File 5) and paleosol carbonate data from four continental North American sites (Supplement Data File 7A). Using the methods outlined below, we have calculated SSTs and error values and compiled a summary of SST estimates for these EECO sites.

We recommend the use of a calibration that covers a wide temperature range and sufficiently replicated analyses of calibration samples from the carbonate type being analyzed, performed with the same analytical and data processing procedures as those employed for the samples (e.g., in the absolute reference frame as defined by Dennis et al., 2011, and done using either gas-standard- or carbonate-standard-based reference frames of Dennis et al., 2011, or the carbonate standard-based reference frame of Bernasconi et al., 2018). Where possible, we recommend that previously published calibrations be converted to these reference frames.

A recent meta-analysis of synthetic carbonate calibrations (Petersen et al., 2019) explores the causes of interlaboratory offsets and makes recommendations on best practice. Best practice for measuring, correcting and reporting Δ₄₇ includes (1) measurement of a large number of gas and/or carbonate standards of different compositions; (2) correction for instrumental nonlinearities (Huntington et al., 2009; Dennis et al., 2011), such as those which arise from secondary electrons (He et al., 2012; Bernasconi et al., 2013), or use of instruments with electron suppression; (3) reporting of data on an absolute reference frame (Dennis et al., 2011); and (4) reporting full methodology, including digestion apparatus, digestion temperature, gas cleaning procedure, mass spectrometer and corrections used, working gas composition, constant sets used for calculations, acid digestion fractionation factor, Δ₄₇ values and errors, temperature calibration, and estimated temperatures. Ideally, the provision of gas and/or carbonate standard results together with the sample data facilitates the broad use of data and future recalculations. We recommend the open archiving of raw intensity values and datasets, with their own digital object identifiers, to enable the reanalysis of data over the long term.

Given the required sample amounts and time-intensive measurements, clumped isotope thermometry is most powerful where other proxies are limited by unknown confounding effects. Rather than providing high-resolution reconstructions, clumped isotopes can be used to ground-truth and improve the accuracy of other proxies, including new constraints on seawater compositions (Evans et al., 2018a). Additionally, systems allowing the repeated measurement of small (∼100–500 µg) sample aliquots give the required replication rates for Δ₄₇, which are averaged across replicates or adjacent samples, and have the potential to provide higher-resolution records of standard δ¹⁸O and δ¹³C analyses. As for δ¹⁸O and Mg∕Ca, samples need to be carefully screened for diagenetic alteration.

4.4.1 Theoretical background of TEX₈₆

The tetraether index of tetraethers consisting of 86 carbon atoms (TEX₈₆) is an organic paleothermometer based on the distribution of isoprenoidal glycerol dialkyl glycerol tetraethers (isoGDGTs) in marine or lake sediments (Schouten et al., 2002)., Within marine environments, isoGDGTs are inferred to be mainly derived from marine Thaumarchaeota (Schouten et al., 2002; Wuchter et al., 2004). Marine Thaumarchaeota occupy much of the epipelagic and mesopelagic zone, but cell numbers are highest in the upper few hundreds of meters of the surface ocean (Church et al., 2010), with TEX₈₆ correlating most strongly with sea surface (SST) or shallow subsurface (sub-T, 50–200 m) temperatures (Tierney and Tingley, 2015). The underlying principal of TEX₈₆ is that the number of cyclopentane rings (moieties) in GDGTs increases with growth temperature in order to alter the fluidity and permeability of the cell membrane (Sinninghe Damsté et al., 2012). Laboratory culture and mesocosm experiments confirm this relationship and indicate that TEX₈₆ values continue to increase with temperature above 30 ^∘C (Wuchter et al., 2004; Schouten et al., 2007; Pitcher et al., 2009; Kim et al., 2010), the upper limit of the modern coretop dataset (Kim et al., 2010; Tierney and Tingley, 2015). Although these studies suggest that a linear relationship between TEX₈₆ and temperature persists at least to 40 ^∘C, the form of the relationship has been shown to vary significantly between different strains of Thaumarchaeota (Elling et al., 2015).

TEX₈₆ has been widely used to reconstruct SST during the Eocene (Brinkhuis et al., 2006; Pearson et al., 2007; Sangiorgi et al., 2008; Bijl et al., 2009, 2010, 2013; Sluijs et al., 2011; Hollis et al., 2009, 2012; Douglas et al., 2014; Frieling et al., 2014; Inglis et al., 2015; Cramwinckel et al., 2018) and across the PETM (Sluijs et al., 2006, 2007a, b, 2011, 2014; Zachos et al., 2006; Schoon et al., 2015; Frieling et al., 2017, 2018).

4.4.2 Strengths of TEX₈₆

TEX₈₆ has several key advantages over other paleothermometers. Isoprenoidal GDGTs are relatively resistant to thermal degradation and diagenesis (Schouten et al., 2004) and can be used to reconstruct SST in settings where well-preserved foraminiferal calcite is absent (e.g., Bijl et al., 2009; Hollis et al., 2012; Pross et al., 2012) and where alkenones (e.g., the $U_{\mathrm{37}}^{K^{\prime }}$ paleothermometer) are either absent or outside their calibration range (Brassell, 2014). Furthermore, TEX₈₆ is not affected by changes in pH and salinity (Wuchter et al., 2004; Elling et al., 2015) and does not require a correction to account for past changes in ocean seawater chemistry (see Sect. 5.1.3 and 5.1.4). TEX₈₆ values are also not substantially impacted by sediment extraction, and processing techniques and interlaboratory variation in TEX₈₆ values are comparable to those of other temperature proxies (Schouten et al., 2013a).

4.4.3 Weaknesses of TEX₈₆

The interpretation of TEX₈₆ values can be complicated by the input of GDGTs from archaea other than marine Thaumarchaeota. Several indices have been developed to screen for such secondary inputs. The branched-to-isoprenoidal tetraether (BIT) index (Hopmans et al., 2004) is used to assess the impact of terrestrial input upon TEX₈₆ values. TEX₈₆-derived SSTs may be warm-biased by 2 ^∘C or more when BIT indices exceed 0.4 (Weijers et al., 2006). However, the impact of terrigenous GDGTs is highly variable, being dependent on the nature and temperature of the source catchment (Bijl et al., 2013; Douglas et al., 2014; Inglis et al., 2015). Sedimentary GDGT production can also affect TEX₈₆ values. Methanotrophic Euryarchaeota can synthesize GDGTs with up to three cyclopentane moieties (Pancost et al., 2001) and may impact TEX₈₆ values in settings in which the anaerobic oxidation of methane (AOM) has occurred. The methane index (MI) can be used to assess AOM, whereby high MIs (>0.5) reflect high rates of AOM and low values (<0.3) suggest normal sedimentary conditions (Zhang et al., 2011). During the Eocene, MI values rarely exceed 0.5, suggesting that most settings are relatively unaffected by diffusive methane flux and associated AOM. Methanogenic archaea can synthesize GDGT-0 as well as smaller quantities of GDGT-1–3 and may bias TEX₈₆ values (Blaga et al., 2009). This can be assessed using the %GDGT-0 index (Inglis et al., 2015). Culturing experiments indicate that a methanogenic source of GDGT-0 is possible when the index is >67 %. Methanogenesis appears to have a minor influence on Eocene records as >90 % of analyses are below this threshold (Inglis et al., 2015).

One of the major uncertainties in TEX₈₆ paleothermometry is the impact of nontemperature factors upon thaumarchaeotal GDGT distributions, including growth phase (Elling et al., 2014; Hurley et al., 2016) and oxygen concentrations (Qin et al., 2015). Although these factors remain difficult to assess in the geological record, the ring index (Zhang et al., 2016) and fCren${}^{\prime }:$ Cren${}^{\prime }+$ Cren (O'Brien et al., 2017) can potentially help to reveal nontemperature controls or exclude the influence of such biases. There is also evidence that sedimentary GDGTs are not exported solely from surface waters but also shallow subsurface waters (Taylor et al., 2013) and that their export may be seasonally biased (Herfort et al., 2006; Bijl et al., 2010). These factors may help to explain some apparent SST anomalies (e.g., Taylor et al., 2018) and discrepancies with other proxies (Eberle et al., 2010; Hollis et al., 2012; Inglis et al., 2015). The suggestion that these discrepancies may be explained by export of GDGTs from the deep ocean (Ho and Laepple, 2016) was refuted by Tierney et al. (2017), who argue that the sedimentary TEX₈₆ signal is predominantly derived from the upper water column (50–300 m) and that temperatures through the upper water column are strongly correlated, therefore minimizing any calibration bias. Export from the upper water column rather than sea surface will impart some uncertainty because TEX₈₆ is calibrated to climatological ocean temperatures, such that the local dominant depth or season of production and export could introduce bias, particularly in settings with relatively deep, sub-thermocline GDGT production (Taylor et al., 2013) or in middle to high latitudes (Hollis et al., 2012). Future calibration studies should provide detailed information about the prevailing oceanographic conditions at each core-top location. In areas with a shallow and steep thermocline, TEX₈₆ may be recording thermocline rather than surface variations (e.g., Huguet et al., 2007). To address this, both SST and sub-T calibrations have been developed (e.g., Schouten et al., 2002; Kim et al., 2008; Tierney and Tingley, 2015) and should be applied with knowledge of the oceanographic characteristics of the site in question.

A final weakness is that many early Eocene sites (e.g., Frieling et al., 2017, 2018; Cramwinckel et al., 2018) are characterized by TEX₈₆ values that exceed the calibration range in modern oceans (0.3–0.8). Whilst some uncertainty may be alleviated by comparison with associated multi-proxy data (e.g., Frieling et al., 2017; Tierney et al., 2017), TEX₈₆ values > 0.8 require extrapolation of the relationship between TEX₈₆ and temperature beyond the modern calibration dataset or incorporation of limited data from laboratory studies (e.g., Pitcher et al., 2009; Elling et al., 2015). Such high values almost certainly indicate hotter temperatures, but the degree of estimated warming depends on assumptions about the mathematical nature of the temperature–TEX₈₆ relationship (see below and discussion in Eley et al., 2019).

4.4.4 Recommended methodologies for TEX₈₆

Here we outline our recommendations for generating TEX₈₆ data from marine sediments, and for converting these data into temperature. We have compiled available TEX₈₆ data from five ODP and IODP sites, five onshore drill holes and nine onshore sections (Supplement Data File 6). Using the methods outlined below, we have calculated SSTs and associated uncertainties, and compiled a summary of SST estimates for LP, PETM and EECO time slices.

Preparation of sediments for TEX₈₆ analysis is straightforward, although large sample volumes (100 g) may be needed where concentrations of organic matter are low. After extraction and purification, the isoprenoidal GDGTs are analyzed by HPLC/MS (see Hopmans et al., 2000). Previously, a Cyano column was used for HPLC separation (Hopmans et al., 2000); however, the use of BEH HILIC UHPLC columns has now been shown to yield better separation while not substantially affecting the TEX₈₆ index (Hopmans et al., 2016). The TEX₈₆ index is calculated as the ratio between isoprenoidal GDGTs as follows:

where GDGT-1, GDGT-2 and GDGT-3 are characterized by one, two and three cyclopentane moieties and cren^′ is the regioisomer of crenarchaeaol.

Several TEX₈₆ calibrations have been developed to estimate sea surface or shallow subsurface temperatures. The original TEX₈₆ core-top calibration was a linear relationship to SST (Schouten et al., 2002). However, the correlation between TEX₈₆ and SST is weak below 5 ^∘C, and samples from the Red Sea deviate from the wider calibration dataset (Kim et al., 2008). To resolve this, TEX₈₆ was recalibrated to exclude samples from the Red Sea and sites where SST <15 ^∘C (e.g., Kim et al., 2008; Naafs and Pancost, 2016; O'Brien et al., 2017).

Various nonlinear relationships have also been proposed (e.g., Liu et al., 2009; Kim et al., 2010). Of these nonlinear calibrations, the most commonly used is TEX${}_{\mathrm{86}}^H$, which assumes an exponential relationship between temperature and GDGT distributions. It uses a log transformation of TEX₈₆ and excludes core-top data from the Red Sea as well as from sites where SST < 5 ^∘C. SST is derived from TEX${}_{\mathrm{86}}^H$ as follows:

However, the TEX${}_{\mathrm{86}}^H$ calibration does not account for observational error and is affected by regression dilution because TEX₈₆ is treated as the independent variable and therefore assumed to be error-free. The regression dilution causes it to systematically underestimate warm SSTs in the modern ocean and, by inference, may also underestimate SSTs in ancient greenhouse climates (Tierney and Tingley, 2014; O'Brien et al., 2017).

Kim et al. (2010) also proposed another index, TEX${}_{\mathrm{86}}^L$, which is also an exponential calibration but uses the entire core-top dataset (with the exception of the Red Sea) and was developed to reconstruct SST across all temperature ranges. However, TEX${}_{\mathrm{86}}^L$ does not strictly reflect the degree of cyclization in GDGTs and consequently lacks obvious biological rationale. It is also particularly sensitive to contributions from other archaea, especially those living in the subsurface (Kim et al., 2015; Taylor et al., 2013), and its application within the geological record has been questioned (Inglis et al., 2015; Taylor et al., 2018).

Available experimental evidence (Pitcher et al., 2009; Schouten et al., 2013b; Elling et al., 2015) suggests a linear relationship between TEX₈₆ and SST, which is consistent with the modern core-top dataset. However, it remains uncertain if TEX₈₆ continues to have a linear relationship with temperature beyond the range of the modern calibration (Cramwinckel et al., 2018). The functional form of the calibration has significant consequences for temperature estimation, especially for values > 0.8 (Fig. 4a), which is particularly relevant for the warm Paleogene, where such values are common.

Figure 4Relationships between TEX₈₆ calibrations and sea surface temperature for low-latitude Paleocene–Eocene sites (data from Cramwinckel et al., 2018). (a) Comparison of the exponential TEX${}_{\mathrm{86}}^H$ calibration (Kim et al., 2010), the linear TEX₈₆ (O'Brien et al., 2017) and BAYSPAR calibrations (Tierney and Tingley, 2014, 2015). Because of the spatially varying regression used for BAYSPAR, only one site (ODP Site 959) is plotted for this calibration. (b) Comparison of SST reconstructions from late Paleocene to early middle Eocene for the three calibrations. The timescale of Gradstein et al. (2012) is used in this figure. Approximate locations of the three targeted time intervals are also shown: latest Paleocene (LP), Paleocene–Eocene thermal maximum (PETM) and early Eocene climatic optimum (EECO).

Download

The Bayesian regression model (BAYSPAR) of Tierney and Tingley (2014, 2015) applies a linear calibration, but the regression terms can vary spatially to accommodate the modern regional variation in TEX₈₆ sensitivity, as in the Mediterranean Sea (Kim et al., 2015). In deep-time settings, BAYSPAR searches the modern core-top dataset for TEX₈₆ values that are similar to the measured TEX₈₆ value within a user-specified tolerance and draws regression parameters from these modern locations. This assumes that modern TEX₈₆ values yield information on the environmental factors that determine the TEX₈₆–SST relationship in the ancient oceans, including factors other than mean annual sea surface temperature that might regionally alter the sensitivity of the proxy. This assumption is difficult to test and has the potential to introduce additional uncertainties or bias.

BAYSPAR can be determined without excluding core-tops from the Red Sea or where SSTs < 15 ^∘C. This is likely important as Red Sea-like conditions cannot be precluded for many Eocene sites (Inglis et al., 2015). The Bayesian model can also be used to estimate uncertainty beyond the modern calibration range. BAYSPAR can be used by downloading the publicly available code from MATLAB (https://github.com/jesstierney/BAYSPAR, last access: 30 June 2019) or Python (https://github.com/brews/baysparpy, last access: 30 June 2019). Both surface (BAYSPAR_SST) and subsurface (BAYSPAR_SubT) calibrations are available, the latter being the gamma-weighted average of the temperature range over 0–200 m water depth, with a maximum probability at ∼50 m (Tierney and Tingley, 2015).

When applied to TEX₈₆ values that exceed the modern calibration range, BAYSPAR typically yields results that lie between the linear and exponential calibrations because it incorporates both “normal tropical marine” and “Red Sea” sensitivities, the latter of which involves a steeper TEX₈₆–SST slope that is closer to the exponential assumption (Fig. 4b). It is thus a good “middle ground” choice for the Paleogene data considered here, and so we have used BAYSPAR_SST in our compilations (Fig. 8, Supplement Data File 6). While an exponential model might be a viable alternative, we do not recommend application of the TEX${}_{\mathrm{86}}^H$ calibration because of the known bias from regression dilution.

Given the ongoing evolution of calibrations, we recommend that researchers include the fractional abundances of all isoprenoidal GDGTs as well as TEX₈₆ values in all future publications as these will facilitate easy recalculation if and when the functional form of the calibration becomes better constrained. If BAYSPAR is used, we recommend that researchers report the parameters used (the prior mean and standard deviation, the search tolerance, and the version of the calibration dataset) so that the results are reproducible. The issue of extrapolation beyond the temperature range of the modern GDGT calibration dataset is highlighted in Eley et al. (2019), which also proposes an alternative calibration methodology.

5.1.1 Theoretical background of branched tetraether (brGDGT) paleothermometry and early work

Branched glycerol dialkyl glycerol tetraethers (brGDGTs) are membrane-spanning lipids produced by bacteria, probably acidobacteria (Weijers et al., 2009; Sinninghe Damsté et al., 2011, 2018). First discovered in a Dutch peat deposit (Sinninghe Damsté et al., 2000), brGDGTs are ubiquitous in terrestrial settings (Schouten et al., 2013b). These compounds have been used increasingly as paleothermometers since a study of a global suite of mineral soil samples (Weijers et al., 2007) found that the degree of methylation in brGDGTs (MBT index) is correlated to LAT and pH, and that the number of cyclopentane rings in the same brGDGTs (CBT index) is related to soil pH, resulting in the MBT–CBT paleotemperature proxy. The proxy has been subsequently revised and recalibrated as MBT'–CBT (Peterse et al., 2012). The following three equations are used to derive LAT from the MBT'–CBT proxy (Eqs. 9–11).

The MBT'–CBT proxy has been used to reconstruct early Paleogene terrestrial temperatures from continental margin marine records from the Southern Hemisphere (Bijl et al., 2013; Pancost et al., 2013; Pross et al., 2012) and the PETM in the Northern Hemisphere (Schoon et al., 2015; Weijers et al., 2007). These studies explicitly assumed that the majority of brGDGTs in marine sediment cores are produced within mineral soils at around sea level close to the core site and transported without alteration into marine sediments of similar age. The resulting LATs are significantly higher than those of the present day. However, despite being in good agreement with other terrestrial proxies for LAT as well as with modeled LAT (Huber and Caballero, 2011), they are ∼5–10 ^∘C lower than SST estimates derived from the same samples. This anomaly is discussed further below.

5.1.2 Strengths of brGDGT paleothermometry

Terrestrial temperatures are a crucial climatic variable, but most of our temperature records during the geological past are from the marine realm. The main strength of the MBT'–CBT terrestrial temperature proxy is that it is one of the few geochemical proxies currently available that can provide high-resolution quantitative estimates of terrestrial temperatures during the Cenozoic. In addition, brGDGTs are ubiquitous in ancient soils (e.g., Chinese loess), peats and lignites from across the Cenozoic (Peterse et al., 2012; Naafs et al., 2018a, b), making it possible to apply this proxy across a range of archives and depositional environments.

5.1.3 New developments in brGDGT paleothermometry

Recent studies suggest that caution should be applied to the interpretation of existing MBT'–CBT data. De Jonge et al. (2013, 2014) used refined analytical techniques to show that the 5-methyl penta- and hexamethylated brGDGTs, which were used to calculate the original CBT and MBT' indices, coelute with newly identified 6-methyl brGDGTs that have a strong pH dependence. Reevaluation of the global mineral soil dataset has led to new indices and calibrations based on 5-methyl brGDGTs only, which are dependent on temperature alone (De Jonge et al., 2014),

as well as a multiple linear regression (MAT_mr) calibration,

Naafs et al. (2017a) extended the calibration dataset and developed a new calibration for LAT as well as one to reconstruct growing degree days above freezing (GDD₀), a measure of cumulative temperature over the growing season. GDD₀ is a guide to warm-season temperatures and helps to constrain seasonal temperature variations, which is particularly important in high latitudes with high seasonality.

5.1.4 Weaknesses of brGDGT paleothermometry

Recent studies have identified challenges in applying brGDGT paleothermometry in marine settings. The assumption that brGDGTs in marine sediments are mainly derived from mineral soils formed at sea level may not be valid in some cases. Recent work has shown that brGDGTs are also produced in situ in the marine water column and in the hosting marine sediments (Sinninghe Damsté, 2016; Weijers et al., 2014), in rivers (Zell et al., 2014) and in lakes (Naeher et al., 2014; Weber et al., 2015). To illustrate this, we plot the relative abundance of the tetra-, penta- and hexamethylated brGDGTs of the published early Paleogene data from five sites together with the modern mineral soil and peat data (Naafs et al., 2017a, b) and early Paleogene data from lignites (Naafs et al., 2018a, b), following the approach of Sinninghe Damsté (2016) (Fig. 5a). The result demonstrates that most of the marine samples analyzed in this dataset have brGDGT distributions that are different from those found in any modern mineral soil or peat with a higher relative contribution of hexamethylated brGDGTs. This difference in brGDGT distribution suggests that these brGDGTs may not originate exclusively from mineral soils or peat. However, it is worth noting that the LATs derived from MBT'–CBT from some of these sites are in close agreement with LATs derived from pollen assemblages using the nearest living relative methods (Pross et al., 2012; Pancost et al., 2013).

Figure 5Distribution of brGDGTs in present-day mineral soils (Naafs et al., 2017a) and peats (Naafs et al., 2017b), compared to that in early Paleogene marine sediments (Weijers et al., 2007; Bijl et al., 2013; Pancost et al., 2013) and terrestrial lignites (Naafs et al., 2018): (a) a ternary plot showing relative proportions of tetra-, penta- and hexamethylated brGRGDTs; (b) a cross-plot showing the average number of cyclopentane moieties (rings) for tetra- and pentamethylated brGDGTs.

Download

Even when derived from mineral soils, the brGDGT assemblage in marine sediment will reflect varying combinations of settings across an entire terrestrial catchment, including different elevations and microclimates (Bendle et al., 2010). Moreover, brGDGT distributions have also been found to be influenced by soil moisture content in some settings, which needs to be considered in paleoclimate reconstruction, especially in arid regions (Dang et al., 2016; Menges et al., 2014).

Although brGDGTs appear resistant to initial degradation in the water column, thermal maturation of brGDGTs on geological timescales does influence the proxy. Schouten et al. (2013b) used artificial maturation experiments to demonstrate that the degree of cyclization of brGDGTs (e.g., CBT-indices) changes at temperatures above 240 ^∘C, and brGDGTs are completed degraded at 300 ^∘C. However, the influence of thermal maturity on the degree of methylation (MBT indices) appears to be limited as long as brGDGTs are not fully degraded.

A further weakness of this proxy is that it fails to capture the super-warm temperatures of the early Paleogene because the calibrations saturate at ∼25 ^∘C (e.g., Peterse et al., 2012; Naafs et al., 2017a). Many of the proxies for terrestrial temperature discussed here suffer from this problem because the modern environments that form the analogues for the calibration do not extend into tropical settings >25 ^∘C. For this reason, datasets should be scrutinized for evidence of LAT plateaus, which suggest the calibration has saturated.

Because the co-elution of newly identified 6-methyl brGDGTs may compromise MBT'–CBT estimates for LAT, earlier datasets should be reanalyzed using the latest analytical methods and calibrations. This problem may not be widespread because the brGDGT distributions reported in several Paleogene studies (e.g., Bijl et al., 2013) are dominated by tetramethylated brGDGTs, which have no 6-methyl homologues.

5.1.5 Recommended methodology for brGDGTs

If the brGDGT-proxy is to be applied to reconstruct terrestrial temperatures during the Paleogene the latest set of calibrations (De Jonge et al., 2014; Naafs et al., 2017a) and analytical techniques that separate the 5- and 6-methyl brGDGTs should be used (e.g., Hopmans et al., 2016). In addition, if applied to marine sediments to reconstruct terrestrial temperatures, supporting evidence should be provided (Fig. 5) to demonstrate a dominant input of brGDGTs from mineral soils or peat.

Another option is to avoid the source uncertainty associated with analyses from marine sediments by directly analyzing terrestrial sediments. For example, Inglis et al. (2017) used the latest analytical techniques and mineral soil calibrations to reconstruct early Eocene terrestrial temperatures from lignites in northern Germany. Although determining accurate ages for lignites can be challenging, these systems are dominated by in situ production of brGDGTs, resulting in reconstructed temperatures with a higher fidelity to local LAT.

Some studies have questioned whether soil-based calibrations should be applied to peat and lignite (Huguet et al., 2013, 2014; Weijers et al., 2011). Consequently, Naafs et al. (2017b) recently developed new peat-specific brGDGT-based temperature and pH calibrations, following Eqs. (16)–(19).

In addition, the relative abundance of H-brGDGTs (brGDGTs with an extra covalent bond between the two alkyl chains) and the degree of methylation of H-brGDGTs is also correlated with LAT (Naafs et al., 2018b). Although the correlation is not strong enough to provide quantitative temperature estimates, a high abundance of H-GDGTs in ancient lignites could provide additional evidence for elevated terrestrial temperatures. Applying the peat-specific brGDGT temperature calibration to the lignite record from northern Germany confirms that LAT was significantly higher than that of the present day, and 2–3 ^∘C higher than those obtained with the mineral soil calibrations (Naafs et al., 2018a). These peat-specific calibrations were recently applied to a range of early Paleogene lignites from across the globe and indicate that LATs were significantly higher than previously thought and are more consistent with the marine-based temperature estimates (Naafs et al., 2018a).

In summary, the application of peat-specific brGDGT-based temperature calibrations to early Paleogene lignites is in the early stages of development but holds great potential to provide new insights into terrestrial temperatures during this key period of the Cenozoic (Naafs et al., 2018a). MBT-based approaches are still judged useful to provide temperature constraints, especially for multi-proxy comparisons in well-dated marine sediments, but must be carefully assessed for calibration saturation at high temperatures, sampling of wide catchment areas, and inputs of brGDGTs that are not sourced from nonmineral soils or peat.

5.2.1 Theoretical background of leaf-morphology-based approaches

Methods of estimating paleotemperatures and other climate parameters from the physiognomy of fossil leaves have been explored since the early 1900s (Bailey and Sinnott, 1915, 1916). They are based on the general premise that plants adapt their physical features to tolerate the climate in which they germinate due to their immobility. Early observations of Southeast Asian floras identified a positive correlation between temperature and the percentage of woody dicotyledonous (dicot) species that had smooth-margined leaves (Bailey and Sinnott, 1915). This correlation was later developed into the widely applied method, leaf margin analysis (LMA; Wolfe 1979; Wing and Greenwood, 1993). Numerous modifications, different calibration datasets and measures of uncertainty related to LMA have since been explored (e.g., Wilf, 1997; Greenwood et al., 2003, 2004; Kowalski and Dilcher, 2003; Hinojosa and Villagran, 2005; Miller et al., 2006; Spicer et al., 2011). Essentially, LMA involves assessment of the margin character (whether smooth – entire – or toothed – nonentire) of all dicot leaf morphotypes in a fossil assemblage. The percentage of entire margins is calculated, and temperature determined using an LMA calibration equation.

Several shortcomings were recognized quite early in the development of the LMA method. For example, the relationship between the percentage of entire margins and temperature was found to vary between geographic regions with no robust means of ascertaining which modern calibration should be used for a fossil leaf flora. In response to these shortcomings and in recognition that other leaf characters may be correlated with climate variables, multivariate methods were explored. CLAMP (Wolfe, 1993; Spicer et al., 2009; Yang et al., 2011, 2015) is the most widely used of these methods. In this method, 31 leaf character variables are measured for each fossil leaf morphotype and estimates of several climate parameters are produced, including temperature, precipitation and humidity. CLAMP relies on comparison of fossil leaf data with large calibration datasets of modern dicot leaf morphology data collected using rigid protocols. The initial calibration dataset (Wolfe, 1993) was dominated by Northern Hemisphere assemblages. Subsequently, calibrations from other geographic areas have been developed (Steart et al., 2010; Jacques et al., 2011; Kennedy et al., 2014; Khan et al., 2014; Yang et al., 2015). Corresponding climate datasets are either based on the closest available long-term climate station measured data, or model-generated gridded datasets (Spicer et al., 2009). Other multivariate methods, such as the nearest neighbor approach (Stranks and England, 1997) and canonical correspondence analysis (Hinojosa and Villagran, 2005; Hinojosa et al., 2006) also use the CLAMP datasets. More information on calibration updates on datasets and methodology, and an online analysis, are available on the CLAMP website (http://clamp.ibcas.ac.cn/, last access: 30 November 2018).

Another multivariate method, digital leaf physiognomy (DiLP), applies multivariate linear regression models based on 92 calibration sites to fossil leaf morphology data captured digitally (Huff et al., 2003; Royer et al., 2005; Peppe et al., 2011). In addition to these leaf-morphology-based methods, leaf size has also been used as a univariate climate proxy, particularly for precipitation (Dolph and Dilcher, 1980; Wilf et al., 1998; Wiemann et al., 1998; Carpenter et al., 2012).

5.2.2 Strengths of leaf-morphology-based approaches

The strengths and weaknesses of leaf-morphology-based approaches to temperature reconstruction have been discussed extensively (e.g., Wolfe, 1993; Wilf, 1997; Uhl et al., 2003; Spicer et al., 2004; Traiser et al., 2005; Spicer et al., 2005, 2011; Green, 2006; Greenwood, 2007; Yang et al., 2007, 2015; Spicer and Yang, 2010; Jacques et al., 2011; Peppe et al., 2011; Kennedy et al., 2014; Khan et al., 2014; Greenwood et al., 2017; Lowe et al., 2018). A key strength is that they provide quantitative estimates that do not depend on correct assignment of fossil leaf forms to modern taxa and their current climate affiliations (e.g., Uhl et al., 2003). Each leaf-morphology-based method also has its own particular strengths. Univariate methods such as LMA are relatively easy to apply (e.g., Yang et al., 2015) whereas a strength of a multivariate method like CLAMP is that the inclusion of several leaf morphology characters alleviates some of the problems associated with taphonomic loss of character information in fossils. An advantage of DiLP is the minimization of any subjectivity or ambiguity in the leaf character scoring process because computer algorithms rather than individuals process the measurements (Peppe et al., 2011).

5.2.3 Weaknesses of leaf-morphology-based approaches

A number of assumptions are made when applying leaf-morphology-based paleoclimate methods, including that calibrations between leaf morphology and climate parameters are stable back through time and that fossil leaf collections are sufficiently representative of the source vegetation. Each method has its own set of weaknesses, many of which are discussed in the references given in Sect. 5.2.2. Univariate methods such as LMA are potentially strongly affected by taphonomic biases because they are reliant on the representative preservation of a single leaf character. LMA is also strongly influenced by geographic variability, with vegetation from different parts of the world yielding significantly different LMA calibrations. CLAMP is less influenced by these weaknesses but is more complicated and time-consuming to apply. As with the MBT'–CBT proxy, the calibration for leaf morphology proxies saturates at ∼25 ^∘C for most calibrations, perhaps at ∼28 ^∘C for CLAMP (Utescher et al., 2014). An exception is LMA using the Kowalski and Dilcher (2003) calibration, where LATs approach 30 ^∘C. This calibration endeavors to account for a cool bias in many other approaches, which is inferred to be due to the over-representation of toothed leaves in the wet environments that typify most settings where fossil leaves are preserved. However, this is still considered a provisional calibration dataset because it includes only 10 sites from a restricted geographic area (Spicer et al., 2011).

5.2.4 Recommended methodology for leaf-morphology-based approaches to estimating paleotemperature

A multi-proxy approach to reconstruction of land-based climate parameters is advocated where possible, including alternative leaf-morphology-based methods and taxonomic (nearest living relative) and geochemical approaches. A statistically assessed consensus approach (Ballantyne et al., 2010; Greenwood et al., 2017; Lowe et al., 2018) is recommended where feasible. For LMA, careful consideration should be given to the choice of calibration dataset, which have significant geographic variation. For CLAMP, the collecting and scoring protocols outlined on the CLAMP website should be followed to ensure the most robust estimates are produced. Updates to the calibration datasets and CLAMP methodology are also available through the CLAMP website (http://clamp.ibcas.ac.cn/, last access: 30 November 2018). It should be noted that both LMA and CLAMP are only applicable to leaf assemblages of woody dicot species and should be applied to fossil assemblages with at least 20 morphotypes (e.g., Wolfe, 1979; Povey et al., 1994; Spicer et al., 2004, 2005).

5.3.1 Theoretical background of nearest living relative (NLR) analysis

Nearest living relative (NLR) approaches are widely used techniques to derive quantitative climate parameters from pre-Quaternary fossil pollen or floral assemblages. In contrast to leaf-morphology-based methods that relate physical attributes of plants to climate, NLR uses a taxonomic approach to derive climate parameters from fossil pollen assemblages based on the climatic ranges of their present-day nearest living relatives. The NLR approach includes several different methods such as the mutual climate range technique (e.g., Pross et al., 2000; Thompson et al., 2012), bioclimatic analysis (e.g., Kershaw and Nix, 1988; Greenwood at al., 2005; Eldrett et al., 2009; Pross et al., 2012; Pancost et al., 2013), climatic amplitude method (e.g., Fauquette et al., 1998; Suan et al., 2017), the coexistence approach (Mosbrugger and Utescher, 1997; Carpenter et al., 2012; Pross et al., 2012; Pancost et al., 2013; Utescher et al., 2014) and variations of the probability density function method (Kühl et al., 2002; Greenwood et al. 2017; Harbert and Nixon, 2015; Hyland et al., 2018).

All these NLR methods apply similar principles. Each method requires a modern vegetation database that can be linked to modern climatology to calculate climate requirements of fossil taxa derived from the modern distribution of their NLR. The climate parameters usually include mean annual air temperature (MAAT or LAT as used here), coldest and warmest month mean temperature, and the mean annual precipitation range. To generate a paleoclimate estimate, the climate ranges in which the NLRs of a given fossil assemblage can coexist must be usually defined for each parameter. The climate estimates are usually shown as climate ranges of equally possible values within which a fossil plant community could have co-existed or as point values with confidence intervals. The latter should be the preferred format for further comparisons with other proxies and climate model outputs.

Concerns about the analytical approach of the coexistence approach (CA), one of the most widely used methods to derive quantitative climate estimates from pre-Quaternary palynological data, have recently been raised in regard to the treatment of outliers, and reliance on single taxa to define upper and lower limits for climate parameters (Grimm and Denk, 2012; Grimm and Potts, 2016). The latter significantly increases the uncertainty in CA-based climate reconstructions of “deep-time” environments that yield many extinct and relict taxa with NLRs whose climate envelope might have changed over time (e.g., Reichgelt et al., 2018). Probability density function (PDF) methods (Chevalier et al., 2014; Kühl et al., 2002; Greenwood et al., 2017) or the coexistence likelihood estimation method (Harbert and Nixon, 2015), which statistically constrain the most likely climatic co-occurrence envelope, offer a statistical solution to reduce the potential impact of wrongly defined climate tolerances.

5.3.2 Strengths of NRL analysis

NLR analysis uses the presence or absence of individual taxa in fossil assemblage rather than relative abundance, which reduces the likelihood of taphonomic biases. This facilitates, to some extent, the reconstruction of past nonmodern analogue climates and environments (Pross et al., 2000; Utescher et al., 2014). Due to the use of presence–absence data, the method also does not require a calibration dataset using modern plant communities. Other methods that assess the direct relationship between pollen abundance and climate (e.g., transfer function techniques) require large calibration datasets from modern environments and their use for deep-time climate reconstruction is therefore very limited (e.g., Bartlein et al., 2011; Garreta et al., 2012). NLR can be applied to a wide range of plant fossils, including seeds, woods, pollen or spores, and leaves. Despite the limitations outlined below, NLR-derived temperatures are generally in good agreement with estimates derived from other paleobotanical and geochemical approaches, including CLAMP and LMA (e.g., Uhl et al., 2003; Roth-Nebelsick et al., 2004; Pound and Salzmann, 2017; Pross et al., 2012).

5.3.3 Weaknesses of the NRL analysis

Quantitative climate estimates from the pre-Quaternary fossil plant record are always accompanied by large uncertainties. The greatest weakness that affects all NLR approaches is the assumption of uniformitarianism, namely that the climate tolerances of modern species can be extended into the past. This assumption inevitably introduces uncertainty that increases with the age of the geological record. Also, the modern distribution of a plant species is not only controlled by its climate envelope, but by a range of other factors such as biogeographic history and competition. By incorporating as many taxa as possible into the NLR analyses, uncertainties relating to such factors can be minimized to some extent (Harbert and Nixon, 2015; Willard et al., 2019).

The accuracy of estimates derived from NLR depends on the quality of the modern climate and plant distribution database, uncertainties in the identification of a fossil taxa and their nearest living relatives, and uncertainties on the total number of taxa within a fossil assemblage (Utescher et al., 2014). Accuracy is greatly reduced in sites with low taxa diversity (such as semideserts), large catchments (such as marine records) or nonmodern analogue vegetation (such as the Arcto-Tertiary relict flora). Also, NLR methods do not account for changes in the relative abundance of fossil taxa or long-distance transport of palynomorphs resulting in “mixed” assemblages from different habitats or vegetation zones. This latter factor places an important caveat on applying NLR to terrestrial assemblages in marine sediments. By its nature, NLR cannot accurately predict environmental conditions beyond the present-day range. For instance, maximum LAT estimates plateau at ∼25 for CA (Utescher et al., 2014).

5.3.4 Recommended methodology for NLR analysis

When deriving climate estimates from NRL analyses, we recommend, whenever possible, a statistically assessed ensemble approach, where multiple terrestrial climate proxies, including geochemical and paleobotanical methods, are combined and assessed using probability density analysis (e.g., Ballantyne et al., 2010; Pross et al., 2012; Greenwood et al., 2017; Lowe et al., 2018). Probability density functions should also be used to overcome the problematic use of absolute limits to describe climate response ranges (Kühl et al., 2002). NLR-derived climate estimates should be presented with confidence ranges including the uncertainty associated with the estimate of climate tolerances (Grimm et al., 2016). For comparisons using larger datasets, each climate estimate should be rated using a qualitative ranking system (i.e., confidence assessment) that addresses major influencing factors, such as taphonomy, diversity and age control (Salzmann et al., 2013). The confidence level should be used in addition to the indicated temperature range as a guide to assess the robustness of a temperature estimate.

The quality of the databases used to derive modern distribution and climate from the NLR is of vital importance for the accuracy of the paleoclimate estimates. Preference should be given to databases with clear policies on data quality and ownership, global coverage and long-term secured funding. We therefore recommend that NLR methods use the GBIF (Global Biodiversity Information Facility) database for sourcing distributions of extant taxa (http://gbif.org, last access: 30 June 2019), and that climate range data be derived from a global climate surface model such as WorldClim2 (Fick and Hijmans, 2017). Most NRL methods require at different stages of their analyses and database queries a pre-selection of modern and fossil data, which is necessary, for example, to exclude introduced garden plants from the modern dataset, redundant occurrences or taxa with doubtful assignments (e.g., Willard et al., 2019). Fossil datasets are often reduced to focus the analyses on a selected set of reliable climate indicators (e.g., Chevalier, 2019). To warrant transparency of the NLR process and reproducibility of the final estimates, it is essential that all steps of the data selection process are accurately described and published.

6.1.1 Theoretical background of boron-based CO₂ reconstructions

Boron isotopes within foraminiferal calcite (δ¹¹B_c) have been used extensively to reconstruct atmospheric CO₂ concentrations on a variety of different timescales (Foster, 2008; Hönisch et al., 2009, 2012; Pearson et al., 2009; Seki et al., 2010; Anagnostou et al., 2016; Gutjahr et al., 2017). The δ¹¹B–pH proxy has a strong grounding in theory (Hönisch et al., 2009; Foster and Rae, 2016) and is based on several foraminiferal calibration studies (e.g., Sanyal et al., 2001; Henehan et al., 2013; Foster and Rae, 2016).

6.1.2 Strengths of boron-based CO₂ reconstructions

Boron-based CO₂ reconstructions spanning the last ∼800 000 years show excellent agreement with ice core CO₂ estimates (Hönisch et al., 2009). The results of boron isotope studies have also been found to be highly reproducible between different laboratories (Foster et al., 2017; Gutjahr et al., 2017), allowing the integration of multiple datasets and providing high confidence in the fundamental measurement of boron isotopes in biogenic carbonates. In common with other marine CO₂ proxies (e.g., alkenones below), age constraints are also typically strong for δ¹¹B_c and it is easy to correlate the resultant CO₂ estimates with marine-based climate proxies, some of which can be measured on the same material (e.g., δ¹⁸O, Δ₄₇ and Mg∕Ca in planktic foraminifera). Finally, detailed studies comparing δ¹¹B_c from contemporaneous diagenetically recrystallized foraminifera and “pristine” glassy foraminifera show no influence of diagenesis on the recorded values (Edgar et al., 2015), indicating that δ¹¹B_c is more diagenetically robust than δ¹⁸O_c; the reasons for this are not yet fully understood (Edgar et al., 2015).

6.1.3 Weaknesses of boron-based CO₂ reconstructions

The complications of the boron isotope proxy can be separated in two categories, based on the two main steps required to estimate CO₂ from δ¹¹B_c values: (i) the calculation of pH from measured δ¹¹B_C values that requires additional constraints on seawater δ¹¹B (δ¹¹B_sw), the vital effects in extinct species, the calcification temperature, the relative depth habitat of foraminifera and the extent of any diagenetic overprint on primary δ¹¹B_c values, and (ii) the translation of the pH estimate into one of atmospheric CO₂, which first requires an estimate of a second carbonate system parameter to infer aqueous CO₂ concentrations at specific seawater concentrations for Mg and Ca, and determination of CO₂ disequilibrium between surface waters and the atmosphere at the study site. The influence of the latter can be reduced by careful site selection targeting oligotrophic settings.

Offsets between the boron isotopic composition of ambient seawater borate ion (δ¹¹B_borate) and that of foraminiferal calcite are generally known as “vital effects” following the usage of Urey et al. (1951) for oxygen isotopes. Many different processes could, in principle, be responsible for such offsets which could also differ between species. Of particular concern is the modification of pH in the diffusive boundary layer around foraminifera by the photosynthesis of algal symbionts (Hönisch et al., 2003). Other possibilities include physiological controls of the local pH during calcification (Rollion-Bard and Erez, 2010) and the incorporation of boron from aqueous boric acid (Vengosh et al., 1991).

6.1.4 Recommended methodology for boron-based CO₂ reconstructions

The recommended procedure for boron isotope reconstructions of atmospheric CO₂ is designed to constrain unknowns and limit uncertainties associated with the proxy method. First, δ¹¹B_sw can be constrained in a number of ways: (i) by modeling the pH profile of the water column using multiple species (Pearson and Palmer, 1999, 2000; Anagnostou et al., 2016), (ii) using pairs of planktic or benthic foraminiferal δ¹¹B_c (Greenop et al., 2017), (iii) utilizing existing records of benthic foraminiferal δ¹¹B_c coupled with assumptions regarding deep water pH evolution (Raitzsch and Hönisch, 2013), (iv) using geochemical box models of the boron cycle (Lemarchand et al., 2000). All of these approaches indicate that the maximum rate of δ¹¹B_sw change is most likely ∼0.1 ‰ Myr⁻¹, which is consistent with the long residence time of boron in the ocean (10–20 Myr), its relatively large oceanic budget, and the modern fluxes of boron in and out of the ocean (Lemarchand et al., 2002). Future studies should reconstruct surface ocean pH depth profiles using foraminiferal δ¹¹B_c to constrain changes in δ¹¹B_sw on shorter timescales, especially during periods when abrupt changes in weathering inputs of boron are possible. It is unlikely that these would be the main control on the evolution of δ¹¹B_sw over longer timescales (e.g., Lemarchand et al., 2002). Alternatively, carefully spaced, short (<1 Myr) time series of relative changes in pH, and hence CO₂ (Foster and Rae, 2016), could bypass the need for accurate δ¹¹B_sw constrains, as they will encompass minimal change in seawater boron composition, typically of a magnitude less than analytical precision (i.e., <0.2 ‰).

Second, the issue of vital effects can be directly investigated within modern species by controlled culturing experiments (Sanyal et al., 1996, 2001; Foster, 2008; Henehan et al., 2013, 2016). For extinct species, vital effects can only be assessed indirectly. Recent work suggests that vital effects of Paleogene foraminifera play a less significant role in determining δ¹¹B_borate than in more recent times (Penman et al., 2014; Anagnostou et al., 2016). Analysis of monospecific (or similarly behaving) foraminiferal species in terms of δ¹⁸O_c, δ¹³C_c and δ¹¹B_c will reduce uncertainties in pH and CO₂ reconstructions.

Third, estimates of salinity and temperature exert a relatively minor influence on estimates of pH (0.005 pH per psu change and 0.01 pH per degrees Celsius) in carbonate system calculations. Nevertheless, seawater temperature estimates are commonly derived from paired analyses of δ¹⁸O_c (Pearson et al., 2007; John et al., 2013; Anagnostou et al., 2016) or Mg∕Ca (Pearson et al., 2009; Foster et al., 2012; Martinez-Boti et al., 2015).

Finally, an additional carbonate parameter is required to calculate CO₂ concentrations at any given seawater salinity and temperature. For pre-Quaternary applications, this is typically provided by assumptions regarding the saturation state of calcium carbonate or alkalinity estimates constrained using geochemical models (Ridgwell, 2005; Pälike et al., 2012; Anagnostou et al., 2016; Gutjahr et al., 2017). Modelling studies indicate that the saturation state of calcium carbonate in surface waters (surface Ω_calc= [Ca]_sw ⋅ [${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$]∕ K_sp) remains essentially constant through the Cenozoic, independent of model boundary conditions (Tyrrell and Zeebe, 2004; Ridgwell and Zeebe, 2005; Hönisch et al., 2012; Pälike et al., 2012). This has recently been confirmed empirically, at least for the last 20 Myr (Sosdian et al., 2018). Estimates of [Ca]_sw are, however, required to derive carbonate ion concentration ([${\mathrm{CO}}_{\mathrm{3}}^{\mathrm{2}-}$]), which is used as the second carbonate system parameter, and these estimates are typically derived from fluid inclusion analysis (e.g., Horita et al., 2002).

At this stage, mean, upper and lower bounds of atmospheric CO₂ are calculated using a Monte Carlo approach to solve the relevant carbonate system equations. All simulations are iterated assuming Gaussian distribution of these parameters within their 2σ uncertainties. For pre-Quaternary timescales, propagated estimates of CO₂ carry an uncertainty of ∼25 % (95 % confidence intervals). Relative CO₂ estimates, however, carry a much-reduced uncertainty of 2 %–10 % (Foster and Rae, 2016).

6.2.1 Theoretical background of alkenone-based CO₂ reconstruction

The saturation state and isotopic composition of alkenone biomarkers, unsaturated ethyl and methyl ketones with carbon chain lengths of C37–C39 are the basis for SST (Brassell et al., 1986; Conte and Eglinton, 1993; Conte et al., 2006) and salinity proxies (e.g., van der Meer et al., 2007; Warden et al., 2016), as well as CO₂ estimation (e.g., Pagani, 2002; Zhang et al., 2013). In the modern oceans, alkenones are produced by a relatively limited number of calcifying haptophyte algae (“coccolithophores”), all within the family Noelaerhabdaceae and including Emiliania huxleyi and Gephyrocapsa species (Volkman et al., 1980; Conte et al., 1998; Pagani, 2014). Where alkenone-bearing sediments predate the first appearance of the modern alkenone-producing species, it is assumed that ancient relatives within the Noelaerhabdaceae, such as the early Eocene to modern genus Reticulofenestra, are the primary alkenone producers (Pagani, 2014). Alkenone-based estimates of CO₂ rely on the relationship between ε_p and [CO_2(aq)] in the surrounding growth medium, which is related to atmospheric CO₂ via Henry's Law (Freeman and Hayes, 1992; Pagani, 2014). Calculation of ε_p requires determination of the difference between the δ¹³C of CO_2(aq) in the growth medium (δ¹³C_CO2(aq)) and the δ¹³C of biosynthesized cellular organic matter (δ¹³C_org). The latter is measured on the haptophyte-specific C37:2 alkenone (δ¹³C_37:2) (Pagani, 2014; Zhang et al., 2013).

As lipids are depleted in ¹³C relative to the carbon isotope composition of the whole cell a correction of +4.2 ‰ is applied when calculating δ¹³C_org from δ¹³C_C37:2 (Schouten et al., 1998; Pagani, 2014). Alkenone-based estimates of ancient CO₂ assume diffusive transport of inorganic carbon across algal cell membranes (discussed below), with the total carbon isotope fractionation (Zhang et al., 2013; Pagani, 2014) being described by the following:

where the carbon isotope fractionation associated with photosynthesis and carbon transport are denoted by ε_f – typically ∼25 ‰ to 28 ‰ (Goericke and Fry, 1994; Farquhar and Richards, 1984) – and ε_t respectively. C_e is the ambient, and C_i the intracellular concentrations of aqueous CO₂. The physiological factors that control the C_i term are complex and poorly known, and so a simplified version of this relationship is typically used (Jasper et al., 1994; Rau et al., 1996):

where the b term integrates a suite of physiological variables, including cellular growth rate and cell geometry, that can influence photosynthetic carbon isotope fractionation (Pagani, 2014, and references therein). Studies reconstructing ancient CO₂ typically assume that the b term is correlated to nutrient levels, represented by phosphate concentrations, with a range of regressions presented by Pagani (2014). Ultimately, CO₂ is calculated from [CO_2(aq)] by applying Henry's Law:

where k_H is defined as

for a salinity of 35 ‰ (Weiss, 1974) and where T is an independently derived estimate of sea surface temperature, best constrained by the $U_{\mathrm{37}}^{K^{\prime }}$ index derived from the same alkenone assemblage, or else another DeepMIP recommended SST proxy methodology if SSTs are beyond the limit of $U_{\mathrm{37}}^{K^{\prime }}$.

6.2.2 Strengths of the alkenone CO₂ method

The relatively high abundance of alkenone biomarkers throughout the Neogene and the relative ease of analysis of δ¹³C_37:2 make alkenone CO₂ estimation suitable for high-temporal-resolution studies of CO₂ dynamics. The alkenone CO₂ method is advantageous over other estimates of ε_p, such as the use of bulk organic carbon or other less specific biomarkers, because of the high specificity of marine alkenone production to one family within the haptophyte algae (Volkman et al., 1980; Conte and Eglinton, 1993). This source specificity also allows for independent constraints on other factors that influence ε_p including cell geometry and physiology, based on the coccolithophore fossil record (Henderiks and Pagani, 2007; Pagani, 2014; Bolton et al., 2016). Finally, as the distribution ratio of alkenones produced by haptophyte algae is closely related to SSTs ($U_K^{\mathrm{37}}$ index; applicable below ∼28 ^∘C) (Conte et al., 2006), similar to δ¹¹B methods, it is possible to directly generate parallel records of two key parameters for CO₂ estimation – SSTs and ε_p37:2.

6.2.3 Weaknesses of the alkenone CO₂ method

There are several uncertainties and assumptions in translating measured δ¹³C_37:2 first into a value for ε_p37:2, and then into an estimate of [CO_2(aq)] and CO₂. In particular, robust proxy estimates of δ¹³C_CO2(aq) and SST are required. These have not always been available in previous alkenone-based estimates of Paleogene CO₂ (Pagani et al., 2005). Inaccurate estimates of SST can cause substantial error in alkenone-based CO₂ estimates (Super et al., 2018), especially at high CO₂ where ε_p approaches the assumed value for ε_f (Bidigare et al., 1997; Riebesell et al., 2000; Pagani et al., 2011; Pagani, 2014). Although the “b” term integrates the effects of a number of physiological variables (e.g., specific growth rate, cell geometry) on net carbon isotope fractionation during photosynthesis (Bidigare et al., 1997; Riebesell et al., 2000; Pagani et al., 2011; Pagani, 2014) it is typically derived from the strong empirical relationship with surface ocean phosphate concentrations [${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$] in the modern oceans (Bidigare et al., 1997). Although there have been attempts to disaggregate cell geometry effects (e.g., Henderiks and Pagani, 2007), most existing records of ancient CO₂ rely on assumed ranges of the b term based on estimates of ancient [${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$] (Pagani, 2014).

Alkenone-based CO₂ estimates also assume the dominance of diffusive CO₂ transport across the algal cell membrane (Zhang et al., 2013). Some studies have shown the operation of carbon concentrating mechanisms (CCMs), including the potential use of bicarbonate carbon (Holtz et al., 2015, 2017), under conditions of CO₂ limitation, resulting in nonlinear relationships between ε_p, growth rate and CO_2aq (Riebesell et al., 2000; Zhang et al., 2013). Such CCM activity could impact the value of ε_p through active transportation of carbon species with varied isotope compositions (Zhang et al., 2013). Ongoing empirical and modeling studies of carbon uptake and partitioning (Bolton and Stoll, 2013; Holtz et al., 2017, 2015; McClelland et al., 2017) and those examining the role of environmental conditions on CCM behavior (Raven and Beardall, 2014; Raven and Hurd, 2012) will hopefully provide constraints on the ubiquity, importance and impact of these mechanisms within the alkenone-producing haptophytes. Finally, changes in growth rate can impact the isotopic offset (Δδ) between alkenone δ¹³C_C37:2 and the δ¹³C_org of the whole cell. Δδ is commonly set at 4.2 ‰, although culture study estimates typically range from 3.1 ‰ to 5.9 ‰ (Pagani, 2014 and references therein), up to a reported maximum of 7 ‰–12 ‰ (Benthien et al., 2007). Where ε_p is high, changes of this order in Δδ can exert a significant impact on reconstructed CO₂ values (Pagani, 2014).

6.2.4 Recommended methodologies for the alkenone CO₂ method

The alkenone CO₂ method requires estimation of several physiological (ε_f, cell size, b) and environmental ([${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$], SST, δ¹³C_CO2(aq)) parameters. Reconstruction of SSTs should follow best practice from one of the established proxies outlined in Section 4 and, where possible, SST uncertainty should be propagated through to final estimates of CO₂. The most reliable estimates of δ¹³C_CO2(aq) are derived from the measurement of test δ¹³C compositions of well-preserved, mixed-layer dwelling foraminifera, within the narrow 250–300 µm size window, which are closest to equilibrium with modern δ¹³C_DIC (Birch et al., 2013). Although very low ε_f values for Emilliania huxleyi of ∼11 ‰ have been published (Boller et al., 2011), these results are yet to be replicated, and we recommend the existing approach of providing propagated CO₂ uncertainty based on a 25 ‰–28 ‰ range of ε_f for algae using C3 pathways (Pagani, 2014).

As with ε_f, establishing robust estimates for the b factor in ancient communities is nontrivial, and we recommend a similar propagation of uncertainty based on a range of b factors drawn from surface water [${\mathrm{PO}}_{\mathrm{4}}^{\mathrm{3}-}$] measurements in analogous modern environments (Zhang et al., 2013). Given the potential importance of cell size within the physiological factors that make up the b term, some studies have sought to use an explicit cell-size “correction” to alkenone CO₂ reconstructions (Bolton et al., 2016). However, the precise application of a cell size correction and its relation to the standard b term remains unclear. Here, we recommend that reticulofenestrid coccolith size is measured along with the determination of δ¹³C_37:2, at least for the purposes of screening out the impact of major shifts in cell size on CO₂ records (Pagani, 2014), and for potential future use to improve quantitative estimates of the b factor.

6.2.5 Other algal substrates for CO₂ estimation

Dinoflagellates. Significant progress has been made towards a CO₂ proxy based on the carbon isotope fractionation of phototrophic dinoflagellates and their preservable cysts (dinocysts) (Hoins et al., 2015, 2016a, b; Wilkes et al., 2017). A strong correlation between ε_pDINO and CO₂ was found in four cultured dinoflagellate species (Hoins et al., 2015). Importantly, there appears to be species-specific effects on ε_pDINO based on uptake and leakage of different carbon species (${\mathrm{HCO}}_{\mathrm{3}}^{-}$ and CO_2(aq)) and light intensity, making single-species records preferable (Hoins et al., 2016a, b). Dinocysts are found in the geological record from the Triassic onwards (MacRae et al., 1996) and several extant species are remarkably long ranging (e.g., Williams et al., 2004). This includes two species used in culture experiments (Hoins et al., 2015), Protoceratium reticulatum and Gonyaulax spinifera, which are preserved in the fossil record as dinocyst species Operculodinium centrocarpum and Spiniferites ramosus, respectively. The geological record for O. centrocarpum and S. ramosus spans the last 60 and 130 Myr, respectively (e.g., Williams et al., 2004). The morphology and ecological affinity of these ubiquitous species is well known and apparently unchanged since the early Paleogene (Zonneveld et al., 2013; Sluijs et al., 2005; Frieling and Sluijs, 2018).

Unlike biomarker-based CO₂ proxies, individual dinocysts of a single species are selected based on morphological characteristics prior to analysis to maximize the similarity in size and shape between analyzed individuals, both within and between samples. A recently developed laser-ablation nanocombustion isotope ratio mass spectrometry (LA-nC-IRMS) method (van Roij et al., 2017) allows for accurate analyses of the carbon isotope composition of single dinocysts (Sluijs et al., 2017).

While results from culture experiments are promising, it should be noted that several, potentially important, factors should be constrained before ε_pDINO can be widely applied as a CO₂ proxy. For example, current efforts aim to resolve the potential offset between the cultured motile cell and the preservable cyst as well as the relation of ε_pDINO to parameters other than CO₂, both in culture and in their natural environment. Although these and perhaps other issues still need to be addressed, the CO₂-dependent ε_pDINO in several dinocyst species may represent an opportunity to reconstruct CO₂ in geological deep time, excluding any morphology, species or genus-related factors in CO₂ dependency that may cause significant uncertainty in CO₂ estimates based on less specific algal compounds.

Phytane and general biomarkers. The carbon isotopic composition of phytoplankton has also been tested in a variety of general biomarkers, compounds which are representative of overall phytoplankton community. Chlorophyll, the primary pigment in all photoautotrophs, is of particular interest for reconstructing secular trends given its spatially and temporally ubiquitous presence in the geologic record. For this reason, studies on both its porphyrin head (Popp et al., 1989) and its side-chain phytol preserved as the diagenetic product phytane (Freeman and Hayes, 1992) are some of the earliest uses of ε_p for CO₂ reconstruction. Here, we will discuss phytane as an example of using general biomarkers for estimating past CO₂.

The overall theoretical background for using general biomarkers, such as phytane, is the same as the alkenone-based reconstructions. Like other lipids, phytane is depleted in ¹³C relative to the carbon isotopic composition of the whole cell; a correction of 3.3±1.3 ‰ standard deviation is applied when calculating δ¹³C_org from phytane, based on the average of 22 cultured species (Witkowski et al., 2018). As seen with the other ε_p-based proxies, the physiological variable b is difficult to constrain. Studies for phytane have used a constant b=170 ‰ kg µM⁻¹ (Bice et al., 2006; Sinninghe Damsté et al., 2008; Naafs et al., 2016), consistent with average b found in modern marine sediments for total organic matter at 168±43 ‰ kg µM⁻¹ (Witkowski et al., 2018) and for alkenone at 165±53 ‰ kg µM⁻¹ (Pagani, 2014). All other parameters are treated as seen in the alkenone-based reconstructions.

Based on these assumptions, a recent study shows that phytane yields a robust and consistent trend over the Phanerozoic, the longest record for CO₂ from a single proxy (Witkowski et al., 2018). This record has the potential to span even further in time given that the δ¹³C of phytane has been measured in 1.73 Ga sediments (Li et al., 2003). This opportunity to observe changes in CO₂ from a single proxy is beneficial to understanding secular trends. In addition to its span in preservation, phytane is one of the most abundant biomarkers in both marine sediments and oils. This availability of phytane is critical during periods where other proxies are rare, extending beyond the alkenone and boron records, most notably used for reconstructing Cretaceous CO₂ (Bice et al., 2006; Sinninghe Damsté et al., 2008; Naafs et al., 2016). The ubiquity of phytane throughout time, location and sample type lends to the possibility for amassing large datasets which can be used for developing robust average estimations for CO₂.

Although the phytane-based reconstructions show remarkably similar trends to other CO₂ proxies over the Phanerozoic, this biomarker requires more investigation. Nearly all studies regarding the mechanisms behind ε_p for CO₂ reconstructions have been species-specific and almost exclusively for alkenone-producing Haptophytes. It is unknown whether integrating the entire phytoplankton community, as opposed to several species, has an effect on the observed signal. Studies testing the different parameters and other potential effects on fractionation are still needed in order to properly constrain the uncertainties on using phytane.

6.3.1 Theoretical background, strengths and weaknesses of leaf stomatal and gas-exchange proxies

Plants interact directly with the atmosphere, making them well suited as climate proxies. Woodward (1987) and subsequent work (e.g., compilation in Royer, 2001) demonstrated that the stomatal density and stomatal index (stomatal density or stomatal + epidermal cell density) in many plant species respond inversely to the partial pressure of atmospheric CO₂; these responses are usually species-specific. If the stomatal response to CO₂ is known (usually from herbaria specimens and elevated CO₂ experiments), it can be applied to fossils of the same species to reconstruct paleo-CO₂. The genetic (Casson and Gray, 2008), signaling (Lake et al., 2002) and functional bases of the stomatal-CO₂ response are known reasonably well: theoretical considerations of the latter predict an inverse power law between stomatal density and CO₂ (Wynn, 2003; Franks et al., 2014; Konrad et al., 2017). The stomatal proxy is popular, having been applied in dozens of studies to produce over 300 estimates of paleo-CO₂ across the Phanerozoic, including the early Paleogene (Royer, 2014). A related proxy called the stomatal ratio relates the ratio of paleo-atmospheric CO₂ to present-day CO₂ to the ratio of stomatal density in a fossil to that in its nearest living equivalent (i.e., not necessarily the same species) via a one-point calibration and an assumed exponent of −1 for the power law. Because of these features, the stomatal ratio is best used for tracking relative (not absolute) changes in CO₂ (McElwain and Chaloner, 1995).

Some limitations of the stomatal proxy are as follows: (1) only species whose stomatal density responds to CO₂ can be used; (2) extant calibrations are laborious to construct, and most fossil plant species are not alive today; (3) the method is empirical, not explicitly mechanistic, meaning that it is difficult to exclude the possibility that the extant relationships differed in the geologic past; and (4) in many species the upper error limit on the CO₂ estimate becomes very large or unbounded when the estimated CO₂ is moderately high or higher. An alternative approach that partly addresses these limitations is a mechanistic model based on the fundamentals of leaf gas exchange (Konrad et al., 2008, 2017; Franks et al., 2014). The core of these models is the well-known relationship between leaf assimilation rate and the product between total leaf conductance (which partly reflects the “supply” side of assimilation) and the CO₂ gradient between the atmosphere and intercellular spaces (the “demand” side) (Farquhar and Sharkey, 1982). The critical measurements needed to estimate atmospheric CO₂ are the leaf δ¹³C, stomatal size and density, and an estimate of the assimilation rate from a nearest living relative. The gas-exchange proxy of Franks et al. (2014) typically has 95 % confidence intervals within −25 % to +35 % of the median CO₂ estimate, on par with other popular paleo-CO₂ proxies.

6.3.2 Recommended methodology of leaf stomatal and gas-exchange proxies

Both the stomatal ratio and leaf-gas-exchange proxies can be applied to most fossil leaves, even if the taxa are extinct today. Because the CO₂ estimates from stomatal ratios are “semiquantitative” (McElwain and Chaloner, 1995), they are best used when complementary CO₂ information is present (e.g., Montañez et al., 2016). The stomatal density and stomatal index proxies work well when calibrated taxa are present as fossils and when estimated CO₂ is low to moderate. At higher CO₂ (above ∼500 ppm for many species; e.g., Smith et al., 2010; Doria et al., 2011), the upper error limit can be unbounded, meaning that only minima can be quantified. One challenge with this proxy is the construction of robust calibrations, especially at the high end of the calibrations, because plants may need multiple growing seasons in experimentally elevated CO₂ concentrations before their stomatal densities and indices fully equilibrate (Royer, 2003; Hincke et al., 2016). The leaf-gas-exchange proxies are comparatively more recent, and so have been tested less. We encourage rigorous testing in nearest living relatives, especially for sensitive model inputs that cannot be measured directly in fossils, for example the assimilation rate at a known CO₂ concentration and the scaling between maximum and operational stomatal conductance to CO₂ (Maxbauer et al., 2014; McElwain et al., 2016a, 2016b, 2017; Franks and Royer, 2017).

In cases where multiple stomatal proxies have been applied to the same fossils, the corresponding CO₂ estimates are broadly similar (e.g., Maxbauer et al., 2014; Montañez et al., 2016; Richey et al., 2018), conferring mutual support for the methods. Generating CO₂ estimates from multiple taxa typically leads to a more accurate consensus estimate outcrop (e.g., Montañez et al., 2016; Kowalczyk et al., 2018; Richey et al., 2018; Steinthorsdottir et al., 2018). Routines for propagating all known uncertainties associated with the leaf-gas-exchange and stomatal density or index proxies have been developed (Beerling et al., 2009; Franks et al., 2014) and should be used.

6.5.1 Theoretical background of C₃ land plant proxy method

Plants are differentiated into three main categories based on the method by which they fix carbon: C₃, C₄ and CAM. The carbon isotope composition (δ¹³C) of C₃ plants, which use the RuBisCO enzyme (ribulose-1,5-bisphosphate carboxylase-oxygenase) to catalyze CO₂ fixation, has long been suspected to be a function of atmospheric CO₂ concentration (Farquhar et al., 1982; Krishnamurthy and Epstein, 1990; Körner et al., 1991; Leavitt and Danzer, 1992; Feng and Epstein, 1995). Subsequent controlled growth chamber experiments conducted across a wide range of CO₂ levels show that carbon isotope discrimination between the δ¹³C value of the atmosphere (δ¹³C_CO2) and the resulting plant tissue (δ¹³C_org) [Δ¹³C = (δ¹³C_atm – δ¹³C_org)/(1 + δ¹³C_org∕1000)] increases with increasing CO₂, according to the following equation (Schubert and Jahren, 2012):

where A=28.6, B=0.22 and C=23.9 (Schubert and Jahren, 2015). These data were consistent with other field and chamber experiments measured across smaller changes in CO₂ and revealed a unifying relationship between changes in CO₂ and changes in the Δ¹³C value (Schubert and Jahren, 2012). Using a large data compilation of terrestrial organic matter spanning the Pleistocene–Holocene transition, Schubert and Jahren (2015) showed how this relationship could be used to reconstruct CO₂ within the fossil record by relating changes in CO₂ to changes in the Δ¹³C value (i.e., Δ(Δ¹³C) = Δ¹³C_(t) – ${\mathrm{\Delta }}^{\mathrm{13}}{\mathrm{C}}_{\left({\mathrm{t}}_{=}\mathrm{0}\right)}$) between a time of interest (t) and a reference time (t=0):

Excellent agreement between the proxy results and ice core data validated changes in the Δ¹³C value as an accurate proxy for CO₂ in the fossil record (Schubert and Jahren, 2015). Much of the application of this work has focused on the early Paleogene greenhouse and accompanying hyperthermals (Schubert and Jahren, 2013; Cui and Schubert, 2016, 2017, 2018). These studies have all revealed background CO₂ across the late Paleocene and early Eocene (prior to the EECO) less than ∼3.5 times preindustrial levels (<1000 ppm), but with significantly higher levels at the hyperthermals.

6.5.2 Strengths of C₃ land plant proxy

The C₃ land plant proxy has the potential for wide applicability in the fossil record because it relies on measurements of the δ¹³C value of terrestrial organic matter, an abundant and widely measured substrate that can integrate carbon contributions across a large number of plants. Notably, the proxy (i) has been validated through comparison with ice core data (Schubert and Jahren, 2015), (ii) includes a robust uncertainty analysis for all input parameters and accounts for potential contemporaneous environmental change (Cui and Schubert, 2016), (iii) can resolve differences in carbon isotope excursion magnitude between terrestrial and marine substrates (Schubert and Jahren, 2013), and (iv) provides a method for producing nearly continuous records of CO₂ at resolution sufficient for resolving early Paleogene CIEs, such as the PETM (Cui and Schubert, 2017); (v) can help identify potential sources of the CIEs (Cui and Schubert, 2018); (vi) is applicable across the entire range of Cenozoic CO₂, including the early Paleogene greenhouse (Cui and Schubert, 2016); and (vii) includes an understanding of the mechanism driving the observed increase in discrimination with increasing CO₂ (Schubert and Jahren, 2018) that fits within the framework for carbon isotope discrimination developed by Farquhar et al. (1982).

6.5.3 Weakness of C₃ land plant proxy

Weaknesses include (i) potential for different responses among different plant groups (Hare et al., 2018; Porter et al., 2017); (ii) the effect of O₂ on Δ¹³C being poorly constrained (Beerling et al., 2002; Berner et al., 2000; Porter et al., 2017); and (iii) the proxy saturating at very high CO₂, yielding positive uncertainties of >120 % above 1000 ppmv (Cui and Schubert, 2016). Large uncertainties in CO₂ can be reduced if changes in other environmental parameters that are known to affect the Δ¹³C value (e.g., water availability, plant community shifts) are constrained (Cui and Schubert, 2016). The proxy also requires estimates of δ¹³C_CO2, which for the Paleocene and Eocene, can be calculated from the δ¹³C value of benthic foraminifera (Tipple et al., 2010). Because the proxy is based on relative changes in the Δ¹³C value, quantitative estimates of CO₂ require either (i) knowledge of Δ¹³C value at a known CO₂ level (Schubert and Jahren, 2015; Eq. 25) or (ii) independent estimates of the change in CO₂ across the two time periods of interest (i.e., ΔCO₂ = CO_2(t) – ${\mathrm{CO}}_{\mathrm{2}\left({\mathrm{t}}_{=\mathrm{0}}\right)}$; as shown in Schubert and Jahren, 2013). Marine carbonate records, which can be used to quantify ΔCO₂, might be affected by dissolution (Cui and Schubert, 2017); changes in the Δ¹³C value of paleosol carbonate mimic changes in the Δ¹³C value of terrestrial organic matter, but might also be affected by diffusion of increased CO₂ into the soil and increased productivity (Schubert and Jahren, 2013). Therefore, careful consideration of the selected substrate and site is required; when information on environmental change is not available, the use of multiple records might lead to better CO₂ estimates (Cui and Schubert, 2017).

6.5.4 Recommended methodology of C₃ land plant proxy

This proxy is ideal for reconstructing the wide range of CO₂ characteristics for the early Eocene and associated hyperthermals. Like many other CO₂ proxies, uncertainties increase with increasing CO₂; positive and negative errors, however, are less than ∼120 % and 40 %, respectively, at the CO₂ levels predicted for the late Paleocene and early Eocene (Cui and Schubert, 2016). The choice of record is important for producing robust CO₂ estimates via Eq. (25) because there is potential for systematic bias in the calculated CO₂ estimate resulting from (i) sampling resolution, (ii) substrate type, (iii) diagenesis and (iv) environmental change (Cui and Schubert, 2017). For example, working in sites with high mean annual precipitation may help to reduce variability in median CO₂ estimates, as Δ¹³C value is less responsive to changes in MAP in ever-wet sites (e.g., Diefendorf et al., 2010). However, Lomax et al. (2019) demonstrated experimentally that plants grown under water stress showed the same relative change in the Δ¹³C value with increasing CO₂ concentration as observed in well-watered experiments (e.g., Schubert and Jahren, 2012), but with lower Δ¹³C values overall. Their experimental data are consistent with results that modeled the predicted effect of water stress on the Δ¹³C value across multiple levels of CO₂ (see Fig. 4a within Schubert and Jahren, 2018). Nevertheless, both chamber (Lomax et al., 2019) and deep-time (Barral et al., 2017) studies have tried to use Eq. (24) to quantify CO₂ and failed to account for the effects of growing conditions and plant species type on the Δ¹³C value. Such attempts highlight the need to use Eq. (25) with robust estimates of Δ(Δ¹³C) (after Schubert and Jahren, 2015) in order to accurately quantify CO₂.