Temperature seasonality in the North American continental interior during the Early Eocene Climatic Optimum

Paleogene greenhouse climate equability has long been a paradox in paleoclimate research. However, recent developments in proxy and modeling methods have suggested that strong seasonality may be a feature of at least some greenhouse Earth periods. Here we present the first multi-proxy record of seasonal temperatures during the Paleogene from paleofloras, paleosol geochemistry, and carbonate clumped isotope thermometry in the Green River Basin (Wyoming, USA). These combined temperature records allow for the reconstruction of past seasonality in the continental interior, which shows that temperatures were warmer in all seasons during the peak Early Eocene Climatic Optimum and that the mean annual range of temperatures was high, similar to the modern value (∼ 26 C). Proxy data and downscaled Eocene regional climate model results suggest amplified seasonality during greenhouse events. Increased seasonality reconstructed for the early Eocene is similar in scope to the higher seasonal range predicted by downscaled climate model ensembles for future high-CO2 emissions scenarios. Overall, these data and model comparisons have substantial implications for understanding greenhouse climates in general, and may be important for predicting future seasonal climate regimes and their impacts in continental regions.


Introduction
The Paleogene was the last major greenhouse period in Earth's history and is characterized by extreme warming events and resultant biological shifts (e.g., Greenwood and Wing, 1995;Wilf, 2000;Zachos et al., 2001Zachos et al., , 2008;;McInerney and Wing, 2011), with prolonged warmth during the Early Eocene Climatic Optimum (EECO) peaking from roughly 52 to 50 Ma (e.g., Zachos et al., 2008;Hyland et al., 2017).The early Eocene in general is thought to represent a warm and "equable" global climate state with high mean annual temperatures (MATs; e.g., Wilf, 2000;Zachos et al., 2008), low mean annual range of temperatures (MARTs; e.g., Wolfe, 1978Wolfe, , 1995;;Greenwood and Wing, 1995), and low pole-to-Equator temperature gradients (LTGs; e.g., Spicer and Parrish, 1990;Greenwood and Wing, 1995;Evans et al., 2018).While high MATs during the Eocene now seems well established, the feasibility of "equable" conditions defined by low MARTs and low LTGs is still in question as a result of increasingly complex global climate models that are unable to reproduce such conditions (e.g., Barron, 1987;Sloan and Barron, 1990;Sloan, 1994;Huber and Caballero, 2011;Lunt et al., 2012).
Recent proxy work on Paleogene warm intervals and hyperthermals such as the Paleocene-Eocene Thermal Maximum (PETM) has suggested that continental interiors may maintain higher or near-modern MARTs during these periods, implying that the "low seasonality" aspect of climate equability may not be reasonable under all greenhouse conditions (e.g., Snell et al., 2013;Eldrett et al., 2014).Despite this suggestion, it remains unclear whether proxy estimates from other basins, regions, and greenhouse periods can be reconciled with the range of feasible conditions provided by climate model studies.Quantitative reconstructions of seasonality (MART) based on precise proxy estimates of MAT, Published by Copernicus Publications on behalf of the European Geosciences Union.E. G. Hyland et al.: Temperature seasonality in the North American continental interior warm month mean temperature (WMMT), and cold month mean temperature (CMMT) could help to resolve some of these model-proxy discrepancies by providing a robust and well constrained set of seasonal observations for comparison to available climate model outputs.Robust proxy reconstructions of seasonality are crucial for understanding this aspect of past greenhouse equability (Lunt et al., 2012;Snell et al., 2013;Peppe, 2013).
Seasonality estimates have previously been made using a variety of proxy paleothermometers in isolation, and can now be made with higher confidence using recently developed methods that target each of these individual temperature parameters: MAT can be estimated using a paleosol geochemistry-based thermometer, WMMT can be estimated using the carbonate clumped isotope ( 47 ) thermometer, and CMMT can be estimated using a nearest living relative (NLR) floral coexistence thermometer.The bulk majorelement geochemistry of modern soils has been used to quantify the effects of weathering processes via a wide range of geochemical indices (see Sheldon and Tabor, 2009).The relationship between modern climate parameters like temperature and indices such as salinization (Sheldon et al., 2002), the paleosol weathering index (Gallagher and Sheldon, 2013), and the paleosol-paleoclimate model (Stinchcomb et al., 2016) has led to the development of climofunctions for MAT that have been used to estimate paleo-MAT during the Cenozoic (e.g., Retallack, 2007;Takeuchi et al., 2007;Bader et al., 2015;Stinchcomb et al., 2016).The clumped isotope ( 47 ) thermometer is based on the temperature-dependent relative enrichment of multiply substituted isotopologues of CaCO 3 ( 13 C 18 O 16 O 2 ) within the solid carbonate phase, which is independent of the isotopic composition of the water in which the carbonate precipitated (e.g., Ghosh et al., 2006;Eiler, 2007).For pedogenic carbonates in temperate regions, this growth temperature is linked to mean warm season soil temperatures (e.g., Quade et al., 2013;Hough et al., 2014), and has been used to estimate paleo-WMMT during the Cenozoic (e.g., Snell et al., 2013;Garzione et al., 2014).The NLR coexistence method has been developed based on the sensitive and highly conserved collective modern cold temperature tolerances of related floras to calculate cold month temperatures (e.g., Wolfe, 1995;Mosbrugger and Utescher, 1997).Those relationships have been refined and used to estimate quantitative paleo-CMMT during the Cenozoic (e.g., Greenwood et al., 2005Greenwood et al., , 2017;;Thompson et al., 2012;Eldrett et al., 2014;Utescher et al., 2014).
Here we employ a multi-proxy approach using paleosol geochemistry, clumped isotope, and floral NLR coexistence thermometry methods from the same localities in order to address seasonality in the past, specifically applying it to the issue of early Eocene greenhouse equability in the North American continental interior.We estimate MAT, WMMT, and CMMT throughout the EECO including both defined peak (∼ 51 Ma) and non-peak conditions (e.g., Hyland et al., 2017), and compare the resultant proxy estimates of temperature seasonality (MART) to the modern climate state of the region, as well as to downscaled climate model predictions of temperature seasonality during the Eocene and for future emissions scenarios.

