Astronomical and tectonic influences on climate and deposition revealed through radioisotopic geochronology and Bayesian age-depth modeling of the early Eocene Green River Formation, Wyoming, USA

Deep-time geologic records of warming climate provide analogs for how Earth systems will respond to rising temperature in the future. The highest temperatures and atmospheric CO₂ levels of the Cenozoic occurred during the early Eocene, including the early Eocene climatic optimum (EECO) (ca. 53–49 Ma; Zachos et al., 2008; Lauretano et al., 2018; Zhu et al., 2019). Superimposed over the long-term warming trends of the EECO are hyperthermal events, transient periods of global warming which are associated with large-scale releases of CO₂ to the ocean and atmosphere and the negative excursions in the carbon isotope record that accompany them (Zeebe et al., 2009; McInerney and Wing, 2011). In the marine record, hyperthermal events have been associated with eccentricity maxima, with the suggestion that eccentricity-modulated seasonality drives increased precipitation, weathering, and sedimentary discharge to the shallow marine environment (Cramer et al., 2003; Zachos et al., 2010; Lauretano et al., 2016, 2018). The dissolution of carbonate material associated with ocean acidification has left offshore strata more condensed than terrestrial or near-shore records, frustrating attempts to identify astrochronological signals in deep marine records (e.g., Westerhold et al., 2017). In contrast, the Wilkins Peak Member (WPM) of the Green River Formation (GRF), Wyoming, USA, records deposition in an alkaline lake (sustained pH >10; Jagniecki et al., 2015), and thus provides a record that is relatively uncompromised by rising atmospheric CO₂, and so is greatly expanded compared to equivalent deep marine sections. Excellent surface exposure of the WPM, numerous drill cores, and a laterally continuous system of distinct bedding surfaces across the basin which permit clear correlation from outcrop to core, provide a high-resolution record of the changing Eocene climate (Culbertson, 1961; Pietras and Carroll, 2006) (Fig. 1). Additionally, the WPM records a portion of the EECO which coincides with the estimated age of a resonance transition in orbital solutions (Laskar et al., 2011; Zeebe and Lourens, 2019). Identifying and radioisotopically constraining the ages of astronomically influenced sedimentary features may provide insights into the validity and accuracy of these astronomical solutions and offer methods to test future solutions.

The WPM is characterized by lithological cycles which have long been interpreted to reflect astronomically forced lake expansion and contraction (e.g., Bradley, 1929; Fischer and Roberts, 1991; Meyers, 2008; Machlus et al., 2008; Aswasereelert et al., 2013). Lake deposits alternate with nine distinctive alluvial intervals that are topped with flooding surfaces, and so were likely deposited during major lowstands (termed marker beds A–I; Culbertson, 1961; Pietras and Carroll, 2006). Prior work has found strong associations between these marker beds and eccentricity (Smith et al., 2010; Aswasereelert et al., 2013). Negative δ¹³C excursions in similar Uinta basin (Wyoming, USA) deposits suggest that the alluvial intervals record deposition during hyperthermals (Smith et al., 2014; Birgenheier et al., 2020). Using the first ²⁰⁶Pb/²³⁸U zircon ages for WPM tuffs, and comparisons of a radioisotopic age model to the astronomical solutions of Laskar et al. (2004) and Smith et al. (2010, 2014) proposed that alluvial sedimentation occurred during minima in short eccentricity. In contrast, Eocene hyperthermals expressed in marine deposits have been interpreted to correspond with short eccentricity maxima (Cramer et al., 2003; Zachos et al., 2010; Piedrahita et al., 2022). To resolve this question and to capitalize on the record of climate change preserved in the WPM, a highly resolved timescale of lacustrine deposition is required.

