The paleosol description and classification in this study are mainly based on the field observations of soil texture and structure, mottled color and mineral accretions. The B-horizon in the soil profile is the key interest of this study, for its important role in reflecting the drainage condition and assisting the pCO₂ reconstruction. Several key features were used to identify the B-horizons during our field examination of paleosols: 1) B-horizons often display reddish, yellowish, or brownish colors due to the accumulation of iron oxides or clay minerals; 2) B-horizons usually exhibit distinctive ped structures, such as wedge-shaped peds or slickensides; and 3) accumulated materials are often present in B-horizons as evidence of illuviation, such as the shiny clay films that form along the wedge-shaped peds or slickensides, the iron oxide coatings, the pedogenic carbonate nodules, etc. (e.g., Mack et al., 1993; Retallack, 1997; 2001; Kraus, 1999; Daniels and Haering, 2006; Kraus et al., 2013).

Carbon isotope analysis in this study includes analysis of bulk organic matter δ¹³C_org, paleosol carbonate δ¹³C_c, and leaf wax n-alkanes δ¹³C _n-alk. Bulk paleosol and mudstone samples were collected every 0.5–1.0 m from the floodplain facies throughout the measured sections, and pedogenic carbonate nodules were collected where available. Fresh bulk paleosol and mudstone samples were collected after removing the surface layers to avoid any possible contamination from recent organic material. Pedogenic carbonate nodules were collected at a depth greater than 30 cm below the paleosol surface in order to reduce the influence of atmospheric CO₂ on δ¹³C_c. Petrographic analysis and X-ray diffraction analysis have been done to these pedogenic carbonate nodules in our earlier study (Gao and Fan, 2018), and the results show that these paleosol carbonates are dominated by low-magnesium calcite, predominantly micritic, suggesting minimal diagenetic influence.

For δ¹³C_org analysis, the samples were examined under a microscope and any obvious modern organic matter was removed by hand picking. Samples were then powdered using a ceramic mortar and pestle and weighted into open silver capsules. Carbonate in the powdered samples was removed using hydrogen chloride fumigation method. Samples were loaded into silver capsules and wetted by adding one drop of deionized water. Open silver capsules were placed in a sealed desiccator containing a beaker of 12 N fuming HCl for at least 4 days, then dried in oven at 75°C for at least 12 h. Each silver capsule was then carefully placed in a tin capsule and tightly packed and loaded into a Costech ECS 4010 Analyzer, which was connected to a Delta V Advantage Mass Spectrometer via ConFlo IV, at the University of Texas at Arlington. The isotope ratios were calibrated using USGS 40, USGS 41, and the precision of δ¹³C_org is better than ±0.2‰ (1σ) based on repeated analysis of standards. The δ¹³C_org values are reported relative to Vienna Pee Dee Belemnite (VPDB).

For δ¹³C_c analysis, the carbonate samples were powdered using a ceramic mortar and pestle. Powdered samples were weighed and loaded into sample vials. After flushing the vials using ultrapure helium, dehydrated phosphoric acid was injected into the vials to react with the samples at 70°C. The δ¹³C values of the resulted CO₂ were analyzed using a Gasbench II and a Delta V Advantage Mass Spectrometer at the University of Texas at Arlington. The isotope ratios were calibrated using NBS-19, NBS-18, and the precision of δ¹³C values is better than ±0.2‰ (1σ) based on repeated analysis of standards. All the isotope values of carbonate were reported relative to VPDB.