Methods
The targeted early Eocene locality is the Green River Basin (GRB) in southwestern Wyoming (USA; Fig. 1).The GRB sequence is composed of a series of terrestrial clastic rocks deposited during the early Eocene and EECO as a result of Laramide synorogenic fluvial and lacustrine sedimentation along the margin of endorheic paleo-lake Gosiute (e.g., Clyde et al., 2001;Smith et al., 2008Smith et al., , 2010Smith et al., , 2015)).Contemporaneous multi-proxy records of peak and non-peak conditions during the EECO are from the interfingering Wasatch Formation, primarily fluvial sandstones and paleosols of the Ramsey Ranch and Cathedral Bluffs members, and Green River Formation, primarily lacustrine shales and carbonates of the Wilkins Peak Member (Fig. 1).The paleosols and pedogenic carbonates were sampled from the Honeycomb Buttes near South Pass, Wyoming (42.24 • N, 108.53 • W; Hyland and Sheldon, 2013), while the floral assemblages were sampled from the Latham coal (41.68  W) outside Rock Springs, Wyoming (Fig. 1; Wilf, 1998Wilf, , 2000)).

Paleosol geochemistry
The bulk major-element geochemistry of modern soils (specifically B horizons) has been used extensively to develop a number of composition-climate relationships, including those predicted by the paleosol-paleoclimate model (PPM 1.0 ), which relates a broad suite of major element compositions to mean annual temperature (among other factors) at the site of soil formation (Stinchcomb et al., 2016).Stinchcomb et al. (2016) developed this nonlinear spline model using the largest available geochemical dataset from 685 modern soils across North America in order to derive proxy relationships between 11 major and minor oxides and MAT.This new proxy is calibrated over a wider range of climatic conditions, soil types, and parent materials than other available proxies (cf.Sheldon et al., 2002;Gallagher and Sheldon, 2013), and has been validated via independent comparisons in both modern climosequences (Stinchcomb et al., 2016) and Miocene paleosols (Driese et al., 2016).Following associated procedures, our bulk paleosol samples from selected upper Bt horizons of defined Alfisols (described in detail by Hyland and Sheldon, 2013) were prepared for majorelement geochemistry by cleaning and grinding to a ho-