The GRF occupies intramontane basins located in modern-day Wyoming, Colorado, and Utah (USA). In the Bridger Basin of southwestern Wyoming this formation was deposited by the ancient Lake Gosiute, and is subdivided into four members: the Luman, Tipton, Wilkins Peak, and Laney. The first and last of these represent overfilled lake conditions, and generally consist of deposits associated with lacustrine high stand and fresh water (Pietras and Carroll, 2006). The middle member, the Wilkins Peak, represents deposition by an underfilled lake, and its deposits are characterized by cyclic emplacement of oil shale, evaporite (i.e., trona and halite), mudstone, and siliciclastic marker beds, which have been interpreted as recurring periods of evaporation and lake expansion (Culbertson, 1961; Fischer and Roberts, 1991; Roehler, 1992; Pietras and Carroll, 2006) (Figs. 1–3). WPM lacustrine strata are intercalated with numerous volcanic tuffs, permitting precise radioisotopic dating (Smith et al., 2008, 2010, 2014; Machlus et al., 2004, 2015). Radioisotopic dating efforts in the WPM have largely focused on seven major tuffs: the Firehole, Second, Boar, Grey, Layered, Main, and Sixth tuffs (Smith et al., 2008; Machlus et al., 2008, 2015). The enclosing system of laterally continuous alluvial marker beds, throughout the basin (e.g., Culbertson, 1961) and the extent to which these tuffs have been characterized and described by prior work, permit the identification of these tuffs in both outcrop and drill cores. This stratigraphic framework allows for detailed inter-site correlation with great confidence across the Bridger Basin (Smith et al., 2008, 2010, 2014, 2015) (Fig. 2). However, efforts to develop a high-resolution chronology for the basin have been hampered by both insufficient precision, and by systematic intercalibration uncertainties between radioisotopic clocks. U-Pb dating by Machlus et al. (2015) approached the precision necessary to link astronomical cycles with lithologic cyclicity, but found significant discrepancies between U-Pb and previous ⁴⁰Ar/³⁹Ar age models. Notably, the Machlus et al. (2015),²⁰⁶Pb/²³⁸U ages for the Sixth, Layered, and Main tuffs are younger than the corresponding Smith et al. (2010),⁴⁰Ar/³⁹Ar ages. Machlus et. al. (2015) interpreted the Smith et al. (2010) multi-crystal dates to have been biased by recycling and inheritance of older crystals (e.g., Wotzlaw et al., 2013; Rivera et al., 2014).

We present a new suite of ²⁰⁶Pb/²³⁸U and ⁴⁰Ar/³⁹Ar dates from seven tuffs in the WPM that are correlated into three drill cores spanning the basin (Roehler, 1992). Our new geochronologic data and age models improve upon prior work in two key ways: (1) an increase in precision and accuracy due to improvements in the chemical abrasion–isotope dilution–thermal ionization mass spectrometry (CA-ID-TIMS) methodology in the case of U-Pb zircon dating, and improved mass spectrometry in the case of ⁴⁰Ar/³⁹Ar dating, and (2) the integration of both chronometers using a Bayesian approach that employs the principle of superposition to further constrain a posterior probabilistic age-depth model (Trayler et al., 2020; Haslett and Parnell, 2008). We leverage these advances to quantitatively evaluate the relationship between eccentricity forcing and sedimentation in the WPM, with implications for terrestrial climate response during Eocene hyperthermal events.

Tuffs were collected from four localities: Sage Creek, Firehole Canyon, Apache Lane, and Tollgate Rock (Fig. 1). At Boise State University (Idaho, USA), zircon grains from seven samples were imaged using cathodoluminescence (CL) and subsequently pre-screened by in situ laser ablation–inductively coupled plasma–mass spectrometry (LA-ICP-MS) dating. Since many of the crystals selected for analysis were small, 2 by 12 by 60 µm trenches transecting the crystal rims and cores were ablated, rather than spots. Zircon crystals that did not yield anomalously old LA-ICP-MS ages, and appeared to be free of older inherited components in CL imaging, were then dated by CA-ID-TIMS. Complete CL imagery and LA-ICP-MS data for each tuff, as well as a summary of LA-ICP-MS screening results (Table S1¹) and full CA-ID-TIMS analytical data (Table S2), may be found in the Supplemental Material. Sanidine crystals from two samples were irradiated for 80 h along with the 28.201 ± 0.046 Ma Fish Canyon sanidine (Kuiper et al., 2008), and dated at the University of Wisconsin–Madison (Wisconsin, USA) WiscAr Laboratory using a Noblesse multi-collector mass spectrometer (Jicha et al., 2016; Supplemental Material, see ^{footnote 1}). The ⁴⁰K total decay constant of 5.463 ± 0.107 × 10⁻¹⁰/yr of Min et al. (2000) was used to calculate ages. Several dozen crystals of Fish Canyon sandine that were either co-located with the Grey and Main tuff crystals or from adjacent wells in a 2.38 cm diameter aluminum disk were measured to determine a mean J value for these samples (Table S3, Supplemental Material). Single crystal fusion experiments were performed on each sample to identify and exclude dates from older inherited or reworked crystals from the weighted mean age determinations.