Based on the trend of δ¹³C_org data and potential PETM, six paleosol samples were selected to analyze δ¹³C _n-alk and δD _n-alk at the Louisiana State University to verify the presence of PETM in our record. These six samples include four within and two outside of the interpreted PETM. Lipids were extracted from 30 to 60 g of powdered rock samples using Soxhlet extractors with a solvent mixture of dichloromethane (DCM) and methanol (MeOH) 2:1 (v/v) for over 48 h. Total lipid extracs (TLE) were evaporated using a TurboVap II evaporator under a stream of purified nitrogen. TLE was separated into compound classes using ∼4.0 g of pre-extracted, activated (@200°C for 2 h) silica gel. Organic compounds in TLE were separated into apolar, intermediate, and polar fractions by using 4 mL hexane, 4 mL DCM, and 4 ml methanol, respectively. n-Alkanes abundances were determined using a Thermal Trace 1,310 Ga Chromatography (GC)-flame ionization detector (FID) fitted with a programmable temperature vaporization (PTV) injector and a fused silica, TG-1MS column (60 m long, 0.25 mm i. d., 0.25 µm film thickness). Samples were carried by helium at a rate of 2 mL/min. GC oven temperature was ramped from 60°C (holding for 1 min) to 320°C at 15°C/min (holding for 20 min).

Measurements of compound-specific δ¹³C _n-alk values were performed using a Thermal Trace 1,310 Ga Chromatography (GC) coupled to a Thermo Delta V Advantage isotope ratio mass spectrometer (IRMS) interfaced with a Thermo Isolink interface. The GC column and carrier flow conditions were identical to above. Compounds were separated on the GC with a temperature program from 60°C (held for 1 min) to 170°C at 14°C/min, to 300°C at 3°C/min, and then to 325 at 14°C/min with an isothermal holding of 10 min. Carbon isotope ratio values were determined relative to a reference, calibrated against Mix A6 (n-C₁₆ to n-C₃₀; Arndt Schimmelmann, Indiana University), and reported relative to VPDB and with two standard deviations.

Measurements of compound-specific δD _n-alk values were performed using a Thermo Trace 1310 GC coupled to a Delta V Advantage isotope ratio mass spectrometer (IRMS) interfaced with an Isolink interface. The GC column and carrier flow conditions were identical to carbon isotope analysis. Compounds were separated on the GC with a temperature program from 60°C (held for 2 min) to 170°C at 14°C/min, to 300°C at 3°C/min, and then to 325°C at 14°C/min with an isothermal holding of 10 min. The H₃ ⁺ factor (Sessions, 2001) was monitored daily prior to δD analysis and was less than eight for the measurement periods. The drift of the instrument was routinely monitored, and individual n-alkane isotope ratios were corrected to n-alkane reference materials (Mix A3, A. Schimmelmann, Indiana University, Bloomington). δD _n-alk values are reported relative to Vienna Standard Mean Ocean Water (VSMOW) and with two standard deviations (Table 1).

We followed Equation 1 that was presented in Cerling and Harris, (1999) for atmospheric pCO₂ reconstruction. p C O 2 p p m v = S z   δ 13 C s − 1.0044 * δ 13 C r − 4.4 δ 13 C a − δ 13 C s (1)

In the equation, pCO₂ is the concentration of atmospheric CO₂, S(z) is the concentration of soil-respired CO₂, δ¹³C_s is the carbon isotope composition of soil CO₂, δ¹³C_r is the carbon isotope composition of soil-respired CO₂, δ¹³C_a is the carbon isotope composition of atmospheric CO₂. δ¹³C_s is determined from the paleosol carbonate δ¹³C_c and paleosol carbonate precipitation temperature using the temperature-sensitive equilibrium isotope fractionation factor between calcite and CO₂ Equation 2 (Romanek et al., 1992): δ 13 C s = δ 13 C c c + 1000 11.98 − 0.12 T 1000 + 1 − 1000 (2)

δ¹³C_r was traditionally represented by paleosol δ¹³C_org values (e.g., Cerling, 1992; Retallack, 2009). However, several studies have shown that bulk δ¹³C_org values in modern soils are comparable to those found in surface litter, but they tend to increase with soil depth (e.g., Torn et al., 2002; Wynn et al., 2005; Wynn and Bird, 2007). Therefore, organic carbon in soil B-horizons, which are typically sampled when working with paleosol carbonate, have δ¹³C_org values higher than the average biomass, and the δ¹³C_r in Equation 1 can be approximated by measured δ¹³C_org values minus 1‰ (Bowen and Beerling, 2004; Breecker, 2013). δ¹³C_a was −6.1‰ during the early Eocene (Tipple et al., 2010). S(z) is commonly assumed as summer minimum soil-respired CO₂ in this equation because soil carbonates are likely to precipitate preferentially during warm seasons (Breecker et al., 2009; Passey et al., 2010; Quade et al., 2013; Hough et al., 2014). We use the value of 2,500 ppm for S(z) as it is commonly used for the subhumid temperature and tropical climates (Breecker et al., 2010).

