Reconstructing the magnitude of Early Toarcian (Jurassic) warming using the reordered clumped isotope compositions of belemnites

The magnitude of temperature changes in the Early Jurassic are not well known. Clumped isotope measurements can potentially be used to provide better constrains, but unfortunately many of the well-studied sedimentary successions that preserve Lower Jurassic fossils experienced burial temperatures above the limits of preservation of D47, which for geological timescales is thought to be between 80–120 C. Samples from these basins are expected to be partially reordered and yield apparent clumped isotope temperatures that are warmer than original values. Here, we explore whether useful paleoclimate information can be recovered from these samples. We test the hypothesis that relative temperature differences are preserved in partially reordered samples when they experience a common burial history. This was done with the use of reordering models and D47 measurements of early Jurassic belemnites from the Aubach section of the SW German Basin, a basin that has a relatively well constrained burial history with maximum burial temperatures above 90 C. We find that even though partial reordering progressively erases the D47 difference between samples, the majority (>50%) of the signal is preserved when samples are buried at temperatures as high of 150 C for up to 200 Ma. Moreover, the models demonstrate that – regardless of burial conditions – partially reordered samples always preserve minimum records of temperature change across climate events. These inferences are supported by the belemnite D47 data that show partially reordered compositions and warming/cooling patterns across the Early Jurassic that closely mimic observations from independent proxies. Model observations are used to interpret a 13 ± 4 C (95% ci) temperature increase that is observed in the belemnite data across the Early Toarcian. The large magnitude of the temperature excursion is explained as a combination of warming and a change in belemnite habitat before and after the Toarcian Ocean Anoxic Event. Our results demonstrate the usefulness of partially reordered samples and further open the use of this proxy in deep time settings. 2020 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/ licenses/by/4.0/).


