Stable isotope evidence for long-term stability of large-scale hydroclimate in the Neogene North American Great Plains

The Great Plains of North America host a stark climatic gradient, separating the humid and well-watered eastern US from the semi-arid and arid western US, and this gradient shapes the region's water availability, its ecosystems, and its economies. This climatic boundary is largely set by the influence of two competing atmospheric circulation systems that meet over the Great Plains – the wintertime westerlies bring dominantly dry air that gives way to moist, southerly air transported by the Great Plains low-level jet in the warmer months. Climate model simulations suggest that, as CO₂ rises, this low-level jet will strengthen, leading to greater precipitation in the spring but less in the summer and, thus, no change in mean annual precipitation. Combined with rising temperatures that will increase potential evapotranspiration, semi-arid conditions will shift eastward, with potentially large consequences for the ecosystems and inhabitants of the Great Plains. We examine how hydroclimate in the Great Plains varied in the past in response to warmer global climate by studying the paleoclimate record within the Ogallala Formation, which underlies nearly the entire Great Plains and provides a spatially resolved record of hydroclimate during the globally warmer late Miocene. We use the stable isotopes of oxygen (δ¹⁸O) as preserved in authigenic carbonates hosted within the abundant paleosol and fluvial successions that comprise the Ogallala Formation as a record of past hydroclimate. Today, and coincident with the modern aridity gradient, there is a sharp meteoric water δ¹⁸O gradient with high (−6 ‰ to 0 ‰) δ¹⁸O in the southern Great Plains and low (−12 ‰ to −18 ‰) δ¹⁸O in the northern plains. We find that the spatial pattern of reconstructed late Miocene precipitation δ¹⁸O is indistinguishable from the spatial pattern of modern meteoric water δ¹⁸O. We use a recently developed vapor transport model to demonstrate that this δ¹⁸O spatial pattern requires air mass mixing over the Great Plains between dry westerly and moist southerly air masses in the late Miocene – consistent with today. Our results suggest that the spatial extents of these two atmospheric circulation systems have been largely unchanged since the late Miocene and any strengthening of the Great Plains low-level jet in response to warming has been isotopically masked by proportional increases in westerly moisture delivery. Our results hold implications for the sensitivity of Great Plains climate to changes in global temperature and CO₂ and also for our understanding of the processes that drove Ogallala Formation deposition in the late Miocene.

1 Introduction

The Great Plains of North America rise slowly in elevation from the wooded lowlands near sea level west of the Mississippi River to the greater than 4 km high “purple mountain majesties above the fruited plain” (Bates, 1911) that comprise the North American Cordillera. Though this region contains some of the flattest landscapes in the United States (Fonstad et al., 2007; Dobson and Campbell, 2014), these plains belie a remarkable climatic setting and geologic history which have conspired to shape the modern-day water resources, ecosystems, and economies of the Plains region. Boreal spring heralds the onset of the Great Plains low-level jet (GPLLJ), a primarily nocturnal, southerly jet responsible for transporting more than 30 % of the water vapor that enters the continental US every year (Helfand and Schubert, 1995). Interactions between the GPLLJ and midlatitude storm systems yield some of the largest and most intense convective systems on the planet (Song et al., 2019). These spring and summertime rains nurture the vast grasslands of the Plains, which return this moisture to the atmosphere via transpiration, often seeding additional precipitation on subsequent days and resulting in one of the tightest couplings between land and atmosphere anywhere on Earth (Koster, 2004). On multi-annual timescales, these same interactions – modified by long-term climatic oscillations such as the El Niño–Southern Oscillation or the North Atlantic Oscillation – are the proximal cause for the extensive floods and deep droughts that frequent this region (Byerle, 2003), perhaps best exemplified by the 1930s Dust Bowl (Schubert, 2004), an event which reshaped American governance and society and has been labeled the worst environmental catastrophe in US history (Egan, 2006).

Climatically, the Great Plains are characterized by a sharp aridity gradient with a more humid climate to the east and a more arid climate to the west (Fig. 1). This aridity gradient spans the 100th meridian and marks a dramatic change in the long-term patterns of precipitation, vegetation, and, consequently, agriculture and human development (Powell, 1879, 1890; Webb, 1931; Seager et al., 2018b). A convenient measure of aridity is the aridity index (AI), which is the ratio of precipitation (P) to potential evapotranspiration (PET): conditions are considered arid – or water-limited – if P $/$ PET is less than 1. The Great Plains today straddle the transition between the wet, eastern US (AI > 1) and the water-limited western US (AI < 1) (Seager et al., 2018b). Consequently, regions to the east are more densely populated and farms rely on rainfed agricultural practices; to the west, settlement is more limited and agriculture relies extensively on groundwater withdrawals or irrigation diversions.

Figure 1 (a) Extent of the continuous Ogallala Formation (thick black outline), with new sites presented in this study (squares) and sites compiled from the PATCH Lab (circles) (Kukla et al., 2022a). Shading is elevation (m). (b) Map of average annual precipitation (mm) from the Global Precipitation Climatology Project (Meyer-Christoffer et al., 2011). (c) Aridity index (AI), which is calculated as P $/$ PET. P as in panel (b); PET is from the Global Land Evaporation Amsterdam Model (GLEAM) (Miralles et al., 2011; Martens et al., 2017). (d) Interpolated distribution of δ¹⁸O of modern meteoric waters, derived from rivers and streams, groundwater, tap water, spring water, and precipitation δ¹⁸O measurements and compiled from Waterisotopes.org (Waterisotopes Database, 2019). Thick black outline in all panels is the areal extent of the continuous Ogallala Formation.

This groundwater – sourced from the Ogallala Aquifer (black polygon outlined in Fig. 1) and replenished by spring and summertime rains – is hosted in the Ogallala Formation, a nearly continuous formation that underlies the Plains from southwest Texas and southeast New Mexico into the southern part of South Dakota, making it one of the largest and most laterally continuous sedimentary formations in North America. Comprised of sediments shed off the Rockies, the Ogallala Formation has been interpreted as a series of alluvial or telescoping megafans (Seni, 1980; Willett et al., 2018; Korus and Joeckel, 2023a) that prograded towards the east and completely buried the pre-existing erosional landscape that, in the southern Great Plains, is cut into Triassic and Permian bedrock, and, in the northern Great Plains, lies on the White River or Arikaree Group. In places, the Ogallala Formation is overlain by alluvial and/or eolian deposits, and, particularly in the southern Great Plains, there is frequently a prominent calcium carbonate caprock that separates the Ogallala from overlying sedimentary units (Gustavson and Winkler, 1988). Deposition began in the middle Miocene and ended in the late Miocene or earliest Pliocene. However, precisely why the Great Plains experienced a prolonged period of deposition, followed by a period of incision that continues to the present day, remains uncertain, with most studies attributing this period of deposition to dynamic topography effects associated with the passage of the Farallon Plate beneath the Great Plains (Moucha et al., 2009; Karlstrom et al., 2011; Willett et al., 2018).

During this time, substantial ecological changes occurred, largely yielding the pre-anthropogenic landscape that characterized the Quaternary Great Plains. Though grasslands were likely present before the middle Miocene, they continued expanding throughout the Plains during the cooling following the peak of mid-Miocene warmth (Jacobs et al., 1999; Stromberg, 2005; Strömberg and McInerney, 2011). These predominantly C₃ grasslands were progressively replaced by C₄ grasslands starting in the latest Miocene and continuing through the Pliocene (Fox and Koch, 2003, 2004). Largely coincident with these changes, large mammal diversity gradually declined on the Plains from its peak during the middle Miocene to a relative low in the Quaternary (Janis et al., 2000; Fritz et al., 2016). Since cessation of Ogallala deposition, a combination of base-level fall and uplift along the Front Range led to incision of the major Plains rivers through the Ogallala Formation (McMillan et al., 2002; Duller et al., 2012), leaving the former Ogallala landscape abandoned and perched above the major Plains valleys (Willett et al., 2018). Today, much of the Ogallala Formation is visible as a prominent escarpment protected by the indurated nature of its many calcic-rich sediments; dotting much of the length of this escarpment are wind turbines powered, to no small extent, by the exceptionally predictable and windy Great Plains low-level jet.

This combination of climate and geology has helped to promote development and agriculture on the Plains, with Ogallala Aquifer water supplying any deficit due to insufficient rainfall. However, the aquifer remains critically overdrawn in many places; further, anthropogenically driven increases in atmospheric CO₂ and associated changes in global climate may affect water availability on the Plains. Indeed, this area appears to be exceptionally sensitive to even small changes in climate and land cover due to the tight coupling between land and atmosphere in this region (Koster, 2004; de Noblet-Ducoudre et al., 2012; Laguë et al., 2019). Thus, small changes in precipitation and/or PET may shift the precise location of the “climatological 100th meridian” (i.e., the approximate location of the boundary between arid and wet ecosystems) with important consequences for ecosystems and agricultural systems. Further, global climate model (GCM) simulations tend to predict only a small change in summer precipitation, though a large shift in the seasonality of that precipitation is driven by dynamical shifts in the westerly jet and the GPLLJ (Cook et al., 2008; Bukovsky et al., 2017; Zhou et al., 2021a). Combined with rising temperatures, this negligible change in summer precipitation implies that the arid conditions characteristic of the western Great Plains may expand eastward with global warming as increases in PET outpace increases in P over the Plains (Seager et al., 2018b; Overpeck and Udall, 2020). However, models still struggle to properly simulate the GPLLJ and its associated precipitation-bearing convective systems. Changes in the interactions between vegetation, soil moisture, and rainfall at higher atmospheric CO₂ remain similarly difficult to model (Bukovsky et al., 2017; Zhou et al., 2021a). Further, paleoclimate data indicate that, in general, warmer periods have actually been wetter and/or greener (Caves et al., 2016; Burls and Fedorov, 2017; Ibarra et al., 2018; Feng et al., 2022), in conflict with many model predictions of drier future conditions (e.g., Scheff, 2018). As a consequence, there is substantial uncertainty regarding how aridity on the Plains – and the interactions between the GPLLJ, precipitation, and vegetation – will change as atmospheric CO₂ rises and global temperature increases.