Bulk paleosol samples from the B-horizons were analyzed for major element concentrations using a Shimadzu EDX-7000 X-ray fluorescence (XRF) spectrometer at the University of Texas at Arlington. The reported weight percentages were calibrated using three USGS standards: SDO-1 (Devonian Ohio Shale), SBC-1 (Brush Creek Shale), and SGR-1 (Green River Shale) (Supplementary Table S2). The weight percentages were then recalculated to molar ratios, following Sheldon et al. (2002). Since the intensity of chemical weathering of paleosols depends in part on precipitation, Sheldon et al. (2002) developed an empirical relationship (Equation (3)) relating the mean annual precipitation (MAP) to the chemical index of alteration without potassium (CIA-K). The error associated with the regression analysis is ±181 mm (Sheldon et al., 2002; Sheldon and Tabor, 2009). MAP millimeters / year = 221 e 0.0197 * CIA ‐ K ;   r 2 = 0.72 (3) CIA ‐ K = Al 2 O 3 /   Al 2 O 3 + Ca O * + Na 2 O × 100 (4)

The chemical index CIA-K, which excludes potassium (K), is derived from Equation 4. K is excluded in calculation because diagenesis can yield elevated K concentrations in paleosols (McLennan, 1993). CIA-K is calculated based on the molar ratio of Al₂O₃, to the sum of Al₂O₃, CaO*, and Na₂O. CaO* is the smaller value between CaO and Na₂O to account for CaO derived from the silicate fraction, excluding the carbonate fraction (McLennan, 1993).

The δ¹³C_org values of the organic-rich pedogenic mudstone that is usually accompanied with coal beds in the Paleocene Fort Union Formation vary between −23.2‰ and −26.5‰, with an average of −25.0‰ (Supplementary Table S1; Figure 4b). The δ¹³C_org values of the paleosol B-horizons in the main body of the lower Eocene Wasatch Formation vary between −24.0‰ and −29.4‰, with an average of −26.3‰ (Supplementary Table S1; Figure 4b).

The δ¹³C_org values of the non-pedogenic mudstone samples have similar variations as the δ¹³C_org values of paleosols. In the Fort Union Formation, the δ¹³C_org values of non-pedogenic mudstone samples vary between −22.8‰ and −26.1‰, with an average of −24.5‰ (Supplementary Table S1; Figure 4b). In the main body of the Wasatch Formation, these δ¹³C_org values vary between −22.5‰ and −27.2‰, with an average value of −25.7‰ (Supplementary Table S1; Figure 4b).

At ∼70 m above the base of the Wasatch Formation, the δ¹³C_org values show a transient drop of ∼4‰, and the CIE spans ∼70 m on the record (Figure 4b). The δ¹³C _n−alk values in this interval show a negative excursion on the order of ∼4‰–5‰ on long−chain n−alkanes with 27, 29, 31, and 33 carbon atoms (Table 1; Figure 5), δD_n−alk values in this interval show increases of 30‰–50‰ (Table 1; Figure 6). Another CIE of ∼2‰ is documented in the δ¹³C_org record after the Paleocene−Eocene boundary, at the 490–590 m level above the base of Wasatch Formation, and is older than ∼53 Ma based on the biostratigraphic age constraint of our studied strata (Figure 4b).

Based on the well calibrated major element data (Table 2; Supplementary Table S2), the reconstructed MAP throughout the Paleocene and the earliest Eocene is generally high, in the range of 800–1,500 mm (Table 2; Figure 4b). An extremely low value of 408 mm occurs at 168.8 m of Wasatch Formation, however, this estimate is not accurate due to its high Na₂O content (5.1 wt%) compared to the other samples (<2.8 wt%). The elevated Na₂O is likely due to diagenetic salt, rather than a contribution from silicate. The other two samples at 405 m and 588 m of Wasatch Formation, where Na₂O content is greater than 2.0 wt%, also show low MAP estimates. The atmospheric pCO₂ during the early Eocene was estimated to be in the range of 600–900 ppm from the six selected paleosol samples in the Wasatch Formation (Table 3).