INTRODUCTION
The Early Jurassic ($201-174 Ma ago) was an interval characterized by dynamic climate when large changes in ocean temperature have been inferred from oxygen isotope https://doi.org/10.1016/j.gca.2020.10.005 0016-7037/Ó 2020 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). ratio measurements of well-preserved calcitic macrofossils (e.g. Dera et al., 2011). Towards the end of the Early Jurassic, temperatures varied from the Late Pliensbachian Cool Event (Korte and Hesselbo, 2011), to the warm Toarcian Oceanic Anoxic event (TOAE) (Jenkyns, 1988;Bailey et al., 2003;Rosales et al., 2004), and to the Early Aalenian cooling Event . One crucial uncertainty that affects inferred paleoenvironmental changes is the relative contributions of temperature and seawater oxygenisotope composition on the d 18 O values of fossil calcite. For instance, when interpreted at face value, fossil d 18 O values suggest seawater temperature changes of up to 20-30°C across some of the major environmental perturbations in the Early Jurassic (e.g., Dera et al., 2011;Korte et al., 2015). Such large temperature changes are hard to reconcile with expectations for mid-latitude sites without invoking large-scale reorganizations of heat transporting ocean currents, or alternatively with large changes in the d 18 O of ocean waters. The result is that the actual magnitudes of temperature change are currently not well constrained.
A potentially useful method to constrain temperatures is the carbonate clumped isotope paleothermometer. Clumped isotopes (expressed with the parameter D 47 ) exploit the preference for 13 C and 18 O to bond with each other in the carbonate mineral lattice, an effect that becomes more prominent with decreasing temperatures (Wang et al., 2004;Ghosh et al., 2006;Schauble et al., 2006). Because D 47 is concerned with isotope exchange reactions within a single mineral phase, it does not depend on the isotopic composition of the fluids from which carbonates form Eiler, 2011). This is an important distinction from the traditional d 18 O paleothermometer, and it is what makes this paleothermometer an attractive method for paleoclimate reconstructions, especially in deep time.
One obstacle for the application of clumped isotopes to the Early Jurassic, and to deep time in general, is the poor preservation that is often observed in carbonate fossils. Carbonates can be altered during burial through a variety of reactions that partially modify or even completely erase the original paleoenvironmental information. Some of these reactions involve the transfer of mass with the surrounding environment and affect the chemical, structural, and isotopic characteristics of carbonates. Because the presence of these types of alterations -such as dissolution and recrystallization -are troublesome, there are a number of mineralogical and geochemical tools that have been developed to screen samples for signs of diagenetic modifications (e.g., Podlaha et al., 1998;Pérez-Huerta et al., 2007;. Some reactions, however, involve no mass transfer and do not affect trace element compositions and d 18 O values (Henkes et al., 2014;Stolper and Eiler, 2015) but they can completely or partially reset clumped isotope compositions. These reactions are called bond reordering reactions and proceed in the solid state when isotopes are exchanged between neighboring carbonate groups as bonds between carbon and oxygen atoms are broken and reformed. Partial or complete reordering occurs when carbonates are exposed to burial temperatures that are different from original conditions and when temperatures are high enough to break the kinetic energy barrier for oxygen and carbon exchange and diffusion within the carbonate mineral lattice, which for calcite and for geological timescales (tens of millions of years) is thought to be between 80 and 120°C (Henkes et al., 2014;Stolper and Eiler, 2015). It is thus possible to encounter samples with pristine textures, trace element concentrations, and stable isotope signatures but with altered D 47 values.
Other than working only in sedimentary successions that were never very deeply buried, there is currently no unambiguous way to detect reordered D 47 compositions. Unfortunately, many of the well-studied sections that preserve Lower Jurassic sediments do not meet this criterion. For instance, sediments in the Cleveland Basin (Yorkshire, England) likely experienced burial depths of more than 4 km (Kemp et al., 2005), and thus temperatures which are above the limits of preservation for D 47 . Similarly, clumped isotope measurements of well-preserved brachiopods from Peniche (the Global Stratotype Section and Point for the Toarcian) show that they are partially reordered (Fantasia et al., 2019). This does not mean that clumped isotopes from these basins cannot provide useful paleoclimate information. On the contrary, we expect relative temperature differences to be preserved when samples experience a common burial history because reordering should shift original D 47 values to higher (or warmer) values by roughly the same amounts (Cummins et al., 2014;Ryb and Eiler, 2018). This hypothesis, however, has not yet been explicitly tested with field measurements of partially reordered fossil calcite, nor have the reordering models been examined in detail to understand what percentage of the signal is lost.
This study has two objectives: the first is to quantify the magnitude of Early Toarcian temperature change in the southwestern Laurasian Seaway (Tethyan epicontinental sea), and the second is to demonstrate the utility of partially reordered clumped isotope compositions to tackle questions in deep time paleoclimate research. We quantify the warming with D 47 measurements of well-preserved belemnites collected from a section in the SW German Basin. Our study site has a relatively well constrained burial history and experienced burial temperatures above the limits of solid-state reordering (Mazurek et al., 2006;Looser et al., 2019). The belemnites show temperatures consistent with reordered clumped isotope compositions and have D 47 values that vary closely with independent proxies of warming and cooling climate across the Early Jurassic. Our clumped isotope data together with the results of reordering models (Passey and Henkes, 2012;Stolper and Eiler, 2015) demonstrate the usefulness of partially reordered samples and further open the use of this proxy in deep time settings.

STUDY SITE
The Aubach section (47. 84437°N, 8.48268°E) is located in the valley of the Aubach stream near the village of Aselfingen in the Wutach area of Southwestern Germany. The section contains shallow marine sediments that were deposited in the SW German Basin of the European Epi-continental Seaway (Laurasian Seaway) during the Early Jurassic (Fig. 1C). The sediments spanning the Pliensbachian to Toarcian stages comprise the Numismalismergel, Amaltheenton and Posidonienschiefer Formations. The stratigraphy of the section is well constrained with ammonite, ostracods, foraminifera and dinoflagellate biostratigraphy (Hahn, 1971;Riegraf et al., 1984;Riegraf, 1985;Schlatter, 1991;Feist-Burkhardt, 2010), carbon isotope stratigraphy (Hougaard et al., in press), and with the presence of several regional marker beds (Riegraf et al., 1984;Rö hl et al., 2001). These beds can be used to correlate the Aubach section to other locations in Southwestern Germany ( Fig. 1; Rö hl et al., 2001;Suan et al., 2015), which includes the well-studied section of Dotternhausen (Fig. 1B), and Northern Switzerland (Kuhn and Etter, 1994;Montero-Serrano et al., 2015;Fantasia et al., 2018b).
The Early Toarcian deposits in the SW German Basin are condensed relative to sections where early Toarcian sediments are preserved; although, a pattern of expansion towards deeper parts of the basin (Dotternhausen) and con-densation towards shallower marginal sites (Swiss Jura) has been documented (Montero-Serrano et al., 2015;Fantasia et al., 2018b). Even though the Early Toarcian is condensed at Aubach, sediments dated to all of the Early Toarcian ammonite zones and subzones identified by Riegraf et al. (1984) can be recognized (Fig. 1A). Thus, this section can be used to constrain the main paleoenvironmental changes that occurred across this time interval.
The regional burial and thermal history was first investigated by Mazurek et al. (2006) who obtained apatite fission track, vitrinite reflectance and maturation-dependent organic biomarker (biomarker isomerisation) data for the Benken borehole, located approx. 25 km southeast of the Aubach section. For Lower Jurassic sediments, they derived a multi-stage burial history including an initial interval of approx. 28 Ma of deep burial during the Cretaceous with temperatures of about 85°C, an interim cooling due to a major erosional event during Paleocene and Eocene, and a second interval of burial with temperatures up to approx. 65°C during the Miocene. The first interval of deep burial was the result of deposition of marine . Correlation between the Pliensbachian-Toarcian successions at Aubach (A) and Dotternhausen (B). Carbon isotope data at Dotternhausen from Rö hl et al. (2001). Lithostratigraphy is from Rö hl et al. (2001) and the abbreviations are: OS= ''Oberer Stein", US= ''Unterer Stein", SG= ''Seegrasschiefer", TF= ''Tafelfleins", SB= ''Spinatum-Bank". Subdivision of the tenuicostatum zone at Dotternhausen after Maisch (2018). (C) Map of Early Jurassic paleogeography and location of the Aubach site, map modified from Korte et al. (2015). Zero meters marks the location of the Pliensbachian-Toarcian boundary in both sections. zoic sediments while the second burial was caused by Alpine Molasse deposits.
More recently Looser et al. (2019) re-examined the thermal history of Northern Switzerland based on combined D 47 data and U-Pb (LA-ICP-MS) ages of various diagenetic carbonates from Late Triassic and Early Jurassic sediments from a site at Frick (Northern Switzerland), which is located approx. 50 km south-west of the Aubach section. Their data confirm the burial history proposed by Mazurek et al. (2006) and also provide a better constraint of 90 ± 3°C for the temperature of burial during the Cretaceous. They also find a population of dolomite veins dated to the late Jurassic and Early Cretaceous with D 47 temperatures of 115 ± 5°C, which raise the possibility of an additional interval of 5-30 Ma of burial at higher temperatures. However, it is not clear if the host rock reached those temperatures because the veins could be the result of hydrothermalism with deep-fluid migration rather than burial. Since both scenarios are possible, we incorporate that uncertainty in our geochemical modeling by considering the effect of two burial scenarios: 1) a scenario that considers only one 90 ± 3°C episode of deep burial (in the following, cold burial) and 2) a scenario where both the 90 ± 3°C and 115 ± 5°C burial events occurred (in the following, warm burial) (Electronic Annex Fig. S1).
There are several reasons to warrant the assumption that burial conditions at Aubach were similar to Frick and Benken. For instance, the thicknesses of the preserved Jurassic formations on top of the Toarcian do not vary greatly between the sites (Naef, 2008), suggesting similar depths at least until the late Jurassic. The sites are positioned at comparable distances to the Alpine frontal thrust, suggesting similar burial since the Oligocene. Additionally, since they are located in the undeformed Tabular Jura, there were likely no variations in burial due to folding and thrusting related to the Jura fold-and-thrust belt. Finally, a clumped isotope temperature of 97 ± 10°C from a calcite cement from the Sinemurian ''Arietenkalk" limestones at Aubach provides direct evidence of burial at similar temperatures as the other two other sites (Electronic Annex dataset).

