Assessing environmental change associated with early Eocene hyperthermals in the Atlantic Coastal Plain, USA

Eocene transient global warming events (hyperthermals) can provide insight into a future warmer world. While much research has focused on the Paleocene–Eocene Thermal Maximum (PETM), hyperthermals of a smaller magnitude can be used to characterize climatic responses over different magnitudes of forcing. This study identifies two events, namely the Eocene Thermal Maximum 2 (ETM2 and H2), in shallow marine sediments of the Eocene-aged Salisbury Embayment of Maryland, based on magnetostratigraphy, calcareous nannofossil, and dinocyst biostratigraphy, as well as the recognition of negative stable carbon isotope excursions (CIEs) in biogenic calcite. We assess local environmental change in the Salisbury Embayment, utilizing clay mineralogy, marine palynology, δ¹⁸O of biogenic calcite, and biomarker paleothermometry (TEX₈₆). Paleotemperature proxies show broad agreement between surface water and bottom water temperature changes. However, the timing of the warming does not correspond to the CIE of the ETM2 as expected from other records, and the highest values are observed during H2, suggesting factors in addition to pCO₂ forcing have influenced temperature changes in the region. The ETM2 interval exhibits a shift in clay mineralogy from smectite-dominated facies to illite-rich facies, suggesting hydroclimatic changes but with a rather dampened weathering response relative to that of the PETM in the same region. Organic walled dinoflagellate cyst assemblages show large fluctuations throughout the studied section, none of which seem systematically related to CIE warming. These observations are contrary to the typical tight correspondence between climate change and assemblages across the PETM, regionally and globally, and ETM2 in the Arctic Ocean. The data do indicate very warm and (seasonally) stratified conditions, likely salinity-driven, across H2. The absence of evidence for strong perturbations in local hydrology and nutrient supply during ETM2 and H2, compared to the PETM, is consistent with the less extreme forcing and the warmer pre-event baseline, as well as the non-linear response in hydroclimates to greenhouse forcing.

1 Introduction

The early Eocene is punctuated by a series of transient phases of carbon input and global warming, which are termed hyperthermals (e.g., Thomas and Zachos, 2000; Cramer et al., 2003; Littler et al., 2014). The well-studied Paleocene–Eocene Thermal Maximum (PETM) may represent the closest approximation to modern-day warming in the geologic record (Zachos et al., 2008). However, debate exists as to the origins of this event with respect to the primary source of carbon (Zeebe and Lourens, 2019; Zachos et al., 2010; Reynolds et al., 2017; Dickens et al., 1997; Frieling et al., 2019; Gutjahr et al., 2017). The mass of C released and rise in pCO₂ are constrained by records of seafloor carbonate dissolution and surface ocean pH (Zachos et al., 2005; Penman et al., 2014; Gutjahr et al., 2017). Along continental margins, for example, in the Salisbury Embayment located within the mid-Atlantic Coastal Plain of the United States, the PETM is marked by sedimentological and biotic evidence of strong hydrological perturbations (Kopp et al., 2009; Sluijs and Brinkhuis, 2009; Stassen et al., 2012b; Self-Trail et al., 2012, 2017b; Rush et al., 2021). The Eocene Thermal Maximum 2 (ETM2) occurred approximately ∼ 2 Myr after the PETM and is characterized by warming and carbon release of approximately one-third that of the PETM (Lourens et al., 2005; Sluijs et al., 2009; Stap et al., 2010; Dunkley Jones et al., 2013; Gutjahr et al., 2017; Harper et al., 2020). ETM2 can provide insight into the climatic impacts of a smaller carbon forcing (pCO₂) when compared to the PETM from geographically similar sites and hence provide insight into sensitivity to greenhouse forcing.

Intensification of the hydrologic cycle is a critical feature of the climatic response to greenhouse forcing (e.g., Carmichael et al., 2017). Evidence for a pronounced mode shift during the PETM includes enhanced sediment fluxes in sections from the mid-Atlantic margin, which suggest an increased frequency of precipitation extremes leading to high-energy fluvial activity and the associated abundances of low-salinity-tolerant biota (John et al., 2008; Self-Trail et al., 2017b; Stassen et al., 2012b; Rush et al., 2021; Sluijs and Brinkhuis, 2009). Similar signals have been found in many regions across various hyperthermals, including the Pyrenees, the margins of the Tethys, Aotearoa / New Zealand, the North Sea basin, and the Arctic (Schmitz and Pujalte, 2007; Jiang et al., 2021; Slotnick et al., 2012; Sluijs et al., 2014, 2007; Jin et al., 2022; Harding et al., 2011). These changes in extreme precipitation patterns are consistent with projections for future change (Pfahl et al., 2017; Swain et al., 2018). Studying the impacts of hydroclimatic changes in response to differing degrees of carbon forcings can inform as to the linearity of these responses, called “tipping point” thresholds, and the background variability in the hydrosphere.

Given the extensive evidence of a mode shift in coastal mid-Atlantic hydroclimate during the PETM, we initiated a regional search for ETM2. We here identify ETM2 and H2 (∼ 54 Ma) in the Knapps Narrows core based on biostratigraphic, paleomagnetic, and chemostratigraphic data. We present proxy records of the regional climatic response, including foraminiferal δ¹⁸O- and TEX₈₆-based temperature reconstructions, clay mineralogy, and organic walled dinoflagellate cyst (dinocyst) assemblages and compare the results to published records of the PETM from this region.

2 Materials and methods

2.1 Material

The Knapps Narrows core was drilled at 38.72129^∘ N, 76.33162^∘ W on the eastern shore of Maryland, in close proximity to drill sites which have been extensively studied for the PETM, South Dover Bridge, and Howards Tract (Fig. 1). The cored target interval lies between 84–102 m in the Nanjemoy Formation, an Eocene glauconitic sand in the upper Pamunkey group deposited in a shallow marine environment (Darton, 1951; Self-Trail et al., 2017b). From the bottom of the section to 88.4 m, the lithology is dominated by very coarse- to fine-grained, poorly sorted, angular to subrounded sand, comprised of ∼ 50 % glauconite, ∼ 30 % quartz, and < 5 % muddy matrix. Between 88.4 m to the top of the section, the core largely consists of very coarse- to medium-grained, subrounded to angular clayey sand with poorly sorted glauconite (∼ 40 %) and quartz (∼ 10 %) in a muddy matrix (20 %). Sand fines are found upsection, and the clay content decreases. The sediments are heavily bioturbated throughout. The Nanjemoy typically overlays Marlboro Clay, which was deposited during the early Eocene PETM, with a significant regional unconformity (Darton, 1951; Stassen et al., 2012a) In this location, due to local faulting, the Marlboro Clay is entirely absent, and an expanded, slightly younger, lower Eocene interval equivalent to the calcareous nannofossil biozone NP11, which contains ETM2, is present (Martini, 1971; Cramer et al., 2003; Lourens et al., 2005). This disconformity and resultant paleobathymetry likely contributed to the enhanced sedimentation rates and expanded section recorded in the core (Appendix Fig. A1). The initial core description noted much coarser grain sizes when compared to the lithological expression of the Nanjemoy in nearby cores, further supporting the interpretation of rapid infilling during this time. All raw data are included in the supplemental datasets that are available through PANGAEA (see the “Data availability” section at the end of the paper).

