# Spatial biases in oxygen-based Phanerozoic seawater temperature reconstructions

Alexandre Pohl, Thomas W Wong Hearing, Arnaud Brayard, Ethan Grossman, Michael M. Joachimski, Guillaume Le Hir, Thomas Letulle, Daniel J Lunt, Mathieu Martinez, Emmanuelle Pucéat, et al.

► To cite this version:

Alexandre Pohl, Thomas W Wong Hearing, Arnaud Brayard, Ethan Grossman, Michael M. Joachimski, et al. Spatial biases in oxygen-based Phanerozoic seawater temperature reconstructions. Earth and Planetary Science Letters, 2025, 663, pp.119418. 10.1016/j.epsl.2025.119418. hal-05061186

HAL Id: hal-05061186

https://hal.science/hal-05061186v1

Submitted on 9 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# Spatial biases in oxygen-based Phanerozoic seawater temperature reconstructions

Alexandre POHL<sup>a,∗</sup> <sup> ID </sup>, Thomas W. WONG <sup> HEARING b </sup>, Arnaud BRAYARD<sup> a </sup>,
Ethan GROSSMAN<sup> c </sup> <sup> ID </sup>, Michael M. JOACHIMSKI<sup> d </sup> <sup> ID </sup>, Guillaume LE HIR<sup> c </sup>, Thomas LETULLE<sup> f,g </sup>,
Daniel J. LUNT<sup> h </sup> <sup> ID </sup>, Mathieu MARTINEZ<sup> i </sup>, Emmanuelle PUCEAT<sup> a </sup>, Guillaume SUAN<sup> g </sup>,
Paul VALDES<sup> h </sup> <sup> ID </sup>, Yannick DONNADIEU<sup> j </sup>

ARTICLE INFO

Editor: Dr Tristan Horner

Keywords:

Stable oxygen isotopes

Salinity effect

Seawater temperature

Paleoclimate

General circulation model

Phanerozoic

A B S T R A C T

Stable oxygen isotopes (δ¹⁸O) are routinely used to reconstruct sea-surface temperatures (SSTs) in the geological past, with mineral δ¹⁸O values reflecting a combination of the temperature and oxygen isotope composition of seawater (δ¹⁸Osw). Temporal variation of mean-ocean δ¹⁸Osw is usually accounted for following estimates of land-ice volume. Spatial variations in δ¹⁸Osw, however, are often neglected or corrected using calibrations derived from the present-day or recent past. Geochemical methods for constraining δ¹⁸Osw and isotope-enabled general circulation model (GCM) simulations are still technically challenging. This lack of constraints on ancient δ¹⁸Osw is a substantial source of uncertainty for SST reconstructions. Here we use the co-variation of δ¹⁸Osw and seawater salinity, together with GCM simulations of ocean salinity, to propose estimations of spatial variability in δ¹⁸Osw over the Phanerozoic. Sensitivity tests of the δ¹⁸Osw-salinity relationship and climate model, and comparison with results of isotope-enabled GCMs, suggest that our calculations are robust at first order. We show that continental configuration exerts a primary control on δ¹⁸Osw spatial variability. Complex ocean basin geometries in periods younger than 66 Ma lead to strong inter-basinal contrasts in δ¹⁸Osw. Latitudinal SST gradients may be steeper than previously suggested during most of the Mesozoic and Cenozoic. This work has limitations, with δ¹⁸Osw-salinity relationships being less reliable in both low-latitude epicontinental settings and high-latitude regions of deep-water formation. Whilst our calculations are limited use in correcting δ¹⁸O measurements for local δ¹⁸Osw, they identify the time slices and paleogeographical regions that should be prioritized for future work using isotope-enabled GCMs.

1. Introduction

Our understanding of ocean temperatures throughout the Phanero-
zoic is primarily derived from stable oxygen isotope (δ¹⁸O)-based tem-
perature reconstructions (Gaskell et al., 2022; Grossman and
Joachimski, 2022; Judd et al., 2024, 2022; Scotese et al., 2021; Veizer

and Prokoph, 2015). Stable oxygen isotope data acquired on carbonate or phosphate biominerals represent half of the entries in a recent Phanerozoic-scale compilation of sea-surface temperature (SST) proxy data and constitute the only proxy extending from the Cambrian to Modern (Judd et al., 2024, 2022). SST reconstructions are pivotal for our understanding of the evolution of Earth's climate, which in turn is

<footnote><sup>a</sup> Biogéosciences UMR 6282, Université Bourgogne Europe, CNRS, F-21000 Dijon, France</footnote>

<footnote><sup>b</sup> School of Geography, Geology and the Environment, University of Leicester, Leicester LE1 7RH, UK</footnote>

<footnote><sup>c</sup> Department of Geology and Geophysics, Texas A&M University, College Station, TX 77843, USA</footnote>

<footnote><sup>d</sup> GeoZentrum Nordbayern, Friedrich-Alexander Universität Erlangen-Nürnberg (FAU), Schlossgarten 5, 91054 Erlangen, Germany</footnote>

<footnote><sup>e</sup> Université de Paris, Institut de physique du globe de Paris, CNRS, F-75005 Paris, France</footnote>

<footnote><sup>f</sup> Department of Geosciences, University of Padua, Padua, Italy</footnote>

<footnote><sup>g</sup> Univ Lyon, UCBL, ENSL, UJM, CNRS, LGL-TPE, F-69622, Villeurbanne, France</footnote>

<footnote><sup>h</sup> School of Geographical Sciences, University of Bristol, Bristol BS8 1SS, UK</footnote>

<footnote><sup>i</sup> Univ Rennes, CNRS, Géosciences Rennes - UMR 6118, F-35000 Rennes, France</footnote>

<footnote><sup>j</sup> Aix-Marseille Université, CNRS, IRD, INRAE, Collège de France, CEREGE, Aix-en-Provence, France</footnote>

<footnote>* Corresponding author.
E-mail address: alexandre.pohl@cnrs.fr (A. POHL).</footnote>
crucial for investigating the co-evolution of marine biodiversity and the
physical environment (Ontiveros et al., 2023; Trotter et al., 2008).
However, major uncertainties persist in the reconstruction of ocean
temperatures during the Phanerozoic, as testified by the range in pub-
lished estimates across the Phanerozoic (Grossman and Joachimski,
2022; their Fig. 5).

Several factors contribute to the uncertainty in δ¹⁸O-based SST re-
constructions. Early meteoric diagenesis can distort or reset primary
environmental isotopic signatures, particularly by shifting δ¹⁸O towards
lighter values. Although the potential impact of early diagenesis can
never be fully excluded, consideration of several visual, textural,
elemental, and isotopic indicators can be used to screen out samples that
are most likely to be diagenetically altered (Brand et al., 2011; Wheeley
et al., 2012). Additional limitations in δ¹⁸O-based temperature re-
constructions arise from uncertainties linked to kinetic effects (Davies
et al., 2023; Huyghe et al., 2022), uncertainties inherent in the oxygen
isotope fractionation equations (Daëron and Gray, 2023; Letulle et al.,
2023; Pucéat et al., 2010), or differences in analytical methods (Wheeley
et al., 2012).

There is also a lack of consensus regarding the long-term evolution of
the δ¹⁸O<sub>sw</sub> value of the global ocean. There is a long-term trend towards
more negative δ¹⁸O in values in fossils and microfossils towards more
ancient periods during the Phanerozoic. The secular trend is interpreted
either as reflecting very warm ocean temperatures (possibly 50 °C at low
latitudes) in the early Paleozoic balanced by a near-constant mean-
ocean δ¹⁸O<sub>sw</sub> value (Grossman and Joachimski, 2022), or reflecting a
modern-like range of SST balanced by a changing δ¹⁸O<sub>sw</sub> value for the
global ocean (Isson and Rauzi, 2024; Veizer and Prokoph, 2015). The
constant global-δ¹⁸O<sub>sw</sub> hypothesis is supported by clumped-isotope ev-
idence detecting no significant change in δ¹⁸O<sub>sw</sub> over the past 500
million years (Bergmann et al., 2018; Henkes et al., 2018); the evolving
δ¹⁸O<sub>sw</sub> hypothesis is supported by temperature estimates based on the
actualistic interpretation of climatic zones reconstructed based on lith-
ological indicators of the climate (Scotese et al., 2021) and pairwise
comparisons of δ¹⁸O measurements from multiple mineral systems with
different temperature sensitivities (Galili et al., 2019; Isson and Rauzi,
2024; Keller et al., 2019; Shields and Veizer, 2002).

Sampling bias is also a problem when interpreting δ¹⁸O-based sea-
temperature reconstructions. Jones and Eichenseer (2021) demon-
strated that temporal variability in spatial sampling bias makes it diffi-
cult to infer changes in global temperatures from available proxy data.
They showed that uneven spatial sampling distorts our vision of changes
in global ocean temperature across the Phanerozoic.