NALMA
Floral collection (Wilf, 2000) Paleosol sampled (Hyland & Sheldon, 2013) Map area  2015) and Hyland and Sheldon (2013).LY is the Lysitean, BF is the Blacksforkian, LU is the Luman Member, NT is the Niland Tongue, TM is the Tipton Member, WPM is the Wilkins Peak Member, LA is the Laney Member, RR is the Ramsey Ranch Member, and CB is the Cathedral Bluffs Member.
mogenous powder.Samples were analyzed using lithium borate fusion preparation and X-ray fluorescence (XRF) measurements at the ALS Chemex Laboratory (Vancouver, BC, Canada), where analytical uncertainty for analyses was maintained at less than 0.1 % for all elements, and replicate analyses had a mean standard deviation of 0.8 % (Table S1 in the Supplement).Resultant major-and minor-element data were not corrected for loss-on-ignition (e.g., Stinchcomb et al., 2016), and were input into the open-access PPM 1.0 model, which produces "low", "best" and "high" MAT estimates; we present the "high" estimates as MAT here (see Sect. 4.1 for explanation; Table S1).Broadly, soil geochemical proxies are consistent with other paleoclimate proxies (e.g., paleobotanical; Sheldon and Tabor, 2009, and references therein), and are more robust to diagenetic alteration under a wide variety of burial conditions (Hyland and Sheldon, 2016).
Pedogenic carbonate nodules from selected Bk horizons (paleosol depths ∼ 20-240 cm) were sliced into thin sections and analyzed under transmitted light and cathodoluminescence to identify primary micritic carbonate (Fig. 2), which was microdrilled and homogenized for clumped isotope ( 47 ) analysis.Extremely shallow (< 50 cm) or deep (> 200 cm) carbonates were analyzed to specifically examine temperature depth profiles in paleosols (Fig. 2), while pedogenic carbonates from commonly sampled depths (50-200 cm; e.g., Cerling, 1984;Koch, 1998;Zamanian et al., 2016) were used for calculating and interpreting paleotemperature records.Powdered samples and carbonate standards were analyzed in replicate at the University of Washington's IsoLab, following methods of Burgener et al. (2016) and Kelson et al. (2017), which are modified after Huntington et al. (2009) and Passey et al. (2010).Briefly, CO 2 is produced from 6 to 8 mg of pure carbonate reacted in a common phosphoric acid bath (∼ 105 % H 3 PO 4 ) at 90 • C. Evolved CO 2 is then cleaned via passage through a series of automated cryogenic traps and a cooled (−20 • C) Poropak Q column using helium carrier gas through a nickel and stainless steel vacuum line, and the purified CO 2 is transferred to Pyrex break seals.Each sample is then analyzed on a Thermo MAT253 mass spectrometer equipped with an automated 10port tube cracker inlet system and configured to measure m/z 44-49, using data acquisition methods and scripts presented by Schauer et al. (2016).
All analyses include an automatically measured pressure baseline (He et al., 2012), are corrected using heated gas (1000 • C; Huntington et al., 2009) and CO 2 -water equilibration (4, 60 • C) lines during the corresponding analysis period (Table S2), and are reported in the absolute reference frame (Dennis et al., 2011).Following recent work (Daëron et al., 2016;Schauer et al., 2016), mass spectrometer data are corrected using the 17 O correction values recommended by Brand et al. (2010).Carbonate standards for these analyses include international standards NBS-19 and ETH-2, as well as internal standards C64 and COR, which are all reported relative to VPDB (δ 13 C, δ 18 O) and absolute refer-ence frame ( 47 ) in Table S2.All samples were analyzed in replicate (3-5) to minimize standard analytical error, and data were reduced following Schauer et al. (2016).Carbonate growth temperatures (T [ 47 ]) were calculated using the most current and extensive inorganic calcite calibration (Kelson et al., 2017), which was produced using the updated 17 O correction values of Brand et al. ( 2010) and is consistent with our analytical methods.Based on preliminary comparisons, the Kelson et al. (2017) calibration produces results that are not significantly different from data calculated using previous calibrations at moderate Earth-surface temperatures (∼ 20-40 • C; Daëron et al., 2016;Cedric John and Matthieu Daëron, personal communication, 2016; Table S2).
Fossil assemblages were selected from the literature (e.g., Wilf, 1998Wilf, , 2000) ) based on temporal fit, floristic diversity, and reliable taxonomy.Fossil taxa were each attributed to a modern taxon based on NLR (e.g., MacGinitie, 1969;Hickey, 1977;Manchester and Dilcher, 1982;Wolfe and Wehr, 1987;Wing, 1998;Wilf, 1998Wilf, , 2000;;Manchester, 2014;SIMNHP, 2015), with unattributed or disputed placements assigned conservatively at higher taxonomic levels (Table S4).Climatic envelopes of modern groups in North America and Asia were retained for the ancient taxa based on environmental niche conservation (e.g., Wang et al., 2010;Fang et al., 2011).Modern taxa distributions (GBIF, 2016) were linked to high-resolution gridded climatic maps (Hijmans et al., 2005) to extract MAT, WMMT and CMMT using the Dismo Package in the R statistical program (R Core Team, 2013).Prior to calculating climatic ranges, plant distribution coordinate files were scrutinized for (1) plants with dubious taxonomic assignments, as not all identifications were rigorous and not all collected specimens were taxonomically assigned by experts (only species-level identifications are included); (2) plants occurring outside of their natural ranges, as many plants occur outside their adapted environment due to agricultural or aesthetic translocation; and (3) redundant occurrences, as many duplicate coordinates or researcher entries exist for the same taxon and their inclusion may skew results toward given localities.
Quantitative paleotemperatures were estimated using a modified bioclimatic analysis approach (e.g., Greenwood et al., 2005Greenwood et al., , 2017;;Thompson et al., 2012;Eldrett et al., 2014).Overlap ranges of climatic tolerances for coexisting species from each assemblage were defined by calculating probability density functions of those climatic envelopes (Fig. 3 and Table S5) consistent with recent work (e.g., Thompson et al., 2012;Harbert and Nixon, 2015;Grimm and Potts, 2016;Greenwood et al., 2017).In order to avoid inclusion of apparent coexistence intervals in which no modern occurrence is recorded, we calculate the collective probability density of taxa co-occurrence for each combination of MAT (x), WMMT (y), and CMMT (z): Calculations are repeated such that the likelihood (f ) is calculated for each climatic combination, for each taxon (t), dependent on the number of taxa (n), using the mean and standard deviation of each taxon (Table S5).Climate input parameters were individual occurrence data points (∼ 32 000) derived from GBIF (2016), excluding combinations unlikely to represent the climatic envelope of the taxa in the assemblage by calculating a maximum likelihood probability density function that defines a precise estimate of temperature parameters with a low standard deviation for each selected assemblage (Fig. 3).