In this contribution, we take advantage of the remarkable spatial extent of Neogene sediments afforded by the Ogallala Formation to understand how changes in global climate impacted the Plains during periods of higher atmospheric CO₂ and warmer global temperatures. Mid-Miocene global temperatures were 7–8 °C warmer than today and atmospheric CO₂ levels were around 500 ppm or higher (Herbert et al., 2016; Steinthorsdottir et al., 2021; CenCO2PIP Consortium, 2023). Since the mid-Miocene, temperature and atmospheric CO₂ have gradually declined, establishing the bipolar glaciation that characterizes Quaternary climate. The Neogene therefore provides an opportunity to answer the question of how aridity has changed on the Plains and whether the climatological 100th meridian shifted eastward in the past as global climate models suggest for the future. The spatial extent of the Ogallala also permits us to examine the mechanisms by which any shifts in aridity may have occurred due to, for example, shifts in the distribution of precipitation as a result of strengthening or weakening of the GPLLJ in a warmer climate. Further, any shifts in climate may point towards the mechanisms that generated the sediment that formed the Ogallala Formation, permitting distinctions between climatically or tectonically controlled generation of the Ogallala Formation (i.e., Molnar and England, 1990; Zhang et al., 2001).

To answer these questions, we rely on the stable isotopes of authigenic carbonates – a material that is particularly abundant within the fluvial and paleosol facies of the Ogallala Formation and which is thought to record precipitation δ¹⁸O. Precipitation oxygen isotopes (δ¹⁸O_p) are sensitive to the moisture source and rainout history of air masses that reach a given location. Because of this sensitivity, spatially resolved datasets of carbonate oxygen isotopes (δ¹⁸O_c) interpreted to track δ¹⁸O_p have been successfully used to constrain how large-scale hydroclimate has changed in response to atmospheric and orographic forcing through time (Fox and Koch, 2004; Mix et al., 2011; McDermott et al., 2011; Kocsis et al., 2014; Caves et al., 2015; Caves Rugenstein and Chamberlain, 2018). The Plains today host a steep δ¹⁸O gradient oriented NE–SW that roughly tracks, but is somewhat oblique to, the aridity gradient at the 100th meridian (Fig. 1c). Given this concurrence, we hypothesize that the Ogallala authigenic carbonate δ¹⁸O_c record over space and time will reflect changes in the large-scale hydroclimate and aridity over the Great Plains. In the following sections, we explain how δ¹⁸O can be used to track hydroclimate and our approach to sampling and to building a spatially extensive dataset. We then investigate the climatic drivers behind the modern δ¹⁸O_p gradient over the Plains and compare this to reconstructed maps of δ¹⁸O_p during Ogallala deposition. Lastly, we apply a recently developed reactive transport model to quantitatively assess our observations of δ¹⁸O and how these data relate to overall hydroclimate on the Plains. We find that large-scale atmospheric circulation and the position of the climatological 100th meridian have been remarkably stable, despite correspondingly large changes in global climate since the late Miocene.

2 Background

2.1 Modern hydroclimate of the North American Great Plains

Spring and summer precipitation in the Great Plains overwhelmingly originates from the Gulf of Mexico, where the Bermuda High drives southeasterly flow while a high-pressure system over the Great Basin effectively blocks Pacific moisture. The combination of these two high-pressure systems creates the conditions for the Great Plains low-level jet, which carries moisture deep into the interior of the North American continent (Helfand and Schubert, 1995). Though the precise mechanisms that generate the GPLLJ are complex (see Shapiro et al., 2016), critically for our purposes, both the existence of the high topography of the North American Cordillera and the east–west-sloping terrain of the Great Plains appear crucial for generating and maintaining a low-level jet over the Plains. For example, in model simulations where cordillera topography is removed, the GPLLJ is weakened or non-existent (Ting and Wang, 2006; Jiang et al., 2007). During winter, the GPLLJ is inactive and precipitation originates largely from the Pacific Ocean due to storms routed by the midlatitude westerly jet. These wintertime storms traverse the wide expanse of topography that comprises the North American Cordillera – including ranges such as the Sierra Nevada, Cascades, and Wasatch – which removes moisture from these midlatitude cyclones. Consequently, precipitation across the Plains typically occurs during interaction with cold Arctic air masses that can penetrate as far south as the southern Great Plains (Nativ and Riggio, 1990; Brubaker et al., 1994). Thus, the combination of topography with atmospheric circulation generates much of the seasonal precipitation pattern that prevails today over the Great Plains.

The position of the 100th meridian aridity boundary is shaped by the relative strength of these two circulation systems (Seager et al., 2018b), yielding a sharp humidity gradient that approximately coincides with the boundary between southerly maritime Gulf of Mexico air and dry continental air (Hoch and Markowski, 2005). On an interannual basis, this boundary can shift depending upon the strength of the midlatitude westerlies relative to the GPLLJ; with increasing westerly wind strength, for example, the 100th meridian shifts east (Hoch and Markowski, 2005). On longer timescales, winter aridity and weaker westerlies have been linked to grassland expansion and forest dieback in the Great Plains in the early Miocene (Kukla et al., 2022b). Besides these two large-scale atmospheric systems, several other factors determine the location and orientation of the sharp aridity gradient that characterizes the Great Plains today (Seager et al., 2018b). First, evapotranspiration of water from the land surface to the atmosphere is critical in determining precipitation. For example, return of water to the atmosphere from the land surface (largely through transpiration) may supply up to 40 % of the precipitation in the Great Plains during spring and summer (Burde et al., 2006; van der Ent et al., 2010). The importance of the land surface may be further enhanced by the widespread grasslands that populate the Plains landscape. Grasses can much more rapidly modify their stomatal conductance and, hence, total transpiration than can trees and shrubs to take advantage of periodic rainstorms (Ferretti et al., 2003; Hetherington and Woodward, 2003). Such rapid water use leads to higher recycling rates of water from the land surface back to the atmosphere. As a consequence, the spread of grasslands onto the Plains during the Miocene has been hypothesized to have fundamentally increased the recycling of water between the land surface and the atmosphere (Mix et al., 2013; Chamberlain et al., 2014).

Model simulations project distinct changes in precipitation and hydroclimate associated with dynamical and thermodynamic responses to warming. GCMs and regional climate models robustly predict that precipitation seasonality will shift from the summer to the spring (Cook et al., 2008). As global temperatures rise, the westerly jet shifts poleward, permitting a stronger and more northerly GPLLJ, producing more precipitation in the late spring (Bukovsky et al., 2017; Zhou et al., 2021a). In contrast, the continued northward shift of the GPLLJ as summer progresses weakens the jet over the Great Plains, leading to enhanced late summer drying (Zhou et al., 2021a). Despite this shift in the timing of the wet season, mean annual precipitation is not expected to change (Bukovsky et al., 2017). With constant mean annual precipitation, the increase in PET owed to rising temperatures will decrease the AI and shift the climatological 100th meridian eastward (Seager et al., 2018a).

2.2 Precipitation oxygen isotopes reflect this hydroclimatic pattern

This annual mixing between the GPLLJ and the westerlies results in a steep spatial gradient in precipitation δ¹⁸O mostly due to the differences in topography traversed by each air mass (Kendall and Coplen, 2001). Mountain ranges tend to increase the net loss of moisture from an air mass, preferentially removing ¹⁸O and decreasing δ¹⁸O_p (Rozanski et al., 1993; Poage and Chamberlain, 2001; Winnick et al., 2014; Mix et al., 2019; Kukla et al., 2019). Thus, precipitation derived from the westerlies has low δ¹⁸O values by the time it reaches the Great Plains. In contrast, GPLLJ moisture, which has not traversed major topographic barriers and is augmented by a high degree of evapotranspiration that replenishes the GPLLJ, is about 10 ‰ higher than equivalent westerly moisture (Mix et al., 2013; Winnick et al., 2014). The aridity gradient, which depends on the relative contributions of westerly vs GPLLJ moisture, is therefore encoded in the value of δ¹⁸O_p across the Great Plains, which can generally be understood as the precipitation-weighted average of the end-member sources. This precipitation δ¹⁸O signal is captured by authigenic carbonates but is further modified by both the temperature of carbonate formation and other potentially spatially variable factors associated with mineral formation, such as differences in precipitation seasonality and evaporation (Breecker et al., 2009; Caves, 2017; Huth et al., 2019; Kelson et al., 2020, 2023).

3 Geological setting and sampling approach