Carbon, oxygen and clumped isotope measurements
We measured the clumped isotope composition of several well-preserved belemnite rostra collected across Late Pliensbachian and Early Toarcian aged sediments. Belemnites were visually examined and screened with minor and trace element data to identify signs of post-depositional recrystallisation. Fresh surfaces of the belemnites were broken and inspected under a binocular microscope. Samples were obtained with a handheld drill only from samples that showed pristine radial structures. Approximately 3-4 mg of material were obtained from each belemnite while carefully avoiding any diagenetic cements, the rims and the apical zones (e.g., McArthur et al., 2000;Li et al., 2013;. These powders were used for both clumped isotope measurements and for measurements of trace element ratios (Mn/Ca, Sr/Ca, Mg/Ca, and Fe/Ca). Trace elements were measured at the University of Exeter (see Ullmann et al., 2020 for detailed methods) and were used as an additional test for potential alteration (e.g., . Clumped isotope measurements were obtained only for samples that were judged to be well preserved with the exception of three diagenetic cements that were sampled from belemnites collected from the falciferum biozone. These were measured to estimate the geochemical signatures of diagenetic alteration products at the site. Clumped isotope and d 13 C and d 18 O measurements were made at the Geological Institute at ETH Zü rich using methods that are described in Meckler et al. (2014) and Bernasconi et al. (2018). Clumped isotope temperatures were estimated with the revised Kele et al. (2015) calibration, which was recalculated with the IUPAC parameters by Bernasconi et al. (2018). Twenty-three belemnite samples were measured on average 11 times each, and uncertainties in D 47 are reported as confidence intervals with margins of error at the 95% confidence level (e.g., Bonifacie et al., 2017;Fernandez et al., 2017). Measurements were obtained over a period of 11 months and have an average external reproducibility, as determined from the four ETH standards and a Carrara Marble, of 0.031 for D 47 , 0.04 for d 13 C, and 0.09 for d 18 O (±1 standard deviation, N = 3777, see Electronic Annex Dataset). Potential isobaric interferences on m/z 47 were monitored according to Davies and John (2017). Oxygen and carbon isotope ratios are reported in the conventional delta notation relative to VPDB, clumped isotope values are presented relative to the carbon dioxide equilibration scale (CDES).
We used the calcite-water oxygen isotope fractionation equation of Daëron et al. (2019) to calculate belemnite d 18 O-based paleotemperatures and assumed a constant seawater d 18 O VSMOW value of À1‰. The same equation and the D 47 temperatures were used to estimate the apparent oxygen isotopic composition of the seawater in which the belemnites lived. We use that equation for belemnitederived estimates because recent work has shown it produces more reasonable seawater d 18 O VSMOW values than alternative formulations (Vickers et al., 2020). On the other hand, the equation of Kim and O'Neil (1997) was used to calculate the d 18 O values of the fluids coexisting with the cements. This formulation produces values that are indistinguishable from the equation of O'Neil et al. (1969), which has been demonstrated to produce accurate d 18 O fluid isotopic compositions of diagenetic cements (Dassié et al., 2018).

Reordering models
We performed several model experiments to determine the extent to which temperature differences between samples are preserved after partial reordering. In the experiments, we used the two types of models that have been developed to describe the changes that occur in D 47 values as a result of burial and exhumation. We considered both models because they are based on different physical reordering mechanisms, and it is not currently known which model more accurately describes changes under geological scenarios. These models are referred to as the transient defect/ equilibrium defect model (Passey and Henkes, 2012;Henkes et al., 2014) (in the following the PH-model) and the exchange/diffusion model (Stolper and Eiler, 2015) (in the following the SE-Model).
To predict changes in D 47 values the models are parameterized with Arrhenian reordering rate constants derived from laboratory experiments where carbonates are heated for varying amounts of time. One important feature observed in these experiments is that carbonates follow non-first order reordering kinetics (e.g., Henkes, 2012: Stolper andEiler, 2015). For instance, at the start of experiments D 47 values change rapidly followed by slower changes at constant rates. Both models are built to account for these observations, but they differ in the mechanisms invoked to explain them. In the PH-model reordering occurs due to the presence of different types of crystallographic defects (e.g., point defects, dislocations, cation substitutions) that provide fast pathways for isotope exchange reactions. The initially fast reordering rates observed in experiments are interpreted as the result of diffusion through a transient defect pool that -with the onset of heating -is progressively annealed away. Once these defects are annihilated, the mobility of C and O are reduced. Subsequent reordering occurs through unannealable defects, referred to as equilibrium defects. These are responsible for the constant reordering rates and the first order behavior observed in latter parts of experiments. The SE-model, on the other hand, considers that reordering occurs through a combination of isotope exchange reactions and diffusion. In this formulation, initial D 47 changes occur as clumped species are quickly broken via isotope exchange between neighboring carbonate groups. This process results in the formation of what are termed pairs in the model, which are defined as two adjacent carbonate groups with either a 13 C or 18 O isotope in them. Because pairs can back react to form clumped species, their presence effectively buffers against further changes in clumped isotopes compositions until the next reordering mechanism takes over. Subsequent reordering occurs when temperatures are higher as pairs are broken up and singletons ( 13 C or 18 O bearing groups) diffuse through the crystal lattice.
In our model experiments, we used the first-order approximation of the PH-model that only considers the reordering that is aided by the equilibrium defect pool (Passey and Henkes, 2012) even though a model that also takes into account the contribution of both defect components is available. The simpler formulation was used because the approximation produces indistinguishable predictions for geological scenarios (Henkes et al., 2014). Changes in D 47 were calculated using the Arrhenian kinetic parameters derived from heating experiments of an optical calcite specimen by Passey and Henkes (2012) using a numerical solution to the integrated first-order approximation model (Equation 13 of Passey and Henkes, 2012). For the experiments with the SE-model, changes in D 47 were modeled using the Python code provided in the supplemental information of Lloyd et al. (2018). The code uses the Arrhenian kinetic parameters for an optical calcite speci-men (Stolper and Eiler, 2015) to solve two differential equations that describe changes in pair concentration and reaction progress with time (Equation 12 and 13 from Lloyd et al., 2018). The solutions to these equations are needed to calculate how clumped isotope compositions evolve.
Our model experiments consisted on two sets of experiments. In the first set, we varied the burial times and temperatures of two calcite samples: one with the initial D 47 value of a sample that formed at 20°C and the other with the value of a sample formed at 30°C. These two calcites were instantaneously subjected to temperatures between 80 and 180°C for 1-200 million years and then immediately exhumed (box-car heating scenarios). Results are reported as apparent temperatures, which are the temperatures that correspond to the reordered D 47 values. Results are also reported as relative temperature differences between the two samples that we considered. These relative temperature differences were calculated using the predicted D 47 values for both samples at the end of the same burial scenarios and the initial temperature difference between the two samples. They are reported in percentage form. For instance, a relative temperature difference of 100% means that the full 10°C difference between the two samples was preserved after partial reordering. In the second set of experiments, we kept the burial temperature constant and varied burial times and the initial D 47 values of samples. These experiments were done to demonstrate that results from first set of experiments can be generalized to samples with other initial clumped isotope compositions. The burial histories consisted of box-car heating scenarios with a burial temperature of 140°C for the PH-model and of 115°C for the SE-model. Different peak temperatures were used because of the different reordering susceptibilities exhibited by the two models (e.g., Stolper and Eiler (2015) and Section 4.4 below).

Trace metal measurements
Belemnite minor and trace element ratios (Mn/Ca, Sr/ Ca, and Fe/Ca) are presented in Fig. 2, where we also plot the compositions of three fracture filling cements that were sampled from the Aubach belemnites. The Figure also includes estimates for upper (Fe/Ca, Mn/Ca) and lower (Sr/Ca) for belemnites that have been judged to be well preserved. Only one sample has Mn/Ca ratios that fall outside of the range. The other samples have element/Ca ratios that fall within the limits of what has been considered excellent preservation, and the vast majority plot within the most conservative limits that have been described in the literature. The cements, on the other hand, plot outside these limits and show much more elevated Fe/Ca and Mn/Ca, as well as low Sr/Ca ratios.

Stable isotopes and clumped isotopes from the Aubach section
The belemnites have d 18 O VPDB values that range between À2.3 and À0.4‰. These values correspond to d 18 O based paleotemperatures of 19 and 28°C, assuming a constant seawater d 18 O value of À1‰ (VSMOW) and the calcite-paleothermometry equation of Daëron et al. (2019). Temperatures approximately 8°C lower are obtained if the more commonly used equation of Kim and O'Neil (1997) is used instead. The most positive d 18 O values (coldest temperatures) occur at the end of the Pliensbachian, and the most negative values (warmest temperatures) occur after the main body of T-OAE (Fig. 3). These data are consistent with a cooling trend from the Early to the Late Pliensbachian, and with warming across the Pl-To boundary. Additionally, our belemnite d 18 O data suggest that the warmest conditions of the Early Toarcian occurred as part of the T-OAE; the warmest temperatures are recorded in the samples measured in the falciferum zone. Belemnites have carbon isotope values ranging from À0.63 to 3.62‰. The lightest d 13 C values were measured in the spinatum zone and the most positive values in the falciferum zone (Fig. 3).
Clumped isotope temperatures range from 30 ± 8 to 54 ± 12°C (95% ci) and are warmer than d 18 O derived paleotemperatures. However, both datasets follow the same trends with the coldest D 47 temperatures at the end of Pliensbachian, and the warmest D 47 temperatures in the samples measured in the falciferum zone. These trends can be seen in Fig. 3 and in Table 1, where we present the average belemnite D 47 temperatures in the different ammonite biozones. We calculated average biozone values because the analytical errors in individual belemnites estimates are relatively high (±7°C 95% ci), and average values -with their much smaller error bars -allow us to examine the main temperature changes across these time periods (Table 1). These values show a long-term cooling trend over the Pliensbachian of 5 ± 4°C (95% ci) and an abrupt temperature increase of 6 ± 4°C (95% ci) across the Pl-To stage boundary. The abrupt warming was followed by an additional warming of 8 ± 4°C (95% ci) in response to the T-OAE, which corresponds to a total temperature increase of 13 ± 4°C (95% ci) for the Early Toarcian (falciferum biozone) relative to the Late Pliensbachian estimate (spinatum biozone).
We used the belemnite d 18 O values and their respective D 47 temperatures to calculate the apparent oxygen isotope composition of Early Jurassic seawater (d 18 O VSMOW ). Seawater estimates were also averaged across ammonite Fig. 2. Element/Ca ratios from belemnite rostra and three fracture filling cements. Element ratios plot within the range of what has been described as pristine compositions. The cements clearly plot outside of these fields. Numbers show the most and least conservative upper (Fe/ Ca and Mn/Ca) and lower (Sr/Ca) for pristine belemnites, as summarized by . Number symbols are: 1) Sr/Ca lower limit from Korte and Hesselbo (2011), 2) Sr/Ca lower limit from Ullmann et al. (2014), 3) Mn/Ca upper limit from Voigt et al. (2003), 4) Mn/ Ca upper limit from Benito and Reolid (2012), 5) Fe/Ca upper limit from Voigt et al. (2003), 6) Fe/Ca upper limit from Jones et al. (1994). biozones because individual estimates have a relatively high uncertainty (average is ±1.3 ‰, 95% ci; Fig. 3). We find an average apparent water d 18 O VSMOW value of 2.4 ± 0.5 ‰ (95% ci) for the Early Jurassic, and some variability between biozones. The variability between biozones, however, is within our analytical uncertainty; a one-way ANOVA shows that the biozones have means that are indistinguishable from each other at the 95% confidence level. The same is true for the individual water d 18 O estimates, which plot within uncertainty of the average water value for the Early Jurassic (Fig. 3).
Clumped, oxygen and carbon isotope data for the three belemnite cements are presented in Table 2. Carbon isotope values are within the range of the belemnites whereas oxy-  temperature of the cements is 51 ± 6°C (95% ci), which is warmer than any of the average biozone temperatures. Similarly, the isotopic compositions of fluids in isotopic equilibrium with the cements are also very different from our estimates of Early Jurassic apparent seawater d 18 O VSMOW values. These values are 4 ± 2 (95% ci) more negative than the fluids in equilibrium with the belemnites.