Modern climate data and model downscaling
The  (Sewall and Sloan, 2006;Thrasher andSloan, 2009, 2010); the specifics of the model configurations can be found in Thrasher and Sloan (2009).Those results were averaged for the final 20 years of the model run at equilibrium and calculated over the same study area (40.5-43 • N by 107-110.5 • W) by integrating data across grid cells monthly for each model year within the above-defined Green River Basin (e.g., Snell et al., 2013).This particular set of Eocene model configurations was chosen because it allows for the highest available resolution over the basin domain using the best available set of boundary conditions (cf., EoMIP; Lunt et al., 2012).All modern climate normals and model downscaling results are reported in Table S7.

Results
PPM 1.0 statistical model results for MAT from these paleosol samples range from 13.5 to 17.6 Uncertainty for these estimates is reported as the root mean squared error of the model fit regression (±2.5 • C).Petrographic observation of carbonate nodules from all depths and selected soils identified dominantly micritic textures with minor components of subangular quartz grains and occasional sparry (> 20 µm) calcite veins and cements; however, we were able to identify and microsample unaltered fine-grained (< 5 µm) calcite material in each of the examined samples (n = 14; Fig. 2).Clumped isotope 47 values for these samples range from 0.582 to 0.631 ‰ (µ = 0.607 ‰; σ = 0.014 ‰), which corresponds to an estimated WMMT range of 18 to 34 Uncertainty for these estimates is reported as propagated error from analytical and equilibrated CO 2 reference frame uncertainty (negligible); replicate standard error  peak conditions into the peak EECO (∼ 51 Ma), after which temperatures decreased back to lower values (Fig. 4).Modern climate normals averaged monthly for the GRB range from −8.4 to 18.1 • C, with a MAT of 4.4 • C (Table S7).Downscaled Eocene climate model results averaged monthly for the GRB range from 4 to 24 • C (LoCO) and 6 to 30 • C (HiCO), with MATs of 13 and 16 • C, respectively (Table S7).Downscaled future climate model results averaged monthly for the GRB range from −5.0 to 20.4 • C (RCP4.5) and −2.9 to 24.7 • C (RCP8.5), with MATs of 7.1 and 10.6 • C, respectively (Table S7).Monthly temperature trends maintain roughly the same shape for modern observational data, future model estimates, and Eocene model estimates.However, the Eocene-modeled cases show substantially higher winter temperatures, and in both modern and Eocene-modeled cases the higher emission or pCO 2 scenario shows an enhanced summer signal relative to the lower emission or pCO 2 scenario from the same time period (Fig. 5).