The sediments of the Ogallala Formation originate from the Rocky Mountains, and eroded material from the Miocene Rockies was transported by braided, high-energy, ephemeral streams and eolian processes across the Plains (Joeckel et al., 2014; Smith and Platt, 2023; Korus and Joeckel, 2023b). The result was the Ogallala Formation, which spans from South Dakota to southern Texas (black outline in Fig. 1). Though the headwaters of these rivers and fans have since been eroded away, except in southern Wyoming, discontinuous remnants of the Ogallala have been mapped nearly up to their sources in eastern New Mexico (Frye et al., 1982). Consisting of gravel, sand, silt, and clay deposits, the Ogallala Formation also contains abundant calcic paleosols and calcic-rich sediments distributed throughout the formation and across the entire N–S extent of the formation (Gustavson, 1996; Joeckel et al., 2014; Smith et al., 2016; Smith and Platt, 2023). The Ogallala Formation in the southern Great Plains unconformably overlies Permian through Cretaceous strata (Gustavson and Holliday, 1999). In the northern Great Plains, the Ogallala Formation is underlain by the late Oligocene to early Miocene Arikaree Group and the White River Group, upper Eocene to Oligocene in age. Chronostratigraphy of the Ogallala Formation is based on fossil vertebrate faunas of Barstovian to Hemphillian North American land mammal ages (NALMAs) and scattered volcanic ash beds and basalt flows (Kitts, 1965; Leonard and Frye, 1978; Frye et al., 1978; Thomasson, 1979; Winkler, 1985; Schultz, 1990; Swisher, 1992; Gustavson, 1996; Tedford, 1999; Tedford et al., 2004; Cepeda and Perkins, 2006; Smith et al., 2016, 2018). Due to the nature of the largely fluvial and eolian deposits, the Ogallala exhibits substantial heterogeneity north to south. Consequently, different workers have classified the Ogallala Formation as a group (in the northern Great Plains) (Tedford et al., 2004) or as a formation, primarily in Kansas and to the south (Gustavson, 1996; Ludvigson et al., 2009). Herein, we refer to the Ogallala exclusively as a formation. However, in the northern Great Plains, there are further distinct formations such as the Valentine, Ash Hollow, and Olcott formations (Joeckel et al., 2014; Smith et al., 2017), each with defined age constraints based upon biostratigraphy and ashes (Tedford et al., 2004). In the southern Great Plains, previous workers have proposed elevating the Ogallala Formation to group status based upon subdividing the Ogallala into the Bridwell and Couch formations (Winkler, 1985; Gustavson and Winkler, 1988). However, in the southern Great Plains, we adopt the terminology of Gustavson (1996), who concluded that these formations are difficult to map and contain little dateable material, suggesting that the Ogallala remain with formation status. The thickness of Ogallala sediments generally varies relative to the underlying topography between 250 m, in regions where it fills paleovalleys, and 10–30 m in the interfluves between paleovalleys (Gustavson and Holliday, 1999). In the southern Great Plains, there is frequently an erosion-resistant caliche or caprock calcrete that separates the Ogallala Formation from the predominantly eolian Plio-Pleistocene Blackwater Draw Formation and, locally, the lacustrine Blanco Formation (Gustavson, 1996; Gustavson and Holliday, 1999). These caliche caprocks are thought to have developed during one or multiple periods of extended landscape stability and likely record a multi-genetic history (Brock and Buck, 2009; Henry, 2017). In contrast, in the northern Great Plains, the Ogallala Formation is overlain by several high-energy deposits, including the Crooked Creek Formation in Kansas and the Broadwater Formation in Nebraska (Swinehart and Diffendal, 1987).

To capture spatial changes in precipitation δ¹⁸O – and hence shifts in the aridity gradient – we collected samples spanning nearly the entire N–S and E–W extent of the Ogallala Formation from paleosol authigenic carbonate material to reconstruct paleo-precipitation δ¹⁸O. We build upon previous work (Fox and Koch, 2003, 2004) that developed a spatially extensive dataset of paleosol carbonate isotopes, collected primarily to understand changes in C₃ and C₄ vegetation during the late Miocene. We build upon these datasets, focusing on filling gaps in the southern and southwestern Great Plains (Texas and New Mexico), while also contributing additional data in the central and northern Great Plains.

Though the Ogallala Formation provides an unparalleled opportunity to collect spatially extensive δ¹⁸O_c data, the precise chronology of Ogallala deposition remains uncertain. In many places, specific formations have been dated using biostratigraphy and radiometrically dated ashes (Tedford et al., 2004); however, in many other places, particularly those Ogallala outcrops to the west and disconnected from the primary exposure of the Ogallala, temporal constraints are provided by lithologic correlations (Frye et al., 1982). Even in the formations which have been dated using biostratigraphy, the age constraints are relatively broad, typically limited by the precision of the North American land mammal ages. Lastly, the relatively thin veneer of Ogallala sedimentation combined with the potentially long time span covered by deposition suggests that there are frequent and temporally extensive unconformities within many sections (Smith and Platt, 2023). As a consequence, it remains difficult to correlate sections across the large expanse of the Great Plains.

For the purposes of our study, we use the published age constraints (Table 1) either from the study that originally analyzed the sampled section in detail or from later publications that provide a more precise age. Because we are only interested in broadly comparing samples spatially, we treat all samples in a given section as having the same age. We bound this age by considering the maximum possible age range for the given formation, which is either the Ogallala Formation or, in the northern Great Plains, one of the formations within the Ogallala Group. For example, in locations where we sample well-defined and dated subdivisions of the Ogallala Group (for example, samples from the Ash Hollow Formation at Lake McConaughy State Recreation Area in Nebraska; Joeckel et al., 2014), we consider our sample ages to be bound by the full age range of the formation. In other cases, we know only that our data lie above or below a certain dated ash (for example, samples from Wildcat Bluff Nature Center in Texas; Cepeda and Perkins, 2006), and therefore one bound on the age is provided by the dated material. In some cases, the sampled unit has only been correlated with the Ogallala Formation and there are no other age constraints to narrow the large possible range of ages. We therefore adopt the full possible range of ages given the identification of the unit as the Ogallala Formation. For example, for many sections in New Mexico, Ogallala outcrops have been identified by lithologic or geomorphic correlation with the contiguous body of the Ogallala Formation to the east, but no dateable material has been recovered (Frye et al., 1982). Despite these broad age constraints, we find that our results are not sensitive to the uncertainty in our correlations or to the precise chronology of Ogallala deposition because δ¹⁸O_c is largely invariant within any given section.

Table 1 Site-averaged data. Values of δ¹³C are reported in ‰ relative to VPDB. Values of δ¹⁸O are reported in ‰ relative to VSMOW. n is the number of samples collected from the site. Mean annual temperature data are retrieved from NARR (Mesinger et al., 2006).

Lastly, we also collected samples from the Miocene-age Tesuque Formation within the Santa Fe Group within the Rio Grande Rift. These samples are the only samples we collected west of the original spatial extent of the Ogallala Formation. Unlike the Ogallala Formation, the Santa Fe Group comprises thick basin fill shed off the Sangre de Cristo Range during ongoing rift extension (Galusha and Blick, 1971; Kuhle and Smith, 2001) and chronological constraints are provided by a combination of magnetostratigraphy (Barghoorn, 1981, 1985), biostratigraphy (Galusha and Blick, 1971; Tedford and Barghoorn, 1993; Aby et al., 2011), and radiometric dates (Mcintosh and Quade, 1996; Izett and Obradovich, 2001; Koning et al., 2013). The Tesuque Formation contains abundant authigenic carbonates, including nodules, root casts, groundwater cements, and laterally extensive Bk horizons (Kuhle and Smith, 2001). To place our samples within an existing stratigraphic framework, we collected samples along the Arroyo de los Martinez section and refer readers to Koning et al. (2013) for a detailed description of this section.

4 Methods

4.1 Stable isotope measurements

We collected carbonates (n=344) from paleosols from 32 distinct sites within the Ogallala Formation or from the Tesuque Formation. In every section, we sampled a wide variety of carbonate types, including rhizoliths, nodules, burrows, carbonate-cemented matrix samples, and caliches in order to test whether these different carbonate types reveal different spatial isotope patterns. Individual sample types for each sample are listed in Table S1 in Manser et al. (2023). During sampling, we first removed the weathered surface layer before selecting samples. Sampling sites are shown in Fig. 1a, sample types are depicted in Fig. 2, and coordinates are reported in Table 1. Samples were powdered for isotopic analysis with a Dremel tool or crushed with a mortar and pestle to obtain a homogeneous powder. Before powdering nodules, caliches, and matrix samples, any outer weathered rind was first discarded. Carbon and oxygen isotope ratios of the carbonate were measured with a ThermoFisher Gas Bench II coupled via a ConFlow IV interface to a Delta V Plus mass spectrometer at ETH Zürich following procedures described in detail in Breitenbach and Bernasconi (2011). Briefly, 140–300 µg of the powdered sample, depending on carbonate content, was reacted with five drops of 104 % phosphoric acid at 70 °C in He-flushed exetainers. Each batch of 79 samples included 16 replicates of the internal standards MS2 (δ¹³C $=+\mathrm{2.13}$ ‰, δ¹⁸O $=-\mathrm{1.81}$ ‰) and ETH-4 (δ¹³C $=-\mathrm{10.19}$ ‰, ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=-\mathrm{18.71}\mathrm{‰}$) interspersed throughout the run. The standards are used for drift corrections and data normalization and are calibrated to the international reference materials NBS 19 (δ¹³C $=+\mathrm{1.95}$ ‰, ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=-\mathrm{2.2}$ ‰) and NBS 18 (${\mathit{\delta }}^{\mathrm{13}}\mathrm{C}=-\mathrm{5.01}$ ‰, ${\mathit{\delta }}^{\mathrm{18}}\mathrm{O}=-\mathrm{23.00}$ ‰; Bernasconi et al., 2018). Analytical reproducibility of the standards was better than 0.1 ‰ (1σ for both δ¹³C and δ¹⁸O). We convert our δ¹⁸O_c data from VPDB to VSMOW using the equations in Brand et al. (2014), and we report all δ¹⁸O_c data relative to VSMOW.