Bayesian age-depth models of the WPM were developed for three separate drill cores using the modified BChron R package of Trayler et al. (2020). The position of tuffs and marker beds in each model represent their height above the base of the WPM (Aswasereelert et al., 2013; Smith et al., 2015). The stratigraphic positions of tuffs in each core are well-constrained due to an enclosing system of laterally continuous alluvial marker beds throughout the basin (e.g., Culbertson, 1961). Crucially, this Bayesian approach does not rely on astronomically tuning the sedimentary record—making assumptions about the sedimentation rate to fit astrochronological targets—and thus permits a more objective test for the association between sedimentary events and astronomical solutions for eccentricity. Moreover, the compound Poisson-gamma accumulation prior used for stochastic interpolation between dated horizons in this approach has been shown to effectively mimic both rapid deposition and near unconformity like events (Haslett and Parnell, 2008). Therefore, the 95% credible interval of the Bayesian model includes and constrains possible hiatuses and periods of highly variable rates of deposition.

Age models were produced for three cores through the WPM (Brownfield et al., 2011): the Blacks Fork 1 core (BF-1, U.S. Geological Survey (USGS) W0080, 41.35647°N, 109.52564°W); White Mountain 1 core (WM, USGS W0081, 41.550856°N, 109.418522°W); and Currant Creek Ridge 1 core (CCR-1, USGS W0420R, 41.181353°N, 109.606624°W). Each model uses the summed probability distribution of all analyses included in the weighted mean for each tuff as inputs for the dated horizons (Tables S2 and S3). ²⁰⁶Pb/²³⁸U and ⁴⁰Ar/³⁹Ar analyses for the Main and Grey tuffs were included as separate probability distributions in each model, as each chronometer represents an independent constraint on the age model at that given depth. All dated horizons were assigned a depth uncertainty of 0.1 m. Age-depth positions were predicted every 0.1 m. Ten-thousand age-depth paths were generated using Markov Chain Monte Carlo methods to produce a median age-depth model, and an associated 95% highest density interval (HDI) or credible interval, the Bayesian analogue to 2σ uncertainty (Trayler et al., 2020).

From five to 13 single-zircon U-Pb CA-ID-TIMS dates were obtained from each tuff. From zircons that yield concordant ²⁰⁶Pb/²³⁸U and ²⁰⁷Pb/²³⁵U dates (Fig. 4), we focus on the ²⁰⁶Pb/²³⁸U dates because they are more precise and accurate (e.g., Mattinson, 2005). Firehole tuff sample 2017-59-FC shows no age dispersion, whereas other samples reveal only minor age dispersion (Table 1). For each, the youngest group of concordant ²⁰⁶Pb/²³⁸U dates that passed the modified Thomson Tau outlier rejection criterion (p = 0.05) were selected for inclusion in a weighted mean that represents the closest approach to the eruption and deposition age of the tuff bed. Crystals returning older dates are excluded as inherited from earlier Paleogene volcanic rocks during magma genesis and/or eruption. Uncertainties for calculated weighted means are reported in the form ± X(Y)[Z] where X is the analytical contribution, Y is the combined analytical and tracer (EARTHTIME 2535; Condon et al., 2015; McLean et al., 2015) uncertainty, and Z combines analytical, tracer, and decay constant uncertainties (Table 1; Fig. 5; Table S2).

Single sanidine ⁴⁰Ar/³⁹Ar dates were obtained from the Main and Grey tuffs (Table 1, Fig. 5), with only one distinct outlier in each set of dates. The large uncertainties in the sanidine dates compared to the U-Pb dates from the same ash beds are due, in part, to the small size of the individual sanidine crystals fused (~75–100 µm). Single-crystal fusions are preferred in this case due to the prior concerns over inheritance and reworking in the multi-crystal analyses of Smith et al. (2008). Uncertainties for calculated weighted means are reported in the form ± X(Y)[Z] where X is the analytical contribution, Y is the combined analytical and J-value uncertainty, and Z combines analytical, J-value, and decay constant uncertainties (2σ; Jicha et al., 2016).

