Carbon isotope stratigraphy and mammal turnover during post- PETM hyperthermals

Paleogene hyperthermals, including the Paleocene-Eocene Thermal Maximum (PETM) and several other smaller 10 events, represent global perturbations to Earth’s climate system and are characterized by warmer temperatures, shifts in floral and faunal communities, and hydrologic changes. These events are identified in the geologic record globally by negative carbon isotope excursions (CIEs), resulting from the input of isotopically light carbon into Earth’s atmosphere. Much about the causes and effects of hyperthermals remains uncertain, including whether all hyperthermals are caused by the same underlying processes, how biotic effects scale with the magnitude of hyperthermals, and why CIEs are larger in paleosol 15 carbonates relative to marine records. Resolving these questions is crucial for their full interpretation and application to future climate scenarios. The Fifteenmile Creek area of the central Bighorn Basin, Wyoming U.S.A., exposes an early Eocene floodplain sedimentary sequence that preserves paleosol carbonates and an extensive fossil mammal collection. Previous analysis of faunal assemblages revealed two pulses of mammal turnover and changes in diversity interpreted to correlate with the ETM2 and H2 hyperthermals that immediately follow the PETM. This was, however, based on long distance correlation 20 of chemostratigraphic records. We present new carbon isotope stratigraphy using micrite δC values from paleosol carbonate nodules preserved in and between richly fossiliferous localities at Fifteenmile Creek to identify the stratigraphic positions of ETM2 and H2. Additionally, we used differential GPS elevations to establish a new stratigraphic framework that assists in correlation and is independent from the biostratigraphy and previous composite lithostratigraphic sections from the area. Carbon isotope results 25 show that the ETM2 and H2 hyperthermals, and possibly the subsequent I1 hyperthermal, are recorded at Fifteenmile Creek. ETM2 and H2 overlap with the two previously recognized pulses of mammal turnover. Comparisons between the new chemostratigraphy and fossil record suggest that the recorded amplitude of these faunal changes may be muted as a result of some stratigraphic averaging of fossils. The CIEs for these hyperthermals are also smaller in magnitude than in more northerly Bighorn Basin records. We suggest that basin-wide differences in soil moisture and/or vegetation could contribute to variable 30 CIE amplitudes in this and other terrestrial records. https://doi.org/10.5194/cp-2021-83 Preprint. Discussion started: 13 July 2021 c © Author(s) 2021. CC BY 4.0 License.

