Clumped isotope evidence for Early Jurassic extreme polar warmth and high climate sensitivity

Periods of high atmospheric CO2 levels during the Cretaceous-Early Paleogene (~140 to 33 My ago) were marked by very high polar temperatures and reduced latitudinal gradients relative to the Holocene. These 15 features represent a challenge for most climate models, implying either higher-than-predicted climate sensitivity to atmospheric CO2, or systematic biases or misinterpretations in proxy data. Here, we present a reconstruction of marine temperatures at polar (>80°) and mid (~40°) paleolatitudes during the Early Jurassic (~180 My ago) based on the clumped isotope (Δ47) and oxygen-isotope (δOc) analyses of mildly buried pristine mollusc shells. Reconstructed calcification temperatures range from ~8 to ~18°C in the Toarcian Arctic and from ~24 to ~28 °C 20 in Pliensbachian mid-paleolatitudes. These polar temperatures were ~10-20°C higher than present along with reduced latitudinal gradients. Reconstructed seawater oxygen isotope values (δOw) of – 1.5 to 0.5‰ VSMOW and of – 5 to – 2.5‰ VSMOW at mid and polar paleolatitudes, respectively, point to a significant freshwater contribution in Arctic regions. This highlight the risk of assuming the same δOsw value for δO-derived temperature from different oceanic regions. These findings provide critical new constraints for model 25 simulations of Jurassic temperatures and δOsw values and suggest that high climate sensitivity is a hallmark of greenhouse climates since at least 180 My.

temperature estimates for ante-Cretaceous periods are mostly derived from the oxygen isotope composition of marine carbonate fossils (δ 18 Oc), with well-known limitation related to uncertainties in the past δ 18 O signature of 40 seawater (δ 18 Ow) (Epstein et al., 1953;Roche et al., 2006;Laugié et al., 2020).
Here, we use carbonate clumped isotope thermometry (Δ47), to simultaneously constrain the calcification temperatures and associated δ 18 Ow values of marine carbonate shells collected from Lower Jurassic sedimentary successions with exceptionally shallow to moderate burial depths spanning subtropical to polar paleolatitudes.
We compare our results to existing Jurassic to Eocene climate proxy data and simulations and discuss their 45 implications for climate sensitivity under greenhouse conditions.