Of 96 zircon crystals analyzed by LA-ICP-MS, one (1%) returned a Cretaceous ²⁰⁶Pb/²³⁸U age and was excluded from subsequent analyses. Eight grains were selected for CA-ID-TIMS analysis. All eight are concordant and have indistinguishable isotope ratios with a weighted mean ²⁰⁶Pb/²³⁸U age of 51.576 ± 0.007 (0.017) [0.058] Ma (mean square weighted deviate [MSWD] = 0.71; n = 8), that we interpret to reflect the time since the eruption and deposition of this tuff.

Of 73 zircon crystals analyzed by LA-ICP-MS, 14 (19%) returned a Cretaceous ²⁰⁶Pb/²³⁸U age and were excluded from subsequent analyses. Five grains were selected for CA-ID-TIMS analysis. The four youngest are concordant and have indistinguishable isotope ratios with a weighted mean ²⁰⁶Pb/²³⁸U age of 51.286 ± 0.029 (0.039) [0.067] Ma (MSWD = 0.23; n = 4), which is interpreted to reflect the time since the eruption and deposition of this tuff.

Of 123 zircon crystals analyzed by LA-ICP-MS, 74 (60%) returned Cretaceous ²⁰⁶Pb/²³⁸U ages and were excluded from subsequent analyses. Eight grains were selected for CA-ID-TIMS analysis. The four youngest zircon dates are concordant and have indistinguishable isotope ratios with a weighted mean ²⁰⁶Pb/²³⁸U age of 51.171 ± 0.009 (0.018) [0.058] Ma (MSWD = 0.65; n = 4) that we interpret to reflect time since the eruption and deposition of this tuff.

Of 190 zircon crystals analyzed by LA-ICP-MS, 79 (42%) returned Cretaceous or older ²⁰⁶Pb/²³⁸U ages and were excluded from subsequent analyses. Initial pilot analysis of eight zircon crystals from the Grey tuff found a bimodal distribution of ages, with a younger population at ca. 50.86 Ma and an older population at ca. 50.90 Ma. An additional six zircons were selected for CA-ID-TIMS analysis, focusing on acicular euhedral crystals with dark outer rims, like those which produced the younger age population in the pilot analyses. Of the 10 crystals analyzed, the five youngest zircons are concordant and have indistinguishable isotope ratios with a weighted mean ²⁰⁶Pb/²³⁸U age of 50.864 ± 0.008 (0.027) [0.061] Ma (MSWD = 0.19; n = 5). We interpret the latter to be the best U-Pb age for eruption and deposition of this tuff.

Single-crystal laser-fusion experiments were performed on 30 sanidine crystals. The youngest population of 29 crystals yields a weighted mean age of 50.912 ± 0.041 (0.052) [0.078] Ma, with an MSWD of 0.87. We interpret this as a best ⁴⁰Ar/³⁹Ar estimate of the age of eruption and deposition of this tuff. For the Grey tuff and Main tuff (discussed next) the ²⁰⁶Pb/²³⁸U and ⁴⁰Ar/³⁹Ar ages are concordant with one another using the Y uncertainties, reflecting in part the large uncertainty in the ⁴⁰K total decay constant of Min et al. (2000). We note that were we to use the more precise ⁴⁰K total decay constant proposed by Renne et al. (2011), the ²⁰⁶Pb/²³⁸U and ⁴⁰Ar/³⁹Ar ages would be discordant for these two tuffs.

Of 69 zircon crystals analyzed by LA-ICP-MS, 19 returned Cretaceous or older ²⁰⁶Pb/²³⁸U ages and were excluded from subsequent analyses. Eight zircon crystals were selected for CA-ID-TIMS analysis. The seven youngest zircon dates are concordant and have indistinguishable isotope ratios with a weighted mean ²⁰⁶Pb/²³⁸U age of 50.114 ± 0.022 (0.034) [0.063] Ma (MSWD = 0.31; n = 7), which is interpreted as the best U-Pb age for eruption and deposition of this tuff.

Single-crystal laser-fusion experiments were performed on 34 sanidine crystals. The youngest population of 33 crystals yield a weighted mean age of 50.133 ± 0.038 (0.048) [0.078] Ma, with an MSWD of 1.10. We interpret this as a best ⁴⁰Ar/³⁹Ar estimate of the age of eruption and deposition of this tuff.