Temperature estimates
Temperature estimates from the PPM 1.0 spline model are based on specifically selected uppermost B horizons of paleosols with comparable parent materials.These horizons were selected based on previous work describing and sampling paleosols from the Cathedral Bluffs Member in the GRB (Fig. 1; Hyland and Sheldon, 2013), and based on the characteristics of soils sampled for the paleosol paleoclimate model dataset (Stinchcomb et al., 2016), in order to generate the most robust input data for the PPM 1.0 spline model.While the PPM 1.0 model produces multiple possible estimates of paleo-MAT, the estimate shown to be most reliable via concurrent comparisons with other paleotemperature methods (paleobotanical and paleosol proxies) is the "high MAT" value we present here (Michel et al., 2014;Stinchcomb et al., 2016;Driese et al., 2016).We further justify our use of the "high" estimate because the PPM 1.0 training dataset heavily samples soils from temperate regions (specifically the conterminous USA), which tend to have lower MAT (≤ 10 • C) and therefore could place excess weight on low values in the model predictive space.This sampling bias likely produces the demonstrated pattern of "best" MAT predictions generally exhibiting positive residuals (Stinchcomb et al., 2016), which means that the PPM 1.0 model would be more likely to skew temperature estimates from paleosol and other modern samples toward lower-than-observed MAT values.The presented MATs appear to coincide with a statistical mean between CMMT and WMMT estimates (Fig. 4), and also agree within uncertainty with independent MAT estimates from other types of paleosol geochemistry (salinization index, δ 18 O; Hyland and Sheldon, 2013) and broadly with updated physiognometric (Table S6; Wilf, 2000) and coexistence analysis paleobotanical estimates from the GRB (Fig. 4).
Based on the assessment of physical and isotopic data, our sampled pedogenic carbonate nodules appear to be primary records of Earth surface temperatures at the time of their formation.All sampled nodules preserve micritic carbonate, and transmitted light and cathodoluminescence images show limited recrystallization or void-filling spar and no evidence of pervasive remineralization (Fig. 2).Clumped isotopic data also suggest primary and uncontaminated carbonate material; 48 values remain low ( 1 ‰; Table S2), indicating a lack of hydrocarbon or sulfide contamination (e.g., Guo and Eiler, 2007;Huntington et al., 2009).Temperature and δ 18 O measurements remain well within the range of reasonable terrestrial values, particularly for continental interior basins with seasonal climates (Table S2; e.g., Quade et al., 2013;Hough et al., 2014).Carbonates forming in temperate regions often exhibit summer or warm-month temperatures due to warm, dry conditions and low soil CO 2 concentrations during those months (e.g., Breecker et al., 2009;Quade et al., 2013).Such conditions are predicted for the GRB during the early Eocene based on regional climate models (Thrasher andSloan, 2009, 2010), and are evident in paleosol features (Clyde et al., 2001;Hyland and Sheldon, 2013) as well as evaporative δ 18 O of source waters from nearby paleo-lakes Gosiute and Uinta (Table S1; e.g., Sarg et al., 2013;Frantz et al., 2014).Further warm biasing of soil temperature with respect to air temperature can be imparted by radiant ground heating, but such effects are likely negligible in shaded forest soils (e.g., Quade et al., 2013;Ringham et al., 2016).Clumped isotope data from two soil depth profiles collected in the GRB agree within uncertainty below ∼ 50 cm (Fig. 2), suggesting that surface heating and depth attenuation of surface temperature variability does not significantly affect the samples used for our MART reconstructions (paleosol depths ∼ 50-200 cm; e.g., Ringham et al., 2016).
These results imply that the temperatures measured from our pedogenic carbonates broadly reflect WMMT of soil as observed in other records (e.g., Peters et al., 2013;Hough et al., 2014;Burgener et al., 2016).Possible exceptions are two samples at the base of the Honeycomb Buttes section (HB-109 and HB-18; Table S2), which appear to correspond to MAT estimates from the same paleosols (PPM 1.0 ; Fig. 4).These lowest temperature estimates from the base of the section may be artificially "cool" as a function of seasonal precipitation regimes spreading carbonate formation across other parts of the year, particularly in soils with deeper Bk horizons like these (e.g., Gallagher and Sheldon, 2016).Because of the likely bias toward MAT in these two samples, we exclude them from calculations of WMMT or MART as indicated in Fig. 4; additionally, this effect means that all of our clumped isotope-based estimates of WMMT may be artificially low, suggesting that our calculated MART values could represent a minimum value.However, our resultant clumped isotope-based temperature estimates are mostly in agreement with both regional climate model predictions of summer month air temperatures (e.g., Thrasher and Sloan, 2009;Snell et al., 2013) and paleobotanical coexistence estimates of WMMTs (Fig. 4).
Paleobotanical coexistence methods have been shown to reconstruct paleo-temperatures robustly, particularly for warm and cold months in well-sampled and taxonomically rich localities such as these (e.g., Thompson et al., 2012;Grimm and Potts, 2016).However, uncertainties may be larger than accounted for by the described statistical methods applied to these assemblages because (1) many fossil classifications within the GRB assemblages are not directly comparable to or identifiable as extant species, and coexistence analyses at a generic or familial level may introduce bias by broadening the temperature tolerance ranges of most groups (e.g., Wang et al., 2010); and (2) evolutionary or climatic preferences of Paleogene fossil taxa may not be fully conserved in extant groups, introducing potential sources of error (e.g., Fang et al., 2011).If we double the estimated error to account for these unquantifiable uncertainties, the collective coexistence probability density functions from these assemblages still produce CMMT, MAT, and WMMT estimates defined by narrow "maximum likelihood" bioclimatic envelopes (< ±3 • C; Fig. 3; Table S5), which suggest that the environmental characteristics of these fossil assemblages are well constrained despite some higherlevel NLR assignments.Additionally, sampling bias from well-sampled temperate regions (e.g., North America) in the modern GBIF (2016) database may place undue weight on the cool end of plant ranges (e.g., Greenwood et al., 2017), constraining paleotemperature estimates to lower values or smaller ranges than is appropriate.This suggests that, similar to clumped isotope-based estimates, our plant-based MART values could also represent a minimum value.Despite this, paleobotanical coexistence CMMT estimates agree with regional climate model predictions of winter month temper-   5.9 atures in the GRB (e.g., Thrasher andSloan, 2009, 2010), MAT estimates agree broadly with multiple paleosol-based proxy estimates (Fig. 4; Hyland and Sheldon, 2013) and with updated paleobotanical physiognomy estimates (Fig. 4; Table S6; Wilf, 2000), and WMMT estimates agree with regional climate model estimates (e.g., Thrasher and Sloan, 2009;Snell et al., 2013) and broadly with clumped isotopebased estimates (Fig. 4).Taken together these proxy results paint a consistent picture of Earth-surface temperatures during the early Eocene, despite uncertainties inherent in each individual method.