The negative carbon isotope excursions (CIEs) associated with the hyperthermals indicate massive injections of isotopically light carbon into the atmosphere, but the source and mechanism of release of these greenhouse gases is unresolved (Higgins 50 and Schrag, 2006). Proposed sources of carbon for the PETM include destabilized methane clathrates (Dickens et al., 1995(Dickens et al., , 1997Katz et al., 1999), volcanism associated with the North Atlantic Igneous Province (NAIP) (Gutjahr et al., 2017), thermogenic methane released from organic-rich sediments during NAIP emplacement (Svensen et al., 2004;Frieling et al., 2016), Antarctic permafrost thaw (DeConto et al., 2010(DeConto et al., , 2012, wildfires burning Paleocene peat deposits (Kurtz et al., 2003;Moore and Kurtz, 2008), and evaporation of epicontinental seas leading to the oxidation of organic matter (Higgins and Schrag, 55 2006), although geological evidence supporting the last two hypotheses is weak. It is also possible that multiple mechanisms and sources may have acted together and that the drivers may not have been the same for the PETM and smaller hyperthermals.
Although the initial trigger for many hyperthermals may be related to orbital cyclicity, the primary cause of warming is likely the associated release of greenhouse gases. Of the smaller hyperthermals, those immediately following the PETM are currently the best studied. These include the Eocene Thermal Maximum 2 (ETM2) at ~54 Ma, also referred to as ELMO (Lourens et al., 60 2005) or H1 (Cramer et al., 2003), as well as the succeeding H2, I1, and I2. Like the PETM, these hyperthermals are associated with negative CIEs, increased CaCO3 dissolution, and foraminiferal turnover (Agnini et al., 2009;Stap et al., 2009Stap et al., , 2010Jennions et al., 2015;Arreguin-Rodriguez and Alegret, 2016;D'Onofrio et al., 2016). The ETM2, H2, and I1 hyperthermals have also been linked to increased radiolarian abundance as well as shifts in calcareous plankton assemblages. Together these suggest upper water column warming, weakening thermal stratification in the upper water column, and increased nutrients in surface waters during the events (D'Onofrio et al., 2016). The high-resolution marine record of the PETM and subsequent hyperthermals has led to a growing understanding of their effects in the ocean. A similar level of detail has not been developed for the terrestrial record of the smaller (non-PETM) hyperthermals. Achieving this level of understanding will require geochemical records from a variety of depositional environments and locations to fully understand spatial heterogeneity and underlying carbon cycle-climate dynamics. Furthermore, understanding the effects of these hyperthermals on the terrestrial 70 ecosystem requires their recognition in richly fossiliferous strata, where both plant and vertebrate fossils occur in abundance.

Terrestrial record of hyperthermals
Terrestrial stratigraphic records of the PETM have been reported from North America (e.g., Koch et al., 1992Koch et al., , 2003Bowen and Bowen, 2008;Baczynski et al., 2013), South America (e.g., Jaramillo et al., 2010), Europe (e.g., Cojan et al., 2000;Schmitz and Pujalte, 2007), Asia (e.g., Bowen et al., 2002;Chen et al., 2014), India (e.g., Samanta et al., 2013), and Australia (e.g., 75 Greenwood et al., 2003). Of these, the Bighorn Basin of Wyoming preserves the most detailed, richly fossiliferous, and best studied terrestrial record of the PETM. Here the PETM corresponds with rapid turnover in both plants and mammals, as the Paleocene flora is replaced by a dry tropical flora (Wing et al., 2005) and the first representatives of several mammalian clades disperse among the Holarctic continents. The "dwarfing" of about 40% of mammalian lineages also occurred during the PETM (Gingerich, 1989;Clyde and Gingerich, 1998;Secord et al., 2012). The hydrologic cycle was also strongly affected (Wing et 80 al., 2005;Foreman et al., 2012;Kraus et al., 2013;Foreman, 2014;Baczynski et al., 2017). A comprehensive summary of terrestrial environmental changes associated with the PETM can be found in McInerney and Wing (2011). In contrast to the PETM, terrestrial records of the post-PETM hyperthermals are limited. They include carbon isotope records from coal exposures in NE China (Chen et al., 2014) and western India (Clementz et al., 2011;Samanta et al., 2013;Agrawal et al., 2017), as well as records from paleosol carbonates preserved in floodplain deposits from the McCullough Peaks region in the 85 central Bighorn Basin (Abels et al., 2012(Abels et al., , 2016D'Ambrosia et al., 2017). Four post-PETM hyperthermals (ETM2, H2, I1 and I2) are preserved in the McCullough Peaks sequence and this area offers one of the best opportunities for constructing a highly detailed terrestrial record, that can also be compared with marine and terrestrial records of the PETM and post-PETM hyperthermals worldwide. The McCullough Peaks fauna, however, is limited compared to other parts of the Bighorn Basin, complicating direct correlation between hyperthermals and faunal change. 90 Records of the PETM show that terrestrial proxies generally record a larger magnitude CIE than marine proxies, leading some to consider the terrestrial record "amplified" (Bowen et al., 2004). In marine carbonate records of the PETM, the CIE generally ranges from ~2-4 ‰ (McInerney and Wing, 2011). In terrestrial records, the PETM CIE ranges from ~3-7 ‰ depending upon the proxy used and can show considerable spatial variability depending on local environmental differences (Bowen and Bowen, 2008). Of the terrestrial δ 13 C proxies, paleosol carbonates record the largest CIE, generally ranging from ~3-7 ‰ (e.g., Koch 95 et al., 2003;Bowen et al., 2004;Bowen and Bowen, 2008), with a mean of ~5.5 ‰ (McInerney and Wing, 2011 Detailed paleosol carbonate records of the post-PETM hyperthermals are currently limited to McCullough Peaks, but also seem to show a larger magnitude CIE relative to the marine record. Marine benthic carbonate records of the ETM2 and H2 100 CIEs are ~1.4-1.5 ‰ and 0.8 ‰, respectively (Stap et al., 2010;Barnet et al., 2019), whereas the terrestrial ETM2 and H2 records from soil carbonates in McCullough Peaks are ~3.8 ‰ and ~2.8 ‰ (Abels et al., 2012(Abels et al., , 2016D'Ambrosia et al., 2017).
No terrestrial organic carbon records exist for these hyperthermals from the Bighorn Basin, although bulk organic carbon records from India and China have an average magnitude of ~2.7-3.5 ‰ for ETM2 and ~2.5 ‰ for H2 (Chen et al., 2014;Agrawal et al., 2017). Moreover, there is debate surrounding the scaling relationship between the amplitude of the PETM and 105 post-PETM CIEs, with some arguing for a consistent linear relationship and therefore common source (e.g., Chen et al., 2014) and others arguing that the post-PETM CIEs appear to scale linearly with each other, but not with the PETM, potentially suggesting a different triggering mechanism for the later hyperthermals (Abels et al., 2016).
The cause of this apparent "amplification" in terrestrial CIEs is not fully understood, although increased soil productivity (Bowen, 2013), increased humidity (Bowen et al., 2004;Bowen and Bowen, 2008), increased carbon isotope discrimination 110 by plants under higher atmospheric pCO2 conditions (Schubert and Jahren, 2013), and the loss of gymnosperms in the PETM (Smith et al., 2007), which have higher δ 13 C values than angiosperms, have been invoked as potential mechanisms.
Alternatively, increased dissolution during hyperthermals could produce a muted CIE signal in some marine records relative to the true atmospheric δ 13 C shifts (Zachos et al., 2005). The paucity of terrestrial records of post-PETM hyperthermals outside of the McCullough Peaks area makes it difficult to understand how these records vary spatially and to compare CIE magnitudes 115 between the terrestrial and marine realms.

Mammal response to hyperthermals
The PETM coincides with one of the most profound phases of mammal turnover during the Cenozoic, marked by the first appearances of several clades including: Perissodactyla, Artiodactyla, Euprimates, and Hyaenodontidae (e.g., Rose, 1981;Gingerich, 1989Gingerich, , 2006Koch et al., 1992;Clyde and Gingerich, 1998;Hooker, 1998;Rose et al., 2012). These taxa 120 disperse among North America, Europe, and Asia during the PETM, but their areas of origin are uncertain (e.g., Bowen et al., 2002;Clyde et al., 2003;Smith et al., 2006;Morse et al., 2019). Several groups also demonstrate transient "dwarfing" during the PETM (Gingerich, 1989;Clyde and Gingerich, 1998;Strait, 2001;Chester et al., 2010), with minimum body sizes corresponding to peak warming in Bighorn Basin equids (Secord et al., 2012), the only clade described in stratigraphic detail through the PETM. 125 Previous work on mammal faunas from post-PETM hyperthermals in McCullough Peaks suggests that mammal body size decreased during ETM2 in at least two lineages, though the fauna did not show major reorganization as is seen during the PETM (Abels et al., 2012;D'Ambrosia et al., 2017). However, PETM reorganization may largely be related to intercontinental and intracontinental immigrants making first appearances, while there is no strong evidence for such immigrants appearing around the time of ETM2 or H2 (Woodburne et al., 2009;Chew, 2015). Until now, the lack of abundant fossil mammals tied 130 closely to the stable isotope record has precluded a detailed assessment of mammal faunal change through these hyperthermals.
https://doi.org/10.5194/cp-2021-83 Preprint. Discussion started: 13 July 2021 c Author(s) 2021. CC BY 4.0 License. Chew (2015) identified two distinct pulses of increased turnover and diversity, known as faunal events B-1 and B-2, in mammal assemblages from the Fifteenmile Creek area in the southern Bighorn Basin. In contrast to the PETM, faunal change during these events was lower in magnitude and driven by increased beta, rather than alpha, diversity. Moreover, the number of first appearances was roughly equal to the number of last appearances and there were no known intercontinental immigrant 135 taxa arriving during the B-1 and B-2 events. Additionally, abundance shifts during these events appear to have favored smaller bodied species (Chew, 2015). Chew (2015) identified the ETM2 and H2 hyperthermals as potential drivers of faunal change during events B-1 and B-2, based on their inferred temporal overlap. Preliminary geochemical work in this region, however, did not identify any clear CIEs in this interval that would confirm this connection (Koch et al., 2003).
Our work presents new, high-resolution carbon isotope records from paleosol carbonates through the fossil-rich Fifteenmile 140 Creek area in the southern Bighorn Basin. This allows, for the first time, the extensive fossil mammal records from this area to be tied directly into the isotope stratigraphy of post-PETM early Eocene hyperthermals. The results offer a direct test of the hypothesis that faunal turnover events B-1 and B-2 are correlated with the ETM2 and H2 hyperthermals (Chew, 2015). This new record also adds to the terrestrial stable isotope record of post-PETM hyperthermals in the Bighorn Basin, allowing for better spatial characterization of these events. 145

Geological setting
The Bighorn Basin is a Laramide structural basin located in northwest Wyoming. Uplift of the Beartooth, Pryor, Bighorn, and Owl Creek mountains resulted in sediment accumulation in the basin during the early Paleogene (Gingerich, 1983;Kraus, 1992). The early Paleogene basin axis was located near the western margin of the basin and sediment thickness generally decreases going from the northwest to southeast within the basin (Parker and Jones, 1986), with regional differences in 150 subsidence potentially due to activation of east-west trending buried faults (Kraus, 1992). The lower Eocene Willwood Formation represents a fluvial depositional system, and it is typified by paleosols of varying maturity, channel sandstones, crevasse-splay deposits, and localized carbonaceous shales (Kraus, 1992(Kraus, , 1997Bown and Kraus, 1993). Near McCullough Peaks, the basin experienced rapid subsidence during this time and the Willwood Formation is characterized by a thick sequence of relatively immature paleosols. 155 The Fifteenmile Creek area of the Bighorn Basin is located ~80 km southeast of McCullough Peaks and, in contrast to the McCullough Peaks, it is characterized by having more mature paleosols, as well as more channel and cut-and-fill deposits associated with slower subsidence and aggradation. The region has produced an extensive collection of nearly 1000 mammal fossil localities that have been tied in to a 700 m thick composite stratigraphic section (Bown et al., 1994). Previous work in the area suggested that the interval containing the ETM2 and H2 hyperthermals was between 380 and 455 meters in the Bown 160 et al., (1994) composite section (Chew, 2015). This interval follows the last occurrence of the mammal Haplomylus speirianus and the first occurrence of Bunophorus estagicus associated with the faunal event Biohorizon B (Schankler, 1980;Bown et al., 1994;Chew, 2009), which occurs just below ETM2 and H2 in sections farther north (Abels et al., 2012;D'Ambrosia et https://doi.org/10.5194/cp-2021-83 Preprint. Discussion started: 13 July 2021 c Author(s) 2021. CC BY 4.0 License. al., 2017). The Chron C24r -C24n.3n geomagnetic polarity reversal provides additional support for this interval containing ETM2 and H2. In McCullough Peaks and other records around the world, the C24r -C24n.3n reversal occurs just above ETM2 165 (Lourens et al., 2005;Abels et al., 2012;D'Ambrosia et al., 2017), meaning ETM2 is closely bracketed by Biohorizon B below and the Chron C24r -C24n.3n polarity reversal above. Magnetostratigraphic work from Elk Creek Rim, about 20 km north of the Fifteenmile Creek sections studied here, placed the C24r -C24n.3n reversal at the ~450 meter level of the Bown et al. (1994) composite section, ~70 meters above Biohorizon B (Clyde et al., 2007; however see Tauxe et al., 1994 for earlier magnetostratigraphic interpretation). Based on these constraints, we constructed detailed carbon isotope datasets from 170 stratigraphic sections along Fifteenmile Creek that span the critical interval between Biohorizon B and the C24r -C24n.3n reversal to identify ETM2 and H2 and determine their correlation to mammal faunal changes.

Field methods
Pedogenic carbonate nodules were collected from five primary stratigraphic sections that span the target interval ( Fig. 1). All 195 samples were collected from fresh, unweathered rock and at least three nodules were analyzed from each sampling site, when possible. The relative stratigraphic positions of sampling sites and fossil localities were measured in local sections using a Jacob's staff and bed tracing. These sections were correlated to the composite stratigraphic section of Bown et al. (1994) by using the composite stratigraphic levels (hereafter referred to as "Bown composite meter levels" or BCM) of fossil localities that were included in the sections and using the measured thickness to infer BCM levels for the rest of the section. Because of 200 the difficulty with tracing individual beds in this low relief area with spatially dispersed outcrops, we employed a second method of correlation using local elevation to develop an independent stratigraphic framework. Elevations were collected using a Trimble GeoXT 6000 handheld differential GPS (dGPS) paired with a Trimble Tornado external dual-frequency antenna. Elevations were collected from sample sites and from important stratigraphic markers (e.g., marker beds and fossil producing horizons) in every section. Elevation is a reasonable proxy for stratigraphic level in the Fifteenmile Creek area since 205 the dip is close to 0°. Checks on local dips were made by tracing well exposed beds using the dGPS. Postprocessing of dGPS data was done using the Trimble Pathfinder Office and Terra Sync Software. After processing, the elevations typically had an accuracy of < 50 cm.

Mammal fossil localities
A total of 18 fossil localities, that include more than 4500 identified fossil mammal specimens, were directly tied to the five 210 primary local stratigraphic sections from which isotope samples were collected ( Fig. 1) and used to correlate these sections to the BCM levels. Three shorter sections were also sampled, one away from the main Basin Draw section (through fossil locality https://doi.org/10.5194/cp-2021-83 Preprint. Discussion started: 13 July 2021 c Author(s) 2021. CC BY 4.0 License. D-1454), and two located near the Kraus Flats and Red Butte sections (through localities D-1207 and D-1532, respectively) (Appendix A).

Laboratory methods 215
Individual nodules were polished using a diamond lap wheel and were visually inspected to identify any primary, secondary, or altered textures. Micritic calcite and sparite were sampled for stable isotope analysis. Micrite is generally interpreted to represent a primary phase formed by precipitation directly from soil CO2, and thus is appropriate for paleoclimate reconstruction (Bowen et al., 2001). Sparite samples were also drilled from a subset of nodules to compare stable isotope values of primary and secondary components (e.g., Bowen et al., 2001;Snell et al., 2013). Sparite in paleosol carbonate nodules 220 generally forms in veins and septarian-style cracks in the nodules. It is thought to precipitate from fluids after deep burial at high temperatures, resulting in lower δ 18 O values than the original soil water δ 18 O (Bowen et al., 2001;Snell et al., 2013). All stable isotope ratios are reported in delta (δ) notation as per mil (‰) deviation relative to the Vienna Pee Dee Belemnite (VPDB) standard, where δ = ([(RSample / RStandard) -1] x 1000), and R is the ratio of the higher mass isotope to the lower mass isotope. 225 Analyses of δ 13 C and δ 18 O values were done at the University of Michigan Stable Isotope Laboratory (UMSIL) using a Thermo Finnegan MAT253 stable isotope ratio mass spectrometer attached to a Kiel IV automated preparation device, and at the University of Colorado Boulder Earth Systems Stable Isotope Laboratory (CUBES-SIL) using a Thermo Delta V continuous flow stable isotope ratio mass spectrometer attached to a GasBench II gas preparation device. Repeated measurements of inhouse standards yield precision of ± 0.1 ‰ or better for both δ 13 C and δ 18 O in both laboratories, although overall uncertainties 230 are likely slightly larger than this in carbon isotopes due to extrapolation of the standard correction lines beyond the lowest standard values for both labs (Appendix B). Data correction and calculations from CUBES-SIL were done using in-house R scripts that utilized Tidyverse and IsoVerse R packages. Powder from a subset of 10 micrite samples was homogenized and replicates were sent to UMSIL and CUBES-SIL to identify if there were differences between the two labs. δ 13 C and δ 18 O values that were measured at UMSIL were then corrected to the CUBES-SIL values based on these replicate samples 235 CIE magnitudes were calculated by detrending the data to account for long-term early Paleogene trends, then taking the difference between the detrended δ 13 C background values and the values from the body of the CIEs, following the method of Abels et al. (2016). δ 13 C values were averaged from each stratigraphic level prior to calculating the magnitude and the mean peak excursion values from each section were also averaged to facilitate comparisons to other records and account for local 240 spatial heterogeneity.

Stable isotopes
A total of 789 carbonate samples from the Fifteenmile Creek area were analyzed for δ 13 C and δ 18 O (Figs. 2 and 3). δ 13 C values for micritic carbonate samples (n = 770) range from -16.4 ‰ to -8.7 ‰ (mean = -11.0 ‰). Micrite δ 18 O values range from -245 16.6 ‰ to -6.5 ‰ (mean = -9.3 ‰). Sparite samples (n = 19) display δ 13 C values between -21.2 ‰ and -9.5 ‰ (mean = -12.6 ‰) and δ 18 O values between -23.4 ‰ and -8.6 ‰ (mean = -14.8 ‰) (Fig. 2). Twelve samples were removed due to relatively low weight percent carbonate (< 50 %). Low weight percent carbonate introduces the possibility that the carbonate in these samples does not reflect pedogenesis, or they could be higher porosity samples that were more susceptible to alteration. Fifty percent is a somewhat conservative cutoff that only retains samples that are majority carbonate. Considering only sampling 250 sites with at least three nodules, the mean within-site range in δ 13 C is 0.7 ‰ but can be up to 3.2 ‰. The mean within-site range in δ 18 O is 0.9 ‰ and up to 5.8 ‰. Together, the within-site range in both δ 13 C and δ 18 O are generally < 1 ‰.
It is possible that small amounts of unrecognized spar or areas of recrystallization were drilled with the micrite, which could contribute to the large spread in δ 18 O towards more negative values (Fig. 2). As such, samples with δ 18 O values greater than 2σ below the mean micrite δ 18 O value (-11.7 ‰) were identified as potentially incorporating a secondary spar phase and are 255 shown as such in Figs. 3 and 4. Secondary recrystallization in a closed system does not affect carbon isotopes to the same extent as oxygen isotopes (Cerling, 1984), suggesting that partially recrystallized samples may still record CIEs reliably. Two stratigraphic levels from the base of the Basin Draw section had anomalously low δ 13 C values (up to -16.4 ‰) that fall outside of the typical range for early Eocene pedogenic carbonate from this region and could incorporate some mixing with sparite, particularly given their relatively low δ 18 O values (up to -12.1 ‰). 260 Replicate samples that were analyzed by the two laboratories (n = 10) show that δ 13 C values differed by up to 0.7 ‰ (mean = 0.4 ‰, standard deviation = 0.3 ‰) while δ 18 O values differed by up to 0.6 ‰ (mean = 0.3 ‰, standard deviation = 0.2 ‰).
Greater differences between the two labs are consistently associated with more negative values, and δ 13 C and δ 18 O values produced in the CUBES-SIL lab are generally more negative than values from the UMSIL. These offsets are likely due to differences in standardization between the labs and this relationship was used to correct the UMSIL data to the CUBES-SIL 265 values (Appendix B). This laboratory offset is within the range of observed within-site variation described above, indicating that the offset is not likely to affect the stratigraphic interpretations relevant to this study.

Carbon isotope excursions
Three stratigraphic intervals have consistently low carbon isotope values (< ~ -12 ‰) relative to background levels (between 275 about -8 ‰ and -12 ‰, mean background δ 13 C = -10.3 ‰), though the precise BCM levels at which these intervals occur varies among the sections (Fig. 4a). The lowest excursion occurs between the 405-430 BCM levels, the middle excursion between the 430-455 BCM levels, and the highest excursion between the 460-500 BCM levels. When the data are plotted against dGPS elevation, the variability in their stratigraphic position is reduced with the lowest excursion occurring between 1320-1340 m in elevation, the middle excursion between 1340-1355 m, and the highest excursion between 1370-1405 m 280 (Fig. 4b). Oxygen isotopes do not show any clear differences between hyperthermal and non-hyperthermal intervals (Appendix C). This is consistent with observations that changes in pedogenic carbonate δ 18 O during the PETM are muted compared to https://doi.org/10.5194/cp-2021-83 Preprint. Discussion started: 13 July 2021 c Author(s) 2021. CC BY 4.0 License.
the δ 13 C record, possibly due to temperature fractionation in the carbonate acting in an opposing direction to temperature fractionation in the atmosphere (e.g., Koch et al., 2003).
The magnitudes of the CIEs also vary, both across the different sections and from one CIE to another. The average magnitude 285 of the stratigraphically lowest CIE is 2.4 ‰, with minimum mean δ 13 C values ranging from -13.3 ‰ to -12.0 ‰ in different sections. The average magnitude of the middle CIE is 1.3 ‰, with minimum δ 13 C values of -12.7 ‰ and -12.2 ‰. The average magnitude of the stratigraphically highest CIE is 1.6 ‰ with minimum δ 13 C values of -12.2 ‰.   dGPS. Each of these methods of correlation has its own inherent uncertainty. For example, early fossil collecting in the region may not have specified the producing level for each fossil specimen. The BCM levels therefore reflect some degree of 310 stratigraphic averaging of fossil horizons. These composite levels also rely on bed tracing, which is difficult to do precisely over long distances in this low-relief terrane. Although the dGPS elevations presented here account for some of this uncertainty in bed tracing, they do not account for changes in local dip. Although dip is generally close to zero in this area, bed tracing combined with dGPS measurements has highlighted some local variation in bedding attitude in the field area. The elevations also may be affected by variations in paleo-topography, including lateral variation in paleosol thickness due to relative position 315 on the floodplain and channel meander. Paleosol and sandstone thickness can vary on a meter scale in this area, which is larger than the vertical error of the dGPS (~50 cm). Vertical dGPS error on the elevations therefore likely contributes very little to the observed variation among sections.
Many aspects that complicate stratigraphic correlation in this region are illustrated with the Basin Draw section. This section is located south of Fifteenmile Creek, across from the main portion of the composite section ( Fig. 1)  Despite uncertainties, stratigraphic levels for each of the three CIE intervals are in better agreement across the five sections using dGPS elevation, than with BCM levels (Fig. 4). Using elevation, the CIE in the Basin Draw section correlates better with the stratigraphically lowest excursion in other sections, particularly when we exclude the outliers at the base of the section that appear to be diagenetically altered and not associated with a well-defined excursion. The local biostratigraphy rules out the possibility that these basal samples could be capturing the PETM and the very negative (< -14 ‰) carbon isotope values 330 are more consistent with diagenesis. It should be noted that the faunal analysis of Chew (2015) places the top of the Basin Draw section within faunal event B-2, and therefore predicts that it would overlap with the middle CIE, which is inconsistent with this correlation.

Identification of post-PETM hyperthermals
To account for the geographic variability in bed thickness described above, the five stratigraphic sections, as well as the three 335 shorter sections through fossil localities, were compiled into a "Fifteenmile Creek Composite Section" (Fig. 5), where the peak negative δ 13 C values from each CIE are aligned and the relative spacing between CIEs is scaled to the average of the spacing from the sections where the two primary excursions were measured together. Using this composite section, the lowest CIE occurs between the ~35 to 55 meter level. The middle CIE occurs between the ~57 to 70 meter level. The highest CIE occurs between the ~85 to 104 meter level, with some relatively low δ 13 C values above this that could possibly represent the onset of 340 a fourth CIE (Fig. 5). We confidently identify the two stratigraphically lowest CIEs as the ETM2 and H2 hyperthermals, based https://doi.org/10.5194/cp-2021-83 Preprint. Discussion started: 13 July 2021 c Author(s) 2021. CC BY 4.0 License. on their position relative to Biohorizon B. Elsewhere in the Bighorn Basin, ETM2 lies just above Biohorizon B (Abels et al., 2012;D'Ambrosia et al., 2017). Biohorizon B occurs at the ~11 meter level in our new composite section, suggesting the lowest CIE that starts at ~35 meters is ETM2 (Fig. 5). Magnetostratigraphic data are consistent with this correlation, as the Chron C24r -C24n.3n polarity reversal is known to lie just above ETM2 and correlates to somewhere in the lower 100 meters 345 of our new composite section (Tauxe et al., 1994;Clyde et al., 2007).
Identification of the stratigraphically highest CIEs is more difficult. One possibility is that they represent I1 and the onset of I2. Another possibility is that they correspond to a less well-defined negative CIE between H2 and I1 near the 260 m level in the McCullough Peaks composite section that was attributed to local environmental factors (Abels et al., 2016) (Fig. 5). Barnet et al. (2019) have recently shown a small CIE at roughly the same level in the marine record, suggesting that this could be a 350 global signal. Biostratigraphy and magnetostratigraphy do little to help resolve these options, due to the small fossil sample sizes in McCullough Peaks, uncertainty on the position of the C24r -C24n.3n reversal in Fifteenmile Creek, and the overall differences in sediment thickness across the basin (Clyde et al., 2007). The stratigraphic spacing between the lower and higher CIEs in our section supports the second scenario (correlation to the poorly defined negative CIE between H2 and I1) supporting the idea that this isotopic signal may be more regional or global in scale (Fig. 5). It is also possible that there is a hiatus in the 355 upper part of the Fifteenmile Creek section or that sediment accumulation rates were much slower in this interval. In either case, the stratigraphically highest excursion could represent I1. Additional sampling and better age control are needed to resolve the precise identification of these upper CIEs. Despite having lower magnitudes than the CIEs in McCullough Peaks, the Fifteenmile Creek CIEs are still amplified relative to the marine benthic record, where the average ETM2 and H2 CIEs are ~1.4 ‰ and 0.8 ‰, respectively (Stap et al., 2010). 375

Comparison to other post-PETM hyperthermal records
The amplification effect, however, seems to be somewhat variable even in similar stratigraphic records that are only ~80 km apart within the same basin. This supports the hypothesis that some (or all) of the observed difference in magnitudes between marine and terrestrial CIEs is due to local environmental or time averaging effects. However, some localities show small offsets from their predicted position relative to the CIE record. Locality D-1311 (442 BCM level), for example, falls within the interval associated with faunal event B-2, but is just above the CIE associated with H2. Locality D-1403 (420 BCM level) is within the ETM2 CIE but fell outside of the interval reported to contain faunal event B-1 (Fig. 6). Error in the stratigraphic correlations used to construct the BCM framework or averaging of specific collecting levels within localities could both contribute to these small offsets. Many of these localities have a relatively small number of 400 specimens and would therefore have little effect on the faunal analysis (e.g., D-1403 has 34 specimens). The misplacement of richer localities, like D-1454 (1151 specimens) and D-1460 (405 specimens), would be more significant. These inconsistencies have likely muted the turnover signal observed by Chew (2015). Moreover, a confounding problem with correlating faunal turnover to many of these localities, is that fossiliferous horizons often occur at multiple levels within a locality. For example, in locality D-1350, mammal fossils have been collected both above and below the onset of the ETM2 405 excursion and most collection records are not precise enough to indicate which specimens came from which levels. Thus, a turnover signal could be muted by mixing taxa from within the excursion with pre-and post-excursion taxa. This problem can be overcome only by recollecting such localities using even higher stratigraphic precision.
More importantly, these results demonstrate that there is a reasonably large miscorrelation in the BCM levels of the localities The δ 13 C record through the locality is consistent with that from the main section and suggests that it is not within ETM2 (Appendix A). It seems likely that the Basin Draw section localities should be shifted downwards by ~20 meters in the Bown et al. (1994) composite section. This would place the productive localities D-1454 and D-1460 within Biohorizon B and would likely strengthen the distinction between Biohorizon B and the subsequent faunal events B-1 and B-2 in a revised faunal 420 analysis should further stratigraphic work demonstrate that this move is warranted.
The isotope results presented here suggest that correlation of key localities should be tested by constructing new localityspecific isotope records that can be directly tied to the mammal collection from the same locality, similar to the sections through D-1454, D-1207, and D-1532 (Appendix A). For example, the locality-specific isotope record for locality D-1207 (448 BCM level) shows that this locality records the onset of the H2 CIE, rather than its end as implied by the stratigraphic 425 range reported for faunal event B-2. Developing similar sections through other individual localities, where the primary fossil producing layer can be reliably constrained and tied into a local dGPS elevation and stable isotope record and then compared to other localities either within or outside of CIEs, would clarify observed faunal turnover and its potential relationship to climate fluctuations. Such precise correlations could also be used to investigate the apparent lack of faunal change near the 486 m BCM level (Chew, 2015), where we do find a small negative CIE (Fig. 4). 430

Effects of landscape heterogeneity on terrestrial carbon isotope records
The degree of isotopic fractionation in terrestrial systems can vary depending on environmental factors, many of which are 440 known to fluctuate during hyperthermals. Together, these factors can complicate efforts to use paleosol carbonate δ 13 C to reconstruct atmospheric δ 13 C fluctuations, which is critical for understanding the source of carbon input during hyperthermals as well as understanding how these changes could be recorded differently across a landscape. These fractionation steps can be broadly grouped into those associated with local vegetation (e.g., fractionation associated with photosynthesis and respiration) and those related to soil processes (e.g., microbial respiration, diffusion through soil, and calcite precipitation) (Koch, 1998). 445 Carbon isotope fractionation in plants is dependent on light, water availability, temperature, soil nutrients, and atmospheric pCO2 (e.g., O'Leary, 1988;Farquhar et al., 1989). The strongest single environmental control on δ 13 C values in C3 vegetation is water availability to plants, which accounts for about half of the variability, with plants discriminating less against 13 C during drier times (Ehleringer, 1993;Stewart et al., 1995;Diefendorf et al., 2010;Kohn, 2010Kohn, , 2016. Position on the landscape can also affect water stress on vegetation by affecting soil drainage. Thus, differences in precipitation and local drainage can have 450 large effects on δ 13 C values (Codron et al., 2005). Floral composition and vegetation density can also affect fractionation of carbon isotopes by plants. For example, gymnosperms have higher δ 13 C values than angiosperms (Smith et al., 2007), so local differences in the relative abundance of these plants could affect soil δ 13 C values. Closed forest canopies tend to hold in moisture and receive lower levels of solar radiation, which result in slower photosynthetic rates in C3 plants (Ehleringer et al., 1986;Stewart et al., 1995). Both of these factors, along with possible recycling of 13 C-depleted CO2 below the canopy, result 455 in lower δ 13 C values in C3 leaves (Van der Merwe and Medina, 1991;Cerling and Quade, 1993;Cerling et al., 2004). In contrast, C3 plants in open areas where plants receive greater solar radiation are enriched in 13 C (Ehleringer et al., 1986).
Increased carbon isotope discrimination by plants under higher atmospheric pCO2 conditions has also been observed and proposed as a significant driver of terrestrial CIE amplification during the PETM Jahren, 2012, 2013). However, more recent studies have suggested this effect is transient, and not relevant to carbon isotope changes over long timescales 460 (Kohn, 2016).
Within soils, respiration from plant roots is traditionally considered not to fractionate carbon significantly (Lin and Ehleringer, 1997;Koch, 1998), although carbon isotope fractionation of up to 6-10 ‰ has been documented during respiration (Ghashghaie et al., 2003). The degree of fractionation during respiration can vary between species and under different environmental conditions. For example, Brüggeman et al. (2011) show a wide range of microbial respiration fractionation 465 values in the presence of C3 plants from -6 ‰ to +8 ‰ with an average of ~ +0.5 ‰. Microbial methanogenesis can also contribute 13 C-depleted methane to soils, thereby decreasing the δ 13 C of the soil gas, although this process is largely restricted to water-logged soils (Wynn et al., 2005), which are not found in our study area. Elevated atmospheric pCO2 has been shown to enhance respiration by both plants and microorganisms (Karberg et al., 2005), potentially contributing to differences in carbon isotope values across the landscape during hyperthermals. Cerling (1984) found that soil CO2, on average, should be ~4.4 ‰ enriched in 13 C relative to soil respired CO2 due to diffusion.
Factors such as grain size or soil moisture may affect the range of expected carbon isotope fractionation due to diffusion within soils, potentially contributing to some of the variation across paleosols in a landscape, though little is known about these potential effects. Carbon isotopes also show a slight temperature dependent fractionation during carbonate precipitation, though these effects are negligible relative to many of the other factors discussed here, and the fractionation associated with 475 carbonate precipitation in soils is generally approximated to be ~ +10.5 ‰ (Koch, 1998). It has also been shown that when atmospheric pCO2 is high or respiration rates are low, atmospheric CO2 can penetrate deeper into the soil profile and soil carbonates will incorporate a larger proportion of atmospheric CO2 (enriched in 13 C) relative to plant respired CO2 (depleted in 13 C) (Cerling, 1984). This process is non-linear with depth, however, and pedogenic carbonate nodules collected below 30 cm depth will continue to incorporate predominantly plant respired CO2. 480 Several of the factors discussed above can reasonably be rejected as being significant contributors to the observed differences in δ 13 C values between McCullough Peaks and Fifteenmile Creek, including atmospheric pCO2 variation and any resulting changes in photosynthesis or penetration depth in the soil, temperature effects on calcite precipitation, plant and microbial respiration in soil, and diffusion. Factors that could play a large role include soil drainage, water stress on vegetation, and floral composition and density. Abels et al. (2016) found that precipitation change during the post-PETM hyperthermals appears to 485 have been negligible, based on mean annual precipitation estimates using the CALMAG proxy. This contrasts with drier conditions interpreted during the PETM from various soil proxies and fossil leaves (Wing et al., 2005;Kraus et al., 2013;Abels et al., 2016). Paleosol analysis across the Bighorn Basin shows that soil drainage varies significantly even across a small geographic area (e.g., Kraus, 1992Kraus, , 1997. Secord et al. (2008) inferred vegetation structure based on carbon isotopes in mammalian tooth enamel from the Fifteenmile Creek area and found that vegetation was relatively dense, but that the forests 490 had an open canopy. The "canopy effect," resulting from highly 13 C-depleted leaves in the understory (e.g., Cerling et al., 2004), also appears to be absent in mammalian enamel from the northern Bighorn Basin (Koch et al., 1995) based on the model in Secord et al. (2008). Given the strong control that water availability has on leaf values, this appears to be the most likely cause of differences in absolute δ 13 C values and in differences in CIE magnitude observed here, although differences in floral composition between the two locations cannot be ruled out as contributing to these differences. 495

Conclusions
The ETM2, H2 and potentially I1 hyperthermals are recognized for the first time in pedogenic carbonate from the southern Bighorn Basin (Fifteenmile Creek area). Correlation of the ETM2 and H2 carbon excursions with the faunal record, supports previous suggestions that some mammalian turnover in the early Wasatchian North American Land Mammal Age, following the PETM, is associated with early Eocene hyperthermals. The new isotope records presented here also highlight some of the 500 complexities of lithological correlation between the low-lying exposures in Fifteenmile Creek area and support the use of an independent means of stratigraphic correlation, such as precise elevation, in conjunction with mammalian biostratigraphy and https://doi.org/10.5194/cp-2021-83 Preprint. Discussion started: 13 July 2021 c Author(s) 2021. CC BY 4.0 License. carbon isotope chemostratigraphy. The magnitudes of the ETM2 and H2 CIEs (2.4 ‰ and 1.3 ‰, respectively) are smaller than what is seen farther north in the McCullough Peaks region of the Bighorn Basin (3.8 ‰ and 2.8 ‰). We suggest that local variation in water availability to plants, and potentially other vegetation and soil processes, likely account for much of the 505 differences in carbon isotope values observed between the two locations.

Carbon isotopes stratigraphy through fossil localities
Three locality-specific stratigraphic sections were collected for stable isotope analysis to identify whether they captured a CIE.
If a CIE fell entirely within the fossil producing levels of a single locality, future fossil collecting could target these areas for 510 more refined faunal analyses that would eliminate uncertainty related to stratigraphic averaging and correlation between These smaller locality-specific sections can then be compared to their nearest main section using the BCM levels and elevation to further assess spatial variability in δ 13 C and confirm the correlations (Fig. A1). The carbon isotope values are relatively 515 consistent over short distances in Fifteenmile Creek (Fig. A1). The Bown et al (1994) composite level (BCM) and dGPS elevation are also both fairly reliable stratigraphic indicators over short traces and correlate the locality sections to the main sections in consistent ways (Fig. A1). An exception to this is correlating the D-1207 locality to the Kraus Flats section, where the BCM level lines up with the H2 excursion, but the dGPS elevations have them offset from each other. D-1207 (448 BCM level) is the only one of these three localities that falls within an excursion (H2). This level is predicted to be within the B-2 520 faunal event in Chew (2015) (435 to 448 m). D-1207 has two taxa that first appear in faunal event B-2 (Hexacodus sp. and Protorohippus venticolum), and this locality therefore represents a place where the carbon isotope stratigraphy and biostratigraphy are consistent and suggest that it falls within H2 and the B-2 faunal event. It also illustrates a need to use both BCM and dGPS for precise correlation of carbon isotope stratigraphy in this area.

Inter-laboratory offsets and correction
Comparison of the carbon and oxygen isotope data from samples run at both CUBES-SIL and UMSIL shows that CUBES-SIL produces generally lower values than UMSIL in this dataset's range of values. This offset to lower values is systematic https://doi.org/10.5194/cp-2021-83 Preprint. Discussion started: 13 July 2021 c Author(s) 2021. CC BY 4.0 License. and is larger at progressively lower values. This issue affects both isotopes but is particularly true for carbon, which suggests 535 that subtle differences in instrumentation and standardization procedure are the causes of this offset. UMSIL uses a ThermoFisher MAT253 attached to a Kiel IV automated preparation device and data are normalized using a best-fit linear regression to NBS-18 and NBS-19 reference materials (δ 13 C = -5.01 ‰ and +1.95 ‰, δ 18 O = -23.2 ‰ and -2.20 ‰, respectively) by regressing measured, 17 O-corrected values against the published values for the standards. Corrections are determined for both carbon and oxygen and these corrections are checked routinely using either or both,540 and tend to be stable over several months.
CUBES-SIL uses a Thermo Delta V gas source, continuous flow isotope ratio mass spectrometer. Samples are measured along with 3-4 standards that either have internationally accepted values or are in-house standards that have values determined relative to these standards (e.g., NBS-18, NBS-19, and LVSEC for light 13 C standards). Measurements of the standards bracket sample measurements in each run and are used to assess behaviors of the instrument that may need to be corrected. Standards 545 covering the full range of signal intensities observed in the samples are measured to assess effects of linearity. Standards are also run intermittently throughout the analytical session to evaluate instrument drift over the course of a run. Lastly, the overall offset of the standards from accepted values is evaluated. These corrections are done independently for carbon and oxygen and are checked using a monitoring standard that is treated as an unknown. In this dataset, linearity and drift were often negligible, and corrections for these effects were only applied when needed. Raw or linearity/drift corrected values were then corrected 550 to final values similarly to UMSIL, by applying a regression between two standards that span a range in values. Runs 1-13 and 15 used NBS-18 as the negative anchor for both δ 13 C and δ 18 O, while runs 14, 16, and 17 used a MERCK carbonate as the negative δ 13 C anchor point (δ 13 C = -35.6 ‰), to better standardize for the range of δ 13 C values found in typical terrestrial carbonates, following recommendations in Coplen et al. (2006).
Both UMSIL and CUBES-SIL used NBS-18 as the negative δ 13 C for most measurements. NBS-18 is ~4-10 ‰ higher than 555 most of the measurements in this dataset, so the linear regression that is used for the scale correction in both labs is extrapolated far beyond that lowest anchor point. As a result, uncertainties in this regression in both labs, while small, can result in differences in sample values that are larger than analytical uncertainty. If this were the cause of the difference in sample values between CUBES-SIL and UMSIL, one would expect to see correlation between the values from both labs, as well as a progressively larger offset between sample values at more negative δ 13 C values. A strong correlation and positive slope are 560 seen between the isotope ratios measured at CUBES-SIL relative to UMSIL in both δ 13 C and δ 18 O values for a small subset of samples that were run in both labs (Fig. B1). The strength of the correlation for δ 13 C suggests that we can correct data run from one lab to be on scale with the other lab (R 2 = 0.88, S = 0.37 ‰). We compared a subset of values produced at CUBES-SIL using either NBS-18 or the lower MERCK standard for the same sample. This showed no systematic difference in the values (Fig. B2), which suggests that CUBES-SIL data are internally consistent, despite the extrapolation of the scale 565 correction for many of the runs. Given that, and that the majority of this dataset was analyzed at CUBES-SIL, we have chosen to correct the UMSIL data to these values, using the regression equations for both carbon and oxygen (Fig. B1). After applying the correction, the UMSIL data become more negative, while still following the same stratigraphic patterns in δ 13 C (Fig. B3).
While we think that this correction scheme is most appropriate for these data, we recognize that it may complicate comparisons to other previously published work that do not use δ 13 C standards that account for values in this range. Further, because the 570 δ 13 C correlation between labs is strong, the magnitudes of the CIE's are not appreciably affected in the combined data after correction. Oxygen isotope values exhibit a weaker correlation between the labs (R 2 = 0.74), indicating that considerably more noise is added by combining the data. However, highly precise δ 18 O values are not important for our interpretations.

Data availability
Raw and processed stable isotope data, as well as R scripts used for data processing and creating figures are available in Open Science Framework (OSF) and can be accessed here: https://osf.io/3z2xc/?view_only=1e9f419b168c43f5b0ae0b8e8f09584a 600