Relationships between geochemical variables
We observe moderate to strong significant correlations (p < 0.05) between Mg/Ca and Sr/Ca ratios, oxygen isotope compositions and clumped isotope temperatures (Fig. 4) Fig. S2).

Reordering models
Results from the first set of model experiments -where we varied the burial times and temperatures of two samples -demonstrate that apparent temperatures increase with higher burial temperatures and, generally, with longer times spent at those temperatures ( Fig. 5A and C). The models, however, do not agree on the temperature at which reordering begins. The SE-model predicts that measurable amounts of reordering occur within 200 million years at temperatures as low as 80°C. In the first-order approximation, on the other hand, similar amounts of reordering do not appear until 120°C. Second, the rates of reordering are very different in the two models. For example, in the SE-model apparent temperatures initially increase very fast (within the first 10 Ma) between 100 and 140°C of burial temperature. Subsequent changes, however, are very subtle and essentially stop. In other words, in the SE-model D 47 values reach a state of partial equilibration and result in the distinctive plateau that is observed in Fig. 5A between 100 and 140°C. This phenomenon occurs because the two reordering mechanisms that operate in the model (exchange and diffusion) procced at different temperatures. At low temperatures the exchange mechanism rapidly destroys clumped species, but higher temperatures are needed for the pairs to be destroyed (i.e., for singly substituted species diffuse away from each other) and for full equilibration to occur (Stolper and Eiler, 2015;Lloyd et al., 2018). Fig. 5B and D show how the temperature difference between the 20°C and 30°C sample changes with increasing amounts of reordering. The figures demonstrate that, in general, longer times and hotter burial temperatures tend to erase D 47 differences between samples. However, both models also show that even though the absolute temperatures are lost, evidence of warming/cooling signals can be preserved. For instance, more than 50% of a signal is preserved after burial for 200 million years at 150°C. This is true for both models. Additionally, the results in these figures closely mirror the results previously presented for the reordering behavior of the 20°C sample. Differences between samples are erased at lower temperatures in the SE-model than in the PH-model. The prominent plateau that we identified in Fig. 6A in apparent temperatures is also present in the relative temperature differences. In this case, the plateau shows that, after an initial partial loss that occurs relatively fast, D 47 differences between samples remain regardless of the time samples spent buried at temperatures of up to 140°C.
Results from the second set of model experiments are presented in Fig. 6. The figure shows that there is a linear relationship between initial and apparent D 47 temperatures when samples are subjected to the same burial history (Fig. 6). Moreover, the slopes of these lines, which are estimates of the fractional amount of D 47 signal that remains after partial reordering (fraction signal remaining), progressively shallow as burial time increases. For example, slopes are close to 1 during the initial stages of burial and gradually change to 0.63 for the SE-model (Fig. 6a) and 0.55 for the PH-model data after 180 Ma of burial (Fig. 6b). Thus, after burial (for 180 Ma) two samples that had an initial temperature difference of 100°C between them would show apparent temperature differences of 63°C and 55°C for the SE-model and the PH-model, respectively (i.e., the product of 100°C and 0.63 equals 63°C).