Temperature seasonality
Because each of these proxies appears to represent different seasonal temperatures robustly, we combine these estimates to produce a new multiply constrained investigation of paleo-MART.By calculating the differences between CMMT from paleobotanical coexistence analysis, MAT from paleosol geochemistry or paleobotanical analyses, and WMMT from 47 composition or paleobotanical coexistence analysis, we can directly estimate MART in the past and compare differences in seasonal temperatures independent of calculation method (cf., Snell et al., 2013).In other words, our approach can define MART as (1) the difference between WMMT and CMMT or (2) twice the difference between MAT and either WMMT or CMMT, assuming that MAT falls half way between those estimates by definition (Table 1).Because our approach can calculate MART using both methods and an average of multiple proxies, this allows for a wide range of independent checks on our estimates, providing the most robust available paleo-MART (Table 1).Each method provides consistent answers that are statistically indistinguishable for a given time period (Student's t test p values = 0.4-0.9),lending confidence to calculations which show that MART ranged from 21-26 • C during the early Eocene (Table 1).
By averaging all data from each population (CMMT, MAT, WMMT) for the peak and non-peak intervals separately, calculated MARTs suggest that seasonality was generally slightly lower than modern across parts of the early Eocene (∼ 21-23 • C, non-peak), but appears to have increased to near-modern ranges during the peak EECO (∼ 26 • C; Fig. 4; Table 1).The calculated uncertainty in the difference between these populations (standard error of a difference) is ∼ 4 • C, which makes the non-peak and peak intervals statistically distinct though nearly overlapping.Overall, this sug-gests that not only is seasonality not reduced during greenhouse periods (e.g., Snell et al., 2013) but it may actually be expanded (Fig. 5).
Estimates from the lower end of our reconstructed MART range are still higher than MART estimates from individual paleobotanical proxies (15-18 • C; e.g., Greenwood and Wing, 1995;Wolfe et al., 1998), but compare favorably to estimates from regional climate models with assumed lacustrine or paludal land cover (20-22 • C; Thrasher and Sloan, 2010).However, estimates from the higher end of the reconstructed MART range compare more favorably to modeled MART values with assumed woodland or forested land cover (24-26 • C; Thrasher and Sloan, 2010;Snell et al., 2013).The transient nature of paleo-lake Gosiute and the variable evolution of environments within the GRB throughout the early Eocene is well documented in stratigraphic archives, indicating that the basin may have been alternately dominated by the paleo-lake or by forested floodplains during this period (Smith et al., 2008(Smith et al., , 2014)).In this context, our results suggest that both lower (though still in excess of any previous paleobotanical estimates) and higher MART states may in fact be reasonable for this region at different points during the early Eocene as the GRB evolved.Moreover, proxy and modeling work does not appear to be contradictory, instead having captured different portions of the range of possible MART values indicated for the peak vs. non-peak EECO in this part of the continental interior (Figs. 4 and 5).Regardless, these results suggest that MART values lower than ∼ 20 • C (e.g., Greenwood and Wing, 1995;Wolfe et al., 1998) may be unreasonable during any part of the EECO, even in the context of variable climate and environmental conditions.This is particularly true because MART estimates using these proxy methods are more likely to underestimate than overestimate seasonality (see Sect. 4.1).

Seasonality implications
Our new proxy data and model comparisons have important implications for continental climates, as they suggest two potential characteristics of seasonality in interior regions during warming events: (1) proxies tend to indicate continental temperatures on the high end of modeled ranges in all seasons and (2) both proxies and regional models indicate that summer temperatures may increase disproportionately, actually broadening MART, at high atmospheric pCO 2 .While proxy and model estimates of paleotemperature generally agree through the early Eocene in the GRB, proxy estimates consistently fall in the top half of all modeled values (Fig. 5).Although these model and proxy results are not statistically distinct, they may suggest that realistic environmental responses could have a skewed distribution within the range of model-predicted climate outcomes, an observation which has been made previously for other regions and time periods (e.g., Roe and Baker, 2007;Diffenbaugh and Field, 2013).
Winter temperatures were generally high during the Eocene (Figs. 4 and 5; e.g., Greenwood and Wing, 1995), but during the peak EECO summer temperatures appear to have increased disproportionally, broadening the range of MART (Figs. 4 and 5).While this apparent trend may be related to the lack of direct CMMT estimates during the peak EECO, the consistency of MART estimates using both reconstruction methods (see Sect. 4.2) suggests the observation is robust.Regional Eocene climate model output for the GRB predicts lower MART (∼ 20 • C) under low pCO 2 conditions (LoCO scenario), and higher MART (∼ 24 • C) under high pCO 2 conditions (HiCO scenario; Fig. 5; Table S7).Therefore, a theoretical transition from lower (≤ 500 ppm) to higher (≥ 1000 ppm) atmospheric pCO 2 during the peak EECO (e.g., Hyland and Sheldon, 2013;Jagniecki et al., 2015) could effectively broaden MART and result in extreme summer temperatures during that period, which would be consistent with both regional model and proxy predictions in the GRB (Fig. 5).Regional model-proxy agreement on the plausibility of variable moderate-to-high MART (20-26 • C) in continental interiors fits with global simulations employing a reasonable set of radiative forcings and climate sensitivities, which project similar seasonality ranges during this and other greenhouse events (Huber and Caballero, 2011;Lunt et al., 2012).These temperature seasonality estimates also corroborate recent work on other regions and warm periods (e.g., Snell et al., 2013;Eldrett et al., 2014), and further support the interpretation that continental interiors were less "equable" than previously thought under greenhouse conditions (Snell et al., 2013;Peppe, 2013).
Increased seasonality and the disproportionate response of summer temperatures during greenhouse climates also has significant implications for predicting future change in continental interiors.Current projections for the next century using downscaled global climate model ensembles (PCDMI, 2014; Table S7) indicate generally increased temperatures and changing seasonality in North America, and GRB temperatures are projected to increase particularly during winter months (Fig. 5).However, for high emissions scenarios that may be closer in character to greenhouse conditions like the peak EECO or the PETM (RCP8.5;e.g., IPCC, 2007;Lunt et al., 2012), summer temperatures in the GRB increase more strongly, broadening MART (Fig. 5; Table S7).This trend in MART from peak EECO proxy data and highemission or pCO 2 model simulations in both the future and Eocene suggests a potential atmospheric pCO 2 threshold for enhanced seasonality, and provides support for models and observations indicating that continental interiors may experience more extreme seasonality in the future under heightened greenhouse conditions (e.g., IPCC, 2007;Diffenbaugh and Field, 2013;Diffenbaugh et al., 2017).The mechanism for producing this increased seasonality remains unclear and requires further study in terms of both proxy applications and model development, although changes in land cover may play a crucial role at least in regional variability (Thrasher and Sloan, 2010;Diffenbaugh and Field, 2013).

Conclusions
Estimates of winter (paleofloral NLR coexistence), mean (paleosol geochemistry), and summer (clumped isotope) temperatures from the early Eocene in the Green River Basin of Wyoming (USA) provide new multi-proxy constraints on seasonality (mean annual range of temperature) in terrestrial settings during greenhouse periods.These records show that MART was variable but near (or above) modern values during the Early Eocene Climatic Optimum, confirming that both seasonality in continental interiors may not remain constant and that EECO conditions likely do not conform to at least the seasonality aspect of greenhouse "equability".Comparisons between proxy data and regional or downscaled climate models further imply that temperature seasonality may respond differently at low vs. high atmospheric pCO 2 .Overall, this suggests that our understanding of past greenhouse climates in continental interiors may be incomplete when it comes to "equability", and proposes the potential for extreme seasonality in these regions during past warming events and in the future, which likely has important implications for natural ecosystems and human infrastructure.

Figure 1 .
Figure 1.Map and stratigraphy of the Green River Basin.(a) Map of the region, showing major sedimentary basins and topographic highs.Stars show proxy record sampling sites (paleosols in yellow, paleoflora in red), and the dashed box is the sampling region for modern climate stations and the downscaling domain for both models.CF is the Cordilleran fold-thrust belt, UU is the Uinta uplift, WR is the Wind River uplift, OC is the Owl Creek uplift, GM is the Granite Mountains, and FR is the Front Range.(b) Simplified stratigraphy of the central to eastern GRB, showing facies for the Green River Formation (GRF) and the equivalent and interfingering Wasatch Formation (WF) based on the work of Smith et al. (2015) andHyland and Sheldon (2013).LY is the Lysitean, BF is the Blacksforkian, LU is the Luman Member, NT is the Niland Tongue, TM is the Tipton Member, WPM is the Wilkins Peak Member, LA is the Laney Member, RR is the Ramsey Ranch Member, and CB is the Cathedral Bluffs Member.
modern temperature dataset was derived from 1981 to 2010 averaged climate normals from the National Oceanic and Atmospheric Administration (NOAA) weather observation stations within the Green River Basin (n = 18; NCDC, 2010), defined as the area 40.5-43 • N by 107-110.5 • W (Fig.1).Future model temperature projection results used a 10-model ensemble from the Coupled Model Intercomparison Project Phase 5 (CMIP5) under standard low (RCP4.5)and high (RCP8.5)emissions scenarios (IPCC, 2014); the specifics of each model and configuration are available from the World Climate Research Programme (2011).Results were averaged monthly for the final 10 years of the model run (2090-2099) and calculated over the same study area using standard bias-correction and spatial downscaling (BCSD) methods developed by PCDMI (2014).Eocene model temperature results used data from a modified threedimensional regional climate model (RegCM3;Sewall and Sloan, 2006;Pal et al., 2007) with established Eocene boundary conditions including low (560 ppm; LoCO) and high (2240 ppm; HiCO) atmospheric pCO 2 scenarios

Figure 3 .
Figure 3. Floral methods description.(a) Probability density functions of hypothetical Taxa A and Taxa B along climatic variable X to form a probability density function representative of the maximum likelihood of co-occurrence.(b) Hypothetical climatic envelope of Taxon Q with climatic variables X and Y , where point R occurs outside the envelope of Taxon Q but within its range of both variables (creating a false inclusion of point R).(c) Probability density function distributions for seasonal temperatures from sampled paleofloral sites, where arrows indicate calculated mean temperatures for each parameter, and n is the number of morphotypes included in each assemblage.
Figure 4. Temperature proxy estimates of CMMT (white plot), MAT (gray plot), and WMMT (black plot) through the early Eocene.Triangles represent paleobotanical coexistence estimates, squares represent paleosol geochemistry estimates, stars represent revised paleobotanical physiognomy estimates, and circles represent clumped isotope estimates.Error bars represent probability density function 2σ (paleobotanical coexistence), root mean squared error (PPM 1.0 paleosol geochemistry), calibration standard error (paleobotanical physiognomy), and propagated analytical and calibration error (clumped isotopes).Shading highlights peak EECO conditions based on previous work (e.g.,Hyland et al., 2017), the long dashed line highlights possible aliasing due to a long sampling interval, and the short dashed line highlights exclusion of two clumped isotope data points (see Discussion).Estimates of peak EECO (51±0.5 Ma) and non-peak EECO MART are defined as described in Table1and the Discussion, with MAT shown by vertical lines.Modern MART and MAT are from averaged climate normals for NOAA weather stations in the GRB (NCDC, 2010).
• N, 107.88 • W), Sourdough coal (41.91 • N, 108.00 • W), Niland Tongue (41.06 • N, 108.77 • W), and Little Mountain quarry (41.28 • N, 109.30 Ma) and non-peak EECO MART are defined as described in Table1and the Discussion, with MAT shown by vertical lines.Modern MART and MAT are from averaged climate normals for NOAA weather stations in the GRB (NCDC, 2010).

Table 1 .
Comparison of Eocene MART estimates using different constraining temperatures and calculation methods.
WMMT − CMMT.MARTC = (MAT − CMMT) × 2. MARTW = (WMMT − MAT) × 2. a Average of all available temperature proxy data across indicated time interval.b Average MART estimate for each calculation method, number in parentheses is the SD of calculation group.