Figure 1 Map showing the location of Knapps Narrows (KN), the South Dover Bridge (SDB), and Howards Tract (HT) coreholes. The Fall Line delimits the crystalline basement rock (gray) from Cretaceous–Cenozoic-age coastal plain sediments (white). Modified from Self-Trail et al. (2017b).

2.2 Methods

2.2.1 Nannofossil biostratigraphy

The calcareous nannofossil biozonation was established, and core descriptions were made in the field and at the U.S. Geological Survey (USGS) Florence Bascom Geoscience Center. A total of 42 samples, taken at approximately 0.15 m intervals for calcareous nannofossil analysis, were extracted from the center of freshly broken core to avoid contamination from drilling fluid. Smear slides for calcareous nannofossil analysis were prepared using the standard technique of Bown and Young (1998), combined with the double-slurry technique of Blair and Watkins (2009) and mounted using Norland Optical Adhesive 61. Slides were examined under cross-polarized light, using a Zeiss Axioplan 2 light microscope (LM) at 1250× magnification. Biostratigraphic zones were based on the NP zonation of Martini (1971) and supplemented by the biohorizons of Agnini et al. (2014). Smear slides are housed in the USGS calcareous nannofossil laboratory in Reston, VA.

2.2.2 Paleomagnetism

Paleomagnetic data were generated at the Paleomagnetism Lab of Lehigh University. A total of ∼ 100 minicores were collected from the Knapps Narrows core (40–110 m), with sampling spacing of 0.2 to 1.5 m guided by preliminary biostratigraphic constraints. Stepwise thermal demagnetization in an ASC Model TD48-SC thermal demagnetizer was applied to the samples, followed by measurement on a three-axis 2G Enterprises 755 superconducting rock magnetometer. Heating involved at least seven treatments, at ∼ 25 ^∘C/50 ^∘C steps from 100 to 350 ^∘C for the weak samples and up to 575–600 ^∘C for stronger ones. Thermal demagnetization generally ceased when the remanent magnetization displayed anomalous surges in magnetization, was too weak for magnetometer precision, or exhibited irregular magnetic directions and/or intensities for two or three consecutive steps. The choice of a thermal rather than an AF (alternating field) demagnetization procedure was due to the potential existence of high-coercivity magnetic minerals (e.g., goethite).