DISCUSSION
We argue that solid state reordering of 13 C-18 O bonds is the most likely explanation for elevated belemnite D 47 temperatures but that despite partial reordering -and the loss Table 2 Stable isotope, clumped isotope from three fracture filling cements in Toarcian belemnites (falciferum biozone), and the apparent d 18 O value of the fluids in equilibrium with the cements. Errors are margins of error at the 95% ci, and n is the number of individual replicate analyses. what is currently understood to be the lethal temperature limit for many marine animals (Brock, 1967;Nguyen et al., 2011;Sunday et al., 2012), and would likely imply impossibly hot seawater temperatures in more tropical latitudes. Belemnite specific D 47 vital effects are unlikely to explain elevated temperatures given previously published belemnite clumped isotope temperatures of 10-28°C (Price and Passey, 2013;Wierzbowski et al., 2018;Vickers et al., 2019Vickers et al., , 2020. In addition, several observations suggest that samples preserve primary trace element concentrations, and that recrystallization reactions during burial cannot explain trends in D 47 temperatures. These observations  are discussed in detail below, where we also discuss the results of our modeling experiments and the paleoclimate implications of our results.

The potential role of recrystallization during burial
We interpret high clumped isotope temperatures as the result of solid state reordering reactions; however, carbonates with high D 47 temperatures can also be produced by dissolution and reprecipitation reactions (Came et al., 2007;Finnegan et al., 2011;Cummins et al., 2014;Bergmann et al., 2018;Stolper et al., 2018). If these reactions occurred in the sediments at larger depths, they could explain the elevated temperatures that we observe. However, several aspects of our dataset are not consistent with recrystallization and point instead to the partial loss of the D 47 signal through the rearrangement of 13 C-18 O bonds.
The minor and trace metal ratios and d 18 O values of our belemnites are inconsistent with the widespread incorporation of diagenetic material. Diagenetic recrystallization generally results in enrichments of Fe and Mn and in depletions of Sr (e.g., Ullmann and Korte, 2015 and references therein). The belemnites, however, have low Fe and Mn compositions that are not correlated with D 47 temperatures (Fig. S2). Moreover, they have Sr/Ca ratios that are similar to what has been observed in specimens considered well preserved (Fig. 2). These patterns are not observed in the three cements that we measured. Instead, cements follow the trends expected for diagenetic carbonates that have higher Fe and Mn concentrations, much lower Sr/Ca ratios and d 18 O values, and higher D 47 temperatures. These observations suggest that -for the style of diagenesis exemplified by the cements -trace elements, d 18 O values, and D 47 temperatures are adequate tools to detect alteration.
Sr/Ca and Mg/Ca are strongly correlated with calcite d 18 O and significantly correlated with D 47 temperatures (Fig. 4). These relationships are also not consistent with the occurrence of diagenetic alteration in the belemnites. Correlations between these elements and d 18 O have been observed in other Early Jurassic belemnites of the Tethyan realm (McArthur et al., 2000;Bailey et al., 2003;Rosales et al., 2004), and have been interpreted as the result of a temperature control in trace element incorporation during the belemnite lifetime. It is tempting to invoke the same explanation at Aubach given that these elements are also correlated with D 47 temperatures. Recent studies, however, have cast doubt into the role of temperature in Sr and Mg incorporation into belemnite calcite (Li et al., 2013;. An alternative explanation is that the relationships are the result of a genus or speciesspecific control in trace element enrichments. This is a possible explanation for our observations given that our belemnites are not differentiated into genus level. Regardless of the reason behind the enrichments, there is a consistent pattern across different sites in the Laurasian Seaway (Yorkshire U.K., southwestern Germany, and northern Spain) with low Sr/Ca and Mg/Ca ratios co-occurring with high d 18 O values during cold periods (i.e., Late Pliensbachian) and vice versa during warm periods (i.e., T-OAE; McArthur et al., 2000;Bailey et al., 2003;Rosales et al., 2004). We argue that these patterns point to primary elemental compositions in the belemnites because otherwise, they are unlikely to be observed across different basins with distinct burial and diagenetic histories.
There is a strong positive relationship between belemnite D 47 temperatures and reconstructed d 18 O waters (Fig. 4F). This pattern has been used in previous studies as a key diagnostic indicator of diagenesis under a closed system with low water to rock ratios (Huntington and Budd, 2011;Cummins et al., 2014;Huntington and Lechler, 2015;Ryb and Eiler, 2018;Veillard et al., 2019). Under this diagenetic regime, recrystallization during burial proceeds at elevated temperatures with fluids that are progressively enriched in 18 O while diagenetic products show little to no change in d 18 O values (Banner and Hanson, 1990;Huntington and Lechler, 2015). For a variety of reasons, however, we do not believe that this mechanism can explain our results. For instance, Henkes et al. (2018) recently showed that these patterns are characteristic of many clumped isotope datasets, and that correlations can be induced even in measurements of replicate standard materials. These trends occur spuriously because errors in D 47 temperatures are nearly perfectly correlated with errors in d 18 O fluid compositions. Thus, relationships have to be interpreted with caution, especially when temperature differences between samples are not much larger than analytical error, and when individual d 18 O fluid estimates are not statistically different from each other, as is the case in the Aubach belemnites (Fig. 3). Moreover, we also observe a negative relationship between D 47 temperatures and belemnite d 18 O values that is not consistent with closed-system recrystallization. The slope of this line is about what we would expect to find if calcite formed at equilibrium with water of a relatively constant isotopic value (0.15 ± 0.05 vs. 0.2°C À1 , Figs. 3 and 4E), and it implies that variability is controlled by open-system reactions, as would occur between primary calcite and seawater.
Even though we find no evidence for large scale secondary recrystallization and believe the data is best explained by reordering, it is possible that small amounts of diagenetic material were included in our measurements. A similar situation was recently encountered by Cummins et al. (2014) who used electron backscatter diffraction (EBSD) to document small amounts (5%) of recrystallized material in a brachiopod with well-preserved fabrics, oxygen and carbon isotope compositions, and minor and trace element concentrations. In the case of belemnites, two ways in which recrystallization might have occurred are: 1) after fluid migration through the apical line and rostrum fractures, and 2) as part of the rostrum wide occlusion of primary porosity originally present in the belemnites. The first of these scenarios very likely happened, but material from fractures and the apical line were largely avoided during sampling, although we cannot completely exclude that some material was incorporated into our measurements. Whether the second scenario occurred is more uncertain because it relates to a still debated aspect of the belemnite rostrum structure. The possibility was raised by recent work that demonstrated the existence of a second ultrastructure in belemnite rostra, which is distinct from the primary tetra-hedral shaped structural units that point away from the apical line (Benito et al., 2016;Hoffmann et al., 2016). In the view of some authors, this second structure is filled during burial with diagenetic cements (Benito et al., 2016;Hoffmann et al., 2016). However, as pointed out by Vickers et al. (2020), that interpretation is not consistent with isotopic and geochemical data that strongly suggests that the second phase also formed during the animal's life.
Given that small amounts of diagenetic material (5-10%) cannot be completely discounted, we estimate their possible impact by mass balance using a diagenetic endmember reflective of the worst-case scenario and assuming linear mixing in D 47 values. Non-linear mixing effects are known to occur when d 13 C and d 18 O values of primary vs. secondary material differ by 10 0 s of per mil (Defliese and Lohmann, 2015), which is not the case for closed system conditions where (by definition) d 13 C and d 18 O values remain constant. Similarly, an assumption of linear mixing is warranted for open system recrystallization products because -as shown by the carbon isotope composition of the cements-the d 13 C values of belemnites and diagenetic products are likely very similar. When d 13 C values do not differ substantially, mixing effects are not significant even with large (>10‰) differences in d 18 O values between original and diagenetic products . Using the maximum burial temperature of 115°C, we can constrain the worst effects at 3-6°C for the incorporation of 5-10% of diagenetic material, respectively. The Maximum burial temperature is reflective of the worst effects because in shallow water settings the temperatures of shallow diagenesis are more similar to where the animals lived . The impacts we predict are smaller than the variability in D 47 temperatures, cannot explain the elevated D 47 temperatures, and we emphasize that they reflect the most conservative scenario because it is unlikely that diagenesis would only proceed in the narrow time window when the section experienced its maximum burial temperatures. Indeed, a more appropriate estimate for the diagenetic end-member might be the actual D 47 temperature of the belemnite cements (51°C), which results in an effect of less than 1°C. For these reasons, we maintain that the presence of small amounts of secondary calcite cannot explain the trends or the temperatures that we observe across the Early Jurassic.