Figure 2 Field photos representing the primary types of authigenic carbonate sampled in this study. (a) Carbonate nodules. (b) Burrows. (c) Root casts. (d) Caprock calcrete, as pictured in SE New Mexico.

4.2 HYSPLIT

Because δ¹⁸O_p is heavily influenced by the pathway that moisture takes to reach a certain site, we use NOAA's Hybrid Single-Particle Lagrangian Trajectory Model (HYSPLIT) (Draxler and Hess, 1998; Stein et al., 2015) to analyze the pathways by which moisture reaches the Great Plains. HYSPLIT is commonly used to understand modern precipitation δ¹⁸O_p data (Sjostrom et al., 2006; Bershaw et al., 2012; Li and Garzione, 2017; Zhu et al., 2018; San Jose et al., 2020) and to yield insights into the controls on reconstructed paleo-precipitation δ¹⁸O (Oster et al., 2012; Lechler and Galewsky, 2013; Caves et al., 2014, 2015; Wheeler and Galewsky, 2017; Zhu et al., 2018). To track air parcels back in time and space from a given location, we use North American Regional Reanalysis (NARR) data, which have a 32×32 km resolution (Mesinger et al., 2006), as the HYSPLIT climatological model input. To resolve spatial variability in the origins and pathways of storms, we simulate the origin and pathway of the air parcels for four selected sites (32° N and −97° E; 32° N and −105° E; 42° N and −97° E; 42° N and −105° E). At each site, we initialize the air parcel at 1000 m above ground. We choose this height as this level encapsulates much of the bulk moisture transported by the GPLLJ, which has maximum wind speeds between 500–1000 m above ground level (Jiang et al., 2007), and also captures moisture transport by the midlatitude westerlies (Lechler and Galewsky, 2013; Wheeler and Galewsky, 2017). The sites were chosen to encapsulate nearly the full latitudinal and longitudinal range represented by the carbonate stable isotope data. At each site, we generate nearly 53 000 back trajectories (i.e., one trajectory every 6 h from 1980–2016) and filter these trajectories for only those that are estimated to produce precipitation within 6 h of reaching the endpoint (e.g., Lechler and Galewsky, 2013; Caves et al., 2015). This results in approximately 4000 to 11 000 trajectories at each site. We further use HYSPLIT's built-in clustering algorithm to calculate the percentage of trajectories that originate from the Gulf of Mexico. Lastly, we note two critical assumptions regarding this HYSPLIT analysis. First, we assume that, by tracking air parcels using HYSPLIT, we are also tracking moisture; however, HYSPLIT does not account for moisture addition by evaporation or removal by precipitation and does not track diffusion of moisture into and out of an air mass. The assumption of moisture transport by advection is sometimes violated in regions of strong air mass mixing, such as in the Great Plains (Draxler and Hess, 1998); however, comparison of HYSPLIT results with the results of a more rigorous moisture tracking model – the Water Accounting Model (WAM-2layers) – generally shows close agreement (Driscoll et al., 2024). Second, our HYSPLIT results are strictly only applicable to understand the modern climate. Nevertheless, we use these results to develop insights into the controls on past precipitation δ¹⁸O.

4.3 Vapor transport model

Given the moisture pathways predicted by HYSPLIT, we use a one-dimensional vapor transport model (Kukla et al., 2019) to predict the isotopic composition of precipitation transported along these pathways. This model links spatial patterns of δ¹⁸O_p to the balance of three moisture fluxes – precipitation (P), evapotranspiration (ET), and transport. It uses energetic and mass balance limits on evapotranspiration to place constraints on the relationship between P and ET, which is a key parameter that controls spatial patterns of δ¹⁸O_p (Salati et al., 1979; Gat and Matsui, 1991; Lee et al., 2007; Winnick et al., 2014; Bailey et al., 2018; Kukla et al., 2019). Together with the dry and moist adiabatic lapse rate (Γ, K m⁻¹) and an assumed environmental lapse rate (γ, K m⁻¹), orographic rainout is incorporated into the model following the work of Smith (1979) and Smith and Barstad (2004). Though there are a variety of topographic parameterizations, for our simulations of westerly-derived moisture we use an idealized topography with a Gaussian-shaped mountain range combined with an orogenic plateau in the lee of the range. In our simulations of GPLLJ δ¹⁸O_p, we assume flat terrain. Even though the Plains gently rise by more than 1000 m from the Gulf Coast to the Front Range, the model of orographic precipitation applies to adiabatic ascent over topography – a process that does not occur over the Plains. Equation (1) is used to calculate the column-integrated precipitable water content (w, kg m⁻²) as a function of advection (u, m s⁻¹; the movement of an air parcel) and eddy diffusion, precipitation (P, kg m⁻²), and evapotranspiration (ET, kg m⁻² s⁻¹) (Kukla et al., 2019).

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial w}{\partial t}}=\mathrm{\nabla }\cdot (D\mathrm{\nabla }\cdot w)-u\mathrm{\nabla }\cdot w+\mathrm{ET}-P,\end{array}$$

where t is time and D is the coefficient for eddy diffusion (m² s⁻¹).

Equation (2) is used to calculate the isotopic ratio of precipitable water (r_w) over space (x) in steady state following Eq. (1). The isotope ratio of precipitation (r_P) is derived by assuming equilibrium fractionation during moisture condensation. The isotope ratio of ET (r_ET) is derived as a function of the transpired fraction of ET ($T/E+T$) and the balance of equilibrium and kinetic isotope fractionation. For further details, we refer the reader to Kukla et al. (2019).

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\mathrm{0}& =D\left({\displaystyle \frac{d^{\mathrm{2}}{r}_{\mathrm{w}}}{\mathrm{d}{x}^{\mathrm{2}}}}+\mathrm{2}{\displaystyle \frac{\mathrm{1}}{w}}{\displaystyle \frac{\mathrm{d}w}{\mathrm{d}x}}{\displaystyle \frac{\mathrm{d}{r}_{\mathrm{w}}}{\mathrm{d}x}}\right)-u{\displaystyle \frac{\mathrm{d}{r}_{\mathrm{w}}}{\mathrm{d}x}}\\ & +{\displaystyle \frac{\mathrm{ET}}{w}}(r_{\mathrm{ET}}-r_{\mathrm{w}})-{\displaystyle \frac{P}{w}}(r_P-r_{\mathrm{w}})\end{array}\end{array}$$

Here, r refers to the ¹⁸O $/$ ¹⁶O ratio and subscripts ET, P, and w refer to evapotranspiration, precipitation, and precipitable water, respectively.

Equations (1) and (2) ensure mass conservation and permit the reduction of precipitable water and of the ¹⁸O of precipitable water by precipitation. Mass conservation relationships between potential evapotranspiration (PET), evapotranspiration (ET), and precipitation (P) are encapsulated within a Budyko hydrologic balance framework that limits moisture recycling based on energy (when PET ≤ P) or water (when P ≤ PET). Equation (3) can be used to calculate a so-called “Budyko curve”, which determines a unique hydroclimate solution for each isotope gradient using a given value of ω – a non-dimensional free parameter that captures land surface characteristics (such as vegetation, bedrock lithology, ruggedness) and modulates the partitioning of precipitation into either runoff or evapotranspiration (Budyko, 1974; Fu, 1981; Zhang et al., 2004; Greve et al., 2015; Kukla et al., 2019).

$$\begin{array}{}\text{(3)}& {\displaystyle }\mathrm{ET}=P\left(\mathrm{1}+{\displaystyle \frac{\mathrm{PET}}{P}}-{\left(\mathrm{1}+{\left({\displaystyle \frac{\mathrm{PET}}{P}}\right)}^{\mathit{\omega }}\right)}^{\frac{\mathrm{1}}{\mathit{\omega }}}\right)\end{array}$$

To populate the parameters in this model, we again use North American Regional Reanalysis (NARR) data (long-term monthly mean for years 1979–2000) (Mesinger et al., 2006). We show input values for model parameters in the Appendix (Table A1). For ω, we use the global mean value of 2.6 (Greve et al., 2015). For the transpired fraction of ET ($T/E+T$) we use a value of 0.64 (Good et al., 2015).

We simplify the longitudinal differences in simulated GPLLJ trajectories (see Fig. 3) as a 1-D storm track that transports moisture from the Gulf of Mexico to the Great Plains. To reflect the general curvilinear trajectories of the GPLLJ, partly caused by the influence of the North American Cordillera (Jiang et al., 2007), we implement a bend in our simplified trajectories at 32° N. From this bend, the simulated trajectory runs along the 101° meridian. This trajectory ends at 43° N, the latitude of our northernmost site. This trajectory passes over the middle of the present-day exposure of the Ogallala Formation and, hence, captures a representation of the atmospheric processes that result in δ¹⁸O_p over the Great Plains.