Of 23 zircon crystals analyzed by LA-ICP-MS, eight (35%) returned Cretaceous or older ²⁰⁶Pb/²³⁸U ages and were excluded from subsequent analyses. Eight zircon crystals were selected for CA-ID-TIMS analysis. The six youngest zircon dates are concordant and have indistinguishable isotope ratios with a weighted mean ²⁰⁶Pb/²³⁸U age of 49.925 ± 0.009 (0.027) [0.060] Ma (MSWD = 1.28; n = 6), which is interpreted to reflect the time since the eruption and deposition of this tuff.

Of 48 zircon crystals analyzed by LA-ICP-MS, nine (19%) returned Cretaceous or older ²⁰⁶Pb/²³⁸U ages and were excluded from subsequent analyses. Fifteen zircon crystals were selected for CA-ID-TIMS analysis. The seven youngest zircon dates are concordant and have indistinguishable isotope ratios with a weighted mean ²⁰⁶Pb/²³⁸U of 49.769 ± 0.016 (0.030) [0.061] Ma (MSWD = 1.12; n = 7), which is interpreted to reflect the time since the eruption and deposition of this tuff.

Prior work by Smith et al. (2008) and Machlus et al. (2015) identified a potential discrepancy between the U-Pb and ⁴⁰Ar/³⁹Ar ages for the Sixth, Layered, and Main tuffs. While sanidine ages tend to skew younger than zircon U-Pb, due to the difference in closure temperature in these mineral systems, Machlus et al. (2015) found U-Pb ages for these three tuffs that were instead younger than the associated ⁴⁰Ar/³⁹Ar ages from Smith et al. (2008). Machlus proposed this discrepancy may be the result of alteration in the sanidine crystals due to deposition in a hypersaline environment, and that the air abrasion methods used by Smith et al. (2008) may not have sufficiently removed altered rims from the sanidine crystals analyzed.

To resolve this discrepancy, we employed several methods to improve the precision of our U-Pb and ⁴⁰Ar/³⁹Ar ages. Unlike prior work by Smith et al. (2008), we analyzed single crystals, rather than multiple crystal aliquots; this approach eliminates possible contamination by inherited crystals, which can skew the resulting ages older. Additionally, we performed ⁴⁰Ar/³⁹Ar analyses on a Nu Noblesse multi-collector mass spectrometer, allowing us to measure multiple isotopes simultaneously. Furthermore, we used cathodoluminescence imaging when selecting zircons for U-Pb analysis to identify crystals free of inclusions, and we screened zircons based on their LA-ICP-MS age to eliminate inherited cores and crystals from analysis. We also employed a double spike tracer to further improve precision relative to the single-spike tracer used by Machlus et al. (2015).

In the case of the Sixth tuff, our U-Pb age of 49.769 ± 0.016 (0.030) [0.061] Ma is intermediate between the Smith et al. (2008),⁴⁰Ar/³⁹Ar and Machlus et al. (2015) U-Pb ages, and agrees with both when total uncertainty is considered (Fig. 5). Our U-Pb ages of 49.925 ± 0.009 (0.027) [0.060] Ma for the Layered tuff and 50.114 ± 0.022 (0.034) [0.063] Ma for the Main tuff coincide more closely with the Machlus et al. (2015) U-Pb ages, as does our ⁴⁰Ar/³⁹Ar age of 50.133 ± 0.038 (0.048) [0.078] Ma for the Main tuff, suggesting the younger ages for these tuffs are likely more accurate. Our ages for the remaining tuffs agree with both the Machlus et al. (2015) and Smith et al. (2008) ages within uncertainty (Fig. 5).

In summary, advances in analytical methods and instrumentation have led to improved isotopic data and we consider the results presented here to represent the most accurate age determinations for these tuffs.

We developed three age-depth models, one each for cores BF-1, CCR-1, and WM (Fig. 6). Each model is congruent and shows an apparent decrease in accumulation rate between the upper and lower portions of the WPM, with the transition occurring somewhere at or above marker bed C, at ca. 51 Ma (Fig. 6). In core CCR-1, located on the southern margin of the basin, the average accumulation rate decreases from 229 mm/kyr between the Firehole and Boar tuffs, to 158 mm/kyr between the Boar and Sixth tuffs; a 31% decrease in accumulation rate. The overall average accumulation rate for the WPM in core CCR-1 is 173 mm/kyr, and agrees remarkably well with a model of astronomical forcing of sedimentation in this same core which supports an accumulation rate of 169.5 mm/kyr (Meyers, 2008). In core BF-1, which is centrally located in the basin, the average accumulation rates for the same sections are 211 and 141 mm/kyr, respectively, a 33% decrease in accumulation rate. These accumulation rates agree with those inferred by Smith et al. (2010) and Machlus et al. (2015) in the BF-1 core. In core WM, near the northern margin of the basin, the accumulation rate decreases 35%, from 183 to 118 mm/kyr.