Reordering models
The Aubach section experienced either one (cold scenario) or two (warm scenario) episodes of deep burial. (Mazurek et al., 2006;Looser et al., 2019;Fig. S1 Electronic Annex). The SE-model predicts reordered D 47 values under both of these burial conditions (Electronic Annex Fig. S3), but when d 18 O temperatures are used to estimate initial D 47 values, only the warm burial scenario predicts D 47 values that are fully in the range of our observations (Electronic Annex Fig. S3). This observation suggests that the Late Jurassic veins identified by Looser et al. (2019) record the temperature of the host rock during burial, and that the warm burial scenario better describes the burial history of the Aubach section. On the other hand, the PH-model does not show appreciable amounts of reordering under either of the two burial scenarios (Electronic Annex Fig. S3). Thus, we find support only in one of the models for our assertion that the belemnites are partially reordered. Given that we have discounted other explanations for the very warm D 47 temperatures, our results suggest that at the low temperature range the PH-model underestimates the effects of burial temperatures.
Both the SE and the PH-model support our hypothesis that relative temperature differences remain when samples are partially reordered. For example, the majority (>50%) of the D 47 difference between two samples remains after burial over a wide range of temperatures (Fig. 5). This is true regardless of the choice of kinetic parameters used to parameterize the models (Fig. S4 Electronic Annex). Moreover, as demonstrated by the results of the second set of model experiments, these results can be generalized to all samples relevant for paleoclimate studies because, for a given burial history, the relationship between reordered D 47 temperatures and initial temperatures is linear (Fig. 6). The important implication behind these observations is that many sections that experienced significant amounts of burial likely still preserve records of relative temperature change. Because of partial equilibration, the SE-model predicts that at temperatures of up to 150°C more than 60% of the relative temperature difference between samples can be preserved even if samples are held at temperature for 100 s of Ma. On the other hand, although the PH-model predicts that D 47 differences are increasingly erased with time, more than 50% of the signal is preserved when samples are buried at temperatures of 150°C for up to 200 Ma. Thus, both models suggest that practically all of the Early Jurassic sections that have been studied can potentially be used to probe paleoclimate questions, with the understanding that reordered records preserve minimum estimates of temperature change across climate events. The only necessary constraint is that different samples must have experienced similar burial conditions. We note that this is the case at Aubach where samples cover a total stratigraphic thickness of 26 m, which translates into a difference in peak burial temperatures of less than 1°C under an average geothermal gradient of 25°C/km. As only part of the temperature signal remains at Aubach, it is useful to closely examine the results of the SE-model to understand how much of the signal could have been erased. With the warm and cold burial scenarios as constraints, the model predicts that after partial reordering 97 to 63% of the temperature difference between samples should remain. The range of percentages correspond to the least and most extreme possible burial conditions that can be expected given the uncertainties in the burial times and temperatures of the cold and warm burial scenarios (Electronic Annex Fig. S1). The upper range (97%) corresponds to the least extreme burial conditions that possibly occurred at Aubach given the ±3°C 95% ci uncertainty in peak burial temperatures under the cold burial scenario (87°C for 28 Ma). Likewise, the 63% corresponds to the most extreme burial conditions that are possible under the warm burial scenario. Taking into account the range of likely burial times (5-30 Ma) and peak burial temperatures (115 ± 5°C; 95% ci), these are 120°C for 30 Ma. The range of percentages is relatively wide, but it can be narrowed down by considering the belemnite D 47 data. For instance, we can discount burial conditions that result in negative temperatures when they are used to calculate original D 47 values. This is the case for many of the long burial times and high burial temperatures under the warm burial scenario (see Electronic Annex Fig. S5). When these conditions are no longer considered, the SE-model predicts that 97 to 73% of the original signal remains at Aubach. We consider that these values are good estimates for the minimum and maximum extent of reordering that occurred, and they can be used to calculate the range of initial D 47 values that existed before burial. It should be emphasized that although we cannot precisely state the amount of the signal that has been lost, the model demonstrates that the temperature changes at Aubach must have been larger than the estimates captured by the reordered belemnites.