Finally, the lack of constraints on regional variations in δ¹⁸Osw can
generate major uncertainties in δ¹⁸O-based temperature re-
constructions, especially in epicontinental settings that represent most
of the Phanerozoic δ¹⁸O database. The impact of the local δ¹⁸Osw
composition is often neglected or approximated using calibrations
established in the modern (O'Brien et al., 2017; Zachos et al., 1994) or
recent geological past (Grossman and Joachimski, 2022; Roberts et al.,
2011). Whilst δ¹⁸Osw can be determined using carbonate clumped
isotope analyses, these analyses cannot currently be performed at scale
and remain vulnerable to uncertainties related to reordering at burial
temperatures, especially in the case of early Paleozoic samples
(Finnegan et al., 2011; Henkes et al., 2018). To a lesser extent, vital
effects may also influence clumped isotope temperature and thus δ¹⁸Osw
determinations (Davies et al., 2023; Huyghe et al., 2022). Alternatively,
regional variations in δ¹⁸Osw can be simulated using isotope-enabled
GCMs (Gaskell et al., 2022; Zhou et al., 2008). This approach is prom-
ising and has proven useful in accounting for local δ¹⁸Osw variations in
δ¹⁸O-based temperature calculations, for instance by Gaskell et al.
(2022) who reconstructed SST evolution over the last 95 million years.
The main limitation (beyond uncertainties in model boundary condi-
tions) is the computational cost that isotope-enabled simulations
represent, which also renders their paleoclimatic application at scale
very challenging. For instance, Gaskell et al. (2022) interpolated

between isotope-enabled simulations conducted at several pCO2 levels
for only two paleogeographical configurations (Eocene, 55 million years
ago [Ma] and Miocene, 15 Ma) when working across the last 95 million
years. Similarly, Valdes et al. (2021) calculated that 18 wall-clock
months would be required to run a fully-equilibrated deep-time
climate simulation using the isotope-enabled version of HadCM3.
Isotope-enabled GCM simulations are limited to only a few intervals in
the geological past (e.g., key Cenozoic time slices including the Last
Glacial Maximum, Miocene, and Eocene, and rare Mesozoic and Paleo-
zoic time slices; Gaskell and Hull, 2023; Macarewich et al., 2021; Zhou
et al., 2008). These technical limitations leave the impact of local δ¹⁸Osw
on Phanerozoic paleo-temperature reconstructions mostly
unconstrained.

Here we use an alternative approach to investigate the magnitude of δ¹⁸Osw variations in space and time during the entire Phanerozoic using salinity, a variable that is commonly calculated in GCM simulations. We make use of the co-variation of seawater salinity and δ¹⁸Osw to estimate spatial and temporal δ¹⁸Osw patterns over the last 541 million years. Salinity and oxygen isotopic composition co-vary in surface seawater because they are controlled by similar processes (Railback et al., 1989; Tiwari et al., 2013). Evaporation both increases seawater salinity and concentrates the heavier isotope (¹⁸O) in seawater, while isotopically light freshwater input directly from atmospheric precipitation or from continental runoff reduces seawater salinity and enriches seawater in ¹⁶O. Recognition of this ‘salinity effect’ permitted the development of δ¹⁸Osw-salinity relationships by which local δ¹⁸Osw can be estimated if local salinity is known (LeGrande and Schmidt, 2006; Railback et al., 1989). Global expressions for δ¹⁸Osw-salinity relationships do not precisely capture spatial variations in δ¹⁸Osw because the relationship is controlled by both oceanic basin geometry and latitude through variations in oceanic circulation (Zhu et al., 2020) and the regional hydrological cycle (LeGrande and Schmidt, 2006; Tiwari et al., 2013). Zhu et al. (2020) notably demonstrated that changes in the sites of deep-water formation have a very significant impact on upper-ocean δ¹⁸Osw. Intensified convection brings ¹⁸O-enriched subsurface water upwards, which can then mix with sea-surface water masses and be advected to lower latitudes through gyre circulations. Nevertheless, δ¹⁸Osw-salinity relationships provide a practical method for estimating local δ¹⁸Osw using results of GCMs that are not isotope-enabled.

This approach allows us to provide a first-order quantification of the
influence of local δ¹⁸Osw on SST reconstructions throughout the Phan-
erozoic. Our results demonstrate that continental configuration exerts a
major control on large-scale δ¹⁸Osw spatial patterns and allow us to
identify which geological intervals and geographical regions are most
vulnerable to SST reconstruction biases.

**2. Methods**

Here, we use the salinity fields from 109 stage-level Phanerozoic
climate simulations conducted using the coupled ocean-atmosphere
Hadley Centre general circulation model HadCM3. These simulations
are largely similar to those of Valdes et al. (2021) but have undergone a
number of developments permitting to better simulate latitudinal tem-
perature gradients (Judd et al., 2024; Supplementary Text). We translate
simulated salinity into local δ¹⁸Osw following the δ¹⁸Osw-salinity rela-
tionship established by Railback et al. (1989; Supplementary Fig. S1).
This provides an estimate of δ¹⁸Osw for each ocean model grid cell over
the Phanerozoic. To illustrate the impact of these spatial and temporal
δ¹⁸Osw patterns on deep-time SST reconstructions, we then calculate a
‘SST interpretation bias’ outlined below that would arise from esti-
mating SST using a globally-uniform δ¹⁸Osw rather than the calculated
local δ¹⁸Osw. All δ¹⁸O values are given relative to Vienna standard mean
ocean water (VSMOW).

For every model ocean grid cell, we first translate the SST simulated
in HadCM3 and the salinity-derived local δ¹⁸Osw value into a biogenic
phosphate δ¹⁸O value (δ¹⁸O<sub>phosphate</sub>) using the oxygen fractionation
equation of Pucéat et al. (2010). We then use this hypothetical δ¹⁸Ophosphate value to reconstruct a new SST value assuming a globally uniform value of δ¹⁸Osw (i.e. neglecting local δ¹⁸Osw effects). The difference between the δ¹⁸Ophosphate-recalculated SST and the model-simulated SST for each grid cell constitutes the ‘SST interpretation bias’ and represents the potential impact of unknown spatial δ¹⁸Osw patterns on proxy-derived paleoclimate reconstructions.

Consider an oxygen isotope fractionation equation between water and biogenic minerals of the form:

$$SST = a(\delta^{18}\text{O}_{\text{mineral}} - \delta^{18}\text{O}_{\text{sw}}) + b \quad (1)$$

where SST is the calculated sea-surface temperature, a and b are equation-specific constants, δ¹⁸Omineral is the measured fossil mineral δ¹⁸O value, and δ¹⁸Osw is the seawater δ¹⁸O value. The SST interpretation bias for each model grid cell can be expressed as:

$$SST \text{ interpretation bias} = a(\delta^{18}\text{O}_{\text{sw-model}} - \delta^{18}\text{O}_{\text{sw-hyp}}) \quad (2)$$

where δ¹⁸Osw-model is the local δ¹⁸Osw calculated from the GCM salinity, and δ¹⁸Osw-hyp is the globally-uniform δ¹⁸Osw hypothetical value. Eq. (2) demonstrates that the calculated SST interpretation bias only depends on the gradient of the oxygen isotope fractionation equation and GCM salinity, and not on the model SST.

We use the phosphate oxygen isotope temperature equation of Pucéat et al. (2010), which has a gradient, a, of -4.22 °C/‰ to illustrate the SST interpretation bias. The narrow range of temperature sensitivities (approximately -4 to -5 °C/‰) of oxygen isotope fractionation equations at marine temperatures means that using other fractionation equations would give similar results (Bemis et al., 1998; Grossman, 2012). Daëron and Gray (2023) further demonstrated that taxonomic differences in foraminiferal (calcite) oxygen isotope temperature equations produce absolute offsets in reconstructed temperatures without affecting the temperature sensitivity (gradient) of the equations, which remains indistinguishable from that of inorganic calcite.

For our main calculations we use the δ¹⁸Osw-salinity relationship of Railsback et al. (1989; Supplementary Fig. S1), which has the advantage of incorporating δ¹⁸Osw-salinity trends from multiple present-day ocean basins with a wide range of salinity values into simple generalized relationships. These general relationships can be adapted to represent ancient conditions by modifying a few key parameters. The δ¹⁸Osw-salinity relationship of Railsback et al. (1989) can be simplified to the general form:

$$\delta^{18}\text{O}_{\text{sw}} = \delta^{18}\text{O}_{0} + \alpha(S - S_{0}) \quad (3)$$

where δ¹⁸Osw is the local value of seawater, δ¹⁸O0 and S0 are the global ocean average δ¹⁸O and salinity, respectively, S is the local seawater salinity value, and α is a regionally-dependent constant, which also varies with the global climatic state (see below). Global-ocean δ¹⁸O (δ¹⁸O0 in Eqn. (3)) was varied during the Phanerozoic to account for variations in land-ice volumes following Grossman and Joachimski (2022) but was not otherwise detrended (Supplementary Fig. S2). The value of α is determined by whether regional salinity (S) is above or below the global average (S0) to account for first-order impacts of changing hydrological and thermal regimes as a function of latitude (Railsback et al., 1989). We followed Railsback et al. (1989) in using α = 0.35 ‰/ppt for regions where salinity exceeds the global average, and α = 21.0/S0 or α = 8.0/S0 where salinity is below the global average respectively for icehouse and greenhouse climates (α = m_e for S > S0 and α = Δfw/S0 for S < S0 in the equations on p. 587 of Railsback et al., 1989). In these calculations, we defined icehouse periods during the Phanerozoic based on the compilation of Cather et al. (2009). These time slices are characterized by global-ocean δ¹⁸O > -1.08 ‰ in Grossman and Joachimski (2022; Supplementary Fig. S2). A value α = 8.0/S0, reflecting more equable climates than the late Cenozoic icehouse state, leads to more ¹⁸O-enriched low-salinity waters than the value α =

21.0/S0 used for icehouse climates.

We evaluated the sensitivity of our calculations to different models of the δ¹⁸Osw-salinity relationship (Gaskell et al., 2022; Tiwari et al., 2013; Supplementary Fig. S1). We also calculated the SST interpretation bias after correcting δ¹⁸Ophosphate for latitude using the formula proposed for the present-day southern Pacific and Atlantic oceans by Zachos et al. (1994), as this correction is widely used in paleoclimate studies (e.g., O'Brien et al., 2017):

$$\delta^{18}\text{O}_{\text{sw}}(\text{‰}, \text{VSMOW}) = 0.576 + 0.041L - 0.0017L^2 + 0.0000135L^3 \quad (4)$$

where L is the paleolatitude.

We further evaluated the impact that the choice of climate model had on our results by repeating our calculations for two time slices with contrasting continental configurations (Late Ordovician, 445 Ma and Late Cretaceous, 95 Ma) using two additional GCMs with different levels of complexity: IPSL-CM5A2 (Laugié et al., 2021) and FOAM v1.5 (Pohl et al., 2014) (Supplementary Text).

# 3. Results and discussion

## 3.1. Validating the approach in the modern

The δ¹⁸Osw-salinity relationship of Railsback et al. (1989) applied to our 0 Ma HadCM3 simulation reasonably captures first-order spatial patterns in present-day sea-surface δ¹⁸Osw but produces substantial deviations in terms of absolute values (Fig. 1). Global-mean calculated sea-surface δ¹⁸Osw is offset by -0.41 ‰ compared with the gridded dataset of LeGrande and Schmidt (2006). This negative bias is relatively uniform in space, with the exception of particularly strong negative biases observed at equatorial latitudes in southeastern Asia and along the western coasts of Africa and central America. The Mediterranean Sea and Russian Arctic are characterized by strong, positive deviations relative to the dataset of LeGrande and Schmidt (2006).

Regional model-data mismatches, both positive and negative, largely arise from anomalies in sea-surface salinity (Supplementary Fig. S3). These regional anomalies, however, do not escalate into any clear systematic bias in the simulated latitudinal SSS gradient (Supplementary Fig. S3c).

Model-data mismatch may further reflect uncertainties in the gridded dataset of LeGrande and Schmidt (2006), who reconstructed modern δ¹⁸Osw using sparse data by means of regional δ¹⁸Osw-salinity relationships. Data coverage is particularly poor in some regions of maximum model-data mismatches, at equatorial latitudes in southeastern Asia and along the western coast of Africa.

The use of a global δ¹⁸Osw-salinity relationship, rather than basin-specific equations, also comes with its own biases. Although evaluating the skill of published δ¹⁸Osw-salinity relationships is beyond the scope on this work, we illustrate the impact of the δ¹⁸Osw-salinity relationship on the bias in calculated δ¹⁸Osw by repeating our calculations using the δ¹⁸Osw-salinity relationship of Gaskell et al. (2022; their supplementary eqn. 6) (Supplementary Fig. S4).

Using the equation of Gaskell et al. (2022) overall improves model-data agreement. The global-mean calculated sea-surface δ¹⁸Osw offset is reduced (from -0.41 ‰ to -0.20 ‰. Spatial patterns of biases in δ¹⁸Osw more closely track sea-surface salinity anomalies. The same regional biases are found at equatorial latitudes in southeastern Asia and along the western coasts of Africa and central America, and in the Mediterranean Sea and Russian Arctic. This model-data comparison shows that global δ¹⁸Osw-salinity relationships permit calculating first-order spatial patterns in δ¹⁸Osw based on salinity fields simulated using a global GCM. Calculated spatial patterns are particularly robust in open-ocean settings, whereas coastal regions at equatorial and polar latitudes, and enclosed basins, are often associated with larger errors arising from local biases in simulated sea-surface salinity.

Therefore, the δ¹⁸Osw-salinity relationship of Gaskell et al. (2022)
<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Biases in calculated sea-surface $\delta^{18}\text{O}_{\text{sw}}$ using the $\delta^{18}\text{O}_{\text{sw}}$-salinity relationship of Railsback et al. (1989). (a) Present-day $\delta^{18}\text{O}_{\text{sw}}$ (LeGrande and Schmidt, 2006), regridded to the HadCM resolution using a bi-linear interpolation. (b) $\delta^{18}\text{O}_{\text{sw}}$ calculated using the 0 Ma HadCM3 simulation, the $\delta^{18}\text{O}_{\text{sw}}$-salinity relationship of Railsback et al. (1989) and a globally-averaged oceanic $\delta^{18}\text{O}_{\text{sw}}$ of 0 ‰. (c) Difference between $\delta^{18}\text{O}_{\text{sw}}$ values calculated using the 0 Ma HadCM3 simulation and (i.e., minus) present-day $\delta^{18}\text{O}_{\text{sw}}$ values of LeGrande and Schmidt (2006). (d) Cross-plot of $\delta^{18}\text{O}_{\text{sw}}$ values calculated using the 0 Ma HadCM3 simulation and the present-day $\delta^{18}\text{O}_{\text{sw}}$ values of LeGrande and Schmidt (2006). The red line represents the 1:1 line that would reflect perfect model-data agreement. (e) As per panel (d) but shown for a narrower $\delta^{18}\text{O}_{\text{sw}}$ range extending from -4 ‰ to 4 ‰. In panels (a)–(c), emerged landmasses are outlined and shaded white; Robinson projections with latitude shown every 30°.</figcaption></figure>

outperforms that of Railsback et al. (1989) for the present-day, icehouse,
climate. However, the $\delta^{18}\text{O}_{\text{sw}}$-salinity relationship of Railsback et al.
(1989) is adaptable to a greenhouse climate state (Supplementary
Fig. S1), unlike that of Gaskell et al. (2022). This climate state adapt-
ability is important for developing a scheme that can be applied
consistently across the Phanerozoic, which was predominantly charac-
terized by greenhouse climates (Judd et al., 2024). Therefore, we use the
equation of Railsback et al. (1989) in our baseline calculations, but we
also show key results using the equation of Gaskell et al. (2022) for
comparison.

3.2. *Spatial patterns in SST interpretation bias*

Estimating SST using a globally-uniform $\delta^{18}\text{O}_{\text{sw}}$ rather than the calculated local $\delta^{18}\text{O}_{\text{sw}}$ generally underestimates SST at low latitudes and overestimates SST at higher latitudes (Fig. 2; Supplementary Fig. S5). This is because evaporation tends to exceed precipitation at low latitudes, leading to higher salinity and hence heavier $\delta^{18}\text{O}_{\text{sw}}$ values (Railsback et al., 1989), while the opposite applies to the mid- and polar

latitudes (Supplementary Figs. S6–S7).

The longitudinal variation of the bias is generally stronger in more
recent, Cenozoic time slices (< 66 Ma) because of the development of
more complex ocean basin geometries which exert a strong impact on
ocean salinity. The ocean basin effect is exemplified by the stronger SST
interpretation bias simulated in the Atlantic Ocean compared to the
Pacific Ocean at 52 Ma, 26 Ma and to a lesser extent at 0 Ma (Fig. 2a–c).
The $\delta^{18}\text{O}_{\text{sw}}$-salinity relationship of Gaskell et al. (2022) produces very
similar spatial patterns but stronger contrasts in calculated $\delta^{18}\text{O}_{\text{sw}}$
(Supplementary Fig. S8).

We simulate large deviations from the general latitudinal pattern in all time slices in coastal areas. SST is overestimated near river mouths due to the massive input of low-salinity and light-$\delta^{18}$O water. This freshwater effect is particularly strong in the vicinity of the intertropical convergence zone due to the strong precipitation rates – in line with isotope-enabled simulations (Macarewich et al., 2021; Zhu et al., 2020). This can be seen at equatorial $\delta^{18}$O sampling locations at, for example, 425 Ma (around Baltica), and 301 Ma and 252 Ma on the western coast of the equatorial landmasses (Fig. 2).
<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Effect of estimating SST using a globally-uniform $\delta^{18}\text{O}_{\text{sw}}$ rather than the calculated local $\delta^{18}\text{O}_{\text{sw}}$. SST over- / under- estimation (positive and negative values respectively) arising from neglecting spatial $\delta^{18}\text{O}_{\text{sw}}$ variations when converting $\delta^{18}\text{O}_{\text{phosphate}}$ measurements to SST. Calculations were conducted using the HadCM3 climatic simulations, the $\delta^{18}\text{O}_{\text{sw}}$-salinity relationship of Railsback et al. (1989) and the temperature-$\delta^{18}\text{O}_{\text{phosphate}}$ relationship of Pucéat et al. (2010). For each time slice, locations of StabisoDB data (Grossman and Joachimski, 2022) are shown with circle markers, the color of which represents the SST bias calculated in the closest oceanic model grid point. Emerged landmasses are outlined and shaded white. Robinson projections with latitude shown every 30°. For readability, only a subset of the 109 stage-level HadCM3 simulations is represented here (every ca. 25 Myrs), although the full series of HadCM3 simulations is used elsewhere. See Supplementary Fig. S5 for zonal averages.</figcaption></figure>

Fig. 3 illustrates the impact of estimating SST using the latitude-based $\delta^{18}\text{O}_{\text{sw}}$ correction of Zachos et al. (1994; Eqn. (4)), instead of a globally-uniform $\delta^{18}\text{O}_{\text{sw}}$ value as per Fig. 2. Although more reliable latitudinal corrections have been proposed over the years (see Gaskell et al., 2022; Grossman and Joachimski, 2022; Roberts et al., 2011), the formula of Zachos et al. (1994) is still very often used in paleoclimate studies (e.g., O'Brien et al., 2017). This approximation, derived from recent ocean values from the southern hemisphere Atlantic and Pacific oceans (Zachos et al., 1994; their p. 362), over-corrects low-latitude $\delta^{18}\text{O}_{\text{sw}}$ according to our method (Supplementary Fig. S9) and, as a result, produces a positive SST interpretation bias over most of the global ocean (Fig. 3). For this reason, this correction will not be used hereafter.

### 3.3. Validating results using geological data and sensitivity testing

Important uncertainties are associated with simulated sea-surface salinity (SSS) values, impacting our calculations of $\delta^{18}\text{O}_{\text{sw}}$, especially in coastal regions. Model SSSs near continental masses are largely dependent on the watersheds imposed as boundary condition in the GCM simulations, which are difficult to reconstruct due to large uncertainties in deep-time land topography. Riverine runoff comprises the input of freshwater, with zero salinity, to coastal ocean grid points. Overestimating runoff, for example due to the inappropriate positioning of an important river mouth in the model, will inappropriately reduce local SSS. This SSS underestimation translates into an underestimation of $\delta^{18}\text{O}_{\text{sw}}$ following the $\delta^{18}\text{O}_{\text{sw}}$-salinity relationship equations

(Supplementary Fig. S1), ultimately resulting in a more positive SST interpretation bias (Eqn. (2)). Applying this locally anomalous positive SST interpretation bias as a correction factor in a $\delta^{18}\text{O}$-temperature calculation will lead to over-correction for local $\delta^{18}\text{O}_{\text{sw}}$ and therefore an excessively cold SST interpretation (see discussion in Section 3.5).

Nevertheless, some SSS spatial patterns are expected on grounds of physical climatology. This is the case for the very saline Paleo-Tethys Ocean during the Late Permian and Early Triassic (Fig. 2k; Supplementary Fig. S6k). Numerical simulations showed that oceanic connection between the Paleo-Tethys and the Panthalassa oceans was limited by the geometry of the Paleo-Tethys basin and regional bathymetry, making it a relatively restricted oceanic basin at tropical latitudes (Wu et al., 2024). Some SSS patterns are further supported by geological evidence. For instance, Mutterlose et al. (2012) interpreted the offset between their Hauterivian to early Barremian TEX86 and $\delta^{18}\text{O}$ data as reflecting SSS values being 3 psu (practical salinity units) higher than global average in northwest Europe. This aligns with the values reaching ca. 37.3 psu (+2.5 psu compared with the model global average) in our GCM simulations (Fig. 2f; Supplementary Fig. S6f). Similarly, remnant Early Cretaceous North Atlantic seawater is characterized by high salinity values (up to twice the modern global average; Sanford et al., 2013), supporting the high model SSS values (up to 45 psu) in the North Atlantic at 149 Ma (Fig. 2g; Supplementary Fig. S6g). During the Permian-Triassic transition, interbedded fluvial and tidal deposits and important terrigenous inputs show that the equatorial northwestern American basin was characterized by important riverine fluxes (e.g., Blakey and Ranney, 2008), a feature that is also consistent
<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Effect of estimating SST using the latitude-based approximation of Zachos et al. (1994) rather than the calculated local δ¹⁸Osw. As per Fig. 2 but correcting δ¹⁸Osw for latitude using the formula of Zachos et al. (1994) instead of assuming a globally-uniform δ¹⁸Osw. The latitude-based correction of Zachos et al. (1994) is valid only within 70° latitude (higher latitudes are left blank).</figcaption></figure>

with model results (see Fig. 2k δ¹⁸O data point on the western coast of equatorial Pangea; Supplementary Fig. S6k).

Although our approach is based on a global δ¹⁸Osw-salinity relationship and therefore does not account for regional or temporal variations, several lines of evidence suggest that our results capture important features of the spatial biases in δ¹⁸O data interpretations. First, spatial patterns and numerical values for our 0 Ma, 26 Ma, and 52 Ma time slices using the latitude-based δ¹⁸Osw correction of Zachos et al. (1994) (Fig. 3) are in relatively good agreement with the results obtained by Gaskell and Hull (2023; their Fig. 1), using the same latitude-based δ¹⁸Osw correction, based on isotope-enabled GCM simulations for the modern and Last Glacial Maximum, Miocene, and Eocene, respectively. The zonally-averaged SST interpretation bias calculated in our Late Cretaceous simulations also aligns with the values simulated by Zhou et al. (2008) for the Cenomanian using an isotope-enabled GCM (Fig. 4c). We tentatively compared δ¹⁸Osw values estimated from carbonate clumped isotope values of Phanerozoic fossil material (Supplementary Text and Supplementary Table S1) with model-derived δ¹⁸Osw values (Supplementary Fig. S10), but the limited number of data points does not permit any statistical assessment. While models and data sometimes align (for instance in the Central Atlantic during the Early Cretaceous; Supplementary Fig. S10b), we also note important mismatches both between models and data and between adjacent data points. Some obvious model-data mismatches arise from the representation of regional paleogeography; this is well illustrated by the lack of Western Interior Seaway in our Late Cretaceous time slice (Supplementary Fig. S10a). This model-data comparison effort is hence of limited use to date due to the small number of data points, but will benefit from additional constraints as new data become available.

Second, Fig. 5 shows that our simulated spatial-temporal patterns in SST interpretation biases stand when alternative δ¹⁸Osw-salinity

relationships are used, including the δ¹⁸Osw-salinity relationship proposed for the Indian and Southern oceans by Tiwari et al. (2013; their Fig. 4b) and the δ¹⁸Osw-salinity relationship proposed for the modern by Gaskell et al. (2022).

Third, results obtained using different GCMs share similarities (Fig. 4). In all three models tested, the primarily latitudinal variation of SST interpretation biases simulated for the Late Ordovician (Fig. 4a,d,g) contrasts with the more complex spatial patterns obtained for the Late Cretaceous (Fig. 4b,e,h), but note that FOAM Late Ordovician results exhibit stronger longitudinal contrasts in SST interpretation bias than the other two models (Fig. 4a,d,g). All Late Cretaceous simulations show a very fresh (low-SSS) Arctic ocean associated with a strongly positive SST bias of ca. +10 °C in IPSL-CM5A2, +20 °C in FOAM and providing the closest match to the results of Zhou et al. (2008) (Fig. 4f) and +25 °C in HadCM3. The magnitude of the SST interpretation bias in the Arctic seems to be directly related to the degree of isolation of the Arctic Ocean from the Tethys and Central Atlantic oceans, which varies across models due to differences in the paleogeographical configurations interpolated at the different model grid resolutions (Fig. 4b,e,h). Specifically, although the three simulations use the paleogeographies of the PALEOMAP project (Scotese and Wright, 2018), the IPSL-CM5A2 simulation is the only one that includes three straits between the Arctic, Tethys, and Central Atlantic, with both HadCM3 and FOAM simulations only including one strait. The stronger mixing of Arctic water masses with the rest of the ocean in the IPSL-CM5A2 simulation probably contributes to the smaller contrast in salinity, hence reduced SST interpretation bias in the Arctic compared to the other two simulations. An SST interpretation bias of a few tens of degrees is the same order of magnitude of bias (> 40 °C) calculated by Gaskell and Hull (2023) for the 66.5–90° N latitudinal band in an Eocene isotope-enabled simulation published by Zhu et al. (2020).
<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Sensitivity to the global climate model. Effect of estimating SST using a globally-uniform $\delta^{18}$O$_{sw}$ rather than the calculated local $\delta^{18}$O$_{sw}$. Calculations were conducted using the $\delta^{18}$O$_{sw}$-salinity relationship of Railsback et al. (1989) and the temperature-$\delta^{18}$O$_{phosphate}$ relationship of Puceat et al. (2010), and 3 different GCMs: (a-c) HadCM3, (d-f) FOAM v1.5 and (g-i) IPSL-CM5A2, for 2 periods of time: the Late Ordovician (left column) and the Late Cretaceous (middle and right columns). Note that all simulations use paleogeographical reconstructions based on or derived from the paleogeographical reconstructions of Scotese and Wright (2018). For the Late Cretaceous (Ordovician), the following simulations are used: 97 Ma – 1426 ppm pCO$_2$ (449 Ma – 782 ppm pCO$_2$) for HadCM3, 90 Ma – 1120 ppm (445 Ma – 2240 ppm) for IPSL-CM5A2, 95 Ma – 1120 ppm (445 Ma – 1960 ppm) for FOAM. Differences in SST interpretation bias arising for each time slice from the use of different pCO$_2$ values are expected to be second-order compared to differences arising from the use of different models (see sensitivity test to pCO$_2$ discussed in Section 3.3). Panels (c), (f) and (i) show the simulated SST (blue line with envelope; left y-axis), the SST reconstructed assuming a globally-uniform $\delta^{18}$O$_{sw}$ (orange line with envelope, left y-axis) and the resulting SST interpretation bias (grey line with envelope, right y-axis). Lines represent zonal means and envelopes extend from the 10th to the 90th percentiles. For comparison, the black line represents the Late Cretaceous SST bias calculated by Zhou et al. (2008) using an isotope-enabled model (right y-axis; to be compared with the gray line). IPSL-CM5A2 fields of temperature and salinity were interpolated from the native model tripolar irregular grid to a 1° by 1° regular grid (using a nearest-neighbor algorithm) before conducting calculations.</figcaption></figure>

Fourth, additional Late Ordovician and Late Cretaceous FOAM simulations run under a wide range of alternative atmospheric pCO₂ levels confirm that palaeogeography has a much greater influence on calculated SST biases than does global climate state (Fig. 6). These results imply that the main patterns and values in SST interpretation bias calculated using the 109 HadCM3 simulations during the Phanerozoic (Fig. 2) would be little impacted by uncertainties in atmospheric CO₂ concentrations.

### 3.4. Temporal trends in SST interpretation bias and implications for proxy data interpretation

In our simulations (Fig. 7a), the time slices and regions that are prone to the strongest SST interpretation biases are the southern middle to high latitudes from the Devonian to the Permian, middle to high northern latitudes (45° N to 90° N) from the Triassic to present-day, and to a lesser extent the southern high latitudes during the Cretaceous and Neogene. The strongest biases correspond to SST overestimations in regions of low salinity at higher latitudes.

Fig. 7c shows that mean SSTs are consistently overestimated in the 45°-60° latitude band in each hemisphere; this latitudinal window represents the high-latitude edge of the region for which most $\delta^{18}$O proxy data are available (Fig. 7a; see also Judd et al., 2022). There are

two clear intervals of substantial biases in these latitudinal bands: the southern hemisphere during the Devonian to Permian, and the northern hemisphere during the Triassic to present-day, with zonally-averaged biases of the order of +5 °C and +2.5 °C, respectively. The strong biases in these two regions arise from the continental arrangement that isolates high-latitude oceanic basins and leads to the development of low-salinity water masses.

To estimate the impact of SST interpretation bias on proxy data interpretations, we evaluated the geographical distribution of the $\delta^{18}$O data available in the StabisoDB database (Grossman and Joachimski, 2022) in Fig. 2 and Fig. 7a. We only used quality-controlled $\delta^{18}$O data of the StabisoDB database, which have been selected from calcite and aragonite $\delta^{18}$O values of samples with metadata documenting preservation of primary isotopic signatures, and removed all incomplete StabisoDB entries, notably those missing geographic coordinates. In our selected StabisoDB dataset, no data are available in the region of strongest bias located at the southern middle to high latitudes during the Devonian and Carboniferous (Fig. 7a). Despite the strong local deviations of $\delta^{18}$O$_{sw}$ in enclosed oceanic basins of the southern middle to high latitudes during this period of time, these local deviations are not expected to impact the interpretation of currently available proxy data (Fig. 2n–p). Some Permian $\delta^{18}$O data points, in contrast, are available at polar latitudes (Fig. 7a). Although they do not come from the enclosed
<figure><img src="image_7.png" /><figcaption>**Fig. 5. Sensitivity to the δ¹⁸Osw-salinity relationship.** Zonally-averaged effect of estimating SST using a globally-uniform δ¹⁸Osw rather than the calculated local δ¹⁸Osw. Calculations were conducted using the HadCM3 climatic simulations, the temperature-δ¹⁸Ophosphate relationship of Puceat et al. (2010) and (a) the δ¹⁸Osw-salinity relationship of Railsback et al. (1989), as per Fig. 2; (b) the δ¹⁸Osw-salinity relationship proposed for the Indian and Southern oceans by Tiwari et al. (2013); and (c) the δ¹⁸Osw-salinity relationship proposed for the modern by Gaskell et al. (2022). C: Cambrian; O: Ordovician; S: Silurian; D: Devonian; C: Carboniferous; P: Permian; T: Triassic; J: Jurassic; K: Cretaceous; P: Paleogene; N: Neogene.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 6. Sensitivity to the global climatic state.** Effect of estimating SST using a globally-uniform δ¹⁸Osw rather than the calculated local δ¹⁸Osw. Calculations were conducted using the δ¹⁸Osw-salinity relationship of Railsback et al. (1989) and the temperature-δ¹⁸Ophosphate relationship of Pucéat et al. (2010), and the FOAM model run at various atmospheric pCO₂ levels (see atmospheric CO₂ concentrations expressed in ppm in the legend) for two time slices: (a) the Late Ordovician and (b) the Late Cretaceous. Lines represent zonal means and envelopes extend from the 10th to the 90th percentiles at each model latitude.</figcaption></figure>

basins on the western side of Pangea (e.g., South America), but from the Neo-Tethys and Panthalassa (e.g., from Australia), they are subject to strong interpretation biases (Fig. 2k–m). Proxy data from the Late Jurassic onwards are also susceptible to SST interpretation biases in the northern middle to high latitudes (Figs. 2a-i, 7a).

SST interpretation biases for specific proxy data locations show that SST values above ca. 45° latitude are systematically overestimated, while the sign of the SST biases is variable at lower latitudes (Figs. 2, 7b; Supplementary Fig. S5). At low latitudes, SSTs are generally

underestimated except near the coast at river mouths. The impact of continental runoff is very visible during the Silurian, Carboniferous, and Permian (Fig. 7b), when δ¹⁸O data mostly come from shallow seas, which in our GCM simulations are under the influence of strong riverine freshwater input (Fig. 2m–o,r).

The middle latitudes represent a transitional domain between the low latitudes, characterized by an overall SST underestimation, and the higher latitudes, characterized by a strong positive SST interpretation bias that increases rapidly with increasing latitude. Geological periods
<figure><img src="image_9.png" /><figcaption>**Fig. 7.** Temporal trends in simulated SST interpretation bias. (a) Zonally-averaged effect of estimating SST using a globally-uniform $\delta^{18}O_{sw}$ rather than the calculated local $\delta^{18}O_{sw}$ (as per Fig. 5a), with paleo-latitudes of StabisoDB data (Grossman and Joachimski, 2022) superimposed as circle markers. (b) As per panel (a) but at StabisoDB locations only. (c) Effect of estimating SST using a globally-uniform $\delta^{18}O_{sw}$ rather than the calculated local $\delta^{18}O_{sw}$ averaged over the northern (orange line) and southern (blue line) mid- to polar latitudes (45 to 60 °N). All calculations were conducted using the HadCM3 climatic simulations, the $\delta^{18}O_{sw}$-salinity relationship of Railsback et al. (1989) and the temperature-$\delta^{18}O_{phosphate}$ relationship of Pucéat et al. (2010). C: Cambrian; O: Ordovician; S: Silurian; D: Devonian; C: Carboniferous; P: Permian; T: Triassic; J: Jurassic; K: Cretaceous; P: Paleogene; N: Neogene.</figcaption></figure>

for which proxy data are available in the 30°–60° latitudinal band may thus be characterized by strong distortions of their reconstructed latitudinal SST gradients. Specifically, an SST underestimation at the low latitudes, combined with an overestimation at the higher latitudes (> 45 °), may artificially flatten latitudinal SST gradients. As previously noted by Zhou et al. (2008), this result is interesting since it may contribute to reconciling proxy data with numerical climate models that have been routinely considered to overestimate the steepness of latitudinal SST gradients during greenhouse periods. Based on StabisoDB data, we estimate that this phenomenon may impact our understanding of latitudinal SST gradients during the Jurassic and Cretaceous (Figs. 2, 7a,b).

### 3.5. Re-examining the phanerozoic low-latitude SST curve using local $\delta^{18}O_{sw}$

In an attempt to quantify the impact of correcting proxy data for local $\delta^{18}O_{sw}$ on long-term SST reconstructions, we revisit in Fig. 8 the recent Phanerozoic tropical SST reconstruction of Grossman and Joachimski (2022). The latter is based on a subset of quality-controlled $\delta^{18}O_{carbonate}$ and $\delta^{18}O_{phosphate}$ data covering the last 500 million years. The black line in Fig. 8 shows the SSTs calculated as per the original publication: global $\delta^{18}O_{sw}$ was varied through time to account for changing land-ice volumes (Supplementary Fig. S2) and local $\delta^{18}O_{sw}$ was corrected for latitude using two different equations representative of either greenhouse or icehouse periods (Grossman and Joachimski, 2022, their eqn. (4) and 5; Roberts et al., 2011). Small differences from the originally published curve (Grossman and Joachimski, 2022; their Fig. 4E) are due to the use of a different smoothing algorithm. In addition, data points that were

not correctly rotated back to their paleo-positions were ignored. The blue line shows the Phanerozoic tropical SST record calculated using the same methods except that, instead of taking a latitude-based $\delta^{18}O_{sw}$ correction as in the original paper, we extracted for each StabisoDB entry the $\delta^{18}O_{sw}$ value simulated at the closest oceanic grid points in the closest of the 109 HadCM3 simulations (see Fig. 7a,b). The offset between the black and blue lines therefore represents the impact of salinity-corrected local $\delta^{18}O_{sw}$ versus a latitude-based $\delta^{18}O_{sw}$ correction derived from modern oceans or the recent geological past. Here we follow the original publication by not detrending for long-term $\delta^{18}O$ changes, which is a broader topic beyond the scope of this study.

Although the long-term temperature trends calculated with and without accounting for local $\delta^{18}O_{sw}$ are mostly similar, we find important deviations on shorter time scales (< 10 million years) during specific time intervals. Taken at face value, our results suggest that the well-known early Paleozoic cooling (Trotter et al., 2008), instead of reaching an acme during the latest Ordovician, may have continued for another 20 million years and culminated during the latest Silurian (but see discussion below). Our results also suggest that the warming documented during the latest Carboniferous, around 300 million years ago, may instead be a sampling artefact with proxy data reflecting low-latitude cooling (see discussion below). Finally, our revised estimates suggest that the background climate during the middle Permian was cooler than previously reported, making the late Permian warming even stronger relative to the earlier part of the Permian. The scarcity of $\delta^{18}O$ data in younger time slices in our StabisoDB dataset, which arises in part from the loss of data points during the calculation of their paleo-positions, does not permit interpreting temporal trends from the Triassic
<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Model-informed Phanerozoic SST reconstruction. (a) Low-latitude (–30:30 °N) SST reconstruction derived from the quality-controlled δ¹⁸Ocarbonate (cross markers) and δ¹⁸Ophosphate (square markers) data of the StabisoDB database as per the original calculations of Grossman and Joachimski (2022) (black markers and line) vs. using for each data point the δ¹⁸Osw value averaged over the closest oceanic grid point and immediately adjacent grid points in the closest of the 109 Phanerozoic HadCM3 time slices considered (blue markers and line). (b) Difference induced by using the δ¹⁸Osw value averaged over the closest oceanic grid point and immediately adjacent grid points in the closest of the 109 Phanerozoic time slices considered, instead of the original calculations of Grossman and Joachimski (2022) (i.e., blue minus black lines of panel (a)). Calculations were conducted using the HadCM3 climatic simulations, the δ¹⁸Osw-salinity relationship of Railsback et al. (1989), the temperature-δ¹⁸Ophosphate relationship of Pucéat et al. (2010) and the temperature-δ¹⁸Ocarbonate relationship of Kim and O'Neil (1997). C: Cambrian; O: Ordovician; S: Silurian; D: Devonian; C: Carboniferous; P: Permian; T: Triassic; J: Jurassic; K: Cretaceous; P: Paleogene; N: Neogene.</figcaption></figure>

onwards with confidence.

The GCM and δ¹⁸Osw-salinity relationships are factors that may strongly impact reconstructed temporal trends. Interestingly, however, very similar results are obtained when using the δ¹⁸Osw-salinity relationship of Gaskell et al. (2022) (Supplementary Fig. S11). The lack of alternative Phanerozoic-scale series of GCM simulations providing access to sea-surface salinity prevents us from quantifying the model-dependence of reconstructed long-term SST trends. Although large-scale spatial patterns in SST interpretation bias share important similarities amongst three models tested on two different time slices, we also note regional deviations (Fig. 4). The middle latitudes constitute a critical transition between the contrasting hydrological regimes of the low and high latitudes, characterized by a steep gradient in SST interpretation bias. This may render the interpretation of SSTs during periods for which oxygen isotope data are found at these latitudes – Jurassic and younger (Fig. 7) – particularly vulnerable to model choice.

Although the δ¹⁸Osw-salinity-SST calculation approach can be reasonably applied through most of the Phanerozoic, the late Carboniferous is a notable aberration during which this method produces strongly negative SSTs (Fig. 8). This is impossible to reconcile with available constraints on the range of variations in Phanerozoic low-latitude SSTs reconstructed using independent approaches based on different proxies (Finnegan et al., 2011; Isson and Rauzi, 2024; O'Brien et al., 2017; Scotese et al., 2021). Clumped isotope temperatures for fossils from the aragonite-bearing Boggy Formation in Oklahoma (305 Ma) yield an average temperature of 24.9 ± 1.1 °C (mean ± 1 SE; Came et al., 2007), similar to clumped isotope temperatures for the older fossils from the Russian Platform (27.0 ± 3.0 °C; Henkes et al., 2018). Other evidence is the presence of widespread late Carboniferous tropical carbonate platforms with abundant warm-water faunas (Kiessling et al., 2003).

In their recent climatic reconstruction obtained by statistically combining proxy data with a large ensemble of HadCM3 simulations, Judd et al. (2024) proposed very small deviations from the SST values of Grossman and Joachimski (2022) during the Carboniferous (their supplementary fig. S20). This contrasts with the very strong, unrealistic cooling reconstructed in this study (Fig. 8). The difference arises, at least in part, from the fact that Judd et al. (2024) used an alternative approach to calculate past δ¹⁸Osw based on paleoclimate simulations. Instead of using a single δ¹⁸Osw-salinity relationship, they calculated δ¹⁸Osw in ancient times by deriving δ¹⁸Osw from analog oceanographic settings found in the present day, based on sea-surface salinity, sea-surface temperature, precipitation, latitude and distance to the coast. This approach, which is largely anchored in the present day and may seem less mechanistic than the one proposed by Railsback et al. (1989), provides more reasonable results in the Carboniferous. This probably reflects the better skill of the approach employed by Judd et al. (2024) in capturing δ¹⁸Osw values in specific oceanographic contexts like low-latitude epicontinental seas and restricted basins where salinity values simulated in climate models and δ¹⁸Osw-salinity relationships are poorly reliable (see discussion below).

The underestimation of low-latitude SSTs during the late Carboniferous reflects in part inaccurate sea-surface salinity values simulated in the GCM. Specifically, sea-surface salinity drops below 20 psu in the restricted equatorial basin on the western side of Pangea at 301 Ma, where most StabisoDB data points are found, as a result of continental runoff in the intertropical convergence zone and poor connection of the basin with the global ocean at the model resolution (Fig. 2m; Supplementary Fig. S6m). Seawater salinity values below 20 psu are unlikely considering that δ¹⁸O was measured on conodonts and brachiopods that are considered to be intolerant of large changes in salinity. In the modern ocean, such seawater salinity values in open-ocean settings are
only found near high-latitude river mouths in the Arctic (Reagan et al., 2024).

The underestimation of low-latitude SSTs also arises from limitations in the δ¹⁸O~sw~·salinity relationship in low-salinity, low-latitude settings. In contrast with the middle to high latitudes where precipitation and hence continental runoff waters are strongly ¹⁸O-depleted (with δ¹⁸O values as low as -20 ‰), water that runs off to the ocean in low-latitude settings is generally characterized by higher (sometimes even positive) δ¹⁸O values (Nan et al., 2019). The higher δ¹⁸O values are due to the combined contributions of less ¹⁸O-depleted precipitation and downstream evaporative enrichment in watersheds extending over arid areas. This induces a local decoupling of salinity and δ¹⁸O~sw~. For instance, surface waters of the present-day Nile river are characterized by mean δ¹⁸O values of +0.8 ‰ and +3.9 ‰ during the wet and dry seasons, respectively (Cockerton et al., 2013). Today, low δ¹⁸O~sw~ values at the low latitudes arise in the specific case of ¹⁸O-depleted runoff (e.g., < -8 ‰ VSMOW) in regions with high-altitude runoff such as the Pacific coast of Panama (Tao et al., 2013).

Similarly, sea-level curves and clumped isotopes show that land-ice volumes reached a maximum during the latest Ordovician (Hirnantian) and decreased afterwards (Finnegan et al., 2011). Although a very similar climatic trend was reconstructed during the Silurian by Judd et al. (2024), featuring a double-peak cooling centered on 30 °C tropical SST (their supplementary fig. S20), aforementioned geological evidence suggests that the drop in Silurian SSTs below Late Ordovician values calculated in our revised SST estimates is unlikely. The same sources of error as for the latest Carboniferous probably apply in the case of the Silurian. We further note that, while the geological record of Silurian glaciations is not sufficient to quantify temporal changes in land-ice volumes (Loydell, 2007), the persistence of substantial land ice during the Silurian, by increasing global-ocean δ¹⁸O~sw~ above the land-ice-free value of -1.08 ‰ used in Fig. 8 (see Supplementary Fig. S2), may be enough to raise Silurian reconstructed SSTs by a few degrees (Pucéat et al., 2010; Eqn. (1)) to Late Ordovician levels.

Fig. 8 demonstrates the overall limitation of our approach in reconstructing δ¹⁸O~sw~ in low-latitude shallow shelves, due to difficulties in accurately simulating sea-surface salinity in these settings in GCMs combined with the potential decoupling of salinity and δ¹⁸O~sw~ that limits the applicability of global δ¹⁸O~sw~·salinity relationships. While these biases apply to all our time slices, their expression becomes prominent during the latest Carboniferous when most StabisoDB data come from such low-latitude, epicontinental settings under the influence of continental runoff. Our δ¹⁸O~sw~·salinity relationship may also inaccurately capture the contribution of upwelling systems to shallow-water δ¹⁸O~sw~ on the western coast of low-latitude continents, while previous isotope-enabled simulations proved that this mechanism significantly influences the δ¹⁸O~sw~ composition in epicontinental seas (Macarewich et al., 2021).

Our results illustrate well the influence of continental runoff on many δ¹⁸O data available in the geological record, which largely come from epicontinental settings especially in pre-Mesozoic time intervals. Previous work demonstrated that today SSTs in these environments are warmer and more seasonal than open-ocean values from the same latitudes, and deviate from the zonal average due to oceanic gyre circulation (Judd et al., 2020). Our analysis demonstrates that, while regional δ¹⁸O~sw~ values largely impact paleo-SST reconstructions, δ¹⁸O~sw~ in epicontinental settings is extremely challenging to constrain, inviting further caution when calculating SSTs based on these coastal data sources. Results also suggest that the large spatial variability in δ¹⁸O~sw~ in epicontinental settings may be an important source of scatter in SST reconstructions deriving from oxygen isotope measurements conducted during the same period of time but at different paleogeographical locations (see Supplementary Fig. S7).

Future work with clumped isotopes (Letulle et al., 2022) and isotope-enabled GCMs may help provide better estimates of δ¹⁸O~sw~ for the calculation of SSTs, although models would also face several of the

issues raised above, notably the difficulty to represent coastal environments and the impact of downstream evaporative ¹⁸O enrichment.

## 4. Concluding remarks

We use the well-documented co-variation of δ¹⁸O~sw~ and ocean salinity, together with GCM simulations, to calculate changes in δ¹⁸O~sw~ and their implications for δ¹⁸O-derived SST reconstructions over the last 541 million years. We demonstrate a first-order impact of continental configuration on δ¹⁸O~sw~ and hence potential biases in SST reconstructions during the Phanerozoic, with the global climatic state having a much smaller impact on this bias. Although our calculations are based on a number of assumptions, large-scale spatial patterns are instructive. Our results for several Cenozoic time intervals and the Late Cretaceous, in particular, are supported by isotope-enabled GCM experiments, and our Late Cretaceous results are relatively robust to alternative climate models and δ¹⁸O~sw~·salinity relationship equations.

Our study identifies the oceanic regions and time intervals that are prone to the strongest δ¹⁸O~sw~·induced SST reconstruction biases, which should constitute a priority for future work using isotope-enabled GCMs. Among the main features of our results, we note that latitudinal SST gradients since the Jurassic may have been significantly underestimated, especially in the northern hemisphere (Fig. 2). Biases in SST gradients could have strong regional imprints, such as over Europe during the Early Cretaceous (Fig. 2f) when SSTs may have been underestimated at 30° N and overestimated at 45–60° N, artificially flattening the northern hemisphere reconstructed SST gradient. We also note that regional SSTs may have been significantly underestimated in the Paleotethys during the Permian and in the European Tethys during the Triassic (Fig. 2i,k). These regional SST interpretation biases, combined with the clustering of available SST proxy data in these same regions, may have led to an overall underestimation of temperatures in the geological past. Warmer background temperatures may have constituted a further challenge for marine organisms facing global climate warming during the Permian-Triassic and end-Triassic mass extinctions (Joachimski et al., 2012; Penn et al., 2018; Sun et al., 2012).

However, our approach fails to accurately calculate δ¹⁸O~sw~ values in low-latitude epicontinental seas under the influence of strong continental runoff. Global δ¹⁸O~sw~·salinity relationships do not permit capturing the complex mechanisms controlling δ¹⁸O~sw~ in these specific hydrographic contexts. In line with previous work that quantified some biases associated with inferring SSTs based on epicontinental-sea data (Judd et al., 2020), our results further invite to remain cautious when dealing with SST proxy data coming from low-latitude shallow-water environments that are strongly influenced by continental runoff.

This type of approach to reconstruct large-scale δ¹⁸O~sw~ spatial patterns during the Phanerozoic based on GCMs (see also Judd et al., 2024) probably represents one of the best options until geochemical analyses permitting to disentangle changes in temperature and δ¹⁸O~sw~ can be readily produced (Finnegan et al., 2011; Henkes et al., 2018; Isson and Rauzi, 2024; Letulle et al., 2022), or increased computational power makes isotope-enabled simulations more widely available.

## Data availability

HadCM3 simulation data are available online (instructions to download: https://www.paleo.bristol.ac.uk/ummodel/scripts/papers/Using_BRIDGE_webpages.pdf). The Late Cretaceous IPSL simulation is available from Laugié et al. (2021). All other climatic fields used in this paper, together with other files required to reproduce the calculations and figures, are hosted on Zenodo (https://zenodo.org/records/13836450). All files are also available on GitHub (https://github.com/alexpohl/d18O_panero); this is the recommended mode of access as it will contain any updates and clarifications. δ¹⁸O~sw~ values calculated using HadCM3 simulations for all Phanerozoic time slices, stored in the form of NetCDF files, are also hosted on Zenodo (https:
//zenodo.org/records/15240528).

### CRediT authorship contribution statement

**Alexandre POHL:** Writing – review & editing, Writing – original draft, Visualization, Validation, Software, Resources, Project administration, Methodology, Investigation, Funding acquisition, Formal analysis, Data curation, Conceptualization. **Thomas W. WONG HEARING:** Writing – review & editing, Visualization, Validation, Resources, Methodology, Conceptualization. **Arnaud BRAYARD:** Writing – review & editing, Validation, Resources. **Ethan GROSSMAN:** Writing – review & editing, Validation, Resources, Investigation, Data curation. **Michael M. JOACHIMSKI:** Writing – review & editing, Validation, Resources, Methodology, Data curation. **Guillaume LE HIR:** Writing – review & editing, Resources. **Thomas LETULLE:** Writing – review & editing, Validation, Resources, Methodology, Formal analysis, Data curation. **Daniel J. LUNT:** Writing – review & editing, Validation, Resources. **Mathieu MARTINEZ:** Writing – review & editing, Validation, Resources. **Emmanuelle PUCEAT:** Writing – review & editing, Validation, Resources. **Guillaume SUAN:** Writing – review & editing, Validation, Resources, Methodology. **Paul VALDES:** Writing – review & editing, Validation, Resources. **Yannick DONNADIEU:** Writing – review & editing, Validation, Resources, Methodology.

### Declaration of competing interest

The authors have no competing interests to declare.

### Acknowledgements

This is a contribution to UNESCO project IGCP 735 "Rocks and the Rise of Ordovician Life" (Rocks n' ROL). AP acknowledges the support of the French Agence Nationale de la Recherche (ANR) under references ANR-22-CE01-0003 (project ECO-BOOST) and ANR-23-CE01-0003 (project CYCLO-SED). AP and MM acknowledge the support of the programme Tellus of the Institut National des Sciences de l'Univers, CNRS (project WAOW). AP and TWW thank the Leverhulme Trust for supporting the project 'Earth System dynamics at the dawn of the animal-rich biosphere' (RPG-2022-233). CM5A2 simulations used in this work have been performed on the HPC resources of the CEA-TGCC using computing hours provided by GENCI under allocations A0130102212 and A0150102212. Calculations were performed using HPC resources from DNUM CCUB (Centre de Calcul de l'Université de Bourgogne). The authors thank Daniel E. Gaskell and an anonymous reviewer for constructive comments that significantly improved the manuscript, and Tristan J. Horner for editorial handling.

### Supplementary materials

Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.epsl.2025.119418.

### References

Bemis, B.E., Spero, H.J., Bijma, J., Lea, D.W., 1998. Reevaluation of the oxygen isotopic composition of planktonic foraminifera: experimental results and revised paleotemperature equations. Paleooceanography 13, 150–160. https://doi.org/10.1029/98PA00070.

Bergmann, K.D., Finnegan, S., Creel, R., Eiler, J.M., Hughes, N.C., Popov, L.E., Fischer, W.W., 2018. A paired apatite and calcite clumped isotope thermometry approach to estimating Cambro-Ordovician seawater temperatures and isotopic composition. Geochim. Cosmochim. Acta 224, 18–41. https://doi.org/10.1016/j.gca.2017.11.015.

Blakey, R.C., Ranney, W., 2008. Ancient Landscapes of the Colorado Plateau. Grand Canyon Association. ed.

Brand, U., Logan, A., Bitner, M.A., Griesshaber, E., Azmy, K., Buhl, D., 2011. What is the ideal proxy of Palaeozoic seawater chemistry?.

Came, R.E., Eiler, J.M., Veizer, J., Azmy, K., Brand, U., Weidman, C.R., 2007. Coupling of surface temperatures and atmospheric CO2 concentrations during the Palaeozoic era. Nature 449, 198–202.

Cather, S.M., Dunbar, N.W., McDowell, F.W., McIntosh, W.C., Scholle, P.A., 2009. Climate forcing by iron fertilization from repeated ignimbrite eruptions: the icehouse-silicic large igneous province (SLIP) hypothesis. Geosphere 5, 315–324.

Cockerton, H.E., Street-Perrott, F.A., Leng, M.J., Barker, P.A., Horstwood, M.S.A., Pashley, V., 2013. Stable-isotope (H, O, and Si) evidence for seasonal variations in hydrology and Si cycling from modern waters in the Nile Basin: implications for interpreting the Quaternary record. Quat. Sci. Rev. 66, 4–21. https://doi.org/10.1016/j.quascirev.2012.12.005.

Daëron, M., Gray, W.R., 2023. Revisiting Oxygen-18 and Clumped Isotopes in Planktic and Benthic Foraminifera. Paleooceanog and Paleoclimatol 38, e2023PA004660. https://doi.org/10.1029/2023PA004660.

Davies, A.J., Brand, U., Tagliavento, M., Bitner, M.A., Bajnai, D., Staudigel, P., Bernecker, M., Fiebig, J., 2023. Isotopic disequilibrium in brachiopods disentangled with dual clumped isotope thermometry. Geochim. Cosmochim. Acta 359, 135–147. https://doi.org/10.1016/j.gca.2023.08.005.

Finnegan, S., Bergmann, K., Eiler, J.M., Jones, D.S., Fike, D.A., Eisenman, I., Hughes, N. C., Tripati, A.K., Fischer, W.W., 2011. The magnitude and duration of late Ordovician-early Silurian glaciation. Science (1979) 331, 903–906. https://doi.org/10.1126/science.1200803.

Galili, N., Shemesh, A., Yam, R., Brailovsky, I., Sela-Adler, M., Schuster, E.M., Collom, C., Bekker, A., Planavsky, N., Macdonald, F.A., Préat, A., Rudmin, M., Trela, W., Stresssson, U., Heikoop, J.M., Aurell, M., Ramajo, J., Halevy, I., 2019. The geologic history of seawater oxygen isotopes from marine iron oxides. Science (1979) 365, 469–473. https://doi.org/10.1126/science.aaw9247.

Gaskell, D.E., Huber, M., O'Brien, C.L., Inglis, G.N., Acosta, R.P., Poulsen, C.J., Hull, P. M., 2022. The latitudinal temperature gradient and its climate dependence as inferred from foraminiferal δ18O over the past 95 million years. Proc. Natl. Acad. Sci. U.S.A. 119, e2111332119. https://doi.org/10.1073/pnas.2111332119.

Gaskell, D.E., Hull, P.M., 2023. Technical note: a new online tool for δ18O–temperature conversions. Clim. Past 19, 1265–1274. https://doi.org/10.5194/cp-19-1265-2023.

Grossman, E.L., 2012. Applying Oxygen Isotope Paleothermometry in Deep Time. The Paleontological Society Papers 18, 39–68.

Grossman, E.L., Joachimski, M.M., 2022. Ocean temperatures through the Phanerozoic reassessed. Sci. Rep. 1–14. https://doi.org/10.1038/s41598-022-11493-1.

Henkes, G.A., Passey, B.H., Grossman, E.L., Shenton, B.J., Yancey, T.E., Pérez-Huerta, A., 2018. Temperature evolution and the oxygen isotope composition of Phanerozoic oceans from carbonate clumped isotope thermometry. Earth Planet. Sci. Lett. 490, 40–50. https://doi.org/10.1016/j.epsl.2018.02.001.

Huyghe, D., Daëron, M., de Rafelis, M., Blamart, D., Sébilo, M., Paulet, Y.-M., Lartaud, F., 2022. Clumped isotopes in modern marine bivalves. Geochim. Cosmochim. Acta. https://doi.org/10.1016/j.gca.2021.09.019.

Isson, T., Rauzi, S., 2024. Oxygen isotope ensemble reveals Earth's seawater, temperature, and carbon cycle history. Science (1979) 383, 666–670. https://doi.org/10.1126/science.adg1366.

Joachimski, M.M., Lai, X., Shen, S., Jiang, H., Luo, G., Chen, B., Chen, J., Sun, Y., 2012. Climate warming in the latest Permian and the Permian-Triassic mass extinction. Geology. 40, 195–198. https://doi.org/10.1130/G32707.1.

Jones, L.A., Eichenseer, K., 2021. Uneven spatial sampling distorts reconstructions of Phanerozoic seawater temperature. Geology. https://doi.org/10.1130/G49132.1.

Judd, E.J., Bhattacharya, T., Ivany, L.C., 2020. A Dynamical Framework for Interpreting Ancient Sea Surface Temperatures. Geophys. Res. Lett. 47. https://doi.org/10.1029/2020GL089044.

Judd, E.J., Tierney, J.E., Huber, B.T., Wing, S.L., Lunt, D.J., Ford, H.L., Inglis, G.N., McClymont, E.L., O'Brien, C.L., Rattanasriampaipong, R., Si, W., Staitis, M.L., Thirumalai, K., Anagnostou, E., Cramwinckel, M.J., Dawson, R.R., Evans, D., Gray, W.R., Grossman, E.L., Henahan, M.J., Hupp, B.N., MacLeod, K.G., O'Connor, L. K., Sánchez Montes, M.L., Song, H., Zhang, Y.G., 2022. The PhanSST global database of Phanerozoic sea surface temperature proxy data. Sci. Data 9, 753. https://doi.org/10.1038/s41597-022-01826-0.

Judd, E.J., Tierney, J.E., Lunt, D.J., Montañez, I.P., Huber, B.T., Wing, S.L., Valdes, P.J., 2024. A 485-million-year history of Earth's surface temperature. Science (1979) 385, eadk3705. https://doi.org/10.1126/science.adk3705.

Keller, C.B., Husson, J.M., Mitchell, R.N., Bottke, W.F., Gernon, T.M., Boehnke, P., Bell, E.A., Swanson-Hysell, N.L., Peters, S.E., 2019. Neoproterozoic glacial origin of the Great Unconformity. Proc. Natl. Acad. Sci. U.S.A. 116, 1136–1145. https://doi.org/10.1073/pnas.1804350116.

Kiessling, W., Flügel, E., Golonka, J., 2003. Patterns of phanerozoic carbonate platform sedimentation. Lethaia. https://doi.org/10.1080/00241160310004648.

Kim, S.-T., O'Neil, J.R., 1997. Equilibrium and nonequilibrium oxygen isotope effects in synthetic carbonates. Geochim. Cosmochim. Acta 61, 3461–3475. https://doi.org/10.1016/S0016-7037(97)00169-5.

Laugié, M., Donnadieu, Y., Ladant, J., Bopp, L., Ethé, C., Raisson, F., 2021. Exploring the impact of Cenomanian paleogeography and marine gateways on oceanic oxygen. Paleooceanogr. Paleoclimatol. 36, e2020PA004202. https://doi.org/10.1029/2020PA004202.

LeGrande, A.N., Schmidt, G.A., 2006. Global gridded data set of the oxygen isotopic composition in seawater. Geophys. Res. Lett. 33, 2006GL026011. https://doi.org/10.1029/2006GL026011.

Letulle, T., Gaspard, D., Daëron, M., Arnaud-Godet, F., Vinçon-Laugier, A., Suan, G., Lécuyer, C., 2023. Multi-proxy assessment of brachiopod shell calcite as a potential archive of seawater temperature and oxygen isotope composition. Biogeosciences. 20, 1381–1403. https://doi.org/10.5194/bg-20-1381-2023.
Letulle, T., Suan, G., Daëron, M., Rogov, M., Lécuyer, C., Vinçon-Laugier, A., Reynard, B., Montagnac, G., Lutikov, O., Schlögl, J., 2022. Clumped isotope evidence for Early Jurassic extreme polar warmth and high climate sensitivity. Clim. Past 18, 435–448. https://doi.org/10.5194/cp-18-435-2022.

Loydell, D.K., 2007. Early Silurian positive δ13C excursions and their relationship to glaciations, sea-level changes and extinction events. Geological Journal 42, 531–546.

Macarewich, S.I., Poulsen, C.J., Montañez, I.P., 2021. Simulation of oxygen isotopes and circulation in a late Carboniferous epicontinental sea with implications for proxy records. Earth Planet. Sci. Lett. 559, 116770. https://doi.org/10.1016/j.epsl.2021.116770.

Mutterlose, J., Malkoc, M., Schouten, S., Sinninghe Damsté, J.S., 2012. Reconstruction of vertical temperature gradients in past oceans — Proxy data from the Hauterivian–early Barremian (Early Cretaceous) of the Boreal Realm. Palaeogeogr. Palaeoclimatol. Palaeococ. 363–364, 135–143. https://doi.org/10.1016/j.palaeo.2012.09.006.

Nan, Y., Tian, F., Hu, H., Wang, L., Zhao, S., 2019. Stable Isotope Composition of River Waters across the World. Water. (Basel) 11, 1760. https://doi.org/10.3390/w11091760.

O'Brien, C.L., Robinson, S.A., Pancost, R.D., Sinninghe Damsté, J.S., Schouten, S., Lunt, D.J., Alsenz, H., Bornemann, A., Bottini, C., Brassell, S.C., Farnsworth, A., Forster, A., Huber, B.T., Inglis, G.N., Jenkyns, H.C., Linnert, C., Littler, K., Markwick, P., McAnena, A., Mutterlose, J., Naafs, B.D.A., Püttmann, W., Sluijs, A., van Helmond, N.A.G.M., Vellekoop, J., Wagner, T., Wrobel, N.E., 2017. Cretaceous sea-surface temperature evolution: constraints from TEX 86 and planktonic foraminiferal oxygen isotopes. Earth. Sci. Rev. 172, 224–247. https://doi.org/10.1016/j.earscirev.2017.07.012.

Ontiveros, D.E., Beaugrand, Grégory, Lefebvre, Bertrand, Marcilly, Chloé M, Servais, Thomas, Pohl, Alexandre, 2023. Impact of global climate cooling on Ordovician marine biodiversity. Nat. Commun. 14, 6098. https://doi.org/10.1038/s41467-023-41685-w.

Penn, J.L., Deutsch, C., Payne, J.L., Sperling, E.A., 2018. Temperature-dependent hypoxia explains biogeography and severity of end-Permian marine mass extinction. Science (1979) 362, eaat1327.

Pohl, A., Donnadieu, Y., Le Hir, G., Buoncristiani, J.F., Vennin, E., 2014. Effect of the Ordovician paleogeography on the (in)stability of the climate. Climate of the Past 10, 2053–2066.

Pucéat, E., Joachimski, M.M., Bouilloux, A., Monna, F., Bonin, A., Motreuil, S., Morinière, P., Hénard, S., Mourin, J., Dera, G., Quesne, D., 2010. Revised phosphate-water fractionation equation reassessing paleotemperatures derived from biogenic apatite. Earth Planet. Sci. Lett. 298, 135–142.

Railsback, L.B., Anderson, T.F., Ackerly, S.C., Cisne, J.L., 1989. Paleogeographic modeling of temperature-salinity profiles from stable isotopic data. Paleogeography 4, 585–591. https://doi.org/10.1029/PA004i005p00585.

Reagan, J.R., Boyer, T.P., Garcia, H.E., Locarnini, R.A., Baranova, O.K., Bouchart, C., Cross, S.L., Mishonov, A.V., Paver, C.R., Seidov, D., Wang, Z., Dukhovskoy, D., 2024. World Ocean Atlas 2023.

Roberts, C.D., LeGrande, A.N., Tripati, A.K., 2011. Sensitivity of seawater oxygen isotopes to climatic and tectonic boundary conditions in an early Paleogene

simulation with GISS ModelE-R. Paleogeography 26. https://doi.org/10.1029/2010PA002025, 2010PA002025.

Sanford, W.E., Doughten, M.W., Coplen, T.B., Hunt, A.G., Bullen, T.D., 2013. Evidence for high salinity of Early Cretaceous sea water from the Chesapeake Bay crater. Nature 503, 252–256. https://doi.org/10.1038/nature12714.

Scotese, C.R., Song, H., Mills, B.J.W., van der Meer, D.G., 2021. Phanerozoic paleotemperatures: the earth's changing climate during the last 540 million years. Earth. Sci. Rev., 103503 https://doi.org/10.1016/j.earscirev.2021.103503.

Scotese, C.R., Wright, N., 2018. PALEOMAP Paleodigital Elevation Models (PaleoDEMS) for the Phanerozoic (PALEOMAP Project, 2018) [WWW Document]. URL https://www.earthbyte.org/paleodem-resource-scotese-and-wright-2018/.

Shields, G., Veizer, J., 2002. Precambrian marine carbonate isotope database: version 1.1. Geochem. Geophys. Geosyst. 3. https://doi.org/10.1029/2001GC000266.

Sun, Y., Joachimski, M.M., Wignall, P.B., Yan, C., Chen, Y., Jiang, H., Wang, L., Lai, X., 2012. Lethally hot temperatures during the early triassic greenhouse. Science (1979) 338, 366–370. https://doi.org/10.1126/science.1224126.

Tao, K., Robbins, J.A., Grossman, E.L., O'Dea, A., 2013. Quantifying Upwelling and Freshening in Nearshore Tropical American Environments Using Stable Isotopes in Modern Gastropods. BMS 89, 815–835. https://doi.org/10.5343/bms.2012.1065.

Tiwari, M., Nagoji, S.S., Kartik, T., Drishya, G., Parvathy, R.K., Rajan, S., 2013. Oxygen isotope-salinity relationships of discrete oceanic regions from India to Antarctica vis-à-vis surface hydrological processes. Journal of Marine Systems 113–114, 88–93. https://doi.org/10.1016/j.jmarsys.2013.01.001.

Trotter, J.A., Williams, I.S., Barnes, C.R., Lécuyer, C., Nicoll, R.S., 2008. Did cooling oceans trigger Ordovician biodiversification? Evidence from conodont thermometry. Science (1979) 321, 550–554. https://doi.org/10.1126/science.1155814.

Valdes, P., Scotese, C., Lunt, D., 2021. Deep Ocean Temperatures through Time. Climate of the Past 17, 1483–1506. https://doi.org/10.5194/cp-17-1483-2021.

Veizer, J., Prokoph, A., 2015. Temperatures and oxygen isotopic composition of Phanerozoic oceans. Earth. Sci. Rev. 146, 92–104. https://doi.org/10.1016/j.earscirev.2015.03.008.

Wheeler, J.R., Smith, M.P., Boomer, I., 2012. Oxygen isotope variability in conodonts: implications for reconstructing Palaeozoic palaeoclimates and palaeoceanography. JGS 169, 239–250. https://doi.org/10.1144/0016-76492011-048.

Wu, Y., Pohl, A., Tian, L., Corso, J.D., Cui, Y., Chu, D., Tong, J., Song, Huyue, Song, Hanchen, Song, Haijun, 2024. Recurrent marine anoxia in the Paleo-Tethys linked to constriction of seaways during the Early Triassic. Earth Planet. Sci. Lett. 643, 118882. https://doi.org/10.1016/j.epsl.2024.118882.

Zachos, J.C., Stott, L.D., Lohmann, K.C., 1994. Evolution of Early Cenozoic marine temperatures. Paleogeography, 9, 353–387. https://doi.org/10.1029/93PA03266.

Zhou, J., Poulsen, C.J., Pollard, D., White, T.S., 2008. Simulation of modern and middle Cretaceous marine δ18O with an ocean-atmosphere general circulation model. Paleogeography 23. https://doi.org/10.1029/2008PA001596.

Zhu, J., Poulsen, C.J., Otto-Bliesner, B.L., Liu, Z., Brady, E.C., Noone, D.C., 2020. Simulation of early Eocene water isotopes using an Earth system model and its implication for past climate reconstruction. Earth Planet. Sci. Lett. 537, 116164. https://doi.org/10.1016/j.epsl.2020.116164.