Figure 3 Map showing the percentage of precipitation-producing storm tracks at four localities on the Great Plains that bound the latitudinal and longitudinal extent of our data. HYSPLIT air parcels at all four sites are initialized at 1000 m above ground. (a) 42° N, −105° E. (b) 42° N, −97° E. (c) 32° N, −105° E. (d) 32° N, −97° E. The resulting mean storm trajectories are demarcated with dashed black lines with an arrow showing the direction of transport. In panels (a)–(d), the extent of the Ogallala is shown in black. (e) An estimate of the percentage of gulf and southeasterly moisture that reaches each location as a function of longitude. Red points represent a transect of HYSPLIT simulations along the 42° N parallel; blue points represent a transect of HYSPLIT simulations along the 32° N parallel.

To test for the influence of westerly moisture in the Great Plains (as indicated by the HYSPLIT results), we initialized the vapor transport model with annual mean NARR data interpolated to a simplified 1-D trajectory representing the westerlies (Fig. 8). The trajectory latitude is chosen to lie between the northern and southern boundary of the Ogallala Formation. To account for the known effects of orographic rainout on westerly moisture (Friedman et al., 2002; Lechler and Galewsky, 2013; Mix et al., 2019), we use a simplified topography, based on the modern observed topographic profile, that follows a Gaussian-shaped mountain with a flat plateau in the lee.

5 Results

5.1 Moisture sources and precipitation trajectories

Contour plots in Fig. 3a–d show the percentage of trajectories that produce precipitation at 1000 m above ground level for each of the plotted locations over the course of a year. Western sites receive a substantial portion of their moisture from westerly or southwesterly trajectories that traverse the high topography of the North American Cordillera. In particular, the northwesternmost site (Fig. 3a) receives little moisture from the gulf. In contrast, the southeasternmost site (Fig. 3d) receives predominantly gulf moisture. These patterns vary somewhat seasonally, with gulf and southeasterly moisture more dominant in the spring and summer months and westerly moisture more dominant in the winter months. Our results from HYSPLIT's clustering algorithm show that northern locations receive less gulf and southeasterly moisture than southern locations; however, regardless of latitude, the percentage of storms sourced from the gulf increases to the east (Fig. 3e).

5.2 δ¹⁸O_c and δ¹³C results

Our new δ¹⁸O_c data range from a minimum of 17.0 ‰ to a maximum of 27.3 ‰ (Table S1 in Manser et al., 2023). When averaged by section, our new δ¹⁸O_c data range from 19.6 ‰ to 26.8 ‰ (Table 1). The compiled data have a much larger range, reflecting their larger latitudinal distribution, with the average section δ¹⁸O_c ranging from 13.6 ‰ to 27.0 ‰ (Appendix Table A2). Within individual sections, there is very little variance: the mean range and standard deviation of δ¹⁸O in individual sections are ±1.8 ‰ and ±0.68 ‰, respectively, indicating that temporal shifts in δ¹⁸O_c in the sampled sections are small (Table 1; Fig. 4b). If we exclude data outside Ogallala Formation, the mean range and standard deviation of δ¹⁸O_c in individual sections are lower at ±1.7 ‰ and ±0.6 ‰, respectively, indicating that carbonates in the Ogallala Formation have very little variance.

Figure 4 Spatial distribution of δ¹⁸O (a) and reconstructed δ¹⁸O_p 1σ (b). (a) Comparison of reconstructed δ¹⁸O_p from δ¹⁸O_c (large squares) calculated using the yearly mean of monthly long-term temperature data retrieved from NARR for each site. Modern meteoric water δ¹⁸O (small circles) is derived from modern groundwater, river water, or stream water retrieved from Waterisotopes.org (Waterisotopes Database, 2019). Data points are colored by their δ¹⁸O values. (b) The 1σ of the reconstructed δ¹⁸O_p derived from the variability in the δ¹⁸O_c values from each stratigraphic section. In both panels, the black polygon marks the extent of modern-day exposure of the Ogallala Formation.

Our new Ogallala δ¹³C data have a mean value of −6.14 ±1.4 (1σ) ‰ (Table 1) (excluding data from the Tesuque Formation). In contrast, our compiled Ogallala Formation data have a mean δ¹³C value of −6.74 ±0.8 (1σ) ‰. This discrepancy is due to relatively high δ¹³C values in the sections in New Mexico identified via lithologic or geomorphic correlation by Frye et al. (1982) to be Ogallala outcrops. When these Frye et al. (1982) data are excluded, the mean δ¹³C value in our new data is −6.69 ±0.9 (1σ) ‰, which is statistically indistinguishable from the compiled data δ¹³C mean.

5.3 Spatial distribution of δ¹⁸O_p

For each site, we calculate a mean reconstructed δ¹⁸O_p value and compare these mean values with modern water δ¹⁸O (retrieved from Waterisotopes.org; Waterisotopes Database, 2019). To reconstruct paleo-precipitation δ¹⁸O, we use the modern 2 m air temperature from NARR to estimate the formation temperature of soil carbonate along with the fractionation factors of Kim and O'Neil (1997). Means, standard deviations, and full δ¹⁸O ranges for each site with new data are presented in Table 1, and the previously published δ¹⁸O data – converted to paleo-precipitation δ¹⁸O_p – are presented in Appendix Table A2. All of the δ¹⁸O_c data are available as supporting material (Table S1) in Manser et al. (2023). The spatial distribution of the reconstructed δ¹⁸O_p shows a nearly identical spatial distribution to that of modern δ¹⁸O (Fig. 4a). From the Gulf Coast inland, paleo-δ¹⁸O follows the modern δ¹⁸O gradient, with lower values to the west and to the north. The maximum change in both paleo-δ¹⁸O and modern δ¹⁸O occurs along a southeast–northwest trend, running roughly from the Gulf Coast in Texas to northern Colorado.

We further test whether decreasing temperatures northward significantly influences our reconstructed δ¹⁸O_p values. The modern mean annual temperature gradient in the Great Plains (from 28.3 to 43° N) is −1.2 °C per 100 km. In June–July–August (JJA), this gradient is reduced to 0.4 °C per 100 km. We use these two temperature gradients to calculate the effect of changing the spatial pattern of temperature on our reconstructed δ¹⁸O_p, as well as testing a hypothetical case with latitudinally constant temperature of 25° C. Using the JJA temperature gradient, the reconstructed δ¹⁸O_p gradient is slightly shallower (red circles in Fig. 5) compared to the modern δ¹⁸O gradient (gray dots) and compared to the reconstructed δ¹⁸O_p gradient (blue circles) using modern yearly mean NARR temperature data (Fig. 5b). The likely maximum difference that a reduced temperature gradient – relative to the modern – can impart on our results is captured in the scenario of spatially uniform temperatures, which represents an extreme end-member scenario where poleward warming in a higher-CO₂ world results in a negligible latitudinal temperature gradient. Even spatially uniform temperatures only result in an increase of 2.9 ‰ in the northernmost sites relative to the assumption of a modern temperature gradient, still within the range of modern δ¹⁸O across the Plains. These results suggest that our reconstructed δ¹⁸O_p values are relatively insensitive to assumptions of carbonate formation temperatures in a warmer world with a potentially reduced meridional temperature gradient (Feng et al., 2016; van Dijk et al., 2020).

Figure 5 (a) Reconstructed δ¹⁸O_p plotted against latitude. Reconstructed δ¹⁸O_p is calculated assuming a spatially uniform temperature of 25 °C (yellow circles), annual mean temperatures (blue circles), or June–July–August (JJA) mean temperatures (red circles). Annual mean and JJA temperatures are taken from monthly long-term mean data retrieved from NARR (Mesinger et al., 2006). Gray circles are modern meteoric water δ¹⁸O, including groundwater and river or stream water, retrieved from Waterisotopes.org (Waterisotopes Database, 2019). (b) The same as (a) but the x axis sets the southernmost data (∼28° N) to zero to emphasize the effect on the spatial gradient in δ¹⁸O.