We identify CIEs in both the δ¹³C _n−alk and δ¹³C_org records near the Paleocene−Eocene boundary in the greater Green River Basin (Figures 4b, 5). The magnitude is on the order of 4–5‰, which is consistent with, or smaller than those of CIEs of the PETM in marine and other continental records. In marine benthic foraminifera δ¹³C_c record, the PETM is expressed as a ∼3‰ CIE (e.g., Zachos et al., 2001). Continental records show that the CIE during the PETM has variable magnitude, depending on the type of samples studied. Such CIE records include a ∼6‰ drop based on studies of bulk organic δ¹³C_org in the Paris Basin, France (Magioncalda et al., 2001), a ∼6‰–8‰ drop based on studies of paleosol carbonates δ¹³C_c in the Bighorn Basin (Bowen et al., 2001; Koch et al., 2003; Abels et al., 2016), a ∼5‰ drop based on a study of mammalian tooth enamel δ¹³C in the Bighorn Basin (Secord et al., 2012), and a 4‰–5‰ drop based on studies of bulk organic δ¹³C_org and leaf wax δ¹³C _n−alk in the Bighorn Basin (Baczynski et al., 2013; 2016). The depositional duration of the ∼4–5‰ CIE in our record seems also to match the duration of PETM in the Bighorn Basin. Our CIE spans ∼70 m, aligning with the stratigraphic interval of PETM in the Bighorn Basin, which varies between 40 m and 100 m depending on the specific study location (Bowen et al., 2001; Koch et al., 2003; Secord et al., 2012; Baczynski et al., 2013; 2016). Therefore, we suggest that the ∼4–5‰ CIE in our record is the PETM. The variations in CIE magnitude in different locations could reflect differences in soil properties, the type of C₃ vegetation (e.g., Cerling and Quade, 1993), the mixing of carbon from different sources (e.g., autochthonous or allochthonous), and different degrees of organic matter degradation (Baczynski et al., 2013; 2016).

The second CIE on our δ¹³C_org record, with a magnitude of ∼2‰, occurred after the PETM but before ∼53 Ma (Figure 4b). This CIE most likely represents the ETM−2 hypothermal event. The ETM−2 is the most pronounced early Eocene hyperthermal event and is expressed as a ∼1‰ negative CIE in marine bulk sediment records (Zachos et al., 2010). In the continental environment, ETM−2 has been documented as a ∼3‰ drop in paleosol carbonate δ¹³C_c in the Bighorn Basin (Abels et al., 2016). Although our δ¹³C_org record does not show an ideal CIE that could precisely represent the ETM−2, the timing and magnitude of this excursion resemble those of the ETM−2, or a combination of ETM−2 and other smaller hyperthermals that happened within the next ∼500 kyr of the ETM−2 (e.g., Cramer et al., 2003; Lourens et al., 2005).

Kraus et al. (2013) reconstructed the MAP in the Bighorn Basin based on a ∼70 m thick interval of paleosols covering the PETM, their results show that the MAP values in the pre− and post−PETM intervals are generally high, in the range of 1,100–1,500 mm, and the values dropped to 800–1,100 mm during the main phase of PETM. However, our estimated MAP remains relatively stable, ranging between 800–1,500 mm during the entire LPEE (Figure 4b). The absence of a decrease in MAP during the PETM in this studied section contrasts with records from other basins and the transient drying observed in leaf wax δD_n−alk and changes in paleosol morphology in this study. For example, a dry PETM was also suggested from a sedimentology study in the Axehandle Basin in central Utah (Bowen and Bowen, 2008) and a leaf physiognomy study in the Bighorn Basin (Wing et al., 2005). Several studies have questioned the accuracy of MAP estimates based on CIA−K, particularly for thinner paleosols (B horizon less than 1 m thick), as they had less time for weathering and may not have a cation distribution representative of precipitation (e.g., Adams et al., 2011). The possibility that the low sampling resolution in our record is responsible for the absence of decreasing MAP cannot be completely ruled out.