Samples from rotary cores may have been shifted from their original declination orientation. The overprint components from drilling-induced magnetization and from normal-polarity, present-day (late Quaternary) field remanence are observed and could often be easily removed during initial demagnetization steps. Instead, normal (N) or reversed (R) polarity could be indicated by natural remnant magnetization (NRM), which pointed downward (N) or upward (R). The present-day magnetic field direction is −11.11^∘ declination and 64.88^∘ inclination (from https://www.ngdc.noaa.gov/geomag-web/#igrfwmm, last access: 10 August 2023, for the modern Knapps Narrows core location), which should be roughly the same as the Eocene, since no major continental movement has occurred (Scotese, 2014).

Polarity zones were assigned to stratigraphic clusters of good-quality characteristic remanent magnetization (ChRM) vectors, while uncertain intervals were assigned to clusters of noisy magnetic behaviors and to significant gaps in the sampling coverage (Przybylski et al., 2010). Detailed paleomagnetic data and additional data analysis protocol are available in the supplementary paleomagnetism Microsoft Excel workbook (see the “Data availability” section).

Given that sediments are largely semi-consolidated, a non-magnetic plastic tube was inserted in the core to extract a sediment plug. The plug was then extruded into aluminum foil pieces and tightly wrapped. Occasionally, paleomagnetic samples were taken by pushing non-magnetic cubes directly into the center of the core halves so that alternate field demagnetization could be used.

The interpreted polarity and characteristic directions of the samples were given a quality rating of “N(R)”, “NP(RP)”, “NPP(RPP)”, “N?(R?)”, or “INT” to samples yielding definite to uncertain ChRM directions, according to a semi-subjective judgment of the behavior of the magnetic vectors during the stepwise demagnetization. This quality rating helps to define a more accurate magnetostratigraphy (Zhang and Ogg, 2003). The N or R ratings were generally assigned to those samples that attained a (semi-)stable endpoint direction (characteristic remanent magnetization or ChRM) during progressive demagnetization. However, only a small subset of our samples displayed this ideal resolution of primary magnetization. We applied a “P” tag of RP or NP when the residual magnetization vector was considered close to attaining an endpoint before losing its residual magnetization or experiencing a surge. The “PP” tag was applied to samples that had a distinct trend toward the polarity hemisphere, which was considered to be indicative of the underlying polarity but also considered to be too far from attaining an endpoint before dying to be used in statistics for computing a mean direction. The “?” qualifier was used to denote possible trends toward an underlying polarity, and INT is for classifying elements that were either entirely uncertain or displayed an endpoint that was intermediate between the N and R poles. Polarity zones were assigned to stratigraphic clusters of ChRM vectors rated as N-NP-NPP or R-RP-RPP. Uncertain intervals were assigned to clusters of N?-INT-R? and to significant gaps in sampling coverage. Summary plots of the lithology of the Knapps Narrows with the interpreted polarity patterns (Fig. A2) were generated with the public TS-Creator database and visualization software (https://timescalecreator.org, last access: 10 August 2023).

2.2.3 Bulk isotope stratigraphy and age model

Bulk carbonate content and bulk carbonate and benthic foraminiferal stable carbon and oxygen isotope analysis was performed in the Earth and Planetary Sciences Department at the University of California, Santa Cruz. CaCO₃ content was determined using a UIC CM140 coulometer. In this system, approximately 20 mg of sediment is dissolved in 2N sulfuric acid. For isotope analyses of bulk carbonate, a mass of ground sediment was weighed to obtain between 35 and 50 µg of carbonate. All samples were analyzed with a Kiel IV carbonate device paired with a Thermo Scientific MAT253 mass spectrometer. The analytical precision of this system, based on replicate analyses of Carrara marble standards, is better than 0.05 ‰ for carbon isotopes.

An age model was constructed for this section on the basis of carbon isotope stratigraphy, specifically the carbon isotope excursion (CIE) minima and maxima (respectively, H1 and H2), correlated to the orbitally tuned deep ocean records of Lauretano et al. (2015). This age model was then used to calculate sedimentation rates.

2.2.4 Spectral analysis

Spectral analysis (eCOCO or evolutionary correlation coefficient) of the gamma ray wireline log of the cored interval was performed via Acycle v2.4.1 to obtain an estimate of the sedimentation rate variations (Li et al., 2019). The eCOCO analysis computes the correlation coefficient (COCO; ρ) between the frequency spectra of target astronomical solutions and the presence of different cycles within studied data, which provides the numerical estimation of the optimal depositional rates (a higher correlation coefficient corresponds to the most likely sedimentation rates) within a sedimentary record (Li et al., 2018, 2019). Monte Carlo simulations are embedded in the COCO method to test the significant level of astronomical forcing assumption.

2.2.5 Palynology

In total, 18 samples were processed at Utrecht University for quantitative palynological analysis. Between 3 and 15 g of freeze-dried, lightly crushed sediment was spiked with a known number of exotic Lycopodium spores and subsequently treated once with 30 % HCl and twice with ∼ 38 %–40 % HF to dissolve carbonates and silicates, respectively. The insoluble residue was sieved using 15 and 250 µm nylon mesh sieves, with ultrasonic bath steps to break up the agglutinated organic matter. The resulting 15–250 µm palynomorph fraction was mounted on glass microscope slides. A general characterization of palynofacies, palynomorph categories (pollen, spores, and aquatic algae), and dinocysts in particular was performed with light microscopy (typically 630× magnification or 1000× magnification, when required), counting at least 200 and typically > 250 elements. For taxonomy, we refer to that cited in Williams et al. (2017), with the exception of wetzelielloid taxa (see Bijl et al., 2017). For paleoecological reconstructions, dinocyst taxa are grouped in complexes of morphologically related species and genera with assumed similar paleoecological affinities, following the empirically based complexes of Frieling and Sluijs (2018).

2.2.6 Clay mineral analysis

Clay mineral analysis was performed in the Earth and Planetary Sciences Department at the University of California, Santa Cruz. Methods were adapted from Kemp et al. (2016), Gibson et al. (2000), and Poppe et al. (2001). Samples were gently crushed and placed in a 1 % sodium hexametaphosphate (Calgon) solution adjusted to a pH of 7–7.5 via the addition of ammonium hydroxide. Samples were disaggregated by placing them on a shaker table at approximately 200 rpm (revolutions per minute) for 48 h. Samples were then wet-washed through a 63 µm sieve to separate the coarse and fine fraction. Samples were then marked 5 cm from the surface of the solution and shaken to resuspend the clay minerals. The suspension was allowed to settle for 4:06 h at 20 ^∘C, in accordance with Stoke's law. The suspension above the 5 cm mark was then extracted via syringe to isolate the < 2 µm size fraction. This suspension was then dried, and approximately 150 mg of this clay-sized fraction was then resuspended in deionized water. The suspension was then passed through a submicron filter attached to a vacuum apparatus to remove the water and orient the clays. The filter paper was removed and placed on the side of a glass beaker before being transferred to a glass slide to produce an oriented mount. For the glycolation step, mounts were placed inside a desiccation chamber with ethylene glycol at 60 ^∘C for a minimum of 4 h.

Oriented mounts were scanned in a Philips 3040/60 X'Pert Pro X-ray diffraction unit from 0–35^∘ 2θ, with a 1^∘ source beam and a $\mathrm{1}/\mathrm{16}$^∘ receiving slit at 45 mA and 40 mV. The results were then printed, and individual sample peaks were weighed to provide a semi-quantitative measurement of clay mineral abundances.

2.2.7 Foraminifer oxygen isotope paleothermometry

Temperatures based on benthic foraminiferal oxygen isotopes ratios were estimated, following the protocol outlined in Hollis et al. (2019). The paleotemperature equation of Marchitto et al. (2014) was used. Then, four to six Anomalinoides acutus specimens were picked from the 180–250 µm size fraction at each interval. All samples were analyzed with a Kiel IV carbonate device paired with a Thermo Scientific MAT253 mass spectrometer. The analytical precision of this system, based on replicate analyses of Carrara marble standards, is better than 0.05 ‰ and 0.10 ‰ for δ¹³C and δ¹⁸O, respectively. In order to calculate an absolute temperature, bottom δ¹⁸O_SW was estimated based on the relationship between the salinity and δ¹⁸O established by Fairbanks (1982) and the modern-day salinity measurements reported in Richaud et al. (2016), with a −1.2 ‰ correction applied to account for the lack of ice during the early Eocene.

2.2.8 TEX₈₆ biomarker paleothermometry

TEX₈₆ analysis was performed on 18 samples at Utrecht University. Glycerol dialkyl glycerol tetraether (GDGT) lipids were extracted from ∼ 10 g of powdered, freeze-dried sediment in 25 mL solvent mixture of dicholormethane (DCM) : methanol (MeOH) (9:1; $v/v$) by a Milestone ETHOS X microwave extraction system, which was set to 70 ^∘C for 50 min. Filtered (using NaSO₄ column) total lipid extracts were dried and separated into apolar, neutral, and polar fractions through AlOx column chromatography, with hexane $/$ DCM (9:1), hexane $/$ DCM (1:1), and 1:1 DCM $/$ MeOH (1:1), respectively, as mobile phases,; thereafter, the extracts were dried again and weighed. Polar fractions, after the addition of a GDGT standard (99 ng of $m/z$ 744) for quantitative analyses, were diluted in hexane $/$ isopropanol (99:1) to a concentration of 2 mg mL⁻¹. A 10 µl aliquot was filtered (0.45 µm polytetrafluoroethylene) and analyzed by high-performance liquid chromatography (HPLC) coupled to an ionization mass spectrometer. The analytical precision was 0.006 TEX₈₆ units, determined using an in-house GDGT standard. We use several indices to constrain potential confounding effects, such as GDGT contributions from methanogenic/methanotrophic microorganisms, terrestrial sources, and deepwater archaea. TEX₈₆ values were calculated from isoprenoid GDGT abundances, following Schouten et al. (2002), and converted to sea surface temperature (SST) using linear and exponential calibrations (e.g., Kim et al., 2010; O'Brien et al., 2017), following recommendations by Hollis et al. (2019). Calibration errors are on the order of ±2.5 ^∘C (1 SD), but any offset from this calibration is assumed to be a constant throughout the study section.

3 Results

3.1 Stratigraphic framework

A distinct magnetic reversal pattern is preserved in the Knapps Narrows core, with the section below 99.21 m (sample KN37) exhibiting a reversed polarity and the upper section (above 91 m, sample KN36) a normal polarity. The reversal boundary, identified as the transition from magnetochron C24r to C24n based on nannoplankton stratigraphic data, thus lies within 99.21–91 m (Fig. 2). The lowermost section is placed in zone NP10, based on the first occurrence (FO) of Rhomboaster spp. and the absence of Tribrachiatus contortus, which has its FO at 96.6 m. The top of zone NP10 is placed at 95.2 m on the last occurrence of T. contortus, and the base of zone NP11 is identified at 95.1 m, based on the FO of T. orthostylus. This reconstructed bio-magnetostratigraphy (Figs. A2 and A3) is generally consistent with the Paleogene timescale 2020 (Speijer et al., 2020) and the timing of ETM2 (Cramer et al., 2003; Lourens et al., 2005; Westerhold et al., 2007; Stassen et al., 2012a).

Figure 2 Knapps Narrows core site (from left to right) depth, nannofossil zones, magnetochrons, lithology, carbon isotopes, percent carbonate, clay mineralogy, and benthic δ¹⁸O temperatures, as well as the TEX₈₆ sea surface temperatures over the ETM2 and H2 interval at Knapps Narrows. TEX₈₆ temperatures use linear calibration from O'Brien et al. (2017). The error in δ¹⁸O represents 1 PSU change in salinity. Note the logarithmic scale for percent carbonate and low-carbonate interval at ETM2, spike in the illite–smectite ratio at ETM2, agreement between benthic and sea surface temperatures, and offsets in timing between carbon isotopes, clay mineralogy, and temperature proxies.

The bulk carbonate δ¹³C record captures a −1.5 ‰ excursion from 98.18–96.41 m before recovering to the baseline, followed by a second −2.5 ‰ excursion between 95.37–93.46 m that peaks at 94.48 m (Fig. 3). This double-peak feature facilitates the recognition of ETM2, as it is documented for several marine and continental sections (e.g., Walvis Ridge and Bighorn Basin in Fig. 3; Cramer et al., 2003; Stap et al., 2010; Abels et al., 2016). Between 93.46 and 88.98 m, the carbon isotope ratios return to near-baseline values before another −1.5 ‰ excursion at 88.51 m, consistent with observations of H2 at other sites (Fig. 3; Cramer et al., 2003; Stap et al., 2010; Abels et al., 2016). Benthic foraminiferal carbon isotope trends largely track those in bulk carbonate. An interval with low-carbonate content (∼ 93–94.5 m; Fig. 2) coincides with the CIE minima at 94.48 m, where the sediments were barren of foraminifera.

Figure 3 Carbonate δ¹³C records of ETM2 and H2 at Walvis Ridge (bulk), Bighorn Basin (pedogenic nodules), Shatsky Rise (benthic), and Knapps Narrows (this study). For Knapps Narrows, bulk carbonate values are in gray, and benthic values are in black. Knapps Narrows at 94.5 m was nearly barren of foraminifera. Walvis Ridge data are from Stap et al. (2010), the age model is from Lauretano et al. (2015), Bighorn Basin data are from Abels et al. (2016), and Shatsky Rise data are from Westerhold et al. (2018).

Figure 4 COCO (a–c) and eCOCO (d–f)analysis of the wireline gamma ray log (42.5–103.5 m). The studied interval lies between 84 and 102 m. In the COCO plots, three clusters of high-ρ sedimentation rates, accompanied with a low-H₀ significance level, and relatively more (6–7 parameters) contributing astronomical parameters are identified, i.e., at ∼ 1, ∼ 5, and ∼ 13 cm kyr⁻¹. For the eCOCO analysis, the sliding window is 18 m, using 0.2 m steps. The target astronomical solution is La2004 at 53 Ma (with seven astronomical cycles at 405, 125, 95, 39.2, 23.1, 21.9, and 18.8 kyr, respectively). All periodograms were analyzed, with the autoregressive (AR1) red noise model removed. The number of Monte Carlo simulations is 2000. The tested sedimentation rates range from 0 to 15 cm kyr⁻¹, with a time step of 0.1 cm kyr⁻¹. Note that from depths of ∼ 82–90 m, which cover the studied interval, a change in the sedimentation rate exists from ∼ 5 to ∼ 13 cm kyr⁻¹. The evolutionary power spectrum (Fig. A5) of this interval also indicates a distortion of cyclicity.

The COCO and eCOCO (Fig. 4) analysis shows a sedimentation rate of ∼ 5 cm kyr⁻¹ for most of the succession from 103.5 to 43.5 m, except for the interval of 90 to ∼ 82 m, where there is a transition to a sedimentation rate of approximately 13 cm kyr⁻¹. Assuming orbital forcing, the ∼ 5 cm kyr⁻¹ and a cycle of ∼ 18 m (Fig. A5) correlates to the 405 kyr long eccentricity period. Given the surge in estimated sedimentation rate, the observed ∼ 18 m cycle should instead correlate to the short eccentricity periodicity. Given the sudden increase in deposition rate, this interval may reflect on other factors; for example, there could be a period of enhanced precipitation, even though this period of elevated sedimentation rates does not correspond to any hyperthermal event, or bathymetric changes related to regional faulting.

3.2 Palynology

Palynological associations are dominated by marine palynomorphs, notably dinocysts, and are common to abundant organic linings of benthic foraminifera. Terrestrial palynomorphs such as pollen and spores are rare. Dinocyst assemblages are well preserved and diverse. They comprise a typical mid–early Eocene shelf assemblage, given the abundant contributions of Spiniferites spp., Operculodinium spp., and species grouped in the Areoligera and Cordosphaeridium complexes (e.g., Pross and Brinkhuis, 2005; Frieling and Sluijs, 2018). The distribution pattern of various chronostratigraphic index taxa such as Biconodinium longissimum and several wetzelielloids (e.g., W. articulata and W. tabulatum) further confirm the assignment to ETM2 and H2 (e.g., Williams et al.,2004; Bijl, 2022). The Senegalinium complex, particularly abundant in the lower part of the section, tolerated low salinities, which is suggestive of strong seasonal river runoff and is in apparent contradiction to the absence of terrestrial pollen and spores (Brinkhuis et al., 2006; Sluijs and Brinkhuis, 2009).

3.3 Clay mineralogy

Clay mineralogy shows moderate variation throughout the studied interval. The clay size fraction is dominated by smectite (> 50 % and up to 95 %) over much of the studied interval, with the remainder being composed largely of illite. There is a marked increase in illite during the body of ETM2, which peaks at a 1.5 illite $/$ smectite ratio, and an additional, smaller, increase in the illite $/$ smectite ratio immediately prior to H2 (Fig. 2). There is no detectable kaolinite within the studied interval. These variations are not correlated with other proxies, and the timing of the mineralogic changes is offset relative to the CIEs of ETM2 and H2; i.e., the shift in mineralogy begins prior to the CIE of ETM2 and after the CIE of H2.

3.4 Oxygen isotope paleothermometry

There is a long-term rise in temperature over the studied interval from 11.1 ^∘C ($-\mathrm{2.8}/+\mathrm{2.9}$) at the base of the core to 13.4 ^∘C ($-\mathrm{2.9}/+\mathrm{3.0}$) at the top (Fig. 2). The largest source of error associated with the δ¹⁸O temperature measurements is the uncertainty associated with local δ¹⁸O_SW. According to the Fairbanks (1982) equation, a change in salinity of ±1 PSU, when propagated across calculations, can result in an error in the temperature on the order of ±2.9 ^∘C, depending on the local freshwater δ¹⁸O. This error cannot be assumed to be constant throughout the core, as hydroclimate is expected to change dramatically during this time period, which could impact local salinity and often seems to be the case with the PETM. Agreement in temperature trends between δ¹⁸O and TEX₈₆ suggests that this effect is minimal. Bottom water temperatures during ETM2 were measured to be 11.6 ^∘C ($+\mathrm{2.9}/-\mathrm{2.8}$), which is 0.9 ^∘C warmer than the baseline value. The greatest warming is seen following the CIE associated with H2 at 14.4 ^∘C ($+\mathrm{3.0}/-\mathrm{2.9}$), which is 2.3 ^∘C warmer than the baseline. Peak warming seems decoupled from minimum δ¹³C values for both ETM2 and H2 (Fig. 2). These temperature variations are likely the minimum estimates given the resolution of the record.

3.5 TEX₈₆ biomarker paleothermometry

SST values based on TEX₈₆ vary by a few degrees, largely mirroring the benthic bottom water temperature trends (Fig. 2). TEX₈₆ records for this site indicate a warming on the order of 3–5 ^∘C over the studied interval. Temperatures appear to increase prior to the body of ETM2 during the pre-onset excursion. Using a linear calibration, SST reached a maximum of ∼ 30 ^∘C, which is ∼ 1.1 ^∘C higher than the baseline. Exponential calibrations result in slightly lower SSTs and slightly muted warming (see the “Data availability” for more information). TEX₈₆ reaches minimum values during the CIE associated with ETM2 at ∼ 28 ^∘C, which 1.0 ^∘C colder than baseline values. The trend reverses and reaches a maximum, following the CIE associated with H2 at ∼ 31^∘C, which is ∼ 1.3 ^∘C warmer than the baseline value.

4 Discussion

4.1 Stratigraphy and sedimentation

Documented sections spanning the ETM2 in marine shelf settings are rare compared to those spanning the PETM. The discovery of the lithologic expression of two smaller hyperthermals, ETM2 and H2, in addition to the PETM, on the mid-Atlantic coast provides a unique opportunity to assess how regional climate responded to two distinct levels of greenhouse forcing. Apart from bio-magnetostratigraphy, the key to identifying ETM2 is the stable carbon isotope stratigraphy, specifically the presence of the paired CIEs of ETM2 and H2, which, at Knapps Narrows, is broadly consistent with trends observed elsewhere (Fig. 1). The magnitude of the CIE here is intermediate compared to records of the ETM2 in open-marine and continental settings but is consistent with observations of the PETM CIE in marginal marine settings (John et al., 2008; Tipple et al., 2011; Sluijs and Dickens, 2012). The marine carbonates of ETM2 document a smaller excursion than those observed in continental sections and in marine organic carbon (Tipple et al., 2011; Stap et al., 2010; Westerhold et al., 2018; Abels et al., 2016; Sluijs and Dickens, 2012). The dampening of the magnitude of change in marine carbonates has been attributed to dissolution and carbonate chemistry changes during the peak of the hyperthermal (e.g., PETM). The enhanced terrestrial signal has been attributed to changes in soil conditions, humidity, plant communities and physiology, and an increase in the organic ¹³C fraction (Bowen et al., 2004; Tipple et al., 2011; Sluijs and Dickens, 2012). The relative lack of dissolution in shallower marine environments and lack of controls that enhance the terrestrial signal likely contribute to this more intermediate CIE.

Although the interval with low-carbonate content prevented the generation of a more detailed foraminiferal isotope record, the covariance between bulk carbonate and benthic foraminiferal δ¹³C throughout, combined with the biostratigraphy and magnetostratigraphy, confirms our interpretation of the carbon isotope stratigraphy. The onset of the PETM in the mid-Atlantic region is characterized by a clay or low-carbonate layer, which has been attributed to both coastal acidification and enhanced siliciclastic sediment fluxes that are climatically driven (Bralower et al., 2018; Gibson et al., 2000; Kopp et al., 2009). High-resolution climate simulations are consistent with enhanced seasonal or peak precipitation, which likely increased the rates of erosion and sedimentation (Carmichael et al., 2017; Rush et al., 2021). Similarly, there are two potential explanations for the low-carbonate interval associated with the body of ETM2, namely coastal acidification and enhanced siliciclastic flux. Coastal acidification would inhibit the calcification and potentially dissolve carbonate directly, leading to the low-carbonate interval. In arguing for increased acidification, Bralower et al. (2018) put forth two possibilities, namely a possible shoaling of the lysocline to the middle shelf and a regional shoaling of the lysocline driven by enhanced eutrophication and microbial activity. However, the increase in pCO₂ levels during ETM2 was much lower than that of the PETM (e.g., Harper et al., 2020), and there is little evidence for a global shoaling of the lysocline to this degree during ETM2. Therefore, it is unlikely that elevated CO₂ levels alone could drive dissolution on the shelf in this region, given the rates of estimated C emissions. It is possible that dissolution could be enhanced by increased oxidation of organic matter for both in situ organic matter and organic matter transported into the shelf (Bralower et al., 2018; Lyons et al., 2019). However, the dinoflagellates in our study interval do not suggest significant changes in productivity in response to the CIEs. Enhanced siliciclastic transport may have contributed to the low-carbonate interval as well, even though this increase appears to be to a much lesser extent relative to the PETM with an estimated 50 %–100 % increase in supply rates within the studied interval, as compared to estimates of 2.8 to 220-fold increases within the PETM. Again, the inflection occurs at ∼ 92 m whereupon the sedimentation rate increase is not correlated to either hyperthermal (Stassen et al., 2012b), although recent studies have suggested much more moderate increases in sedimentation rates across the PETM, which is in good agreement with observations here (Li et al., 2022). Lacking a clear indication of a single driver for the low-carbonate interval, it is likely that a confluence of factors drove the observed changes.

Figure 5 Knapps Narrows core site (from left to right) depth, nannofossil zones, magnetochrons, lithology, carbon isotopes, foraminifera lining counts, relative abundance of wetzelielloids, relative abundance of Goniodomids, and relative abundance of the Senegalinium complex. Foraminifera linings can be used to infer periods of rapid or low siliciclastic influx. Wetzelielloids are associated with periods of high productivity, Goniodomids are associated with high temperatures and strong salinity stratification, and the Senegalinium complex is associated with low salinity (Frieling and Sluijs, 2018). Of note is a long-term increase in salinity from relatively fresh conditions at the base of the core to more saline conditions upsection. Rapid changes in the abundance of wetzelielloids suggest disconformities, periods of non-deposition, or condensed sections.

Figure 6 Sedimentation rates in the Knapps Narrows core based on the correlations of δ¹³C inflection points to the age model of Lauretano et al. (2015) at Walvis Ridge site 391. Sedimentation rates are given in centimeters per thousand years. Note the relative consistency compared to the order of magnitude changes that were documented during the PETM in Stassen et al. (2012b), with moderate increases beginning at 91.8 m, which is consistent with cyclostratigraphy.

Although the coarse resolution of the carbon isotope stratigraphy might not capture a brief increase in sedimentation rates, our age model indicates that the sedimentation rates remained fairly constant across ETM2 and H2, which is in contrast to the lithologic expression of the PETM in this area. Assuming continuous sedimentation, rates range from 4.6 to 10.8 cm kyr⁻¹, with the lower section of the core trending lower and the upper section of the studied interval trending higher, with an inflection point at 91.8 m. The lowest rates occurred prior to ETM2, and the highest occurred during the background state between ETM2 and H2 (Fig. 6). However, sedimentation rates in shelf settings can vary on timescales shorter than the resolution of our age model is able to detect. It is almost certain that the records reflect the amalgamation of periods of enhanced sedimentation and non-deposition (Trampush and Hajek, 2017). Abrupt transitions within the distribution patterns of several dinocyst species suggest minor disconformities, condensed intervals, or depositional hiatuses that would result in highly variable sedimentation rates that we cannot record, based on the resolution of our carbon isotope age model. Therefore, the age model should be seen to be capturing the average over long time periods.

The sedimentation rates predicted based upon the age model are largely consistent with those calculated based on spectral analysis, which found long-term average sedimentation rates on the order of 5 cm kyr⁻¹ in the lower section of the core and 13 cm kyr⁻¹ in the upper section, with an inflection point at the same location as predicted by the age model at ∼ 92 m (Fig. 4). Based on spectral characteristics, the lower section of the core exhibits cyclicity that could be orbitally paced, while the upper section does not. This is unsurprising, given the generally discontinuous or episodic nature of sediment accumulation in such facies. Despite the aforementioned issues surrounding the stochastic nature of shelf deposition, Cyclostratigraphic methods have been used in calculating shelf sedimentation rates for the PETM in the region, and the general agreement between the carbon isotope record and the spectral analysis suggests the general interpretation of a bimodal sedimentation pattern is accurate (Li et al., 2022).

Other factors may have a role in observed sedimentation rates, such as the paleobathymetry and local sea level. Faulting within the region could have led to the development of horst and graben complexes (Fig. A1). These positive and negative structural features could produce localized paleobathymetric highs and lows that could have influenced the observed variation in local sedimentation rates (Sabat, 1977; Self-Trail et al., 2017a).

4.2 Regional hydrology and depositional setting

Evidence supporting a modest intensification of the hydrologic cycle includes the shift in clay mineralogy from smectite-dominated facies to illite-dominated facies during the body of ETM2 and a moderate increase in illite prior to H2, with no detectable levels of kaolinite. Clay mineralogy is used extensively as a proxy for chemical weathering, with a spectrum from smectite and illite to kaolinite being associated with increasingly wet and warm conditions (Singer, 1980). Many locations experience a dramatic increase in kaolinite during the PETM (Robert and Kennett, 1994; Gibson et al., 2000; Kemp et al., 2016; Chen et al., 2016). Previous studies on clay mineralogy over the PETM in the Salisbury Embayment demonstrated that kaolinite increased and that this increase could have resulted from the reworking of older deposits (Gibson, 2000; John et al., 2012). In other regions, the increase in kaolinite is believed to result from an authigenic weathering response to the PETM (Clechenko et al., 2007; Chen et al., 2016). Work on the radioisotopes of the clay minerals, namely strontium and lead, demonstrate little change in the sourcing of sediment in the Atlantic Coastal Plain during the PETM, and work on lithium isotopes seems to suggest rapid changes to global weathering patterns during the PETM (Rush et al., unpublished data; Pogge von Strandmann et al., 2021; Ramos et al., 2022). Since illite is a product of weathering between smectite and kaolinite, alongside recent studies demonstrating an immediate weathering response during the PETM, its presence seems to support the interpretation that the clays formed in situ during ETM2. As the warming of ETM2 is intermediate to that of the background state and the PETM, the environmental perturbations during ETM2 would likely illicit a weathering response intermediate to the background state and that of the PETM.

The palynological assemblages imply a complicated story. The assemblages are typical of those seen in neritic settings (Pross and Brinkhuis, 2005; Frieling and Sluijs, 2018). They show a large degree of variability in freshwater runoff, stratification, and productivity, but the majority of these changes show no relation to δ¹³C or temperature change, which is atypical for dinocysts (Fig. 5). Wetzelielloids, which include the Apectodinium complex and Apectodinium spp., show (quasi-)global abundance peaks across the PETM but demonstrate a small relative increase prior to ETM2, with their highest relative abundances occurring several meters above both ETM2 and H2 (Crouch et al., 2003; Frieling and Sluijs, 2018; Fig. 5). The only anomaly that clearly corresponds to a hyperthermal is the peak abundance of Goniodomideae during H2 (Fig. 5). Goniodomideae represent a thermophile group thriving in lagoonal environments (variable salinity) and further offshore under strongly stratified conditions (e.g., Reichart and Brinkhuis, 2003; Frieling and Sluijs, 2018). Taking into account the decreasing abundance of the Senegalinium complex, this likely points to an increase in salinization during this interval, due to decreased runoff or increased evaporation. This interval also yields anomalous abundances of benthic foraminifer linings, which generally indicate a temporary reduction in sedimentation, further supporting the interpretation of decreased runoff during this interval.

An overall shallowing trend in the upper section of the core may be inferred, based on a long-term decline in the abundance of the low-salinity-tolerant Senegalinium complex and a rise of the Goniodomids, which typically occur in massive (harmful) blooms in (seasonally) high-salinity, warm, lagoonal conditions in modern and ancient settings (Zonneveld et al., 2013; Frieling and Sluijs, 2018; Sluijs et al., 2018). This shallowing trend also indicates a shift from brackish conditions to more saline conditions upsection that do not necessarily correspond to hyperthermal events and may indicate an overall longer-term trend.

Increased abundances of the low-salinity-tolerant Senegalinium complex and the supply of terrestrial palynomorphs suggest increased river runoff in shelf sections across hyperthermal events (e.g., Crouch et al., 2003; Sluijs and Brinkhuis, 2009; Sluijs et al., 2009). These abundances in Knapps Narrows show no correspondence to the CIEs or temperature anomalies. The complex was abundant before, during, and after ETM2 but was not present during H2 after its gradual decline, which does not support a strong hydrological response to this event. Terrestrial palynomorphs are rare throughout the record, reminiscent of the upper Paleocene and lower Eocene of the New Jersey Shelf (e.g., Zachos et al., 2006). This suggests changes in transport of terrestrial palynomorphs or that the hinterland was sparsely vegetated in the region, also explaining the lack of response during the PETM, ETM2, and H2 (Self-Trail et al., 2017b). This may be consistent with the paleolatitude of about 33^∘ N (http://www.paleolatitude.org/, version 2.1, last access: 10 August 2023; van Hinsbergen et al., 2015), which arguably would represent an arid zone. Quantitative reconstructions of salinity could aid in this interpretation. However, at present, we are unable to reconstruct surface salinities due to a paucity of surface-dwelling foraminifera species within the studied interval and the qualitative nature of dinoflagellate-based paleoenvironmental reconstructions.

Goniodomids reach a maximum at H2 that also corresponds to the highest reconstructed temperatures. A rise in salinity seems to precede this warming, but it should be noted that dinoflagellate blooms are strongly seasonal and therefore should not necessarily correspond to our salinity estimates. The correspondence of Goniodomid abundances with transient global warming events has been previously recorded, for example, in brief intervals of the PETM on the New Jersey Shelf, the Nigerian Shelf, and across the middle Eocene climatic optimum in the eastern equatorial Atlantic (Sluijs and Brinkhuis, 2009; Frieling et al., 2017; Cramwinckel et al., 2019). Very high stable carbon isotope ratios suggest that these Goniodomid abundances indicate that at least some of these abundances represent an increase in the geographical range and frequency of harmful blooms (Frieling et al., 2017; Sluijs et al., 2018).

In contrast to the records here, the Lomonosov Ridge palynological assemblages demonstrate a marked shift during ETM2, with an increase in dinoflagellate markers that are indicative of freshening and eutrophication and the appearance of palm pollen that indicates high winter temperatures (Sluijs et al., 2009; Willard et al., 2019). Further evidence of hydrologic changes is noted in the hydrogen isotope records of leaf wax n-alkanes, which demonstrate a large rise during the PETM, a rise prior to ETM2, and a drop during the body of the event, as well as during H2 (Pagani et al., 2006; Krishnan et al., 2014).

Looking at other trackers of hydrologic changes across hyperthermals, records of enhanced siliciclastic flux were archived at Mead and Dee streams in Aotearoa / New Zealand across several hyperthermal events (Nicolo et al., 2007; Slotnick et al., 2012). While focused more on biotic and ocean chemistry conditions, the Nile Basin is also interpreted as having a similar response to the PETM but on a smaller scale (Stassen et al., 2012a). Such evidence of enhanced hydrology driving siliciclastic fluxes is also found in hemipelagic sections such as the Terche section of northeastern Italy (D'Onofrio et al., 2016).

The records from continental sections suggest that changes in the regional precipitation patterns were not entirely uniform. An increase in siliciclastic sediment flux in a terrestrial section at Ellesmere Island is interpreted as evidence of an intensified hydrologic system during ETM2 (Reinhardt et al., 2022). The Bighorn Basin record suggests contrasting regional climatic changes between the PETM and later hyperthermal events (Abels et al., 2016), which is in line with observations from this study.

The overall hydrologic response to ETM2 appears much more modest than that of the PETM and somewhat muted, given the scale of global warming and assuming a roughly linear response of precipitation patterns to global warming. Given that the baseline temperature prior to ETM2 is a few degrees higher than prior to the PETM, it is possible that sensitivity to an equivalent temperature perturbation would be less extreme than for the PETM, i.e., a state-dependent response. This also suggests variability in the hydroclimate responses to greenhouse forcing in different regions. While the Arctic records demonstrate a strong hydrological response across ETM2, here the interpretation is less straightforward (Sluijs et al., 2009). Depending on the region, the local hydroclimate seems to respond in a somewhat non-linear fashion to global warming, characterized more by abrupt mode shifts (e.g., frequency of extremes or duration of wet vs. dry seasons), as opposed to gradational changes (e.g., mean annual precipitation; Carmichael et al., 2017; Rush et al., 2021). In the mid-Atlantic region, our palynological record appears to confirm this overall mild response to the ETM2. In effect, to record only minor changes is rather unique, given that other transient climatic events in the Cenozoic are typically characterized by massive concomitant dinocyst assemblage changes and more significant changes in sedimentation.

4.3 Temperature

GDGT distributions indicate a dominant pelagic thaumarchaeotal source based on low abundances of terrestrial GDGTs and of lipids derived from methanotrophs or methanogenetic microorganisms (Fig. A3). GDGT2 $/$ 3 ratios between 1.6 and 2.1 (Fig. A4) indicate a dominant lipid origin from the upper ∼ 150 m of the water column (Hurley et al., 2018; van der Weijst et al., 2022). Although this may imply contributions from sub-thermocline Thaumarchaeota, potentially compromising direct quantitative SST assessment based on the surface sediment calibration dataset, we present the results as SST reconstructions, following various linear and non-linear calibrations following convention (van der Weijst et al., 2022; Hollis et al., 2019). Variability in the TEX₈₆ record should reflect SST variability very well (Ho and Laepple, 2016).

The ETM2 has already been documented in numerous pelagic sections, which uniformly demonstrate 2–4 ^∘C of warming during the event (Harper et al., 2018). The uniformity observed in open-ocean sections is mirrored in previously documented shelf sections. The central Lomonosov Ridge has been argued to have experienced an increase in warm, wet conditions during ETM2, similar to reconstructions for the PETM (Sluijs et al., 2009, 2020).

There is broad correspondence between temperature variability derived from TEX₈₆ and δ¹⁸O, suggesting the trends are accurate. TEX₈₆ records indicate minor variations, with a rise of 1 ^∘C prior to ETM2, a drop of 1 ^∘C relative to the baseline during the body, and an increase of 1.3 ^∘C during H2. The magnitude of these shifts is much lower than those of the PETM at nearby sites, with paleotemperature reconstructions from the PETM derived from TEX₈₆ and δ¹⁸O at other sections in the region suggesting 5–8 ^∘C surface warming without much scatter (Sluijs et al., 2007; Zachos et al., 2006; Babila et al., 2022). The elevated temperatures prior to the CIE associated with ETM2 are followed by a drop in temperatures during the event. There are multiple potential explanations for the apparent lack of correspondence between the CIEs and temperature variability. The first is that we were unable to measure δ¹⁸O during the peak of ETM2 due to the low-carbonate interval, and we may be missing some degree of internal variability within the event. However, the TEX₈₆ record had the lowest temperature recorded within this interval (Fig. 2). In shelf settings, deposition is generally discontinuous and variable, particularly during periods of dramatic climatic changes such as Eocene hyperthermals (Trampush and Hajek, 2017). This has the potential to introduce preservation and sampling biases within the sedimentary record and, given the coarse resolution of the temperature records, the possibility of not capturing a warming anomaly that lasted only a few tens of thousands of years. Another possibility is the contribution of regional-scale changes resulting from changing ocean and atmospheric circulation patterns in a modestly warmer world. Another factor to consider is orbital forcing. The Eocene hyperthermals appear to coincide with 400 kyr eccentricity maxima (Lourens et al., 2005; Zachos et al., 2010). Although the lower-frequency cyclical changes in CO₂ appear to be paced by eccentricity (Zeebe et al., 2017; Vervoort et al., 2021), depending on the phase of precession, the local temperature, particularly on or near land, could deviate from global patterns (Kiehl et al., 2018). This effect would be further exaggerated if local sediment deposition tended to be biased toward one phase. The variation in responses between various proxy systems and muted responses relative to the PETM highlights the non-linearity of the climate system to CO₂ forcing and demonstrates the importance of considering additional forcings.

5 Conclusions

This study provides the first documentation of ETM2 and H2 in the Salisbury Embayment of the mid-Atlantic Coastal Plain. Clay mineralogy, δ¹⁸O, and TEX₈₆ records portray climate perturbations across this interval. The shift in clay mineralogy from smectite to illite, taken in the context of kaolinite-dominated facies during the preceding PETM, suggests that the clay mineralogy changes are authigenic and represent an intermediate stage of weathering relative to the PETM and the background state resulting from an enhanced hydrologic cycle. Temperature proxies show minor increases in temperature that appear asynchronous with changes to the carbon cycle and may represent alternative forcings, such as changes in ocean and atmospheric circulation or orbital forcing or a potential sampling artifact. Dinocyst assemblages demonstrate little response to the short-lived hyperthermals and point more towards being dominated by long-term trends. Similarly, sedimentation rates also appear to be dominated by long-term trends occurring in a bimodal fashion, with the transition between the states having no correlation to either of the hyperthermals. When compared to the PETM and other records of ETM2, these contrasting records demonstrate the complexities of climatic changes and the regional variability associated with rapid warming events, thus highlighting the need for extensive spatial coverage for interpreting climatic changes. In essence, the local changes in precipitation did not scale linearly with the changes in temperature. The relatively muted signal in this region compared to the PETM may be due, in part, to the warmer baseline state prior to ETM2 compared to the PETM. Future research into high-resolution climate models of these smaller-scale hyperthermal events may assist in interpretation of this record and records in different regions.

Appendix A

Figure A1 Regional calcareous nannofossil biostratigraphy. Locations are from Fig. 1. Note that the disconformity at the base of Knapps Narrows has a truncated lower-zone NP10, as indicated by the absence of the Marlboro Clay. The scale shows thickness of the sediments in each core (in meters), as measured from the base of the PETM.

Figure A2 (a) The magnetostratigraphy and calcareous nannofossil zones (nanno-zone) of Knapps Narrows core. Note that the C24r $/$ C24n reversal boundary was not well determined due to a sampling gap but was well constrained to the interval, where the question mark is placed, based on nannofossils. (b) The latest Paleocene to earliest Eocene timescale from GTS2020 (Gradstein et al., 2020), showing the geomagnetic polarity scale, the calcareous nannofossil zones, and the composite δ¹³C isotope curve (major δ¹³C excursions and events are also marked). The excursion of ETM2 aligns with the C24r $/$ C24n reversal.

Figure A3 Ring index (Zhang et al., 2016) plotted against TEX₈₆ values for Knapps Narrows. Values generally fall within the acceptable range, indicating little non-thermal overprint. The temperature reconstructions obtained from the outliers were in line with independent measurements from δ¹⁸O. Moreover, the cutoff for identifying anomalous values is rather conservative (Yi Ge Zhang, written communication, 2022). Therefore, we included these data points in the figures.

Figure A4 Analysis of potential confounding factors for TEX₈₆ paleothermometry for all samples. The red line indicates the cutoff value; red crosses indicate samples exceeding the cutoff (a) TEX₈₆, with the blue line indicating the maximum modern core-top value (∼ 0.72). (b) Low branched and isoprenoid tetraether (BIT) index values indicate the insignificant supply of terrestrial GDGTs. (c) The fractional abundance of crenarchaeol regioisomer relative to total crenarchaeol (fcren) shows no non-thermal contribution of crenarchaeol isomer. (d) Methane index values imply no contribution by methane-metabolizing archaea. (e) The anaerobic oxidation of methane (AOM) ratio values show no contribution by anaerobic methane oxidizers. (f) GDGT2 $/$ 3 values show no contributions of deep-dwelling archaea. (g) Low GDGT0 $/$ cren (crenarchaeol) values show no contribution by methanogenic archaea.

Figure A5 Spectral analysis of the detrended gamma data of the Knapps Narrows core (43.5–103.5 m). From left to right are the raw gamma data (solid black line) with the ca. 21 m LOWESS (locally weighted scatterplot smoothing) trend (solid blue line), the detrended gamma data overlain by the ∼ 18 m filter (Gaussian bandwidth at 0.055±0.025), and a spectral analysis of the detrended series. An interpolation of 0.03 m is applied before the spectral analysis. The evolutionary power spectrum (fast Fourier transform; Kodama and Hinnov, 2015) uses an 18 m window. The 2π multi-taper method (MTM) power spectrum shows that the dominant cycles have wavelengths of ∼ 18, ∼ 5.5, 3, ∼ 1, and 0.6 m, respectively.