Controls on reordering rates
There is an important caveat to our hypothesis that relative temperature differences are preserved after partial reordering. The models that we used do not consider differences in reordering kinetics between samples. If differences exist, they can erase paleoclimate information. Several possible factors have been hypothesized to affect rates. These include mineralogy (Stolper and Eiler, 2015;del Real et al., 2016;Lloyd et al., 2017Lloyd et al., , 2018Ryb et al., 2017), minor and trace element substitution (Passey and Henkes, 2012;Henkes et al., 2014), deformation (Siman-tov et al., 2016;Ryb et al., 2017), the presence/absence of water (Passey and Henkes, 2012;Brenner et al., 2018), transformation from aragonite to calcite (Staudigel and Swart, 2016;Ritter et al., 2017;Chen et al., 2019), and overburden pressure (Brenner et al., 2018). Since the rostra of belemnites only consist or low Mg calcite (Benito et al., 2016;Hoffmann et al., 2016) and samples experienced similar overburden pressures and no deformation, the relevant variables that could control kinetics in this case are minor and trace metal substitutions and changes in water fugacity. In this section, we argue that these two variables likely did not impart a measurable control on the reordering kinetics of the different belemnites.
Cation substitutions are hypothesized to affect rates mostly because they impact the amounts of crystallographic defects (Passey and Henkes, 2012). Thus they are more relevant for the PH-model than for the SEmodel, although it is possible that they could also affect rates of the latter given the strong dependence of oxygen self-diffusion in calcites on Mn concentrations (e.g., Kronenberg et al., 1984). There are currently no controlled heating experiments that systematically evaluate the role of substitutions, but Arrhenian kinetic parameters are available from calcites with different Mn and Mg contents. These can be used to evaluate the effects of these two elements. For instance, Henkes et al. (2014) compared the optical calcite of Passey and Henkes (2012) with a brachiopod calcite. They found that these two calcites have parameters for the PH-model that are indistinguishable even though the brachiopod has approximately 10 times more Mg than the optical calcite. Similarly, Stolper and Eiler (2015) compared reaction rate constants for the SE-model of the same two calcites plus the constants derived from an additional optical calcite. They observed that the rate constants for all of these materials follow (within uncertainty) the same trend in an Arrhenius plot, suggesting that regardless of starting material a single set of parameters can be used to describe their kinetics.
These results indicate that Mg concentrations do not control reordering rates. They also hint the same for Mn because the optical calcite of Stolper and Eiler (2015) has aproximately two to three times more Mn than the other two calcites. However, there is another set of experiments that can be interpreted to suggest the opposite. A spar calcite reported by Passey and Henkes (2012) shows reordering rates that are much faster and more complex than any of the other calcites. This result was interpreted by Passey and Henkes (2012) as evidence for a trace element control in reaction rates because the spar has an elevated Mn contented (up to 60 times more than the optical calcites). An alternative explanation, as discussed by Stolper and Eiler (2015), is that the knetics are similar to the other calcites, but they are observed to be different because the material has a different starting concentration (as would occur after partial reordering) of pairs than predicted by the SE-model.
The presence of water during high temperature burial is another variable that could affect reordering rates because oxygen diffusion in calcites is orders of magnitude faster under wet conditions (Farver, 1994). However, in controlled heating experiments water has only a very modest effect on reordering rates (Brenner et al., 2018).
All this shows that it is unlikely that the reordering rates of our belemnites were affected by variations in minor and trace element substitutions and/or water fugacity, but ultimately the issue of variable kinetics can be resolved with new experimental data. This is the case because the SEmodel envisions that kinetics are a universal property of a mineral, while in the PH-model kinetics can be thought of as a distinct property of each carbonate grain, controlled by grain-specific defect pools (Passey and Henkes, 2012;Lloyd et al., 2018). Recently, Chen et al. (2019) published results of low-pressure heating experiments of aragonite and found that during the transformation of aragonite to calcite D 47 values increased, paradoxically, by 0.05-0.150‰. They were able to explain their results with the SE-model and showed that D 47 values increased as a direct consequence of a high concentration of pairs (immediately adjacent singly substituted carbonate ion isotopologues) during the transition. Their results, on the other hand, cannot be explained with the PH-model because the presence of crystallographic defects cannot result in D 47 values that increase during heating against the direction of thermodynamic equilibrium. These observations strongly suggest that the SE-model describes the mechanisms of solid state 13 C-18 O reordering correctly, which in turn is consistent with similar reordering kinetics for the different belemnite samples at Aubach.