The observed reduction in accumulation rate after ca. 51 Ma in all three cores coincides with the cessation of faulting along the northern flank of the Uinta uplift (Fig. 7). Uplift along the Uinta fault would have caused flexural subsidence of a southern basin trough (Roehler, 1992), resulting in increased accommodation in the Bridger Basin and relatively rapid accumulation of sediment, including bedded evaporite. As uplift ceased, the reduction in flexural subsidence reduced accommodation space, accumulation rates decreased, and shallowing of the lake prevented the deposition of bedded evaporites (Smith et al., 2008).

To test the association between the alluvial marker beds and eccentricity observed by Smith et al. (2010), we compared the predicted ages for the tops of the marker beds (Table S5) with the ages of peak eccentricity from the orbital solutions of Laskar et al. (2011) and Zeebe and Lourens (2019). The tops of the marker beds were chosen for this comparison because these transitions represent relatively abrupt flooding surfaces (Aswasereelert et al., 2013; Smith et al., 2014) (Fig. 8), suggesting that they are more likely to be isochronous across the basin than the bases of the marker beds, which commonly exhibit evidence of scour (Smith et al., 2015). The upper flooding surfaces of the marker beds are therefore interpreted as optimal targets for comparison with astronomical solutions and the marine record.

Of the nine marker beds, the upper six display a close association between the median predicted age of the top of the marker bed and eccentricity maxima in both the Laskar et al. (2011) and Zeebe and Lourens (2019) astronomical solutions (Fig. 9). Of note, all six marker bed flooding surfaces that coincide with eccentricity maxima occur in succession and were deposited after the observed shift in sedimentation rate at ca. 51 Ma. To quantify this association, we employed the Rayleigh test for phase directionality, using the rayleigh.test function of the R package “circular” (Agostinelli and Lund, 2017) and the hilbert and resample functions of the R package “astrochron” (Meyers, 2014). We adopt as a formal null hypothesis (H₀) that alluvial-lacustrine transitions are independent of eccentricity maxima, and the alternative hypothesis (H₁) that alluvial-lacustrine transitions are paced by eccentricity maxima. The phase of eccentricity associated with these paleoenvironmental transitions is determined using the Bayesian age-depth model of core BF-1 to determine the median age for the top of each marker bed, which is then used to sample the phase of eccentricity in the astronomical solution, determined via Hilbert transform. The results for the La10b eccentricity solution (Laskar et al., 2011) are shown in Figure 10, since it was used to tune the benthic N. truempyi δ¹³C and δ¹⁸O records of Ocean Drilling Program Site 1263 of Lauretano et al. (2018), and therefore permits direct comparison with the marine record. Other astronomical solutions are also evaluated (Fig. 11).

The observed eccentricity phase angles produce a Rayleigh test value R = 0.6835, with a p-value of 0.0012, thus we reject H₀ (Fig. 10). This strong correlation between an untuned sedimentary record and a key theoretical eccentricity solution convincingly supports the findings of Smith et al. (2010, 2014), Aswasereelert et al. (2013), and Machlus et al. (2015) that alluvial marker bed deposition was paced by short eccentricity. When this test is applied only to the upper six consecutive marker beds, the Rayleigh test value improves to R = 0.8098 (p-value of 0.0013), suggesting that cessation of uplift along the Uinta fault and the accompanying reduction of accommodation space rendered the basin more susceptible to astronomical forcing, or more capable of preserving such signals (Fig. 10).

Our model also indicates a relationship specifically between marker bed deposition and rising eccentricity. Eccentricity primarily influences Earth’s climate system through amplitude modulation of precession, potentially driving variations in seasonality, precipitation, and physical and chemical weathering across different phases of eccentricity (e.g., Zachos et al., 2010; Birgenheier et al., 2020). Recent findings in early Eocene marine strata suggest that the carbon isotope excursion associated with the Paleocene-Eocene thermal maximum and earlier carbon cycle perturbations were triggered by climate changes and light carbon release to the oceans that reflect the rise to maxima in short eccentricity (Piedrahita et al., 2022).