The 30‰–50‰ increase in leaf wax δD _n−alk values in our record show transient drying during the PETM (Figure 6). Studies of terrestrial plants and soils have suggested that leaf waxes are refractory, and their δD _n−alk values primarily reflect meteoric water δD values (e.g., Sauer et al., 2001; Smith and Freeman, 2006). The increase in δD _n−alk values during the PETM most likely reflects a combination of increasing meteoric water δD due to decreasing precipitation amount, increasing soil evaporation, and reducing hydrogen isotope fractionation between water and leaf wax during climate drying. It has been suggested that soil evaporation enriches deuterium in soil water and decrease the apparent fractionation factor, especially in arid climate (e.g., Smith and Freeman, 2006). The development of different types of paleosols throughout our measured section further supports the inferred transient drying. During the Paleocene, the floodplains were generally water−logged and had histosols (coal beds) developed, indicating that the climate was wet. During the early Eocene, paleosol types within the PETM are different compared to those in the rest of the section. Within the PETM, poorly drained red paleosols and calcisols developed, indicating dry environment. Outside of the hyperthermal, the floodplain deposition was dominated by well drained green paleosols, indicating wet climate. Soil type representing wet paleoclimate changed from histosols during the Paleocene to green paleosols during the early Eocene, possibly reflecting a change in ecosystem when climate became very hot during the early Eocene.

The transient drying during the PETM may be explained by the mid−latitude bypass mechanism, which was proposed to explain the elevated runoff and deuterium−enriched precipitation in the Arctic during the PETM (Pagani et al., 2006). Decreased meridional temperature gradients during the PETM might have reduced rainout of water vapor by synoptic eddies in the mid−latitudes, and increased moisture delivery to the high latitudes, because a decreased temperature gradient provides a smaller energy supply for synoptic eddies (Pierrehumbert, 2002; Pagani et al., 2006). This mechanism predicts temporally drier conditions during the PETM than other times during the early Paleogene in the mid−latitudes, and spatially drier conditions in the southerly sites than in the northerly sites in Wyoming during the PETM (Bowen and Bowen, 2008). Bowen and Bowen (2008) found that the CIE in central Utah was ∼2‰ smaller than the ∼6‰–8‰ drop in the Bighorn Basin during the PETM, both determined from paleosol carbonates δ¹³C_c, and suggested that the smaller CIE could have been caused by drier conditions. The speculative correlation of the PETM with a ∼15−m−thick evaporite sequence deposited in the Flagstaff lake system in central Utah also supports a significant increase in regional aridity during the PETM at this location compared to the Bighorn Basin (Bowen and Bowen, 2008). The CIE magnitude in the greater Green River Basin is consistent with that based on leaf wax δ¹³C _n−alk in the Bighorn Basin (Baczynski et al., 2013; 2016), indicating no evidence of southward aridification.

The general wet climate during the LPEE is consistent with climate model predictions. However, the model also predicts that the PETM, with higher pCO2 levels compared to the rest of the LPEE, should be wetter. By using the Community Climate System Model (CCSM3) with atmospheric CO₂ concentrations of 4×, 8×, and 16× the preindustrial value, Winguth et al. (2010) demonstrated that precipitation during the PETM slightly increases, and evaporation minus precipitation slightly decreases in mid to high latitudes with increased surface temperatures, and attributed the changes to responses to CO₂ radiative forcing (also seen in Sloan and Rea, 1996; Shellito et al., 2003). This prediction is not supported by the observed transient drying. The discrepancy between model predictions and proxy data could arise from several uncertainties. These include, but are not limited to, that model simulations may overestimate CO₂ radiative forcing, as our records show that the pCO₂ levels during the PETM were only 1–2 times higher, not at least 4 times higher, than preindustrial levels; model simulations may not fully account for the topography development in the western U.S.A. during the early Paleogene; and seasonal distribution of precipitation and potential biases in leaf wax records towards growing season precipitation.