New constraints on Early Jurassic Climate
The main changes in seawater temperatures that have been previously documented across the Late Pliensbachian and the Early Toarcian are captured by the apparent D 47 temperatures (Fig. 7). For instance, we observe the coldest D 47 temperatures in the Late Pliensbachian where a cold climate has also been inferred from oxygen isotope measurements of carbonate fossils (e.g., Dera et al., 2011;Korte et al., 2015;Ruebsam et al., 2019b), and from the occurrence of glendonites and glacial deposits (Price, 1999;Rogov and Zakharov, 2010;Suan et al., 2011). Similarly, warming across the Early Toarcian is well supported by oxygen isotope measurements (e.g., Bailey et al., 2003;Dera et al., 2011;Ullmann et al., 2014;Korte et al., 2015), large perturbations of the global carbon cycle (e.g., Hesselbo et al., 2000Hesselbo et al., , 2007Hesselbo et al., , 2009Them et al., 2017;Fantasia et al., 2018a), higher abundances of clay minerals typical of warm and humid conditions (e.g., Dera et al., 2009;Fantasia et al., 2018c), and the latitudinal migrations of marine faunas (e.g., van de Schootbrugge et al., 2005;Zakharov et al., 2006). These cooling and warming patterns are also reflected in the reordered belemnites, which show cooling towards the end of the Pliensbachian and warming in response to both the Pl-To boundary event and the T-OAE (Figs. 3 and 7). The observation that our D 47 measurements closely mimic trends that were documented with independent proxies lends additional weight to the validity of these climate patterns, and to our hypothesis that meaningful paleoclimate data can be recovered from partially reordered samples. Moreover, our measurements not only confirm previous observations, but they additionally provide new constraints on the magnitude of temperature changes.
In contrast to large scale climatic trends, the magnitude of temperature changes from the Late Pliensbachian to the Early Toarcian are not well understood. Across these intervals, the oxygen isotope compositions of well-preserved belemnites, brachiopods, and bivalves vary by as much as 6‰ in subtropical latitudes (Hebrides and Cleveland Basins), and by 3‰ in slightly lower latitude basins (Lusitanian Basin, Basque-Cantabrian, Asturian Basins) (Fig. 7). If these data are interpreted in terms of temperature alone, they imply a maximum Early Toarcian warming of approx- Error bars in the corrected temperatures reflect the range in possible burial scenarios (see text and electronic annex for details). (C) Oxygen isotope composition of waters relative to the spinatum zone average for both the apparent and corrected scenarios. Grey shading is the isotope signal that would be observed if a 30-90 m glacio-eustatic sea level rise had occurred during this time period. The stratigraphic location of data was plotted at the biosubzone level using the biostratigraphic zonation of Korte et al. (2015). Samples were equally distributed when there was more than one sample per subzone.
imately 28 or 13°C, respectively. These magnitudes are thought to be too large to represent realistic changes in average seawater temperatures for mid-latitude sites (Suan et al., 2008(Suan et al., , 2010Ullmann et al., 2014Ullmann et al., , 2020, and they are much larger than 3-5°C suggested by numerical models of the T-OAE (Dera and Donnadieu, 2012). For these reasons, previous workers have argued that not all of the signal in the d 18 O data is due to temperature. For example, the more negative carbonate d 18 O values observed at higher latitudes could indicate a pronounced north-south salinity gradient during the T-OAE that resulted from enhanced rainfall/runoff (e.g., Bailey et al., 2003) or from the southward flow of Artic low salinity waters (e.g., Hesselbo et al., 2000;Bjerrum et al., 2001). Similarly, the d 18 O of seawater could have evolved towards lighter values after the decay of proposed Late Pliensbachian high-latitude ice sheets, a hypothesis that is tentatively supported by a close correspondence between temperature and sea-level curves (Suan et al., 2010;Ruebsam et al., 2019a).
Other workers, on the other hand, have noted that that a large change in water oxygen isotope compositions is not necessarily needed to explain the belemnite d 18 O data. For instance, Ullmann et al. (2014) argued that moderate warming in combination with a change in belemnite habitat could result in the 6‰ drop in belemnite d 18 O values that is observed in the Cleveland Basin. They showed that, in response to T-OAE seafloor anoxia, belemnites likely adapted to a nektonic lifestyle relative to the nektobenthic niche that they previously occupied. Thus, d 18 O temperatures can be interpreted in the warmest part of the record as upper water column temperatures, and as bottom water temperatures in the coldest parts. It is currently not known which of the different scenarios outlined above can best explain the carbonate d 18 O data. Moreover, these scenarios are not mutually exclusive and could have operated in concert to produce the d 18 O paleotemperature curve.
We observe a temperature increase between the Late Pliensbachian spinatum and the Early Toarcian falciferum biozones of 13 ± 4°C (95% ci) (Fig. 6), a value that should be considered a minimum estimate because the belemnites are partially reordered. The temperatures before burial (corrected temperatures) can be estimated with the SEmodel and the burial scenarios that correspond to the maximum (73%) and minimum (97%) likely extents of reordering (see Electronic Annex for details). When this is done, the temperature difference between the spinatum and falciferum biozones is 16 ± 7°C (Fig. 7). As noted above, a change of this magnitude is regarded as likely too large for a subtropical site. However, because this estimate is independent of seawater d 18 O, the first two scenarios (seawater d 18 O changes from salinity or ice sheet growth/decay) cannot be invoked. A combination of warming and a change in belemnite habitats, on the other hand, can potentially explain its large magnitude. We currently cannot disentangle the relative contribution of these two processes to our belemnite D 47 data, but we note that this can be done in the future when records of bottom water temperatures from benthic organisms become available. Moreover, it is also important to emphasize that our data cannot be used to infer the ocean temperature changes that occurred at other sites of the Laurasian Seaway, as magnitudes likely differed across latitudes and different oceanographic conditions. Clumped isotope data from other basins are needed to understand the regional patterns of Toarcian warming.
Both of our temperature change estimates (corrected and apparent) show that warming between the spinatum and the falciferum biozones was larger than the amount implied by the Aubach d 18 O VPDB belemnite data. For example, across the same time period belemnite d 18 O VPDB values vary by 1.4 ± 0.5 in biozone averages, which corresponds to 6.3 ± 2.3°C under constant water d 18 O values (Fig. 3). Given that the D 47 data shows more warming, the only way to reconcile the two proxies is to envision a change in seawater d 18 O values. In fact, the biozone averaged seawater d 18 O VSMOW estimates show that the value of the falciferum biozone is 1.9 ± 1.9 (±95% ci) more positive than the value of the spinatum biozone under the corrected temperature scenario, or 1 ± 1.2 (±95% ci) more positive under the apparent scenario (Fig. 7).
These d 18 O water values are not consistent with the sudden decay of large high-latitude ice sheets. Globally distributed transgressive facies during the Early Toarcian provide strong evidence for eustatic sea-level rise (e.g., Hallam, 1997;Nikitenko et al., 2008;Hermoso et al., 2013;Han et al., 2016;Krencker et al., 2019). There is currently no consensus, though, on the mechanism driving the sea-level change. Hallam (1997) used sediments from the Cleveland Basin to infer a magnitude in the order of 30-90 m and suggested a tectonic origin for the event. More recently, Ruebsam et al. (2019a) argued for a 300 kyrs duration for the main phase of sea-level rise, which implies rates that are more consistent with changes in continental icevolume . Because a glacio-eustatic mechanism would result in more negative Early Toarcian ocean d 18 O values relative to the Late Pliensbachian, our data can be used to directly test this hypothesis. Based on the Quaternary relationship between ocean d 18 O and sea-level (approx. 0.1‰ per 10 m; Miller et al., 2005), a 30-90 m sealevel rise would drive global ocean d 18 O values in the falciferum biozone to values 0.3 to 0.9 lighter than the spinatum biozone. However, both the corrected and apparent d 18 O water estimates show a trend that is in the opposite direction. Moreover, the error bars of these estimates are far off from the expected isotope signals (Fig. 7). Although our data is not consistent with glacio-eustasy, we express some caution in this interpretation because we are not able to rule out local or surface expressions in water d 18 O values. Additional data from other basins/bottom waters are necessary to verify the reproducibility of this signal (i.e., to allow the differentiation between regional vs. local effects). Given these caveats we speculate that the shift to more positive d 18 O values during the falciferum biozone may be due to a change in the water mass balance at the site with an increase in 18 O-enriched southern Tethys-sourced water masses relative to the more 18 O depleted water from the Boreal Sea. In any case, our data highlight the utility of reordered clumped isotopes to tackle these kinds of problems in 'deep-time' paleoclimate questions.

CONCLUSIONS
One of the major drawbacks of the clumped isotope proxy for deep-time climate reconstructions, is its relative sensitivity to alteration by solid state reordering. Here we demonstrate that climate patterns can still be reconstructed in sections where samples are partially reordered. The reordering models of Passey and Henkes (2012) and Stolper and Eiler (2015) were used to determine the extent to which initial temperature differences between samples are preserved when samples experience a common burial history. The models show that even though long burial times and high burial temperatures (>80°C) erase absolute temperature information and decrease the D 47 differences between samples, warming and cooling signals are preserved. For instance, the models show that more than 50% of a signal is preserved after burial at 150°C for 200 million years.
Model inferences are supported by D 47 measurements of belemnites that have partially reordered D 47 temperatures (35-48 ± 3°C) and trends that closely mimic climate patterns previously reported. The reordered belemnites record an apparent 13 ± 4°C (95% ci) temperature increase between the Late Pliensbachian the Early Toarcian. Using the burial history of the basin and the SE-reordering model we constrained the actual magnitude of change to be 16 ± 7°C. These results provide direct evidence for large temperature changes in the southwestern Tethyan epicontinental sea (Laurasian Seaway) and demonstrate the utility of partially reordered samples for paleoclimate reconstructions.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.