In the Bridger Basin, our model suggests that increasing eccentricity enhanced delivery of lowstand alluvial sediment onto the desiccated lake floor, possibly due to increased seasonality, “monsoon-like” precipitation, and increased physical and chemical weathering (e.g., Zachos et al., 2010; Birgenheier et al., 2020). Alluvial sediments dominated deposition in the basin under rising eccentricity until reaching peak eccentricity, during which further intensification of seasonality appears to have shifted the hydrologic budget of the Bridger Basin, altering accommodation by moving the locus of alluvial sedimentation toward the basin margins and drowning alluvial deposits nearer to the basin center. This interpretation contrasts with Smith et al. (2014), who argue that lacustrine deposits of the WPM, including bedded evaporites that record expansion and contraction of the lake, reflect greater seasonality and precipitation fluctuations expected during periods of high eccentricity. Accordingly, Smith et al. (2014) propose that climatically drier and more stable periods associated with eccentricity minima led to deposition of the alluvial marker beds.

Our findings and interpretations differ in detail from those of Smith et al. (2014) in large part because new U-Pb zircon ages are significantly more precise than previous geochronology in Machlus et al. (2015) and Smith et al. (2008, 2010), and the ages are integrated using Bayesian age-depth modeling. However, despite these advances, further work is needed to refine the pacing of alluvial marker bed deposition during maxima in orbital eccentricity. This need for further refinement is apparent in the observation that, whereas the median predicted ages of the alluvial marker bed tops in our Bayesian model—the statistically most likely ages—show a strong association with eccentricity maxima, Monte Carlo simulations using the full 95% HDI envelope of the Bayesian model produce R values closer to ~0.3 with reduced statistical significance (Fig. 11). From this degraded significance, we conclude that the conservative error envelope resulting from the extrapolation between dated tuffs in our Bayesian age-depth modeling algorithm is likely over-estimated. U-Pb ages from tuffs that have yet to be dated, coupled with additional, more precise, ⁴⁰Ar/³⁹Ar dates from several tuffs, are needed to reduce the uncertainty envelope of the Bayesian model, and thus better constrain the tempo and genesis of the alluvial marker beds and their relationship with the phase of eccentricity.

²⁰⁶Pb/²³⁸U zircon and ⁴⁰Ar/³⁹Ar sanidine ages from tuffs in the WPM are concordant and sufficiently precise to resolve discrepancies between the two chronometers noted in previous work. The Bayesian approach integrates these clocks into a single age-depth model, with a two- to ten-fold increase in temporal resolution compared to previous radioisotopic age models for the WPM. This model reveals basin-scale changes in sedimentation rate coincident with cessation of faulting along the Uinta uplift at ca. 51 Ma, suggesting a uniform and rapid reduction in flexural subsidence as a likely driver for reduced sediment accumulation. Our new age model also permits a direct comparison between an untuned terrestrial sedimentary record and theoretical astronomical solutions, which suggests a close association between the deposition of nine alluvial marker beds and the phasing of orbital eccentricity. In our proposed model, flooding surfaces atop six of these alluvial beds correspond with short eccentricity maxima. Confirming this terrestrial response of climate and the hydrologic cycle to increasing eccentricity, during the Eocene greenhouse, is within reach of Bayesian age-depth models that integrate additional radioisotopic ages of tuffs. Our findings thus motivate further efforts to correlate between terrestrial and marine records of responses to astronomically driven climate change.

This contribution represents Bruck’s MSc thesis at the University of Wisconsin–Madison. We thank James Crowley, Debbie Pierce, and Bryan Wathen for assistance in the labs, Robin Trayler for guidance regarding the modified Bchron package, and Mike Smith for his comments on, and discussion of, the manuscript. Comments of two anonymous reviewers propelled us to clarify many points and are greatly appreciated. This research was supported by National Science Foundation–Division of Earth Sciences grants 1813278 (Carroll, Singer, Meyers) and 1813088 (Schmitz), an ExxonMobil/Geological Society of America graduate research grant (Walters), an American Association of Petroleum Geologists Foundation grant-in-aid (Walters), and the Department of Geoscience at the University of Wisconsin–Madison.