There is not only a pronounced N–S δ¹⁸O gradient in the paleo-data, but also substantial W–E variation. Figure 6 shows modern and reconstructed δ¹⁸O_p plotted against longitude using the NARR mean annual temperature to constrain the formation temperature of carbonate. Again, as with the comparison against latitude (Fig. 5), all reconstructed δ¹⁸O_p estimates fall within the same range as the modern δ¹⁸O when projected to longitude. Though again this conclusion is tempered by our assumption of carbonate formation temperatures, we suggest, as above, that this is likely to be a negligible effect, particularly given the relative insensitivity in the ¹⁸O fractionation factor to temperature (∼0.2 ‰ K⁻¹) (Kim and O'Neil, 1997).

Figure 6 Modern meteoric water δ¹⁸O (small, gray circles) and reconstructed δ¹⁸O_p (large, colored circles) versus longitude. Reconstructed δ¹⁸O_p is colored with respect to latitude. Modern meteoric water data include δ¹⁸O data from groundwater, river water, and stream water retrieved from Waterisotopes.org (Waterisotopes Database, 2019). Vertical error bars are 1σ of the mean section δ¹⁸O data.

5.4 Reactive transport modeling of δ¹⁸O_p

The vapor transport model (Kukla et al., 2019) allows us to predict δ¹⁸O_p along a given storm track, reflecting the balance between transport, precipitation, and evapotranspiration. Our model of GPLLJ δ¹⁸O_p simulates higher δ¹⁸O at the end of the trajectory compared to both modern and reconstructed δ¹⁸O_p values (Fig. 7a). Consistent with higher δ¹⁸O values, the model also predicts higher atmospheric moisture content and relative humidity than indicated by NARR, suggesting it is underpredicting rainout from air masses that reach the northern plains or that the single air mass assumption implicit in this model is not valid in this region (Fig. 7b).

Figure 7 Vapor transport model simulation of δ¹⁸O_p along an idealized Great Plains low-level jet (GPLLJ) trajectory. (a) The red line is the modeled δ¹⁸O_p along the simulated storm track (exact location of the storm track is shown as a red line in the inset). Green points are reconstructed δ¹⁸O_p within ±4° longitude of the trajectory. Gray circles are data greater than ±4° longitude from the simulated trajectory. Vertical error bars are 1σ of the mean section reconstructed δ¹⁸O_p. (b) Simulated precipitable water (solid red) and P $/$ ET (solid black) plotted along the simulated storm track. The dashed line shows the actual annual mean precipitable water along this storm track retrieved from NARR.

Figure 8 Vapor transport model simulation of δ¹⁸O_p along an idealized westerly trajectory. (a) The red line is the modeled δ¹⁸O_p along the simulated storm track (exact location of the storm track is shown as a red line in the inset). Green points are reconstructed δ¹⁸O_p data within ±4° latitude from the storm track. Large, unfilled gray circles are data greater than ±4° latitude from the simulated trajectory. Vertical error bars are 1σ of the mean section reconstructed δ¹⁸O_p. Small gray circles are modern meteoric water δ¹⁸O, including stream water, river water, and groundwater, retrieved from Waterisotopes.org (Waterisotopes Database, 2019). The gray line is a kernel-smoothed average of the meteoric water δ¹⁸O against distance along the storm track. (b) Simulated precipitable water (solid red) and P $/$ ET (solid black) plotted along the simulated storm track. The dashed line shows the actual annual mean precipitable water along this storm track retrieved from NARR.

Our model of westerly moisture δ¹⁸O_p suffers from a similar mismatch between simulated δ¹⁸O_p and both modern δ¹⁸O and reconstructed δ¹⁸O_p (Fig. 8a). Though our westerly trajectory uses an idealized topography to simulate the topography of the North American Cordillera, we are able to approximately reproduce the decrease in δ¹⁸O_p that occurs close to the coast of North America. Variations in simulated δ¹⁸O_p after the initial orographic rainout are primarily driven by temperature that varies with elevation in the NARR data. This leads to minor discrepancies where elevation decreases δ¹⁸O_p in modern waters but increases δ¹⁸O_p in the model due to colder temperatures raising δ¹⁸O_p. Nevertheless, simulated δ¹⁸O_p approximates meteoric water δ¹⁸O along the westerly trajectory until this trajectory reaches the Great Plains at a distance of around 1550 km. At distances larger than this, the model substantially underestimates meteoric water δ¹⁸O and reconstructed δ¹⁸O_p on the Great Plains. Again examining climatological measures related to δ¹⁸O_p, there is a distinct discrepancy in atmospheric moisture content between the modeled output and the NARR data over the Great Plains at a distance of around 1550 km from the starting point of the simulated westerly trajectory (Fig. 8b).

6 Discussion

Our new spatially resolved data provide insight into moisture transport to the Great Plains during the late Neogene. In brief, the overall constancy of δ¹⁸O between the late Miocene and the present suggests that the features of atmospheric circulation most responsible for moisture delivery to the Great Plains today – notably the Great Plains low-level jet and the wintertime westerlies – have likely been the dominant features since at least the late Miocene. Further, that the reconstructed δ¹⁸O_p gradient is indistinguishable from today's δ¹⁸O gradient indicates that the balance of GPLLJ and westerly moisture has hardly changed since the late Miocene. Either these circulation features have undergone no substantial change in net rainout, or, if one has, its effect was masked by countervailing changes in the other. Overall, these results bolster earlier findings by Fox and Koch (2004), who found a consistent south-to-north gradient in δ¹⁸O_c in the Great Plains that they attributed to a similar latitudinal temperature gradient in the Miocene.

Below, we discuss our findings in more detail, place these data in the context of previous work in the region, and discuss how our new data permit a re-interpretation of the controls on long-term climate in the Great Plains. We also discuss several important caveats in our data that point towards the need for future work to study more nuanced changes in climate than can be resolved with this dataset.

6.1 Constancy of δ¹⁸O_p in relation to climate change

The Miocene was warmer than the present day, particularly during the Miocene Climate Optimum (MCO) (Westerhold et al., 2020; Steinthorsdottir et al., 2021), and was characterized by a long-term cooling trend after the MCO (Herbert et al., 2016). Even after this cooling interval, the lack of extensive Northern Hemisphere ice sheets and likely a smaller Antarctic ice sheet indicate that the late Miocene was warmer than today, with a reduced latitudinal temperature gradient (LaRiviere et al., 2012; Feng et al., 2016). The warmer climate and shallower temperature gradient could impact δ¹⁸O_p, yet the decrease in reconstructed δ¹⁸O_p with latitude appears identical to today (Fig. 4). This result is perhaps not surprising, since temperature appears to have only a secondary effect on the latitudinal gradient of δ¹⁸O_p (Fig. 5a). Instead, such constancy in δ¹⁸O_p through time indicates that the mixing of westerly and southerly moisture – the primary control on the latitudinal δ¹⁸O gradient – was similar to today.

The results of our vapor transport modeling support the contention that mixing between dry, low-δ¹⁸O westerly air masses and moist, high-δ¹⁸O southerly air masses is required to explain the long-standing presence of this steep latitudinal gradient in δ¹⁸O_p. The vapor transport model adequately predicts δ¹⁸O_p in westerly- or GPLLJ-dominated regions, but it performs poorly where these air masses meet (Figs. 7 and 8). Air mass mixing is neglected in the 1-D vapor transport model, so the model's poor performance in these regions points to mixing as an important control on the latitudinal δ¹⁸O_p gradient. Notably, this region of air mass mixing is closely tied to the east–west aridity gradient because it tracks the trade-off between drier westerly air and the wetter GPLLJ. Thus, the surprisingly static spatial δ¹⁸O_p pattern across the Great Plains since the Miocene could be interpreted to reflect no change in the relative mixing of dry and moist air masses over time, maintaining the spatial pattern of the modern aridity gradient. The presence of the North American Cordillera – likely high since the Eocene (Chamberlain et al., 2012; Mix et al., 2013) – suggests that orographic rainout on the western margin of North America is also a long-standing feature, resulting in low-δ¹⁸O moisture that is advected to the Great Plains from the west. What is more surprising is that the strength of the GPLLJ appears to be similar to today even in the warmer Miocene, given that this jet and its moisture transport are likely to be sensitive to global climate change (Cook et al., 2008; Zhou et al., 2021a).

Overall, our vapor transport modeling indicates that a purely southerly source of moisture would result in δ¹⁸O_p values that are too high in the northern plains (Fig. 7), whereas a purely westerly source brings δ¹⁸O_p values that are too low (Fig. 8). To be clear, a static δ¹⁸O_p gradient does not require static precipitation and evaporation rates over time. In contrast, precipitable water, P, and E likely all increased in the warmer Miocene (Held and Soden, 2006; Siler et al., 2018), but in such a way that net rainout (and thus the δ¹⁸O_p fingerprint of hydrological change) did not vary substantially.

6.2 Additional factors that influence δ¹⁸O_c

There are a number of additional factors, both climatic and non-climatic, that may affect the δ¹⁸O_c data, potentially decoupling δ¹⁸O_c values from δ¹⁸O_p and conflicting with our interpretation above. Such factors include changes in soil temperatures and evaporation and also imprecision in the chronologies of the sections we sampled. Below, we discuss these factors and why our interpretations are robust to assumptions regarding these factors.

6.2.1 Effect of temperature and evaporation on δ¹⁸O_c

Given that temperature affects the fractionation of ¹⁸O between water and calcite (Kim and O'Neil, 1997), spatial changes in temperature may alter our reconstructed δ¹⁸O_p gradient. However, we found that even extreme scenarios for changes in the spatial pattern of temperature over the Great Plains (ranging from the modern annual average latitudinal temperature gradient to a null latitudinal temperature gradient) had only a small effect on the overall latitudinal gradient of reconstructed δ¹⁸O_p (Fig. 5). We therefore conclude that, while temperature change since the Miocene has likely impacted the absolute value of δ¹⁸O_c, it does not substantially alter our conclusion that the reconstructed δ¹⁸O_p spatial pattern is similar to today.

While the spatial pattern of temperature change in the Plains does not appear to have had a major influence on our reconstructed δ¹⁸O_p gradient, for any reasonable temperature scenario (annual, JJA, or spatially uniform 25 °C), reconstructed δ¹⁸O_p absolute values are typically higher than modern δ¹⁸O_p by 2 ‰–3 ‰. Fox and Koch (2004) also found that δ¹⁸O_c was typically higher than predicted, particularly in the southern Great Plains, and attributed this observation to both higher temperatures and δ¹⁸O_p. We suggest that, while the gradient of δ¹⁸O_p does not appear to have substantially changed since the Miocene, the starting value of marine vapor δ¹⁸O may have been slightly higher in the warmer Miocene due to elevated precipitable water, elevating δ¹⁸O_p across the Great Plains (i.e., van Dijk et al., 2020).

Recent work has indicated that evaporative conditions (i.e., low AI values) may elevate δ¹⁸O_c, thereby leading to an overestimate of δ¹⁸O_p (Kelson et al., 2023). We suggest that our overall conclusion of a similar spatial pattern of δ¹⁸O_p in the late Miocene is likely not overly influenced by changes in the spatial pattern of evaporation. Today, the latitudinal gradient in AI is negligible (Seager et al., 2018b); instead, the strongest gradient in AI is west to east (Fig. 1). Consequently, the stark latitudinal gradient in reconstructed δ¹⁸O_p is unlikely to be driven by differential evaporative conditions. While higher marine vapor δ¹⁸O may explain the elevated reconstructed absolute δ¹⁸O_p values (2 ‰–3 ‰ above modern meteoric water δ¹⁸O), another explanation may be pervasive evaporative enrichment of pedogenic carbonates across the Great Plains, though this evaporative enrichment does not appear to vary latitudinally.

6.2.2 Changes in precipitation seasonality

Soil carbonates are thought to form mostly seasonally, particularly in the warm season and perhaps as soils dry and plants senesce, causing soil CO₂ to decline and carbonates to form (Breecker et al., 2009; Huth et al., 2019). However, the timing of soil carbonate formation may shift as the timing of precipitation and plant productivity varies due to climate change, and such shifts have been invoked to explain a wide range of soil carbonate records in the western US and elsewhere (Caves, 2017; Kukla et al., 2022b; Rugenstein et al., 2022). GCMs indicate that precipitation seasonality over the Great Plains will change substantially in a warmer world, with precipitation shifting towards the spring and away from the late summer (Cook et al., 2008; Bukovsky et al., 2017; Zhou et al., 2021a). Though it is difficult to assess the timing of even modern soil carbonate formation in the Great Plains, our reconstructed Miocene δ¹⁸O_p data agree most closely with modern water δ¹⁸O when using mean annual temperatures. Though we have no a priori reason to suspect that Miocene carbonates may record mean annual temperatures, this close agreement suggests that Miocene carbonates may have formed in the shoulder seasons (i.e., spring and fall), thereby recording intermediate temperatures and, likely, a mixture of summertime GPLLJ moisture and wintertime westerly moisture.

That Miocene soil carbonates formed in the shoulder seasons would help to reconcile our findings of a largely invariant spatial pattern of δ¹⁸O_p since the Miocene with previously published records that show a large (2 ‰–4 ‰) increase in δ¹⁸O_c during the Pliocene–Quaternary in the Great Plains (Fox et al., 2012; Mix et al., 2013; Chamberlain et al., 2014). We suggest that this discrepancy (i.e., a static late Miocene δ¹⁸O_p spatial pattern but increasing δ¹⁸O_c in the Pliocene–Quaternary) likely arises due to changes in the seasonality of soil carbonate formation between the late Miocene and Quaternary. For example, as the world cooled into the Quaternary, the southward shift of the westerly jet would suppress the springtime GPLLJ, while enhancing GPLLJ precipitation over the Great Plains during the summer – the opposite response than observed in GCMs as CO₂ rises (Cook et al., 2008; Bukovsky et al., 2017; Zhou et al., 2021a). Such a seasonality shift might shift the timing of soil carbonate formation towards the summer, elevating δ¹⁸O_c as it records higher summertime δ¹⁸O_p (Liu et al., 2010) and thereby elevating δ¹⁸O_c relative to the late Miocene.

6.2.3 Imprecision in chronologies and timing of carbonate formation

Our interpretation hinges upon correlating sections across the vast expanse of the Great Plains; however, ages for many Ogallala sections are only poorly constrained, typically via biostratigraphy (Kitts, 1965; Skinner et al., 1977; Thomasson, 1979; Winkler, 1987, 1990; Tedford et al., 2004) or, less frequently, via radiometric dates (Swisher, 1992; Cepeda and Perkins, 2006; Henry, 2017; Smith et al., 2018). Where biostratigraphy has been used, age resolution is typically limited to the scale of the North American land mammal ages. In many places in eastern New Mexico, the Ogallala Formation is recognized by its clay lithology and its topographic position relative to the laterally extensive caprock to the east (Frye et al., 1978, 1982). Further, the character of Ogallala deposition – large alluvial megafans that gradually filled in pre-existing valleys and covered intervening interfluves combined with its relative thinness – suggests that deposition was often sporadic with frequent and cryptic unconformities associated with many of the paleosols distributed throughout the formation. Consequently, correlations across the entirety of the Ogallala Formation are likely to be imprecise and our treatment of the data likely mixes data from the late middle Miocene and the latest Miocene.

Nevertheless, this imprecision in correlation likely does not affect our conclusions. We note that nearly all of our sections have low δ¹⁸O_c variability (Fig. 4b; Table 1), suggesting that, throughout deposition of the Ogallala, there was very little change in δ¹⁸O_p. Ludvigson et al. (2016), analyzing two cores drilled through the Ogallala Formation in western Kansas, found similarly very low variability in δ¹⁸O_c across the 50–60 m of sampled core material. This low variability, which contrasts with sections located to the west within the North American Cordillera, indicates that even large discrepancies in correlated sections likely do not affect our overall conclusion that the south-to-north δ¹⁸O_p gradient has remained similar to today since the late Miocene.

Further, previously published Neogene δ¹⁸O_c records show almost no change in δ¹⁸O_c during the Miocene (Fox and Koch, 2004; Mix et al., 2013; Chamberlain et al., 2014), again suggesting that imprecision in our correlations is not likely to impact our conclusions. In contrast, modeling studies that have examined changes in δ¹⁸O_p in the Great Plains in the Miocene have suggested that the south-to-north δ¹⁸O_p gradient should increase by several per mille as a consequence of higher CO₂ and a shallower Equator-to-pole temperature gradient (Feng et al., 2016; Lee, 2019). These studies used boundary conditions thought to approximate the middle Miocene, with atmospheric CO₂ of 560 ppm. Disagreement between our results and these model results may arise for a variety of reasons. One such reason may revolve around uncertainty in dating of the Ogallala Formation. While these modeling studies used middle Miocene boundary conditions and 560 ppm CO₂, our data likely come from the late Miocene, following substantial cooling after the peak of Neogene warmth during the middle Miocene (Herbert et al., 2016). Consequently, the environment in which much of the Ogallala Formation was deposited may differ from that simulated by Feng et al. (2016) and Lee (2019). Alternatively, given the difficulty in simulating the land–atmosphere coupling over the Great Plains in modern-day simulations (Bukovsky et al., 2017; Laguë et al., 2019; Zhou et al., 2021a) and the importance of the land surface in modulating δ¹⁸O_p, minor misrepresentations in the land surface parameterizations within these models may yield outsize impacts on the simulated δ¹⁸O_p gradient. While a thorough review comparing our results to simulated Miocene δ¹⁸O_p over the Great Plains is outside the scope of this paper, our new spatially resolved dataset provides an opportunity to test model predictions of late Miocene climate simulation skill over central North America.

We additionally assume that the authigenic carbonates that we sampled formed in the late Miocene. Several studies (Joeckel et al., 2014; Smith and Platt, 2023) have noted that some caliche units in the northern and central Great Plains may reflect post-depositional case hardening rather than syn-depositional carbonate formation based upon lateral discontinuities in these units. Given the low variance in our Ogallala samples (typically much less than 1 ‰), we suggest that all of the carbonates in our study formed from the same meteoric water. Further, the δ¹³C values of our carbonate samples are relatively low ($\sim -\mathrm{7}$ ‰). Given the expansion of C₄ plants during the Plio-Pleistocene, the δ¹³C values of authigenic carbonates formed in the Plio-Pleistocene are substantially higher ($\sim -\mathrm{2}$ ‰) (Fox and Koch, 2003). Thus, the low δ¹³C values of our carbonate samples indicate that they formed during the late Miocene under a dominantly C₃ grassland environment. The only exception to these low δ¹³C values in our sections is from carbonates sampled in eastern New Mexico from outcrops identified as Ogallala by Frye et al. (1982). In many of these sections, there are authigenic carbonate samples with δ¹³C values $>-\mathrm{5}$ ‰ and frequently the caprock caliche samples have δ¹³C values approaching 0 ‰ (Table S1 in Manser et al., 2023). This suggests that these outcrops either may not be correlative with the Ogallala (Henry, 2017) or that authigenic carbonate formation in these sections occurred substantially later than Ogallala deposition. However, we have no independent age constraints with which to better constrain the ages of these units, and we therefore include them in our study. We do note that the reconstructed δ¹⁸O_p from these sections is indistinguishable from modern δ¹⁸O_p in eastern New Mexico, further suggesting that the long-term pattern of δ¹⁸O_p in the Great Plains has remained invariant. Thus, the combination of the low δ¹³C values and the low variance in our δ¹⁸O_c data suggests that our assumption of late Miocene carbonate formation is robust.

7 Implications

That atmospheric circulation over the Great Plains has remained relatively constant since the late Miocene provides important context for understanding how Great Plains hydroclimate and environments may respond to higher atmospheric CO₂. Further, our data provide insight into the climatic and tectonic controls that may have driven deposition of the Ogallala Formation. In both cases, our data places critical constraints on our understanding of both past environments and the future evolution of climate in the Great Plains.

7.1 Response of Great Plains hydroclimate to CO₂

Of concern as atmospheric CO₂ rises is whether and how the climatological 100th meridian, which demarcates the semi-arid west from the humid east, will shift. The current position of this climatological 100th meridian is partly set today by the boundary between dry westerly air masses that have lost much of their moisture from passage over the North American Cordillera and moist southerly masses transported by the GPLLJ. If dynamical relative shifts in the strength of these two predominant circulation systems were to occur as CO₂ rises, one might expect the boundary between low δ¹⁸O_p and high δ¹⁸O_p to shift. For example, model simulations indicate that the poleward shift of the westerly jet and North Atlantic Subtropical High with warming should enhance GPLLJ moisture transport in the spring but suppress it in the summer (Zhou et al., 2021a, b). If the overall tendency were a weaker GPLLJ with warming, we might expect less northward moisture transport and thereby a shift of this δ¹⁸O_p boundary southward and eastward. Instead, we find that the spatial pattern of δ¹⁸O_p has remained unchanged since the late Miocene, suggesting that the relative strength of these two circulation systems was not substantially different in a warmer, higher-CO₂ world. Such findings bolster GCM results that, overall, mean wet-season precipitation does not change substantially and that therefore the boundary between the GPLLJ and the westerlies is not overly sensitive to warming.

Though dynamical shifts in circulation may affect hydroclimate over the Great Plains, additional thermodynamic mechanisms have been invoked to explain shifts in the climatological 100th meridian with warming. Because PET will rise faster than P across the Plains, aridity will increase (Seager et al., 2018a, b). However, this will be complicated by how actual ET – and therefore the partitioning of moisture returned to the atmosphere versus to runoff – responds to warming. How ET changes with warming depends on a variety of factors including plant water use efficiency (Swann et al., 2016; Lemordant et al., 2018) and the temporal distribution of precipitation (Scheff et al., 2022). Our results imply that either (1) the effect of increasing PET relative to P had a small effect on net rainout – and hence δ¹⁸O_p – in the Miocene, or (2) decreases in ET efficiency with warming, perhaps driven by more efficient plant water use (Lemordant et al., 2018) and which would reduce ¹⁸O return to the atmosphere, counteracted the expected increase in δ¹⁸O_p due to increasing PET relative to P.

Though the late Miocene was likely warmer than today with higher than pre-industrial CO₂ (Herbert et al., 2016; Mejía et al., 2017; Sosdian et al., 2018; CenCO2PIP Consortium, 2023), the precise global climate that our new data reflect is uncertain. This uncertainty arises not only because of the chronological uncertainty associated with the Ogallala Formation, but also due to uncertainty in both global temperature and CO₂ reconstructions. The late Miocene may have seen global temperatures elevated by 5° higher relative to pre-industrial (Westerhold et al., 2020; Ring et al., 2022) and atmospheric CO₂ between 350 and 500 ppm (Tanner et al., 2020; Brown et al., 2022; CenCO2PIP Consortium, 2023). However, these estimates – particularly for CO₂ – remain imprecise, and it is therefore difficult to determine just how insensitive the Great Plains hydroclimate is to warming and higher CO₂. Further, while the land surface plays a critical role in modifying Great Plains hydroclimate, the late Miocene Great Plains were already dominated by grassland ecosystems (Stromberg, 2005). Though C₄ plants have a greater water use efficiency than C₃ plants (Osborne and Sack, 2012), the invariant δ¹⁸O_p spatial pattern suggests that overall ecosystem water use (i.e., transpiration to the atmosphere) was likely similar to today and overall hydroclimate in the late Miocene resembled the hydroclimate today in the Great Plains.

7.2 Implications for understanding the origin of the Ogallala Formation

Why much of the Great Plains shifted from a dominantly erosive landscape to a depositional landscape for several million years, before returning to the dominantly erosive landscape of today, has remained a major outstanding question in North American geology (Heller et al., 2003). As the sediments of the Ogallala Formation are sourced in the Rocky Mountains, this question has been intimately linked to what process drove this major late Cenozoic erosional event, producing the alluvial fans and megafans that blanketed the Great Plains in the late Miocene (Trimble, 1980; Molnar and England, 1990). Paleoaltimetry data have suggested that much of the North American Cordillera has been high since the Eocene (Forest et al., 1995; Mix et al., 2011). If mean elevation has not changed for much of the Cenozoic, then changes in climate – perhaps driven by increasingly variable, orbitally controlled glacial cycles – may have driven the increase in erosion (Molnar and England, 1990; Gregory and Chase, 1994; Zhang et al., 2001). In contrast, multiple studies that have analyzed paleo-slopes of the Ogallala Formation fluvial sediments conclude that the Ogallala Formation must have been tilted post-deposition due to a long-wavelength uplift (McMillan et al., 2002; Heller et al., 2003), perhaps associated with dynamic topography effects driven by continued subduction of the Farallon Plate (Willett et al., 2018).

Our data suggest that hydroclimate was not substantially different in the late Miocene and was therefore unlikely to be the proximal cause of widespread erosion of the Rocky Mountains and deposition of the Ogallala Formation. Instead, we suggest that, because overall hydroclimate – and, in particular, runoff – was similar to today, deposition of the Ogallala Formation must have been driven by long-wavelength tilting of the North American plate that uplifted portions of the Rocky Mountains. Estimates of the uplift necessary to explain post-depositional tilting are less than 1 km (McMillan et al., 2002; Leonard, 2002). Such uplift was likely not substantial enough to either (1) be detectable in paleo-altimetry datasets or (2) substantially modify the strength of the GPLLJ (Jiang et al., 2007). Uplift might be expected to result in lower δ¹⁸O_p in late Cenozoic basins in the Rocky Mountains due to orographic forcing of precipitation. However, in this region, orographic forcing is limited to GPLLJ storms that reach the Front Range, and continued mixing with low-δ¹⁸O_p westerly moisture would make it very difficult to discern such a signal. We thus conclude that hydroclimate was not substantially different in the late Miocene in the Great Plains and that long-wavelength tilting must be a necessary component of the formation of the Ogallala Formation.

8 Conclusion

Our new spatially resolved late Miocene δ¹⁸O_c data from across the Great Plains reveal that the spatial pattern of reconstructed δ¹⁸O_p is indistinguishable from the spatial pattern of modern meteoric water δ¹⁸O. Despite changes in global climate, we suggest that this static δ¹⁸O spatial pattern reflects an atmospheric circulation system over the Great Plains that in the late Miocene was largely identical to today: wintertime westerly air masses, dried by transit over the North American Cordillera, delivered low-δ¹⁸O moisture to the northern Great Plains, while southerly moisture, transported by the Great Plains low-level jet, brought high-δ¹⁸O moisture to the southern Great Plains. Thus, on the timescales preserved within the Ogallala Formation, the mixing zone between these two circulation systems was similar to today. Given that these two atmospheric circulation systems are responsible for a nearly 15 ‰ δ¹⁸O gradient north to south along the Great Plains, these results appear to be robust to assumptions regarding temperature, evaporation, and precipitation seasonality changes since the late Miocene and also insensitive to uncertainties in our correlations across the expanse of the Ogallala Formation.

These results suggest that, on geological timescales, large-scale hydroclimate in the Great Plains is relatively insensitive to changes in global temperature and atmospheric CO₂. Model projections tend to indicate that mean annual precipitation will not increase substantially as CO₂ rises (Bukovsky et al., 2017; Zhou et al., 2021a), and our results support this contention, in that a static δ¹⁸O spatial pattern suggests that net rainout in the late Miocene across the Great Plains was likely not substantially different than today. In contrast, our results do not support projections that aridity will increase over the Plains (Seager et al., 2018a). While increases in PET are robust in models (Scheff et al., 2017), our results suggest that countervailing decreases in actual ET efficiency likely offset any change in PET, yielding an approximately constant spatial pattern of P $/$ ET. Further, our results indicate that large-scale changes in climate (Zhang et al., 2001) are likely not the cause of the widespread erosional event responsible for deposition of Ogallala sediments. Instead, the fact that Great Plains hydroclimate in the late Miocene was similar to today supports the notion that long-wavelength uplift (Moucha et al., 2009; Karlstrom et al., 2011; Willett et al., 2018) shifted the Great Plains from a dominantly erosive landscape to an aggradational one.

The Great Plains lie at a unique climatic and geologic intersection, and this confluence has inspired more than a century of work to understand the relationship between climate, geology, and ecosystems over geologic time in the Plains (Powell, 1879, 1890; Seni, 1980; Fox and Koch, 2003; Jiang et al., 2007; Seager et al., 2018b; Willett et al., 2018). However, the sharp geologic and climatic gradients that characterize this landscape complicate efforts to understand past changes. Nevertheless, our results suggest that spatially resolved datasets provide a powerful means to constrain the position and shape of these gradients in the past. Additional work that utilizes newly developed proxy systems and/or pairs well-established proxies with our continental-scale data will provide an even sharper picture of the climatic and geologic forces that have shaped this remarkable landscape.

Appendix A

Table A1 Input parameters for the reactive transport model. All inputs are retrieved from NARR data except ω, which is set to the global average of 2.6 (Greve et al., 2015), and the transpired fraction of ET, which is set to 0.64 (Good et al., 2015).

Table A2 Site averages for previously published data, accessed via the PATCH Lab (Kukla et al., 2022a). Values of δ¹³C are reported in ‰ relative to VPDB. Values of δ¹⁸O are reported in ‰ relative to VSMOW. No standard deviation or range is reported for sites that only include one sample. Note that n/a means not applicable.