50
The Polovinnaya River section is located in northern Siberia (72°36'05" N, 107°58'52.2" E), and was located near the north pole during the early Jurassic (Fig. 1). Our bivalve samples come from between 0 and 14m in the section and belong to the Toarcian (Suan et al., 2011). This interval consists of silty shale slightly enriched in organic carbon (TOC ~0.5%). Except for localized carbonate concretions, there is no carbonate fraction in the sediment. The studied interval has been previously correlated to the lower Toarcian Serpentinum ammonite zone 55 based on biostratigraphy of foraminifera and dinoflagellate cyst, and lithostratigraphic correlation with other sections of the basin (Suan et al., 2011). This section records very abundant Dacryomya bivalve shells (Fig. S1), an opportunistic suspension-feeder genus tolerant to poorly oxygenated waters, which preferred conditions with weak hydrodynamics (Zakharov and Shurygin, 1978). Few belemnite rostra were also recorded as well isolated fish scales and teeth. The section undergone exceptionally low burial as suggested by the low values of Rock-60 Eval Pyrolysis Tmax (mean = 420°C) previously measured in the host sediments (Suan et al., 2011). Regional stratigraphy from the more distal Anabar area suggest local overburden not exceeding 1000m: a total overburden (Lower Toarcian to Valanginian) of about 380m is recorded in the Anabar River area (Nikitenko et al., 2013) located 200 km East of the Polovinnaya section, which may be extended to about 1000m when adding Valanginian-Cenomanian overburden from the more distal Bol'shoi Begichev islands. Modern local geothermal  Dera et al. (2009). (b) Arctic map modified from Nikitenko and Mickey (2004). (c) Tethyan map 70 modified from Thierry, (2000). Localities: PR = Polovinnaya River; Wq = Warcq.

Warcq
Samples from north-eastern Paris Basin were collected in 2014 from a temporary road cutting located near Warcq, Ardennes, France (49°45'21.6"N 4°39'28.8"E). They consist of grey silty claystone with lenses of packed 75 carbonated shell fragments, mainly from a variety of bivalves (Grammatodon, Malletia, Limea, Oxytoma) and few ammonoids (Beaniceras, Aegoceras?, Dactyloceratidae) (Fig. 2). The lithology, fossil preservation and assemblages of the sampled beds is similar to those described by Thuy et al. (2011), from a nearby site of Sedan and dated from the Pliensbachian Davoei zone. The sampled levels are therefore tentatively attributed to the lower Pliensbachian Davoei ammonite zone. Mean Tmax values of 425°C and maximum burial temperatures near 80 60°C have been reported for Pliensbachian sediments of the NE Paris Basin (Disnar et al., 1996;Blaise et al., 2014). These burial temperatures should be regarded as an upper limit, as the very proximal sampling area near Warcq was repeatedly emerged during the Mesozoic and hence shows a much-thinner Mesozoic cover than these more distal sites (Waterlot et al., 1960).  Dactylioceratidae indet., specimen ARD-07.

Sampled material
The two studied sites present exceptionally rare records of aragonite preservation for the Lower Jurassic interval.
Dacryomya shells are the most abundant macrofossil and the unique bivalve genera to occur in the Polovinnaya River section. They are very abundant in the lower part of the section (0 to 8m). They are mainly represented by adult shells, while juveniles are common in only few levels. They appear as ~1cm distinct individual or detached 100 valves, sometimes close to each other (Fig. 2). The carbonate shells, often flattened and partially to entirely preserved, are a few millimetres thick but brittle and detached easily from their inner and outer mould. Their cream to white colour contrasts with the dark aspect of the sediment, and few thicker individuals are iridescent.
Molluscs shells from Warcq, clearly show a more energetic environment as they mostly appear as packed shell fragments showing a higher taxonomic diversity relative to the other site. Few complete individuals and 105 separated valves can be observed among the debris with their associated mould in or around the remaining shell.
Shells are cream to clear white, with some showing iridescence.
The microstructural preservation state and mineralogy of the analysed bivalve and ammonite shells were investigated using a Phenom Pure G2 scanning electron microscope in backscatter mode (BSEM G2) and Raman spectroscopy using an XploRA Raman microscope in Laboratoire de Géologie de Lyon (LGLTPE).

Geochemical analysis and data processing
The remnants of carbonate shells were sampled as a whole using dental tools under a binocular-microscope. A carbonate vein and matrix from the carbonate nodule Pol 29 were also sampled to constrain the geochemistry of 115 this potential diagenetic phase.
The Δ47 and δ 18 O values of 13 samples were measured (1 to 5 replicates each) using methods described by (Daëron et al., 2016). Carbonate samples were converted to CO2 by phosphoric acid reaction at 90 °C in a common, stirred acid bath for 15 minutes. Initial phosphoric acid concentration was 103 % (1.91 g/cm3) and each batch of acid was used for 7 days. After cryogenic removal of water, the evolved CO2 was helium-flushed 120 at 25 mL/mn through a purification column packed with Porapak Q (50/80 mesh, 1 m length, 2.1 mm ID) and held at −20 °C, then quantitatively recollected by cryogenic trapping and transferred into an Isoprime 100 dualinlet mass spectrometer equipped with six Faraday collectors (m/z 44-49). Each analysis took about 2.5 hours, during which analyte gas and working reference gas were allowed to flow from matching, 10 mL reservoirs into the source through deactivated fused silica capillaries (65 cm length, 110 μm ID). Every 20 minutes, gas 125 pressures were adjusted to achieve m/z = 44 current of 80 nA, with differences between analyte gas and working gas generally below 0.1 nA. Pressure-dependent background current corrections were measured 12 times for each analysis. All background measurements from a given session are then used to determine a mass-specific relationship linking background intensity (Zm), total m/z = 44 intensity (I44), and time (t): Zm = a + bI44 + ct + dt2. Background-corrected ion current ratios (δ45 to δ49) were converted to δ13C, δ18O, and "raw" Δ47 values 130 https://doi.org/10.5194/cp-2021-79 Preprint. Discussion started: 29 July 2021 c Author(s) 2021. CC BY 4.0 License. as described by Daëron et al., (2016), using the IUPAC oxygen-17 correction parameters. The isotopic composition (δ13C, δ18O) of our working reference gas was computed based on the nominal isotopic composition of carbonate standard ETH-3 (Bernasconi et al., 2018) and an oxygen-18 acid fractionation factor of 1.00813 (Kim et al., 2007). Raw Δ47 values were then converted to the I-CDES Δ47 reference frame by comparison with four "ETH" carbonate standards (Bernasconi et al., 2021) using a pooled regression approach 135 (Daëron, 2021). Full analytical errors are derived from the external reproducibility of unknowns and standards (Nf = 89) and conservatively account for the uncertainties in raw Δ47 measurements as well as those associated with the conversion to the "absolute" Δ47 reference frame.
Complementary δ 13 C and δ 18 O analyses of the smallest Arctic shells were performed at LGLTPE, using a Multiprep TM automated sampler coupled to a dual-inlet GV Isoprime TM mass spectrometer. Samples were

145
Clumped isotope temperatures were computed based on the I-CDES calibration of (Anderson et al., 2021).
Temperature uncertainties correspond to the fully-propagated 95% confidence intervals from Δ47 measurements of each sample (Daëron, 2021), neglecting the much smaller uncertainties in the calibration. The δ 18 O values from aragonite samples were adjusted considering the different phosphoric acid fractionation factors for calcite and aragonite (Kim et al., 2007). The δ 18 Ow values relative to VSMOW was estimated using Δ47-derived 150 temperatures and the equations of Grossman and Ku (1986) and Kim and O'Neil (1997) for mollusc shells and calcite vein respectively.
Paleolatitude of the studied sites was computed using the online paleolatitude calculator paleolatitude.org (van Hinsbergen et al., 2015) computed with the model of Torsvik et al. (2012).

Results
The SEM observations of shell fragments of Dacryomya jacutica revealed well preserved sheet nacreous microstructures underlying a prismatic layer we interpret as the outer shell layer (Fig. 3). All Raman spectra gathered from Dacryomya jacutica shells confirm that the original aragonite mineralogy is preserved. Mollusk shells from Warcq showing an aragonite mineralogy revealed microstructures similar to those observed in 160 Dacryomya jacutica, the main differences being that sheet nacreous structures of the studied ammonite shell (ARD-05) shows thinner tablets than those of bivalve shells (Fig. 3). Both SEM and Raman data indicate that the sample ARD-03 (bivalve fragment) is in calcite, showing a darker colour and no iridescence, with a much simpler and massive structure observed in SEM (Fig. 3

Sample preservation
The SEM and Raman observations reveal that the analysed mollusc shells from both sites retain pristine 185 aragonite mineralogy and microstructures with no evidence of recrystallization (Fig. 3). Despite their aragonite mineralogy, the Dacryomya shells from sample Pol-29 record unusually low δ 13 C that are ~8‰ lower than the other Dacryomya shells analysed from the same succession. The carbonate matrix of the nodule where these shells are embedded also records a very low δ 13 C value (-21.43‰) but a δ 18 O value within the range of the bivalve shells. We therefore attribute the extremely low δ 13 C values of bivalves shells of this level to an early carbonate minerals was recently suggested as an alternative process that may alter the clumped isotope signature of biogenic carbonates without substantially affecting the stable isotope signature of the shell nor its mineralogy 200 (Nooitgedacht et al., 2021). In their heating experiments, these exchanges resulted in a significant decrease of the Δ47 value of the bivalve shells compared to the original shell, and a minor (~0.1‰) decrease in δ 18 O of the heated shell. We cannot exclude that this process has altered the fossils studied here even at low temperature, nor do we have evidence that it occurred. The Δ47 temperature of 31.1 ± 4.8 °C for the fracture-infilling calcite vein in Arctic Russia is significantly higher than those inferred from bivalves and is consistent with a formation depth 205 <1 km assuming a geothermal gradient of 25°C/km. The reconstructed δ 18 Ow value of -10.7± 0.9 ‰ for this calcite vein is also substantially lower than those inferred from associated bivalves, consistent with a late-phase meteoric source for the mineralizing fluid.

Event
Bivalve shells record a marked rise in δ 13 C along the section up to ~5‰ that parallels that recorded by organic carbon δ 13 C data ( Bivalve shell growth can be highly variable during the animal life (Schöne, 2008), making any paleoenvironmental record derived from bivalve shell either incomplete (because of growth cessation) or at least biased towards the period of maximum growth rate. Shell growth rate can be controlled by both environmental parameters (Temperature, salinity, food availability …), biological processes such as spawning and ontogeny 225 (Schöne, 2008). One major aspect of shell growth that may bias the geochemical signal data is seasonal shell growth-cessation. In modern high latitude bivalves, seasonal shell growth-cessation generally occurs during the winter, triggered by low temperatures or low food availability (Peck et al., 2000;Vihtakari et al., 2016;Killam and Clapham, 2018).
It appears likely that food availability declined markedly during the Early Jurassic polar night, which would have 230 certainly led to winter growth cessation in the analysed Dacryomya shells. In the present-day Nucula annulata, an aragonite bivalve with similar ecology to the analysed Dacryomya jacutica, growth cessation occurs in winter and during spawning at peak local temperatures, its average δ 18 Oc hence records late spring to early fall SST (Craig, 1994). By contrast, growth band δ 18 Oc offers evidence for summertime-only growth cessation in highlatitude Eocene bivalves from Antarctica, with inferred winter SST of 11.1±0.6 and summer SST of 17.6±1.3°C 235 (Buick and Ivany, 2004;Douglas et al., 2014) very close to our maximal and minimum SST estimates of the Toarcian Arctic. A comparable seasonal δ 18 Oc record could not be generated from our Russian Arctic material owing to the very small size of the available Dacryomya shells (1 to 2cm). In any case, the temperate data from https://doi.org/10.5194/cp-2021-79 Preprint. Discussion started: 29 July 2021 c Author(s) 2021. CC BY 4.0 License. NE France should be minimally affected by seasonal biases as shell precipitation occurs more continuously throughout the year in modern temperate molluscs (Killam and Clapham, 2018). Besides, both sites were 240 deposited in near shore environments at very shallow depths likely not exceeding a few tens of meters (Thuy et al., 2011;Suan et al., 2011). Although bivalves from both sections record temperatures near the sea bottom that were likely slightly cooler than the sea surface, the difference should not exceed a few degrees. We therefore conservatively interpret the reconstructed temperatures as reflecting polar warm-season SST (summer; SSTPWS) in Arctic Russia and low latitude annual SST in NE France. These SSTPWS for the T-OAE are still 10-20°C 245 higher than present-day SSTPWS (Fig. 5).
The reconstructed polar δ 18 Ow values ranging from -4.9 ±1.2 to -2.5 ±0.8‰ during the T-OAE are significantly lower than the value of -1‰ expected for an ice-free world mean open ocean (Shackleton and Kennett, 1975).
These results imply a substantial freshwater contribution to the studied basin during the T-OAE, probably resulting from coastal runoff at this relatively proximal site, as evidenced from abundant terrestrial organic 250 matter (Suan et al., 2011). High temperatures and reduced salinity are in broad agreement with paleontological evidence for warm and humid temperate conditions during the T-OAE interval in Arctic Siberia (Rogov et al., 2019). Arctic coast salinity can be inferred from Δ47-derived δ 18 Ow with reasonable assumptions of local δ 18 O of precipitations and runoff (δ 18 Op). Higher polar temperatures should have produced higher δ 18 Op than those prevailing today (Rozanski et al., 1992), a hypothesis supported by terrestrial plants n-alkanes hydrogen isotopes 255 and paleosol siderite Δ47 data indicating Early Eocene δ 18 Op of -10 to -15‰ van Dijk et al., 2020). Assuming a similar range of δ 18 Op values in the Early Jurassic Arctic, and a global mean ocean with a salinity of 34.5‰ and a δ 18 Ow of -1‰, mass balance calculations indicate mean salinity of 23.9±2.9‰ (1σ, n=8) and 27.7±1.8‰ (1σ, n=8) with δ 18 Op of -10 and -15‰, respectively (see Supplementary Data). These estimates point towards brackish waters at Polovinnaya River during the Toarcian, consistent with the fossil assemblages 260 of the succession that includes abundant terrestrial organic matter and wood debris, marine to brackish elements such as abundant dinoflagellate cysts, benthic foraminifera (preserved as organic linings and agglutinate forms) but rarer typically marine elements that are only represented by a few belemnite rostra and unidentifiable ammonite internal moulds (Suan et al., 2011). Such values are also comparable to the salinity of 28‰ estimated using a fully coupled ocean-atmosphere model for the Toarcian (Dera and Donnadieu, 2012), although the 265 Arctic temperature obtained by the same model are in strong disagreement with our data (Fig. 5). Such observations should be replicated around the Arctic realm to test whether the brackish environment evidenced here is of local or more regional nature.  δ 18 Ow was calculated using the oxygen isotope fractionation equation of Grossman and Ku (1986).

275
The mid paleolatitude SST reconstructed by our new clumped isotope data are in good agreement with recent Sinemurian-Pliensbachian and Toarcian TEX86 H data pointing to summer SST~20-30°C at slightly lower paleolatitudes (Robinson et al., 2017;Ruebsam et al., 2020). Considering the scarcity of other Early Jurassic temperature proxy-data, model-based SST and δ 18 Ow estimates, we extend the comparison to SST and δ 18 Ow estimates based on various proxy-data and published Earth system simulations for other Jurassic to Eocene 280 intervals ( Fig. 5; Supplement). This compilation shows that Δ47 SST from NE France agree with most previous TEX86 H and Δ47 SST for the Jurassic-Eocene interval, with values > 5°C higher than present-day SST. This comparison also shows that the Maastrichtian was characterized by substantially lower Δ47 SST, in line with independent evidence for global cooling during this interval (O'Brien et al., 2017;Pucéat et al., 2003), although many Δ47 records from this stage come from the Western Interior Seaway (WIS) where temperatures were likely 285 influenced by specific regional patterns, such as southward influence of arctic water or significant freshwater contribution in the basin (Coulson et al., 2011;Petersen et al., 2016). This compilation also reveals lower SSTs near 50°N in the Callovo-Oxfordian. Russian Arctic SST are very close to Campanian-Maastrichtian and Early Eocene TEX86 H -derived SST (Jenkyns et al., 2004;Sluijs et al., 2020) and Early Eocene bioclimatic SST (Suan et al., 2017) from the Arctic (Fig. 5). The Δ47 data presented herein suggest a decrease in mean SST of 290 0.26±0.05°C per ° of latitude, i.e., a reduction of the latitudinal SST gradient of 32±10% relative to present, consistent with the most conservative Early Eocene estimates (Evans et al., 2018).
Our Δ47 SST for the Lower-Jurassic can be compared to published results from Earth System models that simulate past intervals of global warmth to discuss model-data discrepancies. Proxy-data indicate an atmospheric pCO2 of 1000±500 ppmv during the Early Jurassic, with maximum values of 1750±500 ppmv, i.e., 6x pre-295 industrial levels (PIL), during the T-OAE (McElwain et al., 2005;Li et al., 2020). Earth system models run at 6x PIL for Early Jurassic (Dera and Donnadieu, 2012) or Cretaceous-Eocene paleogeography almost invariably https://doi.org/10.5194/cp-2021-79 Preprint. Discussion started: 29 July 2021 c Author(s) 2021. CC BY 4.0 License. produce lower SST than those inferred from our Δ47 data, with a maximum model-data discrepancy of >15°C at high latitudes (Fig. 5). To achieve such polar warmth, the Eocene CCSM3 simulations require 16x PIL, more than twice that indicated by Lower Jurassic and Eocene proxy data (Huber and Caballero, 2011 (Laugié et al., 2020).
The hypothesis of shell growth restricted to warmest month in the analysed Toarcian Arctic bivalves, however, remains questionable given the evidence for summertime-only growth cessation in Eocene bivalves from Antarctica (Buick and Ivany, 2004). Finally, Arctic SST as high as 15-20°C are successfully achieved in the 305 Eocene CESM1.2 CAM5 at 6 to 9x PIL (Zhu et al., 2020), in which climate sensitivity increases with rising CO2 due to low-altitude cloud albedo feedbacks and improved radiative parameterization. As this model produces an increase in climate sensitivity with CO2 in both Eocene and modern conditions, our results thus support the growing body of evidence that the amplitude of the future anthropogenic warming may be underestimated by conventional state-of-the-art models.

310
The reconstruction of δ 18 Ow values using proxy data provides a complementary aspect to assess model capabilities, as this indicator is sensitive to both climate parameters (moisture, humidity and temperatures) and paleogeography. Our mid-latitude δ 18 Ow are broadly similar to those reconstructed using marine turtle bones δ 18 OPO4 and Δ47 data from Jurassic to Eocene bivalves, ammonites, foraminifera, as well as belemnites (Fig. 5).
As for temperatures, regional processes in the Western Interior Seaway explains the mid latitude very low δ 18 Ow 315 indicated by most late Cretaceous Δ47 data (Petersen et al., 2016).
We are aware of only few Earth-system δ 18 Ow simulations for the interval considered here (Zhou et al., 2008;Tindall et al., 2010;Zhu et al., 2020), hence limiting model-data comparisons. The higher freshwater contribution near high-latitude landmasses of the Northern Hemisphere in both models produced lower δ 18 Ow that are broadly consistent with previous and our proxy data (Fig. 5). This good agreement, however, might be 320 partly fortuitous, as proxy data suggest SST much higher than those produced by both models (Fig. 5). As mentioned above, such higher-than-predicted polar warmth would have substantially increased high-latitude δ 18 Op, so that higher runoff would be required to reproduce the magnitude of the poleward drop in δ 18 Ow indicated by proxy data. This highlights the usefulness, in future models of past greenhouse climates, to systematically provide δ 18 Ow predictions so that δ 18 Ow estimates derived from Δ47 data may serve as a constraint 325 on Earth system models.

330
(e, f). Marker colour in proxy data shows sample age (see key). Results of Earth system simulations are shown as annual averages (bold lines) and summer and winter seasonal averages (colour shading). See supplement for references used in this compilation.

335
The clumped isotope compositions of pristine, minimally buried, marine molluscs shells yield SST >25°C at mid-latitudes during the early Pliensbachian and SST >10°C at polar paleolatitudes during the T-OAE. The reconstructed δ 18 Ow values point to higher freshwater contribution toward Arctic regions, illustrating the dangers of assuming a fixed global δ 18 Ow value for δ 18 O-derived temperature reconstructions. Although further work should clarify the influence of seasonal changes in the recorded SST values at polar sites, these results strengthen 340 a growing body of evidence for higher climate sensitivity under high atmospheric CO2 conditions and suggest that this higher sensitivity is a general feature of greenhouse climates since at least 180 Ma.

Data availability
Detailed data supporting this study are available in the supplementary material. Raw data are available on 345 request to the author.

Author contribution
TL and GS designed the study and led the writing in close cooperation with CL, MR and MD. MR and GS participated to the field work and collected the samples. TL prepared and sampled the shell material for 350 geochemistry and performed the SEM observations. MR, JS and OL identified the fossils. MD and TL performed the clumped isotope analyses and data processing. AV-L and TL performed the stable isotopes https://doi.org/10.5194/cp-2021-79 Preprint. Discussion started: 29 July 2021 c Author(s) 2021. CC BY 4.0 License.