Impact of Mexican Border rift structural inheritance on Laramide rivers of the Tornillo basin, west Texas (USA): Insights from detrital zircon provenance

Large rivers (>10³ km long, >10⁵ km² catchment area) are the main conduits for transfer of sediment from erosional hinterland source areas to coastal depositional sinks (e.g., Schumm, 1977; Allen, 2008; Romans et al., 2016). Therefore, controls on continental-scale drainage organization and trunk river routing have fundamental impacts on the position and evolution of sedimentary basins and continental margin sedimentation. Plate-tectonic setting has long been posited as the first-order control on the development and organization of river systems (Potter, 1978; Miall, 2006; Tandon and Sinha, 2007). In actively deforming orogenic hinterlands, rivers are sensitive to tectonically driven surface uplift and topography, and drainage patterns and trunk river routing can be a direct expression of orogenic structure (Brookfield, 1998; Hallet and Molnar, 2001; Clark et al., 2004; Castelltort et al., 2012; Seagren et al., 2022). For example, in active fold-and-thrust belts, drainage basins are commonly dominated by arrays of smaller transverse tributaries and distributive fluvial systems that flow parallel to the regional slope of the orogenic wedge and feed axially routed trunk rivers within the associated flexural foreland basin (e.g., Burbank, 1992; Hovius, 1996; Horton and DeCelles, 2001; Weissmann et al., 2010). Alternatively, structural inheritance of preexisting mechanical and rheological weaknesses in convergent orogens can cause out-of-sequence deformation and the development of complex orogenic topography in a “broken foreland,” leading to more tortuous drainage patterns that respond to local uplift, topography, and basin fragmentation (Hain et al., 2011; Strecker et al., 2011; Seagren et al., 2022). In tectonically quiescent regions, tectonic inheritance, such as thickened or thinned crust and paleorift fault architecture, can also impart a discernible control on drainage organization and trunk river position (Potter, 1978; Cox, 1989, 1994). Hence, structural and tectonic inheritance likely play a key role in long-term continental drainage organization from the orogenic hinterland to the passive continental margin and, ultimately, can impart a primary control on the long-term position of terminal sedimentary sinks.

During the Late Cretaceous–Paleogene, numerous large rivers drained the southern half of North America toward the Gulf of Mexico (e.g., Lawton et al., 2009; Galloway et al., 2011; Sharman et al., 2017; Blum et al., 2017). However, the location and persistence of these large rivers and tectonic controls on sediment routing remain unclear. The Late Cretaceous onset of Laramide basement–involved deformation segmented the Sevier foreland into a complex landscape of basement-cored uplifts and intermontane basins (Fig. 1; e.g., Dickinson and Snyder, 1978; Dickinson et al., 1988; Lawton, 2019). Sediment routing through the evolving Laramide landscape likely comprised intra- and extra-orogen rivers and intermontane lakes (Dickinson et al., 1988; Lawton, 2019). Multiple recent syntheses have attempted to reconstruct detailed sediment dispersal pathways from the Laramide orogen to the Gulf of Mexico with objectives ranging from predicting deep-water depositional trends to studying tectonic and climatic controls on continental-scale routing systems within a linked source-to-sink framework (e.g., Galloway et al., 2011; Sharman et al., 2017; Blum et al., 2017; Hessler et al., 2017; Zhang et al., 2018). Despite these efforts, key uncertainties remain surrounding Laramide sediment dispersal to the Gulf of Mexico, including the capture and integration of Laramide intra-orogen fluvial systems and internally drained Laramide hinterland basins by extra-orogen drainages (e.g., Sharman et al., 2017, 2022; Blum et al., 2017; Lawton, 2019). Resolving these questions will be fundamental to understanding the tectonogeomorphic evolution of broken forelands and tectonic controls on the development of large drainages.

While many sedimentary provenance and sediment-routing studies have focused on the Rocky Mountain part of the Laramide orogen, the nature and pathways of Late Cretaceous to Eocene large rivers in the southwest United States and along the U.S.-Mexico border region remain poorly understood. The tectonic evolution of this corridor is complicated by interactions among the thin-skinned Mexican orogen, Laramide basement uplifts, and various inverted early Mesozoic rift basins, which created a complex broken-foreland basin system (Fig. 1; e.g., Fitz-Díaz et al., 2018). Late Cretaceous–Eocene rivers draining the Laramide-Mexican broken foreland had to navigate a transitional tectonogeomorphic landscape from the Laramide hinterland through the inverted Mexican Border rift–Mexican orogenic belt toward the Gulf of Mexico. These paleorivers have been the subject of limited detrital zircon studies (e.g., Lawton et al., 2009; Mackey et al., 2012; Clinkscales and Lawton, 2015; Bataille et al., 2019; Craddock et al., 2021), and a more detailed provenance analysis is warranted to create a robust sediment dispersal reconstruction that delineates U.S.-Mexico border region sediment-routing pathways and considers the role of inherited tectonogeomorphic topography. Here, we present new depth-profile detrital zircon U-Pb data and a provenance analysis for Upper Cretaceous–Eocene fluvial strata of the Tornillo basin in west Texas and integrate them with other provenance studies to reconstruct regional sediment dispersal pathways from the Laramide hinterland. The Tornillo basin was situated in a complex tectonogeomorphic setting created by the interaction of the Laramide and Mexican orogens with inherited Mexican Border rift architecture and therefore offers a unique opportunity to reconstruct how the inversion of inherited tectonic and structural features impacted the drainage evolution of Laramide rivers through the U.S.-Mexico border region (Fig. 1). Furthermore, these strata represent an underappreciated link in recent attempts to create comprehensive sediment-routing models from the U.S. and Mexico hinterlands to the Gulf of Mexico (e.g., Galloway et al., 2011; Sharman et al., 2017; Blum et al., 2017), and delineation of their provenance will help to connect fluvial systems of the northwest Gulf of Mexico margin to their hinterland sources and may have cascading implications for the hinterland sources of other major Paleogene fluvial sedimentation axes along the Gulf of Mexico.

The Campanian–Lutetian Tornillo basin occupies a position along the former southwestern margin of North America in west Texas, sitting atop the late Paleozoic Marathon fold-and-thrust belt, and it is bounded by the Chihuahua fold belt to the southwest and the Marathon uplift to the northeast (Figs. 1 and 2; Item S1 in the Supplemental Material¹; Lehman, 1991; Hennings, 1994). The Tornillo basin has been variably described as a flexural foreland or intermontane basin (Lehman, 1991; Lehman et al., 2018) that formed in response to the Laramide (e.g., Muehlberger, 1980; Cobb and Poth, 1980; DeCamp, 1985) and Mexican orogens (sensu Fitz-Díaz et al., 2018). Whereas sediment accumulation rates (~15–45 m/m.y.) would place the Tornillo basin in a distal foreland basin setting (Lehman, 1986, 1991; Schiebout et al., 1987; Bataille et al., 2016, 2019), the tectonic setting of the Tornillo basin is more in line with a broken Mexican foreland (sensu Ingersoll, 2011) fragmented by various uplifts and inverted Mexican Border rift basins in the transition region of the Laramide-Mexican orogen.

During the early Mesozoic, the continental Mexican Border rift system, which extended from northeastern Mexico to southeastern California along the former Paleozoic continental margin of North America, formed a series of northwest-southeast–trending extensional rift basins and horst blocks (Figs. 1 and 2; Lawton and McMillan, 1999; Dickinson and Lawton, 2001a, 2001b; Stern and Dickinson, 2010). In the western part, the Mogollon-Burro rift shoulder exposed Yavapai-Mazatzal (ca. 1800–1620 Ma) basement in central Arizona and southern New Mexico, forming the Mogollon highlands (Bilodeau, 1986; Reynolds et al., 1989; Dickinson et al., 1989). The corresponding McCoy and Bisbee rift basins to the south (Bilodeau, 1979, 1982; Dickinson et al., 1986; Lawton and McMillan, 1999) accumulated syntectonic Upper Jurassic–Lower Cretaceous clastic sediments derived from adjacent exposed eroding Precambrian basement, Paleozoic sedimentary successions, and Jurassic volcanic assemblages (e.g., Dickinson and Lawton, 2001b; Barth et al., 2004; Spencer et al., 2011; Lawton et al., 2020). The Mexican Border rift axis trended southeast into west Texas and northeastern Mexico as the Chihuahua and Sabinas troughs, which received Upper Jurassic–Lower Cretaceous clastic detritus from flanking uplifts and platforms (Dickinson and Lawton, 2001b; Eguiluz de Antuñano, 2001). By Aptian–Albian time, widespread subsidence and high sea level resulted in extensive carbonate deposition in the greater Gulf of Mexico region (McFarlan and Menes, 1991).

The Late Cretaceous saw a major change in tectonic regime in the U.S.-Mexico border region as the Laramide flat slab impinged on the Cordilleran margin ca. 95–85 Ma (Livaccari et al., 1981; Liu et al., 2010). This event resulted in regional uplift and retro-arc shortening extending from southeastern California into southern Arizona, New Mexico, and northern Sonora before migrating into the interior of the Colorado Plateau and Rocky Mountains (Lawton, 2019, and references therein). The McCoy basin was overthrust by Laramide reverse faults likely reactivating rift basin–bounding faults (Hoisch et al., 1988; Tosdal, 1990; Dickinson and Lawton, 2001b). Similarly, the Bisbee basin was disrupted by Late Cretaceous–Eocene Laramide reverse faults that flanked basement-cored uplifts and resulted in syntectonic Late Cretaceous–Eocene strata deposited disconformably or angularly over the Bisbee Group (e.g., Davis, 1979; Seager and Mack, 1986; Dickinson and Lawton, 2001b; Clinkscales and Lawton, 2015, 2017). The orientation and spacing of Laramide structures in southern Arizona, New Mexico, and northern Sonora appear to have been in part controlled by inherited normal faults associated with the Mexican Border rift, and Laramide inversion of closely spaced normal faults created a border region landscape characterized by significant topographic relief dissected by a mixture of coarse clastic, lacustrine, and basin-parallel fluvial depositional systems (Basabilvazo, 2000; Bayona and Lawton, 2003; Clinkscales and Lawton, 2015, 2017).

Contemporaneously, the Late Cretaceous Mexican orogen, which extends for ~2000 km along the length of Mexico (Fig. 1; Fitz-Díaz et al., 2018, and references therein), transitioned from a thin-skinned, salt-detached fold-and-thrust belt and foreland basin system to a thick-skinned orogen that reactivated early Mesozoic extensional structures of the Mexican Border rift in north-northeast Mexico (McKee et al., 1990; Lawton, 2000; Eguiluz de Antuñano, 2001; Haenggi, 2001, 2002; Chávez-Cabello et al., 2005; Mauel et al., 2011; Fitz-Díaz et al., 2018). However, reactivation of inherited Mexican Border rift structures by Mexican orogenic deformation varied spatially due to the nature of shortening and inversion experienced by different segments of the broken foreland. For instance, the salt-dominated Parras and La Popa Basins preserve largely intact sedimentary successions, and the Coahuila block experienced only limited deformation and remained relatively high standing through the Eocene (Weidie and Murray, 1967; McBride et al., 1974; Lawton et al., 2009). In contrast, the Chihuahua trough was uplifted and eroded (Haenggi, 2002) between the Maastrichtian and the Eocene, and its synrift basin fill was severely deformed, inverted, and thrust eastward as the Chihuahua fold belt (Hennings, 1994; McDowell and Mauger, 1994; Haenggi, 2002), which ultimately formed the southwestern margin of the Tornillo basin (Lehman, 1991). Similarly, the Sabinas basin in Coahuila was deformed and eroded during structural inversion (e.g., Eguiluz de Antuñano, 2001).

From the late Paleozoic to early Cenozoic, abundant zircon grains were produced in intermediate-to-felsic magmatic arcs along the western margin of North America, including the voluminous Mesozoic Cordilleran magmatic arc (Figs. 1 and 2). Permian–Triassic arc magmatism was located in the East Mexico arc from 287 to 232 Ma and encompassed felsic plutons in eastern Mexico as well as isolated plutons that extended northwest across north-central Mexico (Figs. 1 and 2; Torres et al., 1999; Coombs et al., 2020). By latest Permian–Triassic time, the nascent Cordilleran arc was established in northern Sonora (ca. 270–260 Ma) and southern California (ca. 275–210 Ma; Miller et al., 1995; Barth et al., 1997, 2011; Barth and Wooden, 2006; Arvizu et al., 2009; Cecil et al., 2019; Riggs et al., 2020; Dobbs et al., 2021). It remained active into the Jurassic, extending the length of eastern California/western Nevada and doglegging southeast into southeastern California, southern Arizona, and northern Sonora (Busby-Spera, 1988; Tosdal et al., 1989; Busby-Spera et al., 1990; Saleeby and Busby-Spera, 1992; Haxel et al., 2005). By the Early Jurassic, arc magmatism and diffuse back-arc magmatism, referred to as the Nazas arc, extended from Arizona and northern Sonora through north-central Mexico into eastern and southern Mexico (Fig. 1; e.g., Bartolini, 1998; Dickinson and Lawton, 2001b; Bartolini et al., 2003; Lawton and Molina Garza, 2014). It is best represented in Mexico by volcanic, volcaniclastic, and sedimentary rocks in north-central and northeastern Mexico (Jones et al., 1995; Bartolini, 1998; Bartolini et al., 2003; Barboza-Gudiño et al., 2008), as well as plutonic rocks and caldera complexes in southern Arizona and northern Sonora (Dickinson and Lawton, 2001b; Anderson et al., 2005). The Nazas arc is prominently represented in detrital zircon ages from late Mesozoic–early Cenozoic Mexican foreland basin strata as a ca. 180–150 Ma age group with high-amplitude modes from ca. 175 to 165 Ma (Lawton et al., 2009; Lawton and Molina Garza, 2014; Fitz-Díaz et al., 2018; Juárez-Arriaga et al., 2019).

During the Late Jurassic, contemporaneous with initiation of the Mexican Border rift, magmatism in the Nazas arc ceased (Dickinson and Lawton, 2001b; Stern and Dickinson, 2010), although Late Jurassic and earliest Cretaceous magmatism persisted along the Cordilleran margin of Mexico (e.g., Schaaf et al., 2020). Subsequently, the Cordilleran arc of the western United States experienced a major lull in magmatism during the Early Cretaceous from ca. 140 to 125 Ma (Armstrong and Ward, 1993; Christiansen et al., 1994; Barton, 1996; Ducea, 2001; Ducea and Barton, 2007; Barton et al., 2011; Paterson et al., 2011). In contrast, arc magmatism persisted in the form of the Early Cretaceous Alisitos arc and continental arcs within the Guerrero terrane in western Mexico (Fig. 1; Fitz-Díaz et al., 2018; Schaaf et al., 2020). The detrital zircon signature of the Alisitos arc is distinctive because it consists of an age mode centered at ca. 132 Ma that ranges from ca. 156 to 117 Ma (Martini et al., 2012; Fitz-Díaz et al., 2018). Voluminous continental arc magmatism resumed in the middle Cretaceous along the margin of the western United States and western Mexico. In detail, Cordilleran arc–related plutonic and volcanic rocks were emplaced in the southwestern United States and western Mexico as a series of batholiths, including the Sierra Nevada, the Mojave Desert–Salinia, the Peninsular Ranges, and the Sonora and Sinaloa batholiths (Fig. 1; Yonkee and Weil, 2015; Fitz-Díaz et al., 2018, and references therein). The Sierra Nevada arc was active from ca. 125 to 80 Ma and experienced particularly high rates of magmatic flux from ca. 100 to 85 Ma (Chen and Moore, 1982; Saleeby et al., 1987; Coleman and Glazner, 1997; Ducea, 2001; Ducea and Barton, 2007). The Mojave Desert–Salinian batholith is represented by granitoid ages that range from ca. 100 to 85 Ma in the Salinian block and from ca. 85 to 70 Ma in the western Mojave Desert region (Jacobson et al., 2011; Sharman et al., 2015). The Peninsular Ranges batholith ranges from ca. 140 to 80 Ma in age (Silver and Chappell, 1988; Kimbrough et al., 2001; Ortega-Rivera, 2003), but the La Posta suite, an especially voluminous exposure of plutonic rocks, represents unusually high magma production from ca. 99 to 92 Ma (Walawender et al., 1990; Kimbrough et al., 2001; Grove et al., 2003; Todd et al., 2003).

Voluminous continental arc magmatism along the western margin of the United States and Mexico was disrupted by Laramide flat-slab subduction under southern California by ca. 95–85 Ma (Livaccari et al., 1981; Liu et al., 2010). In the Sierra Nevada, arc magmatism ceased at ca. 80 Ma (Dumitru et al., 1991; Saleeby, 2003), and shallow subduction resulted in the uplift and exhumation of the Mojave Desert–Salinia segment of the Cordilleran arc and much of the retro-arc region (e.g., Malin et al., 1995; Saleeby, 2003). Laramide-aged magmatism progressively swept inboard (Figs. 1 and 2; Coney and Reynolds, 1977), starting with ca. 85–70 Ma magmatic rocks in the western Mojave region (Jacobson et al., 2011; Sharman et al., 2015) before migrating into eastern Sonora, southern Arizona, and southwestern New Mexico from ca. 80 to 50 Ma (the Laramide porphyry copper province; Valencia-Moreno et al., 2007; Greig and Barton, 2019), and culminating in ca. 75–43 Ma plutons in Colorado (Colorado Mineral Belt; Chapin, 2012). These inboard Laramide plutons are characterized by significant crustal contamination, as exhibited by abundant inherited zircon cores with U-Pb ages at 1700 Ma, 1650–1630 Ma, and 1450–1400 Ma as well as cores at 166–141 Ma and 124–112 Ma (Fig. 2; Valencia et al., 2005, 2006; Del Rio-Salas et al., 2013, 2017; Barra and Valencia, 2014; Clinkscales and Lawton, 2015; Amato et al., 2017; Mizer, 2018; Seedorff et al., 2019). Similarly, arc magmatism migrated eastward in northwest Mexico from the Peninsular Ranges to eastern Sonora between ca. 90 and 50 Ma, where it is known as the Laramide magmatic arc of Sonora (Figs. 1 and 2; Damon et al., 1983; McDowell et al., 2001; González-León et al., 2011). U-Pb zircon data of the Sonoran batholith constrain the timing of major Laramide magmatism between ca. 80 and 50 Ma, with evidence for Proterozoic inheritance from zircon cores (McDowell et al., 2001; González-León et al., 2011).

Strata of the Tornillo basin in Big Bend National Park in west Texas are among the few and best-preserved records of Upper Cretaceous–Eocene sedimentation between Arizona and the Gulf of Mexico (Fig. 3; Lehman, 1986, 1991). The Tornillo Group comprises an entirely nonmarine fluvial stratigraphic section, including the Maastrichtian Javelina, the Maastrichtian–Ypresian Black Peaks, and the Ypresian Hannold Hill Formations (Fig. 4). The basal Javelina Formation overlies shallow-marine strata of the Campanian–Maastrichtian Aguja Formation, which caps a stratigraphic succession that records the withdrawal of the Western Interior Seaway and migration of the Gulf of Mexico shoreline to the southeast (Maxwell et al., 1967; Schiebout et al., 1987; Lehman, 1991; Lehman et al., 2018). The upper Tornillo Group is unconformably overlain by fluvial strata of the Ypresian–Lutetian Canoe Formation, and the contact is characterized by highly variable erosional relief (Fig. 4). The most complete section of the Tornillo Group, preserving the Javelina, Black Peaks, and Hannold Hill Formations capped by the Canoe Formation, is found in the eastern part of Tornillo Flat (Fig. 3B). Carbon (δ¹⁸C) isotope measurements from pedogenic carbonate nodules and youngest detrital zircon U-Pb ages suggest that the Black Peaks and Hannold Hill Formations at Tornillo Flat were deposited between 67 and 52 Ma without a major hiatus (Bataille et al., 2016, 2019). A carbon isotope excursion associated with the Paleocene–Eocene thermal maximum has not been identified, but its predicted stratigraphic position based on magnetostratigraphic and carbon isotope tie points coincides with a sandstone (the “South Wall” sandstone) in the upper Black Peaks Formation (Fig. 4; Schiebout et al., 1987; Bataille et al., 2016, 2019). The Black Peaks and Hannold Hill Formations are unconformably overlain by the “Big Yellow” sandstone of the lowermost Canoe Formation at Tornillo Flat. However, in other parts of the Big Bend region, the Canoe Formation unconformably overlies strata as old as the Aguja Formation, and its basal contact is interpreted as a regional angular unconformity (Lehman et al., 2018). At Dawson Creek, only the Javelina Formation and an abbreviated Black Peaks Formation are preserved beneath the overlying Canoe Formation (Figs. 3C and 4). A combination of magnetostratigraphy, detrital sanidine ⁴⁰Ar/³⁹Ar dating, and biostratigraphy at Dawson Creek indicated that the upper Aguja through Black Peaks interval was deposited there from >74 Ma to 62.5 Ma in four 10⁵ yr intervals separated by multi-million-year depositional hiatuses (Leslie et al., 2018).

Variable pre-Canoe erosion was likely a primary consequence of the syntectonic setting of the Tornillo basin. The basin abuts several Laramide monoclines to the west, southwest, and northeast, all of which currently expose Lower Cretaceous carbonate rocks (Fig. 3; Cobb and Poth, 1980; Lehman, 1986, 1991). Intraformational thickness variations and angular stratal relationships indicate that these structures were active during deposition of the Tornillo Group in the latest Cretaceous through early Eocene. Increased uplift within the basin in the early Eocene resulted in erosional removal of parts of the upper Black Peaks Formation at Dawson Creek and local deposition of the Hannold Hill Formation at Tornillo Flat (Lehman et al., 2018). In contrast, the Canoe Formation shows no signs of syndepositional deformation (Lehman et al., 2018). The tectonic significance of the Canoe Formation has been debated, with some studies arguing that the lowermost Big Yellow sandstone of the Canoe Formation marks the cessation of Laramide-Mexican deformation in the Big Bend region (Maxwell et al., 1967; Runkel, 1988; Lehman, 1991; Lehman and Busbey, 2007; Lehman et al., 2018), whereas others envision the Big Yellow sandstone in the context of continued regional deformation (Schiebout, 1970; Hartnell, 1980; Rigsby, 1982, 1986; Rapp et al., 1983; Schiebout et al., 1987). Regardless, deposition in the Tornillo basin was largely disrupted by the advent of local magmatic activity in the middle Eocene (Maxwell et al., 1967; Lehman and Busbey, 2007; Lehman, 1991).

In total, 28 new sandstone samples were collected from Upper Cretaceous to Lower Eocene strata for detrital zircon laser ablation–inductively coupled plasma–mass spectrometry (LA-ICP-MS) U-Pb dating from three different localities in Big Bend National Park, west Texas (Fig. 3). Sandstone samples were stratigraphically correlated to published measured sections of Atchley et al. (2004), Bataille et al. (2016, 2019), Leslie et al. (2018), and Lehman et al. (2018) to refine the chronostratigraphic framework as well as to constrain sedimentary provenance (Fig. 4). Seventeen samples were collected from western and southern Tornillo Flat, spanning the uppermost Javelina, Black Peaks, Hannold Hill, and Canoe Formations. Nine samples were analyzed from Dawson Creek, including samples from the upper Aguja, Javelina, lower Black Peaks, and Canoe Formations, and two additional samples were collected from the Black Peaks Formation at Glenn Springs Draw (Figs. 3 and 4). Refer to Item S2 for sandstone sample descriptions (^{footnote 1}).

All sandstone samples were subjected to standard mineral separation techniques to isolate zircon, including crushing, water table heavy mineral concentration, heavy liquid density separation by bromoform and methylene iodide, and magnetic separation. For LA-ICP-MS depth-profile U-Pb analysis, unpolished zircons from each sample were then sprinkled onto 1-inch-diameter (~2.5-cm-diameter) acrylic disks with double-sided adhesive tape. From each sample, between 130 and 150 zircon grains were randomly selected for analysis to ensure a >95% confidence level that any age fraction >5% of the total population would be represented (Vermeesch, 2004). The zircons were analyzed by depth-profile LA-ICP-MS employing a Photon Machines Analyte G.2 excimer laser using a 30 μm spot size at 10 Hz in a large-volume Helex cell coupled to a Thermo Element 2 ICP-MS. This technique determined unique U-Pb ages from multiple zircon age domains with ablation depth into the zircon crystal ultimately resulting in core and rim U-Pb ages for single zircon crystals (Marsh and Stockli, 2015). GJ1 was used as the primary reference standard (601.7 ± 1.3 Ma; Jackson et al., 2004), while Pak1 (43.03 ± 0.01 Ma, in-house) or Plešovice (337.13 ± 0.37 Ma; Sláma et al., 2008) standards were used as secondary reference standards. The data were reduced using Iolite data reduction software and VizualAge (Paton et al., 2011; Petrus and Kamber, 2012). The reported date uncertainties (α; sensu Horstwood et al., 2016) are the internal precision for each analysis plus a propagated precision that added the uncertainty for excess scatter in the primary reference material during the analytical session (Paton et al., 2010). See the Maximum Depositional Age Calculation section for additional uncertainty handling protocols. For U-Pb dates younger than 850 Ma, the ²⁰⁶Pb/²³⁸U dates are reported as the U-Pb best age, and for dates older than 850 Ma, the ²⁰⁷Pb/²⁰⁶Pb date is reported. A ±15% ²⁰⁶Pb/²³⁸U versus ²⁰⁷Pb/²³⁵U discordance filter and 10% 2σ α uncertainty filter were applied to ²⁰⁶Pb/²³⁸U ages. A ±30% ²⁰⁶Pb/²³⁸U versus ²⁰⁷Pb/²⁰⁶Pb discordance filter was applied to ²⁰⁷Pb/²⁰⁶Pb ages. All mineral separation and LA-ICP-MS analyses were completed at the UTChron facilities at the Jackson School of Geosciences at the University of Texas at Austin. Refer to Item S3 in the Supplemental Material (^{footnote 1}) for detrital zircon U-Pb analytical results reported at 2σ α uncertainty.

Detrital zircon U-Pb age data are displayed as kernel density estimate (KDE) plots in Figures 5, 6, and 7, and detrital zircon age groups (also referred to as components) are presented in order of abundance in the Results section. The bin boundary between Hauterivian–Coniacian (130–87 Ma) and Coniacian–Ypresian (87–52 Ma) age groups was chosen based on the most prominent data gap in the KDE curves. Detrital zircon age distributions are displayed using KDE plots with a 1σ bandwidth of 2 m.y. for ages from 300 to 0 Ma and a 1σ bandwidth of 8 m.y. for ages from 2000 to 300 Ma. Colors under the KDE curves in Figures 5–7 match both the generalized basement age domains shown in Figure 2 and the detrital zircon age groups in the text—with the sole exception that age groups between 300 and 87 Ma are broken out in further detail in the text. In addition, the narrower color bars in Figures 5–7 correspond to specific age ranges of North American tectonic-magmatic events discussed in the text and shown in Figure 1, labeled at the top of each plot.

Depth-profile U-Pb dating systematically resolves multiple age domains within individual zircons. Core-rim age pairs that are resolved by depth profiling provide important insights into the magmatic and metamorphic history of detrital zircon source regions, and the unbiased random sample of the core-rim parent population generated by depth profiling gives more weight to the interpretation of their bivariate distribution. Hence, in addition to plotting these pairs as scatter plots, we also generated a bivariate KDE plot to guide interpretation (Fig. 8). KDE plots were made, and age modes were selected, using detritalPy (Sharman et al., 2018). Scatter plots, bivariate KDE plots, and ternary plots were generated using seaborn, matplotlib, and python-ternary libraries in Python (Hunter, 2007; Harper et al., 2019; Waskom, 2021).

The primary goal of detrital zircon maximum depositional age (MDA) calculation is to establish the true crystallization age of the youngest detrital zircon age component to provide a maximum age constraint for the true depositional age of a sample. The benefits and challenges of various routinely employed detrital zircon MDA methods are related to the (often contrasting) effects of the youngest age component sampling density and several sources of analytical, systematic, and geologic uncertainty (e.g., Coutts et al., 2019; Herriott et al., 2019; Copeland, 2020; Sharman and Malkowski, 2020; Vermeesch, 2021; Schwartz et al., 2022). The chronostratigraphic significance of an (accurate) MDA increases as it approaches stratal age, and the likelihood an MDA coincides with the true depositional age is foremost a function of the abundance of near-depositional-age zircons within the sampled unit and the total number, n, of zircon grains sampled (Coutts et al., 2019; Sharman and Malkowski, 2020). However, well-sampled, high-abundance, near-depositional-age zircons also increase the likelihood of obtaining zircon dates that are less than the true zircon crystallization age, commonly resulting in single-grain or low-n MDA calculations (e.g., youngest single grain or youngest three grains) that are younger than the true depositional age (Coutts et al., 2019; Sharman and Malkowski, 2020). Zircon dates that are too young can result from various factors, including analytical dispersion, laser-induced matrix effects, Pb loss, and/or sample contamination (Dickinson and Gehrels, 2009b; Allen and Campbell, 2012; Spencer et al., 2016; Herriott et al., 2019). More-conservative MDA methods (e.g., youngest cluster overlapping at 2σ uncertainty or youngest statistical population) that use multiple zircon ages are less prone to include or be heavily weighted by too-young zircon dates and thus are less likely to underestimate true depositional age (Coutts et al., 2019; Sharman and Malkowski, 2020). Conversely, when detrital zircon samples have a low abundance of near-depositional-age zircon grains, less-conservative MDA methods are more likely to coincide with true depositional age (Sharman and Malkowski, 2020). We calculated both the less-conservative youngest single grain age (YSG; Dickinson and Gehrels, 2009b) and the more-conservative weighted mean of the youngest statistical population (YSP; Coutts et al., 2019) of three or more ages that had a mean square of weighted deviation (MSWD) closest to 1 for each sample.

Detrital zircon samples from the Aguja, Javelina, and Black Peaks Formations yielded relatively high-abundance youngest age clusters such that the YSG ages typically overlapped within 2σ uncertainty of the YSP weighted mean ages (e.g., 16TO-20 and 16TO-01 in Figs. 9 and 10). Therefore, we report the more-conservative YSP age. Conversely, detrital zircon samples from the Hannold Hill and Canoe Formations have a low abundance of scattered youngest ages such that the YSP ages are significantly older than the YSG ages (e.g., 16TO-11 in Fig. 11). Therefore, we report the YSG age as the MDA, although we caution that the YSG ages exhibit characteristics of too-young biases. Both YSG and YSP ages for each sample are reported in Table 1. Prior to MDA calculation, detrital zircon analyses were filtered to <5% discordance to exclude possibly problematic analyses and then sorted by their age plus 2σ α uncertainty (Sharman et al., 2018; Sharman and Malkowski, 2020). Weighted mean ages were calculated using the date uncertainties (α). Following weighted mean calculation, external uncertainties, including the long-term excess variation (2.3%) in the secondary reference material (Plešovice), the age uncertainty of the primary reference material (GJ1; Jackson et al., 2004), and the decay constant uncertainty (Jaffey et al., 1971), were propagated in quadrature into the MDA uncertainty (β; sensu Horstwood et al., 2016). Both α and β uncertainties are reported at 2σ for each MDA using the convention MDA ± α/β. For better comparison, we recalculated previously published detrital zircon MDAs from Bataille et al. (2019) (see also Colliver, 2017) and Lehman et al. (2022) using the same discordance filters and MDA algorithms as the new detrital zircon MDAs. For the recalculated MDAs, we approximated the long-term excess variation at 2% (Horstwood et al., 2016). Detrital zircon MDA calculations were generated using detritalPy (Sharman et al., 2018).

This study presents a total of 3768 new detrital zircon U-Pb analyses from 28 sandstone samples (~135 analyses per sample), of which 3679 (~98%) analyses passed concordance filters (Item S3 [^{footnote 1}]). In summary, the detrital zircon U-Pb data from Campanian–Lutetian strata of the Tornillo basin consist of dominant Hauterivian–Coniacian (130–87 Ma; ~20%) and Coniacian–Ypresian (87–52 Ma; ~30%, age fraction [a] = 1118/3679) age groups with modes at 92 Ma and 78 Ma, respectively (Fig. 5A). Subordinate age components are Paleoproterozoic–Mesoproterozoic (1900–1300 Ma; ~12%) with modes at 1674 Ma and 1425 Ma, Ectasian–Tonian (1300–900 Ma; ~8%), Tonian–Pennsylvanian (900–300 Ma; ~10%) with a mode at 451 Ma, Permian–Triassic (300–200 Ma; ~8%) with a mode at 256 Ma, and Jurassic–Early Cretaceous (200–130 Ma; ~11%) with modes at 166 and 151 Ma. The remaining ages are Paleoarchean–Paleoproterozoic (3300–1900 Ma; ~2%) and Ypresian–Lutetian (52–44 Ma; ~0.2%).

Detrital zircon U-Pb ages from three samples (16TO-18, 16TO-17, and 16TO-19) from the Aguja Formation at Dawson Creek (Figs. 3C and 4) range from Archean to Maastrichtian. The ages are dominated by a Coniacian–Ypresian age component (~25%, a = 102/411) with a mode at 83 Ma (Figs. 5B and 6A). Subordinate age groups are Paleoproterozoic–Mesoproterozoic (~15%) with modes at 1676 and 1421 Ma, Ectasian–Tonian (~9%), Tonian–Pennsylvanian (~16%) with a mode at 447 Ma, Jurassic–Early Cretaceous (~10%) with a mode at 164 Ma, and Hauterivian–Coniacian (~17%) with a mode at 99 Ma. The remaining minor age components are Paleoarchean–Paleoproterozoic (~2%) and Permian–Triassic (~6%).

Detrital zircons in two samples (16TO-20 and 18TO-02) collected from the Javelina Formation at Dawson Creek and one sample (16TO-01) collected at Tornillo Flat (Figs. 3 and 4) show a range of ages from Archean to Paleocene. The dominant age group is Coniacian–Ypresian (~42%, a = 166/399) with a mode at 70 Ma (Figs. 5B, 6A, and 7). Subordinate age groups are Paleoproterozoic–Mesoproterozoic (~12%) and Hauterivian–Coniacian (~18%) with a mode at 93 Ma. The remaining age groups (<10% each) are Paleoarchean–Paleoproterozoic (~2%), Ectasian–Tonian (~6%), Tonian–Pennsylvanian (~7%), Permian–Triassic (~6%), and Jurassic–Early Cretaceous (~8%).

Lower Black Peaks detrital zircons for one sample (18TO-06) from Dawson Creek, two samples (16TO-21 and 16TO-22) from Glenn Springs Draw, and three samples (16TO-02, 16TO-03, and 16TO-04) from west Tornillo Flat (Figs. 3 and 4) exhibit ages ranging from Archean to Paleocene. The two dominant age components are an Hauterivian–Coniacian (~22%) age group with a mode at 92 Ma and a Coniacian–Ypresian (~32%, a = 245/773) age group with modes at 77 and 64 Ma (Figs. 5B, 6A, 6B, and 7). Subordinate age components are Paleoproterozoic–Mesoproterozoic (~12%) and Jurassic–Early Cretaceous (~10%) with modes at 181, 165, and 150 Ma. The remaining age groups (<10% each) are Paleoarchean–Paleoproterozoic (~1%), Ectasian–Tonian (~7%), Tonian–Pennsylvanian (~8%), and Permian–Triassic (~9%) with a mode at 258 Ma.

Detrital zircon ages of six samples (16TO-05, 16TO-06, 16TO-07, 16TO-09, 16TO-10, and 16TO-08) from the upper Black Peaks Formation at Tornillo Flat (Figs. 3B and 4) range from Archean to Paleocene. The dominant Coniacian–Ypresian (~28%, a = 208/743) age group has modes at 79 and 67 Ma (Figs. 5B and 7). Subordinate age groups are Paleoproterozoic–Mesoproterozoic (~14%) with modes at 1657 and 1410 Ma, Tonian–Pennsylvanian (~9%), Jurassic–Early Cretaceous (~12%) with modes at 172 and 152 Ma, and Hauterivian–Coniacian (~18%) with a mode at 94 Ma. The remaining minor age groups (<10% each) are Paleoarchean–Paleoproterozoic (~3%), Ectasian–Tonian (~8%), and Permian–Triassic (~7%) with a mode at 265 Ma.

Detrital zircons of four samples (16TO-11, 16TO-15, 18TO-01, and 16TO-12) from the Hannold Hill Formation at Tornillo Flat (Figs. 3B and 4) show ages ranging from Archean to Eocene. The dominant age groups are Hauterivian–Coniacian (~23%) and Coniacian–Ypresian (~31%, a = 160/520) with modes at 89 Ma and 80 Ma, respectively (Figs. 5B and 7). Subordinate age groups include Paleoproterozoic–Mesoproterozoic (~9%), Tonian–Pennsylvanian (~10%) with a mode at 450 Ma, Permian–Triassic (~10%) with a mode at 265 Ma, and Jurassic–Early Cretaceous (~11%) with modes at 166 and 155 Ma. The remaining age groups are Paleoarchean–Paleoproterozoic (~1%) and Ectasian–Tonian (~6%).

Detrital zircon ages of three samples (18TO-04, 18TO-03, and 18TO-05) from the lower Canoe Formation at Dawson Creek and three samples (16TO-13, 16TO-14, and 16TO-16) from the Big Yellow sandstone of the Canoe Formation at Tornillo Flat (Figs. 3 and 4) range from Archean to Eocene. The dominant Coniacian–Ypresian (~28%, a = 237/833) age group has a mode at 78 Ma (Figs. 5B, 6A, and 7). Various subordinate age groups are Paleoproterozoic–Mesoproterozoic (~12%) with modes at 1673 and 1431 Ma, Ectasian–Tonian (~8%), Tonian–Pennsylvanian (~9%) with a mode at 453 Ma, Permian–Triassic (~10%) with a mode at 255 Ma, Jurassic–Early Cretaceous (~11%) with modes at 167 and 150 Ma, and Hauterivian–Coniacian (~19%) with a mode at 93 Ma. The remaining very minor age groups are Paleoarchean–Paleoproterozoic (~1%) and Ypresian–Lutetian (~1%).

Out of 3768 detrital zircon analyses, 216 analyses (~6%) represented 108 concordant core-rim age pairs from individual zircon grains (Fig. 8). Considering the relative paucity of core-rim age pairs from each sample and the general similarity of detrital zircon age distributions throughout the Campanian–Lutetian section, we combined the age pairs into an aggregate data set. The large sample size of the aggregate data set reduced uncertainty and allowed for more robust interpretation. The core-rim age pairs have two dominant age-pair clusters comprising Ectasian–Tonian cores with Paleozoic rims (~19%, a = 21/108) and Jurassic cores with middle to Late Cretaceous rims (~18%) (Fig. 8). Subordinate age-pair clusters include Paleoproterozoic–Mesoproterozoic cores with Late Cretaceous–Paleocene (~9%), Permian–Triassic (~7%), and Ectasian–Tonian rims (~7%) and Cretaceous cores with Cretaceous rims (~6%). Various minor age-pair clusters include: Paleoarchean–Paleoproterozoic cores with Paleoproterozoic–Mesoproterozoic rims (~3%), Paleoproterozoic–Mesoproterozoic cores with Paleoproterozoic–Mesoproterozoic (~5%), Paleozoic (~3%), and Jurassic rims (~2%); Ectasian–Tonian cores with Ectasian–Tonian rims (~3%); Ectasian–Tonian (~5%) and Tonian–Pennsylvanian (~5%) cores with Late Cretaceous–Paleocene rims; and Neoproterozoic cores with Tonian–Pennsylvanian rims (~4%) (Fig. 8).

This study provides new detrital zircon U-Pb MDA constraints on the chronostratigraphic framework of the Tornillo basin (Table 1; Figs. 9–13; Items S4–S7 [^{footnote 1}]). The detrital zircon data from the Aguja, Javelina, and Black Peaks Formations has relatively high-abundance youngest age clusters such that the YSG and YSP ages typically overlap within 2σ uncertainty (Figs. 9–11). Furthermore, MDAs from the Aguja, Javelina, and Black Peaks Formations get progressively younger up section, suggesting that they provide robust Late Cretaceous–Paleocene near-depositional-age constraints (Figs. 12 and 13). This notion is corroborated by the close match between the MDAs for these formations and existing depositional age constraints from magnetostratigraphy, detrital sanidine ⁴⁰Ar/³⁹Ar dating, and biostratigraphy (Figs. 12 and 13; Leslie et al., 2018; Bataille et al., 2019). Therefore, we consider that near-depositional-age zircons from the Aguja, Javelina, and Black Peaks Formations were well sampled. In contrast, the youngest detrital zircon U-Pb ages for the Hannold Hill and Canoe Formations are scattered and low in abundance such that the YSP ages are significantly older than the YSG ages (Fig. 11). All Hannold Hill and Canoe Formation YSP ages are not chronostratigraphically relevant, and the YSG ages, although sometimes consistent with preexisting age constraints, do not steadily get younger up section, indicating that near-depositional-age zircons were sparsely sampled (Figs. 12 and 13). Therefore, the Hannold Hill and Canoe Formation MDAs were of limited use in constraining the timing of Eocene deposition in the Tornillo basin, and the YSG ages that did appear to coincide with near-depositional-age constraints were interpreted with caution toward too-young biases.

Detrital zircon U-Pb MDAs from the Aguja, Javelina, and Black Peaks Formations at Dawson Creek provide support for existing age constraints and models that indicate the Dawson Creek stratigraphic section was typified by punctuated deposition separated by multi-million-year depositional hiatuses (Leslie et al., 2018; Lehman et al., 2022), as well as providing new constraints bracketing the location of the Cretaceous-Paleogene boundary. Upper Aguja Formation samples 16TO-17 and 16TO-19 at Dawson Creek yielded Campanian–Maastrichtian 78.6 ± 0.9/2.0 Ma and 72.0 ± 1.1/2.0 Ma YSP ages consistent with Judithian North American Land Mammal ages (NALMA; Cifelli et al., 2004) and 76.9 ± 1.2 Ma (Befus et al., 2008) and 72.6 ± 1.5 Ma (Breyer et al., 2007) volcanic deposits from other Aguja Formation locales (Table 1; Figs. 9 and 12). At Dawson Creek, Leslie et al. (2018) correlated the local polarity stratigraphy from the upper Aguja through the lower Javelina interval to C31r and C31n of the geomagnetic polarity time scale (GPTS), indicating that the interval was deposited between 71.5 and 68.4 Ma, which is reasonably corroborated by the 72.0 ± 1.1/2.0 Ma YSP age from sample 16TO-19 (Fig. 12). Leslie et al. (2018) correlated the local polarity stratigraphy for the upper Javelina Formation to C29r, indicating that the interval was deposited between 66.4 and 65.7 Ma and stipulating an unconformity representing a >2 m.y. hiatus between the lower and upper Javelina Formation. Sample 16TO-20, collected from the middle of the Javelina Formation near the base of the upper Javelina interval, yielded a YSP age of 66.3 ± 0.5/1.6 Ma, supporting the local polarity correlation of Leslie et al. (2018) and the presence of a multi-million-year hiatus within the middle Javelina Formation (Table 1; Figs. 9 and 12).

The local polarity stratigraphy of the uppermost Javelina–lower Black Peaks contact at Dawson Creek has been magnetically overprinted, complicating the time scale correlation of the interval (Leslie et al., 2018). Below the magnetically overprinted interval, the upper Javelina Formation contains Lancian dinosaur fossils, indicative of a latest Cretaceous age (Fig. 12). Wiest et al. (2018) located the Cretaceous-Paleogene boundary in the upper Javelina Formation just below the magnetically overprinted interval at 172 m based on an abrupt decrease in burrow diameter that may represent a postextinction recovery community (Fig. 12). Above the magnetically overprinted interval, the abbreviated section of Black Peaks Formation contains mammal fossils of the late Danian Torrejonian NALMA (Leslie et al., 2018; Lehman et al., 2022), and Leslie et al. (2018) correlated the Black Peaks Formation’s local polarity stratigraphy to C27r (63.5–62.5 Ma) of the GPTS. This correlation implies the presence of a >2.2 m.y. hiatus within the magnetically overprinted interval.

Sample JC 8, from Lehman et al. (2022), yielded a robust (n = 9) recalculated Danian YSP age of 62.5 ± 0.2/1.3 Ma out of a lithologically unique conglomerate (the “odd conglomerate”) at nearby Rough Run Creek (Table 2; Fig. 9). This conglomerate crops out at the base of the magnetically overprinted interval at Dawson Creek, indicating that the overlying overprinted interval is Torrejonian, and that a disconformity representing an ~3 m.y. hiatus is located within the uppermost Javelina Formation at or just below the odd conglomerate (Figs. 4 and 12). Sample 18TO-02, collected from the uppermost sandstone of the Javelina Formation just above the odd conglomerate, yielded a YSP age of 61.2 ± 0.6/1.5 Ma, corroborating a Torrejonian depositional age for the magnetically overprinted interval and the placement by Lehman et al. (2022) of an ~3 m.y. hiatus at the base of the interval within the uppermost Javelina Formation at Dawson Creek (Table 1; Figs. 9 and 12).

These results bracket the Cretaceous-Paleogene boundary at Dawson Creek below the magnetically overprinted interval and above the last in situ dinosaur fossils in the upper Javelina Formation (Fig. 12). Lehman et al. (2022) noted the presence of in situ dinosaur fossils in a mudstone interval below the odd conglomerate at nearby Rough Run Creek and suggested that the Cretaceous-Paleogene boundary at Dawson Creek may be located at the disconformity surface. The high abundance of near-depositional detrital zircon ages in the sandstones across the Cretaceous-Paleogene boundary indicates that chemical abrasion–isotope dilution–thermal ionization mass spectrometry (CA-ID-TIMS) U-Pb dating for detrital zircon MDAs could be used to refine the chronostratigraphy of the upper Javelina–lower Black Peaks section at Dawson Creek (e.g., Herriott et al., 2019).

Detrital zircon MDAs from the Javelina and Black Peaks Formations at Tornillo Flat provide new age constraints on the location of the Cretaceous-Paleogene boundary and the Paleocene–Eocene thermal maximum at Tornillo Flat and lend support to existing age models of the Javelina–Black Peaks interval (Bataille et al., 2019). Sample 16TO-01, collected from the uppermost Javelina Formation, yielded a robust (n = 21) Maastrichtian YSP age of 67.8 ± 0.4/1.6 Ma; this is consistent with Lancian fossils reported from this interval, and it corroborates the carbon isotope tie point of Bataille et al. (2019), which indicated the top Javelina was deposited at 68 Ma (Table 1; Figs. 10 and 13). In situ Lancian fossils in the lowermost Black Peaks Formation at nearby Grapevine Hills indicate that the Cretaceous-Paleogene boundary is likely within the lower Black Peaks Formation at Tornillo Flat (Lehman et al., 2018, 2022). Lower Black Peaks sample 113015LC-02 yielded a robust (n = 13) Danian YSP age of 63.4 ± 0.4/1.3 Ma (recalculated from Bataille et al., 2019) that provides a firm upper bracket on the location of the Cretaceous-Paleogene boundary (Table 2; Figs. 10 and 13). Conversely, sample 16TO-02, collected from the lowermost Black Peaks Formation, yielded a similarly robust (n = 7) Maastrichtian YSP age of 67.6 ± 0.5/1.7 Ma (Fig. 10). Although the exact correlation of 16TO-02 into the stratigraphic section of Bataille et al. (2019) is uncertain, it is likely that the sample was collected between 10 and 30 m below sample 113015LC-02. Given the robust n values of the YSP ages for both samples, we consider that their YSP ages are closely representative of the true depositional age, and, therefore, samples 16TO-02 and 113015LC-02 bracket the Cretaceous-Paleogene boundary at Tornillo Flat (Fig. 13). Furthermore, both the proximity of samples 16TO-02 and 113015LC-02 to each other in the stratigraphic section and the ~4 m.y. difference in their YSP ages support the presence of a major disconformity within the lowermost Black Peaks Formation, consistent with the apparent absence of early Danian Puercan fauna assemblages (Fig. 13; Bataille et al., 2019; Lehman et al., 2022). The abundant near-depositional-age zircons in both 16TO-02 and 113015LC-02 suggest that high-resolution sampling of sandstones and detrital zircon CA-ID-TIMS U-Pb dating could further constrain the location of the Cretaceous-Paleogene boundary disconformity at Tornillo Flat.

The “log jam sandstone” is a regionally correlative marker interval within the middle of the Black Peaks Formation (identifiable by the presence of petrified tree logs up to 10–15 m in length) that was herein used to informally split the formation into lower and upper units (Fig. 4; Lehman et al., 2018). The log jam sandstone interval contains Torrejonian–Tiffanian NALMAs, placing the interval within the Danian–Selandian (Lehman et al., 2018, 2022; Bataille et al., 2019), and a carbon isotope tie point correlates the interval to 60.6 Ma (Fig. 13; Bataille et al., 2019). Sample 16TO-03, collected from a sandstone just below the log jam sandstone, yielded a YSP age of 62.9 ± 1.2/1.9 Ma, corroborating the Danian–Selandian age for the interval (Table 1; Figs. 10 and 13).

The upper Black Peaks Formation from the log jam sandstone interval to the base of the Paleocene–Eocene thermal maximum interval has yielded Selandian–Thanetian Tiffanian NALMAs (Lehman et al., 2018; Bataille et al., 2019), and Bataille et al. (2019) placed carbon isotope tie points at 58.7, 57.1, and 56.6 Ma through the interval. The 58.7 and 57.1 Ma tie points were further constrained by two thin intervals of geomagnetic normal polarity correlated to top C26n and top C25n (59 and 57.1 Ma) in the GPTS (Schiebout et al., 1987; Bataille et al., 2016). Sample 16TO-05, collected just above the log jam sandstone, sample 16TO-06, collected midway up the interval, and sample 16TO-10, collected just below the Paleocene–Eocene thermal maximum interval, yielded upward-younging Danian–Selandian YSP ages of 62.7 ± 1.0/1.7 Ma, 60.5 ± 0.9/1.6 Ma, and 59.6 ± 0.6/1.5 Ma, respectively (Figs. 11 and 13). Although these YSP ages are generally consistent with a Selandian–Thanetian age for the interval, they are notably 3–4 m.y. older than the corresponding carbon isotope tie points of Bataille et al. (2019). It is noteworthy that the YSP age n values for these three samples are relatively low (n = 4, 6, and 4, respectively), and that three other samples collected from that interval (16TO-07, 16TO-09, and 16TO-08) returned MDAs older than stratigraphically underlying MDAs (Table 1; Fig. 13). These older YSP ages and low n values indicate that near-depositional-age zircons were sparsely sampled or absent through the upper Black Peaks interval, limiting the use of detrital zircon MDAs as age constraints for true depositional age. We suggest that ultrahigh-density (n = 300–1000) sampling of detrital zircon samples from the upper Black Peaks interval could yield more near-depositional-age zircon to further resolve and better evaluate the upper Black Peaks age model.

A combination of biostratigraphy (NALMAs), magnetostratigraphy, and carbon isotope chemostratigraphy has located the Paleocene–Eocene thermal maximum horizon (56 Ma; Westerhold et al., 2020) within the Black Peaks Formation at the “south wall sandstone” (Figs. 4 and 13; Bataille et al., 2016, 2019). However, an unequivocal location of the Paleocene–Eocene thermal maximum horizon has not been possible due to the lack of absolute age dating constraints and the absence of isotopic excursions within the Paleocene–Eocene thermal maximum interval that would explicitly mark the Paleocene–Eocene thermal maximum and early Eocene hyperthermal events (Lehman et al., 2018; Bataille et al., 2019). Bataille et al. (2019) reported a youngest detrital zircon age of 54.5 ± 1.0 Ma from the south wall sandstone (sample 120115LC-04), in support of locating the Paleocene–Eocene thermal maximum horizon there. We recalculated the MDA for sample 120115LC-04 and report a robust (n = 7) YSP age of 55.2 ± 0.4/1.2 Ma, strongly corroborating the location of the Paleocene–Eocene thermal maximum horizon at the south wall sandstone or in a proximally underlying mudstone (Table 2; Figs. 11 and 13). Given that YSP MDAs are less prone to too-young biases than a single detrital zircon age, and that the YSP age for sample 120115LC-04 is younger than the Paleocene–Eocene thermal maximum (but overlaps within β uncertainty), it would be difficult to locate the Paleocene–Eocene thermal maximum much higher in the section than the south wall sandstone. CA-ID-TIMS U-Pb dating of the youngest zircons from the south wall sandstone could confidently resolve whether the south wall sandstone slightly postdates or coincides with the Paleocene–Eocene thermal maximum.

Unfortunately, older MDAs from the Hannold Hill and Canoe Formations provide limited new constraints on the age models of the uppermost part of the Tornillo Group (Lehman et al., 2018; Bataille et al., 2019). For example, the basal Exhibit Sandstone Member of the Hannold Hill Formation yielded in situ Wasatchian fauna, placing the formation firmly within the Ypresian (Wilson and Runkel, 1989). However, detrital zircon samples collected from the Hannold Hill Formation all yielded Maastrichtian or older YSP ages (Table 1; Figs. 11 and 13). Similarly, early Uintan vertebrate fauna suggests the basal Big Yellow Sandstone Member of the Canoe Formation was deposited during the Lutetian (Robinson et al., 2004; Lehman et al., 2018), but detrital zircon samples all yielded Selandian or older YSP ages (Table 1; Figs. 11 and 13). Even the YSG ages from the Hannold Hill and Canoe Formations are frequently older than the robust 55.2 ± 0.4/1.2 Ma YSP age from the south wall sandstone in the upper Black Peaks Formation, effectively making most of the uppermost Tornillo Group YSG ages irrelevant as age constraints (Fig. 13). We attribute these older MDAs to limited-to-scarce near-depositional-age zircon within the Hannold Hill and Canoe Formations (e.g., Coutts et al., 2019).

The low abundance of near-depositional-age zircon necessitates a careful assessment of the potential effects of too-young biases, such as Pb loss, in the few Eocene Hannold Hill and Canoe Formation YSG ages that have been obtained. For example, sample 16TO-16 yielded a YSG age of 49.5 ± 0.8/1.4 Ma from the Big Yellow Sandstone Member (Table 1; Figs. 11 and 13). This YSG age falls within the late Ypresian age interval (51.6–48.1 Ma) estimated by Bataille et al. (2019) for the Big Yellow Sandstone Member (this age interval is also based on two youngest detrital zircon ages), but it does not overlap with the Lutetian as indicated by the NALMAs (Robinson et al., 2004). However, the YSG age for 16TO-16 is the youngest date of a long tail of decreasing dates from ca. 70 to 50 Ma (Fig. 11). For U-Pb ages younger than 100 Ma, zircon can experience significant Pb loss, and the analyses will remain nominally concordant but form a negatively skewed tail of too-young dates (Spencer et al., 2016). Most Hannold Hill and Canoe samples show this date pattern, with a major cluster of dates between ca. 80 and 70 Ma and then a long tail of dates between ca. 70 and 50 Ma (e.g., samples 16TO-11 and 16TO-16; Fig. 11; Items S4 and S7). Conversely, many Aguja, Javelina, and Black Peaks samples show distinct youngest age clusters that do not have prominent negatively skewed tails (e.g., samples 16TO-01, 16TO-02, 16TO-06, and 16TO-10; Figs. 9–11; Items S4–S5). Based on these date patterns, it is likely that Pb loss contributed to the youngest dates in the Hannold Hill and Canoe samples. In a high Pb-loss scenario, it is possible that the youngest detrital zircon age component in the Hannold Hill and Canoe Formations has a true age between ca. 80 and 70 Ma. Therefore, we suggest that the YSG ages do not provide useful constraints on the depositional age of the Hannold Hill and Canoe Formations, although CA-ID-TIMS U-Pb dating of the youngest detrital zircon grains analyzed by LA-ICP-MS could be used to test this Pb-loss hypothesis. Refer to Items S4–S8 in the Supplemental Material (^{footnote 1}) for ranked date plots for samples not discussed in the text.

Detrital zircon MDAs from the upper Javelina and lower Black Peaks Formations at Dawson Creek and Tornillo Flat provide further evidence for a basinwide ~3–4 m.y. hiatus in the early Danian (Figs. 12 and 13). Faunal assemblages from the upper Javelina–lower Black Peaks interval transition from Lancian (Maastrichtian) to middle Torrejonian (late Danian) NALMAs, indicating that the entire early Danian Puercan faunal assemblage is missing from the stratigraphy (Leslie et al., 2018; Lehman et al., 2022). In western stratigraphic sections of Big Bend National Park, Lancian fossils are separated from Torrejonian fossils by the “odd conglomerate” (Lehman et al., 2022). At Rough Run Creek, the odd conglomerate (sample JC 8) yielded a recalculated YSP age of 62.5 ± 0.2/1.3 Ma, which, juxtaposed overtop an interval containing Lancian fossils, indicates that the conglomerate represents a disconformity and an ~3–4 m.y. hiatus (Lehman et al., 2022). At Dawson Creek, a 61.2 ± 0.6/1.5 Ma YSP age (sample 18TO-02) from just above the odd conglomerate, which again overlies a Lancian fossil interval, supports the inference that the conglomerate represents a disconformity having a >3–4 m.y. hiatus at that location as well. The odd conglomerate has not been identified in eastern exposures of the stratigraphy, but based on the absence of Puercan fauna and a lower Black Peaks Formation YSP age of 63.4 ± 0.4/1.3 Ma (sample 113015LC-02), Bataille et al. (2019) suggested that a disconformity is present in the lowermost Black Peaks Formation at Tornillo Flat. A Maastrichtian YSP age of 67.6 ± 0.5/1.7 Ma (sample 16TO-02) from 10 to 30 m below sample 113015LC-02 potentially brackets the location of a disconformity separating Maastrichtian and upper Danian strata and representing an ~3–4 m.y. hiatus within the lowermost Black Peaks Formation. Altogether, the NALMA faunal assemblages and detrital zircon MDAs from multiple locales of the upper Javelina–lower Black Peaks interval indicate the presence of a basinwide ~3–4 m.y. hiatus representing early Danian time. Both eustatically and tectonically driven changes in fluvial base level have been suggested for the cause of the early Danian basinwide hiatus in sedimentation (Lehman et al., 2022). In particular, uplift along the Terlingua monocline was responsible for the erosional removal of the middle and upper parts of the Black Peaks Formation at Dawson Creek during the early to middle Eocene (Lehman et al., 2018). An earlier Late Cretaceous–Paleocene period of uplift of the monocline may explain the abbreviated Javelina–Black Peaks section at Dawson Creek, the restriction of the odd conglomerate to the western part of the basin, and the basinwide early Danian depositional hiatus (Lehman et al., 2022).

In recent mapping, the Javelina–Black Peaks formational contact was placed between the upper “sandstone-dominated” interval of the Javelina Formation and the lower “mudstone-dominated” interval of the lower Black Peaks Formation (Turner et al., 2011; Lehman et al., 2018). This choice made practical sense for mapping purposes, as the contact is generally recognizable, and it had the added benefit of making the Javelina Formation mostly Upper Cretaceous and the Black Peaks Formation mostly Paleocene in age. However, it was recognized that this contact is gradational, and, as Javelina Formation sandstone bodies are generally lenticular, the stratigraphic level of the top Javelina sandstone varies laterally (Lehman et al., 2018). The updated chronostratigraphic framework presented herein recognizes a basinwide lower Danian disconformity within the lowermost Black Peaks Formation everywhere except at Dawson Creek, where the disconformity is located in the uppermost Javelina Formation (Figs. 12 and 13). This result indicates that the Javelina–Black Peaks lithostratigraphic formational contact is significantly time-transgressive: upper Danian–Selandian at Dawson Creek and Maastrichtian at other well age-constrained locales in the basin. Furthermore, the new age constraints indicate that the uppermost ~15 m section of the Javelina Formation at Dawson Creek, starting from the odd conglomerate and including sample 18TO-02, should be considered cogenetic with lower Black Peaks strata.

Previous provenance work indicated that strata of the Tornillo basin were sourced from both proximal and distal sources in western North America (Lehman, 1991; Colliver, 2017; Bataille et al., 2019). Paleocurrent indicators show generally southeastward-directed paleoflow that persisted through deposition of the upper Aguja Formation, the Tornillo Group, and the Canoe Formation (Schiebout, 1970, 1974; Lehman, 1986; Rigsby, 1986; Haenggi, 2001). Aguja Formation sandstone compositions, dominated by volcanic rock fragments, fresh plagioclase feldspar, and quartz, led Lehman (1991) to hypothesize that the Aguja sediment was derived from a Cretaceous magmatic arc in western Mexico (Fig. 14A). Sandstone compositions of the Javelina and Black Peaks Formations are also rich in volcanic rock fragments, plagioclase feldspar, and quartz, but, importantly, they also include subordinate sedimentary rock fragments such as carbonate rock fragments, chert, and reworked fossils (Figs. 14A and 14B; Lehman, 1991). These compositional data suggest that Tornillo Group strata were sourced from distant volcanic and plutonic rocks, such as the Late Cretaceous magmatic rocks of New Mexico and Arizona, distant Lower Cretaceous chert-bearing carbonates, and nearby Upper Cretaceous marine shales and carbonates (Lehman, 1987, 1991). Importantly, the inclusion of local sedimentary sources is consistent with the onset of Laramide-Mexican deformation in the Big Bend region.

Lehman (1991) documented a compositional shift in the early Eocene as Hannold Hill and Canoe Formation conglomerates contain abundant limestone, sandstone, and chert clasts. Similarly, sandstone compositions are dominated by carbonate rock fragments, chert, quartz, and orthoclase and plagioclase feldspar, signaling a reduced contribution from distant magmatic sources and a dominant contribution from nearby sedimentary sources, consistent with peak early Eocene Laramide-Mexican shortening in the area (Figs. 14A and 14B). Conglomerates and sandstones in the Big Yellow sandstone also contain pebbles of Paleozoic chert and metamorphic rock fragments, respectively, that are not found in the underlying Hannold Hill, as well as rounded limestones and reworked petrified wood likely recycled out of underlying Tornillo Group strata (Lehman, 1991). Overall, Tornillo basin sandstone and conglomerate compositions indicate a sediment provenance shift from predominantly distal volcanic, plutonic, and sedimentary sources in the Paleocene to local syntectonic sedimentary sources in the Eocene. The first detrital zircon U-Pb age data from Tornillo basin strata in the Big Bend region were published by Bataille et al. (2019), from which they extracted the youngest ages to refine their Tornillo Flat chronostratigraphic age model. They noted that the predominant detrital zircon U-Pb age groups are consistent with derivation from western U.S. and Mexican Mesozoic magmatic sources, including Laramide, Cordilleran, and East Mexico arcs, as well as from Laramide basement uplifts of the southwestern United States, and they interpreted that the Tornillo basin drained portions of northwestern Mexico, west Texas, southern New Mexico, and southern Arizona (Colliver, 2017; Bataille et al., 2019).

Detrital zircon U-Pb data in this study from Campanian–Lutetian strata of the Tornillo basin show a dominance of ages younger than 300 Ma (~69%, n = 2529/3679; Fig. 5A), which prominently reflect magmatic sources in the western North American Cordilleran arc, including the nascent Permian–Triassic Cordilleran, the Jurassic Cordilleran-Nazas, the middle Cretaceous Cordilleran, and the Late Cretaceous–Eocene Laramide magmatic arcs (Figs. 1 and 2). Subordinate Paleoproterozoic–Mesoproterozoic ages with modes at 1674 and 1425 Ma (~12%; Fig. 5A) were likely derived from Yavapai-Mazatzal basement and ca. 1400 Ma granite-rhyolites of the midcontinent-to-western United States and northwestern Mexico (Fig. 2). The remaining Ectasian–Tonian and Tonian–Pennsylvanian (mode at 451 Ma) age groups (~17%) likely represent zircons recycled from Phanerozoic sedimentary units of the western United States and northern Mexico. These provenance interpretations are in broad agreement with previous detrital zircon U-Pb provenance interpretations for the Tornillo basin (Colliver, 2017; Bataille et al., 2019). However, this extensive data set combined with new information from core-rim age pairs (Fig. 8), which can fingerprint zircons by the spatial overlap of magmatic arcs and basement terranes, allowed for improved provenance interpretations that isolated specific source terranes and source regions that supplied sediment to the Tornillo basin and refined Late Cretaceous–Eocene sediment-routing reconstructions for the southwest United States and northern Mexico.

The Paleoproterozoic–Mesoproterozoic age component (~12%, n = 457/3679) is represented in Tornillo basin strata by two distinct subcomponents with age modes at 1674 and 1425 Ma, respectively (Figs. 5–7). The Statherian subcomponent reflects zircons from the Yavapai-Mazatzal basement province (ca. 1800–1620 Ma), whereas Calymmian–Ectasian zircons are diagnostic of the midcontinent Granite-Rhyolite Province (ca. 1480–1340 Ma), which largely intruded Yavapai-Mazatzal basement (Fig. 2). These two subcomponents are genetically linked and occur together in the Phanerozoic detrital zircon record and are subsequently referred to as Yavapai-Mazatzal. Possible first-cycle source regions for Yavapai-Mazatzal zircons in strata of the Tornillo basin include the relict Mogollon highlands, Laramide basement–cored uplifts in the Arizona–New Mexico–Sonora border region, and/or basement-involved uplifts in the southern Laramide province of northern New Mexico and Colorado (Figs. 1 and 2). The Mogollon highlands are composed of Yavapai-Mazatzal basement that was exhumed as part of the northern rift shoulder to the Bisbee basin during Mexican Border rifting in the Late Jurassic (Bilodeau, 1986; Dickinson and Lawton, 2001b). Basement horst blocks within the Mexican Border rift contributed Yavapai-Mazatzal zircons into the McCoy and Bisbee basins into the Cretaceous (Barth et al., 2004; Spencer et al., 2011; Lawton et al., 2020), while the Mogollon highlands contributed detrital zircons northward into the Sevier foreland (e.g., Dickinson et al., 2012; Szwarc et al., 2015; Pecha et al., 2018; Pujols and Stockli, 2021). Mogollon topography may have become subdued by the Late Cretaceous as marine transgressions covered the region (Lawton et al., 2020), but it was renewed during Laramide uplift and Mexican Border rift reactivation (e.g., Favorito and Seedorff, 2021), as evidenced by Yavapai-Mazatzal zircons reaching the Cordilleran forearc (Jacobson et al., 2011; Sharman et al., 2015) and the appearance of quartzite clasts of the Proterozoic Mazatzal Group in the Eocene Mogollon Rim Formation (Potochnik, 1989).

In the southern Laramide province of northern New Mexico, detrital zircons from the Cretaceous–Paleogene San Juan, Galisteo–El Rito, and Raton basins show increasing proportions of Paleoproterozoic–Mesoproterozoic ages as Yavapai-Mazatzal basement was progressively exposed by Laramide basement–involved uplifts starting in the Maastrichtian and continuing through the Eocene (Figs. 1, 15, and 16; Bush et al., 2016; Pecha et al., 2018; Smith et al., 2020). Although Laramide rivers draining these basins have been interpreted to drain to the south-southeast toward Texas (e.g., Lawton, 2019), the dominance of Yavapai-Mazatzal ages in detrital zircon data from these basins is incompatible with detrital zircon ages of Tornillo basin strata. It is therefore unlikely that rivers from the San Juan, Galisteo–El Rito, and/or Raton basins were routed to the Tornillo basin, as this would require significant mixing with rivers from other source regions to dramatically dilute the Yavapai-Mazatzal zircon signal (Figs. 15 and 16). Furthermore, the presence of Paleoproterozoic–Mesoproterozoic zircon cores with Permian–Triassic and Late Cretaceous–Paleocene rims in Tornillo basin strata provides strong evidence for a Mogollon highlands or border region basement-cored Laramide uplift origin for the subordinate Yavapai-Mazatzal zircon component (Fig. 8). These core-rim age pairs represent zircons from the nascent Permian–Triassic Cordilleran and Laramide arcs that intruded into Yavapai-Mazatzal country rock in the Arizona-Sonora region (Fig. 2). We suggest that the Paleoproterozoic–Mesoproterozoic age component represents a combination of first-cycle zircons from Yavapai-Mazatzal basement exposed by Laramide uplifts and Mexican Border rift reactivation, inherited zircons from Laramide magmatic rocks, and recycled zircons from the McCoy and Bisbee basins.

Detrital zircons with ages between Ectasian–Tonian (~8%, n = 277/3679) and Tonian–Pennsylvanian (~10%, n = 350/3679, peak at ca. 451 Ma) tend to comprise a minor proportion of Tornillo basin strata age distributions (Figs. 5–7). Various basement terranes currently situated along the eastern margin of North America yield Mesoproterozoic–Paleozoic zircons, including the Grenville and Oaxaquia (ca. 1250–950 Ma) and various Neoproterozoic (ca. 800–500 Ma) and early–middle Paleozoic (ca. 500–300 Ma) arcs and terranes (Fig. 2; e.g., Thomas et al., 2019). Regionally, however, these Mesoproterozoic–Paleozoic ages are ubiquitous in Paleozoic–Mesozoic strata of the western United States and Mexico, including Paleozoic–Mesozoic strata of the Colorado Plateau and Sevier thrust belt (e.g., Dickinson and Gehrels, 2003, 2008, 2009a, 2010; Laskowski et al., 2013; Pujols et al., 2020), Lower Cretaceous Mexican Border rift strata of the McCoy and Bisbee basins (Dickinson et al., 2009), Upper Triassic–Middle Jurassic strata in central Chihuahua (Lawton et al., 2018), Upper Jurassic–Lower Cretaceous strata in Coahuila (Thomas et al., 2019), and Pennsylvanian–Permian strata of the Late Cretaceous–Eocene Marathon Uplift (Thomas et al., 2019). Whereas Neoproterozoic–Paleozoic age components are regionally represented by dispersed, multimodal age distributions, the Ordovician–Silurian mode in Tornillo basin strata is relatively well defined. In fact, Ectasian–Tonian cores with Paleozoic rims (of which 16/21 are Ordovician–Silurian) constitute a prominent bivariate core-rim age cluster in Tornillo basin strata despite being part of subordinate univariate Ectasian–Tonian and Paleozoic components individually (Figs. 5 and 8).

Ectasian–Tonian cores with Paleozoic rims have been found in detrital zircons from strata of the southern Appalachian foreland basin, where the zircons were derived from the Appalachian orogen (Liu et al., 2022). Therefore, since Appalachian-derived zircons were recycled from Mesozoic strata of the Colorado Plateau into the Bisbee basin following Border rifting (Dickinson et al., 2009), Bisbee basin strata could have provided a proximal recycled source of Appalachian-derived Ectasian–Tonian and Paleozoic zircons and core-rim age pairs to the Tornillo basin. Alternatively, Upper Ordovician magmatic rocks that intruded 1250–950 Ma basement rocks of the Oaxaquia terrane exposed in southwest Tamaulipas have yielded zircon ages and core-rim age pairs that closely resemble the Ordovician–Silurian mode and core-rim age pairs of Tornillo detrital zircons (Alemán-Gallardo et al., 2019). Hence, the Ectasian–Tonian and Tonian–Pennsylvanian zircons in Tornillo basin strata may have originated in the Appalachian orogen and been transported west to Mesozoic sedimentary strata of the Colorado Plateau, or these zircons may have originated from the Coahuila terrane, where Oaxaquia basement rocks were intruded by Paleozoic magmatic rocks (Fig. 2). Regardless, rocks from the Colorado Plateau and/or the Coahuila terrane were eroded and shed into Mexican Border rift basins, including the Bisbee basin and the Chihuahua trough, and then subsequently recycled into the Tornillo basin during Laramide-Mexican inversion.

Two possible magmatic sources can be invoked for the Permian–Triassic age component centered at ca. 256 Ma (~8%, n = 304/3679; Fig. 5A), namely, the Permian–Triassic East Mexico arc in eastern Mexico (Torres et al., 1999; Coombs et al., 2020) and the latest Permian–Triassic nascent Cordilleran arc in the southwest United States and Sonora (Figs. 1 and 2; Arvizu et al., 2009; Riggs et al., 2013, 2020). The Permian–Triassic component is dominated by latest Permian ages, which, with the appearance of a Late Triassic mode (ca. 240–200 Ma), support derivation from the nascent Cordilleran arc. This interpretation was further confirmed by core-rim age pairs, as eight of 11 Permian–Triassic rims occurred on Yavapai-Mazatzal (Paleoproterozoic–Mesoproterozoic) cores (Fig. 8). These observations support derivation of Permian–Triassic zircons from the nascent Cordilleran arc in southeastern California, southern Arizona, and northern Sonora, where Permian–Triassic magmatic rocks intruded Yavapai-Mazatzal basement, and effectively preclude derivation from the East Mexican arc, which intruded into Mesoproterozoic Oaxaquian basement (Figs. 1 and 2).

The Jurassic–Early Cretaceous age component (~11%, n = 388/3679) is bimodal with modes at ca. 166 Ma and ca. 151 Ma (Figs. 5–7). The ca. 166 Ma mode likely corresponds to the Jurassic Nazas arc, whereas the ca. 151 Ma mode overlaps with Mexican Border rift–related magmatism. In terms of core-rim relationships, Jurassic cores exhibit both middle Cretaceous Cordilleran and Late Cretaceous Laramide arc rims, indicative of a source in the southwestern U.S.-Mexico border region, where Laramide, Cordilleran, Mexican Border rift–related, and Nazas arc magmatism spatially intersected (Fig. 1). Detrital zircons from Lower and Upper Cretaceous strata in the border region (Barth et al., 2004; Spencer et al., 2011; Clinkscales and Lawton, 2015; Lawton et al., 2020) contain Jurassic U-Pb age modes consistent with derivation from Jurassic magmatic rocks during Mexican Border rifting in the Late Jurassic–Early Cretaceous and Laramide deformation in the Late Cretaceous–Eocene (Dickinson and Lawton, 2001b). Laramide uplifts in southern Arizona, New Mexico, and northern Sonora likely exhumed Jurassic magmatic rocks that provided first-cycle zircons and inverted Bisbee Group synrift strata that provided a prominent source of recycled Jurassic zircons to the Tornillo basin (Fig. 1; Clinkscales and Lawton, 2015; Lawton et al., 2020).

Possible sources for Hauterivian–Coniacian Cordilleran arc zircons (130–87 Ma; ~20%, n = 719/3679) are the Sierra Nevada, the Mojave Desert–Salinia, and the Peninsular Ranges batholiths (Fig. 1). Middle Cretaceous detrital zircon age modes from the Tornillo strata range between 100 and 90 Ma, which coincides with a period of high magmatic flux in the Sierra Nevada (100–85 Ma), the Salinian block (100–85 Ma), and the supervoluminous La Posta suite (99–92 Ma) in the Peninsular Ranges (Figs. 5–7). Alternatively, Cordilleran zircons could also have been recycled locally from ash-fall tuffs in middle Cretaceous sedimentary strata of the southwestern United States and northern Mexico during Laramide tectonism. Zircons from distal tephras in the Mojado, Dakota, and Mancos formations have been documented in the U.S.-Mexico border region (Clinkscales and Lawton, 2015; Lawton et al., 2020). Similarly, the Boquillas Formation in northern Mexico and the Eagle Ford Formation in Texas contain abundant tephra with ages between 100 and 85 Ma that resemble age modes in the Tornillo strata (e.g., Pierce et al., 2016; Deluca, 2016). Considering the large proportion of Cenomanian–Coniacian ages in the Tornillo detrital zircon data, we interpret them to be primarily derived from the Cordilleran arc batholiths, although the degree of contribution of recycled ash-fall zircons remains unresolved. Differentiating between the different Cordilleran batholiths that were tapped by the Tornillo drainage is difficult, but the occurrence of Jurassic cores with middle Cretaceous rims is most consistent with source areas where the eastern zone of the Cordilleran arc overlapped with the Jurassic Cordilleran-Nazas arc (Figs. 1 and 8). Hence, the southern Sierra Nevada or the Mojave Desert–Salinia batholiths represent the more likely sources for these Cordilleran zircons (Fig. 1). However, this does not preclude source contributions from other Cordilleran batholiths, such as the northern Peninsular Ranges.

Another likely source for Cordilleran zircons in the Tornillo basin is recycling of Upper Cretaceous strata in the retroarc of the Sierra Nevada–Mojave Desert–Salinia batholith complex. Exposures of Upper Cretaceous retroarc strata, for example, the McCoy Mountains Formation in western Arizona, exhibit a prominent 110–80 Ma age component and otherwise highly variable detrital zircon age distribution. This Cretaceous age component is typically accompanied by proportionally subequal-to-larger Jurassic and/or Yavapai-Mazatzal components and reflects a combination of zircon sources from the Cordilleran arc and proximal exposures of Jurassic magmatic rocks and Yavapai-Mazatzal basement unroofed during Laramide deformation (Barth et al., 2004; Spencer et al., 2011). Although containing similar age components, the age proportions in Tornillo detrital zircon ages are distinctly different, exhibiting a dominant 100–87 Ma (along with the 87–52 Ma) age component with subordinate Jurassic and Yavapai-Mazatzal components (Figs. 5–7). Although it is possible that Upper Cretaceous retroarc strata contributed recycled zircons to the Tornillo basin, recycling of these strata alone would have been insufficient to produce the proportions of Cordilleran zircons observed in the Tornillo basin. Direct first-cycle bedrock-derived input from the Sierra Nevada–Mojave Desert–Salinia batholith complex to the Tornillo basin is likely required to explain the observed dominant 100–87 Ma and more subdued Jurassic and Proterozoic age components.

Possible sources of Coniacian–Ypresian Laramide arc (87–52 Ma) detrital zircons (~30%, n = 1118/3679) include the Colorado Mineral Belt, the Laramide porphyry copper province of southern Arizona, southwest New Mexico, and northern Sonora, and the Laramide Sonoran batholith (Figs. 1 and 2). Whereas proximity to the Tornillo basin alone would favor sourcing these Laramide zircons from the porphyry copper province, this interpretation is also supported by two prominent zircon core-rim age-pair relationships. First, Paleoproterozoic–Mesoproterozoic cores with Late Cretaceous–Paleocene rims represent instances where Laramide magmas intruded Yavapai-Mazatzal basement (Fig. 8). However, Laramide magmatism intruded into Yavapai-Mazatzal basement in the Colorado Mineral Belt, the porphyry copper province, and the northern part of the Sonoran batholith (Fig. 2; McDowell et al., 2001; Valencia et al., 2005; González-León et al., 2011; Amato et al., 2017; Seedorff et al., 2019). A second prominent cluster of core-rim age pairs of Jurassic cores and Late Cretaceous rims (Fig. 8) appears to correspond to Laramide magmas of the porphyry copper province that intruded magmatic rocks of the Jurassic Cordilleran-Nazas arc in southern Arizona, southern New Mexico, and northern Sonora (Fig. 1). This relationship is supported by zircon U-Pb data from porphyry copper province magmatic rocks in New Mexico and northern Sonora characterized by Jurassic inheritance (Valencia et al., 2005; Amato et al., 2017). Therefore, the presence of Jurassic–Late Cretaceous core-rim age pairs in the Tornillo detrital zircons demonstrates input of Laramide zircons from the porphyry copper province and, possibly, the Sonoran batholith to the Tornillo basin. Hence, abundant volcanic rock fragments and fresh plagioclase feldspar documented in Tornillo basin sandstones (Lehman, 1991) were likely derived from volcanic and plutonic rocks of the porphyry copper province.

Overall, the detrital zircon U-Pb age distributions of the Tornillo basin strata display remarkable temporal stability from the upper Aguja Formation through the lower Canoe Formation (Figs. 5B, 6, and 7). There are no large shifts in age group proportions or sudden appearances or disappearances of components that would indicate drastic provenance changes. However, there are subtle temporal trends already partially alluded to in the MDA section. For the Paleocene–Eocene, the Laramide magmatic component (87–52 Ma) comprises a diverse range of age modes that can generally be lumped into two principal groups: a dominant 87–70 Ma subcomponent and a subordinate 70–52 Ma subcomponent (Figs. 5–7). In particular, the data from the Paleocene–Eocene Black Peaks, Hannold Hill, and Canoe Formations reproducibly show multimodality in the Laramide magmatic component, comprising multiple major age modes between 87 and 70 Ma coupled with smaller age modes between 70 and 52 Ma (Figs. 6 and 7). In contrast, the Laramide arc detrital zircon component in the Upper Cretaceous Aguja and Javelina Formations tends to be unimodal and systematically tracks the depositional age from 84 to 69 Ma (Figs. 6A and 7). This shift from unimodal Laramide arc age groups in the Javelina Formation to more diverse, multimodal Laramide arc age components in the Black Peaks Formation is accompanied by an ~10% overall decrease (from ~42% down to ~32%) in the Laramide arc age component (Fig. 5B). Both trends are broadly consistent with Laramide magmatism becoming less voluminous and migrating eastward from the western Mojave region to eastern Sonora, southern Arizona, and southwestern New Mexico from ca. 85 to 60 Ma. Late Cretaceous unimodal Laramide arc age modes represent the dominant zircon contribution from voluminous Laramide magmatic centers that provided a high abundance of syndepositional Laramide zircons to the Tornillo basin. As Laramide magmatism waned into the Paleocene, contributing a lower proportion of syndepositional zircons, older Laramide U-Pb age modes from uplifted and eroding older Laramide magmatic centers as well as older U-Pb age components from basement uplifts significantly increased again.

Compositional data from Tornillo basin strata indicate that there was an apparent shift in provenance around the Paleocene-Eocene boundary based on a change in sandstone compositions from those rich in volcanic rock fragments, plagioclase feldspar, and quartz in the Javelina and Black Peaks Formations to compositions dominated by carbonate rock fragments, chert, quartz, and orthoclase and plagioclase feldspar in the Hannold Hill and Canoe Formations (Figs. 14A and 14B; Lehman, 1991). Lehman (1991) attributed this shift to a reduction in sediment contribution from more distant volcanic and sedimentary rocks and an increase in nearby sedimentary sources involved in local Laramide-Mexican deformation. This shift, however, is not evident in the detrital zircon data, which show similar detrital zircon age distributions in the Black Peaks, Hannold Hill, and lower Canoe Formations, and which do not indicate a major provenance shift (Figs. 5 and 14B). For example, sample 18TO-01, collected from a sand lens in a Hannold Hill conglomerate dominated by limestone and chert clasts, does not exhibit a change in detrital zircon signal (Fig. 7). The disconnect between compositional and detrital zircon data is likely best explained by high input of carbonate rock fragments and chert from Lower Cretaceous carbonates that were progressively unroofed in Laramide monoclines along the Tornillo basin flanks. Given the presumably low zircon fertility of these carbonates, their detrital contributions would significantly impact the sandstone composition but not be reflected in the detrital zircon record.

In summary, Campanian to Lutetian strata of the Tornillo basin were sourced from a diverse array of volcanic, plutonic, and sedimentary rocks in the southwest U.S.-Mexico border region that encompassed southern New Mexico, Arizona, and northern Sonora, plutonic rocks of the Cordilleran arc in southeast California, and recycled Mexican Border rift sedimentary rocks of northern Sonora and Chihuahua (Figs. 14C and 17). The dominant component of Laramide-aged volcanic and plutonic zircon was eroded from a border region landscape characterized by mixed volcanic and plutonic complexes of the porphyry copper province. The southwest border region was also involved in Laramide deformation and Mexican Border rift inversion, resulting in uplift and unroofing of Jurassic and Permian magmatic rocks, Yavapai-Mazatzal basement, and various Mexican Border rift and Sevier foreland sedimentary units that supplied many of the subordinate zircon components. In total, the complex structural inheritance, numerous basement domains, and sedimentary basins resulted in a diverse source region landscape consisting of volcanic and plutonic complexes and a mix of basement and sedimentary cover units exhumed by Laramide tectonic structures.

The reconstructed >1400-km-long, ~250,000 km² Tornillo drainage basin suggests that complex structural inheritance and Laramide inversion of the Mexican Border rift were primary controls on the axial-trunk river pathways through the border region (Fig. 14C). A dominant Cordilleran arc zircon component indicates that headwaters of the Tornillo river extended as far west as the easternmost reaches of either the Peninsular Ranges, Mojave Desert–Salinia, and/or southern Sierra Nevada batholiths of the Cordilleran arc in southeastern California or at least into proximal retroarc Upper Cretaceous sedimentary units. From there, the Tornillo river and its tributaries threaded their way through a structurally complex landscape of Laramide volcanic complexes and Bisbee rift basin fault blocks and basins that were being reactivated and inverted as Laramide basement–involved uplifts in southern Arizona and New Mexico and northern Sonora. This tectonic landscape contributed the non-Cordilleran arc–related dominant detrital zircon components that make up the detrital zircon ages observed in the strata of the Tornillo basin. Upon entering west Texas, the Tornillo river was flanked to the south by the emerging topography of the inverted Chihuahua trough. From west Texas, the Tornillo river flowed toward the northwest Gulf of Mexico margin of northern Mexico or southern Texas, where it likely contributed sediment to major Late Cretaceous–Eocene coastal depocenters, including Wilcox Group deltas (e.g., Galloway et al., 2011). This large (>10⁵ km²) drainage basin, with its headwaters near the Pacific Ocean, was characterized by a trunk river that flowed along a narrow inversion-flank drainage corridor adjacent to topography created by inverted early Mesozoic Mexican Border rift features, including Laramide basement–cored uplifts in the Bisbee basin and the inverted Chihuahua trough. We posit that inherited Mexican Border rift architecture was the first-order control on the sediment routing to the Tornillo basin, funneling fluvial systems parallel to the former rift axis (Fig. 14C). For this reason, the axial-trunk Tornillo river essentially followed the U.S.-Mexico border corridor for the entirety of its drainage (Fig. 17). The interpretation of a major, axially routed river to the Tornillo basin from at least the Maastrichtian to the earliest Ypresian is supported by the abundance of volcanic rock fragments and plagioclase feldspar in the Javelina and Black Peaks Formations (Lehman, 1991), suggesting that this axial-trunk river was responsible for the bulk of the sediment supply to the Tornillo basin (Figs. 14A and 14B).

This axial-trunk river was likely supplemented by transverse tributaries feeding sediment from the inverted Bisbee basin and the Chihuahua fold belt, the inverted Chihuahua trough, northward (Fig. 14C). The Chihuahua trough is composed of sedimentary strata, including early Mesozoic siliciclastic rocks and thick sequences of Cretaceous carbonates. Despite the low zircon fertility of Chihuahua trough rocks, several subordinate zircon age components are likely attributable to recycling of early Chihuahua trough strata and Bisbee basin strata during Laramide-Mexican inversion of the Mexican Border rift. In particular, Ectasian–Tonian and Tonian–Pennsylvanian zircons were likely recycled from the Bisbee basin and Chihuahua trough, as underscored by core-rim age pairs linking zircons to Mesozoic strata of the Colorado Plateau and/or the Coahuila terrane, both of which were sources to the Mexican Border rift (Figs. 2 and 8). Additionally, Jurassic and Permian zircon age components, originally supplied to the Bisbee basin and Chihuahua trough during Mesozoic rifting, were likely recycled during rift inversion. Importantly, the emerging topography of Laramide uplifts in Arizona, Sonora, and New Mexico and the Chihuahua fold belt also appears to have created a regional drainage divide forming the southern margin of the Tornillo drainage basin, as evidenced by a lack of detrital zircons derived from western Mexico (Fig. 17). In contrast to the Parras–La Popa basin, the Tornillo detrital zircon data lack ca. 132 Ma zircons derived from the Alisitos arc and Guerrero terrane in western Mexico (Fig. 15; e.g., Lawton et al., 2009; Schaaf et al., 2020).

This provenance scenario is supported by available sandstone petrography (Figs. 14A and 14B; Lehman, 1991). Abundant volcanic rock fragments and plagioclase feldspar support an axially dominant river system sourced from distal Laramide magmatic and basement sources with subsidiary transverse sediment input to the Tornillo basin via tributaries from the Maastrichtian to Ypresian. However, a dramatic increase in abundance of carbonate rock fragments and chert in the Hannold Hill and Canoe Formations points to an increase in transverse sediment input from local emerging Laramide monoclines and uplifts during the early Eocene that is not readily evident in the detrital zircon data (Fig. 14B). This is likely due to the low zircon fertility of the transverse source regions, and, thus, the detrital zircon signal likely continued to predominantly reflect the axially derived sediment input component through the Ypresian–Lutetian. Alternatively, the detrital zircons of the Hannold Hill and Canoe Formations could reflect erosional recycling of the Black Peaks Formation in these emerging local monoclines. This is at least partially supported by the similarity in the detrital zircon signatures, evidence for contemporaneous erosion of parts of the upper and lower Black Peaks Formation, now unconformably overlain by the Canoe Formation (Fig. 4; Lehman et al., 2018), and evidence for recycling of the Black Peaks Formation into the Big Yellow sandstone, including recycled silicified logs (Lehman, 1991). Overall, this alternate scenario would imply a local tectonic reorganization across the Black Peaks–Hannold Hill contact and a switch to predominantly transverse sediment input from recycled sources during the early Eocene. Analysis of sand body thicknesses and fluvial paleocurrent data, however, does not suggest a major decrease in discharge or change in paleoflow direction that would point to a major disruption and abandonment of the axial drainage delivering sediment to Tornillo basin from the Paleocene to the Eocene. This implies that there was not a major decrease in the contribution of axial sediment delivery from distal sources during the early Eocene, just the addition of transverse sediment input from carbonate-dominated local Laramide monoclines, uplifts, and thrust sheets (e.g., Hennings, 1994). By the middle Eocene, however, the onset of Trans-Pecos volcanism in the region severely disrupted sedimentation and the Tornillo drainage basin (Lehman, 1991).

Early provenance interpretation based on sandstone petrography favored a scenario in which the Tornillo basin was sourced from distant volcanic highlands in the western United States and Mexico (Fig. 14A; Lehman, 1991). However, more recent reconstructions commonly interpreted the Campanian–Lutetian Tornillo basin as a locally sourced basin predominantly draining the Chihuahua fold belt and west Texas or incorporated it as part of a northern Mexican drainage basin (e.g., Galloway et al., 2011; Blum and Pecha, 2014; Sharman et al., 2017; Blum et al., 2017). Most recently, Colliver (2017) interpreted detrital zircon data to indicate that the Tornillo basin was likely sourced from more distant magmatic sources, in agreement with the compositional data (Lehman, 1991), and implied a source region encompassing west Texas, southern New Mexico, southern Arizona, and northwest Mexico (Bataille et al., 2019). Depth-profiled detrital zircon data presented here confirm that the Campanian–Lutetian Tornillo basin was not locally sourced, but rather received sediment from a large drainage basin that extended into southern New Mexico, Arizona, northern Sonora, and possibly southeastern California (Fig. 17). Aside from simple consideration of the U-Pb age components, this source region for the Tornillo basin is strongly supported by detrital zircon core-rim age relationships that match zircon inheritance documented in the border region (Fig. 8).

Part of the challenge of expanding the Campanian–Lutetian Tornillo drainage basin to encompass a large part of the border region is that it creates source area competition with four other published models of Laramide and Mexican river drainages (Fig. 17). The Tornillo drainage basin was bordered to the northwest by north-northeast–flowing rivers of the Colorado Plateau, including the California river, which flowed to the Uinta basin (Davis et al., 2010; Dickinson et al., 2012; Foreman and Rasmussen, 2016), to the northeast by east-southeast–flowing rivers of the Four Corners region (Donahue, 2016; Pecha et al., 2018; Smith et al., 2020), to the southwest by west-flowing rivers that navigated the Salinian embayment and exported material from northwest Sonora to the Transverse Ranges forearc (Jacobson et al., 2011; Sharman et al., 2015), and to the south by the Parras–La Popa drainage basin that drained much of northern Mexico (Lawton et al., 2009). Therefore, the Campanian–Lutetian Laramide hinterland of southeastern California, western Arizona, and northern Sonora appears to have served as a common headwater region for the north-flowing California river, the east-flowing Tornillo river, and southwest-flowing rivers of the Salinian embayment—all tapping the same source terranes (Figs. 15–17). This inference is supported by the similarities in the detrital zircon age distributions from the Transverse Ranges forearc basin and the Tornillo and Uinta basins (Figs. 15 and 16). There is a marked similarity between the Tornillo and Transverse Ranges detrital zircon age distributions, such as the shared codominant middle Cretaceous (Cordilleran arc) and Late Cretaceous–Paleocene (Laramide arc) age components, which suggest that the Tornillo and Salinian embayment rivers accessed the same source areas in southeastern California, western Arizona, and northern Sonora. Furthermore, both the Tornillo and California rivers may have accessed the same source region in southeastern California based on the dominant presence of Cordilleran arc–derived zircons younger than 300 Ma in the Tornillo and Uinta basins. Full resolution of the Laramide drainage basin boundaries between the various paleorivers that tapped the Laramide hinterland of California, Arizona, and Sonora will require further provenance testing and integration of preexisting data sets from the southwestern United States and Mexico.

Several provenance studies and synthesis papers have interpreted north-to-northeast drainage from central Arizona to the Four Corners region during Late Cretaceous–Eocene time (e.g., Galloway et al., 2011; Donahue, 2016; Sharman et al., 2017; Blum et al., 2017; Pecha et al., 2018). For example, Pecha et al. (2018) interpreted that Laramide-aged zircons in Campanian strata of the San Juan basin were transported from the porphyry copper province in southern Arizona and southwest New Mexico along north-flowing rivers. Beginning in the Maastrichtian, however, there was a major provenance shift as rivers in the San Juan basin switched to predominantly south-southeast–directed paleoflow, and Laramide-aged zircons were instead derived from the Colorado Mineral Belt to the north (Figs. 15–17; Pecha et al., 2018; Cather et al., 2019). This late Campanian switch in San Juan basin provenance coincided with deposition of the Aguja Formation in the Tornillo basin, which has clear porphyry copper province provenance that continued through the Maastrichtian–Ypresian Javelina and Black Peaks Formations. Therefore, we suggest that rivers of the Campanian–Lutetian Tornillo drainage basin likely captured much or all of the sediment derived from the porphyry copper province from the north-flowing rivers that fed the San Juan basin during the Campanian. Despite the shift in provenance in the San Juan basin, regional provenance syntheses have tended to maintain a hypothetical river that flowed north and east from the porphyry copper province into west-central New Mexico throughout the Campanian–Lutetian (e.g., Galloway et al., 2011; Sharman et al., 2017; Blum et al., 2017; Pecha et al., 2018). However, strata of the Baca basin in west-central New Mexico are middle Eocene in age (e.g., Cather, 2004), making the existence of a northeast-flowing river through central Arizona and New Mexico prior to the Lutetian difficult to evaluate. Regardless, by the onset of Baca basin deposition, there is evidence for an east-flowing river system through west-central New Mexico (Potochnik, 1989; Elston and Young, 1991; Lawton, 2019). The development of this middle Eocene river through New Mexico coincided with the apparent deflection of the Tornillo inversion-flank river system away from the Tornillo basin by Trans-Pecos volcanism.

The Tornillo river and its tributaries were bordered to the south by the inverted Chihuahua trough, which separated the Tornillo river drainage from the Parras–La Popa drainage basin, a continental-scale routing system that ran axially along the Mexican orogen thrust front in northern Mexico and terminated into the Parras and La Popa basins of northeast Mexico (Fig. 17). Although Tornillo and Parras–La Popa detrital zircons exhibit similar overall age distributions, the Parras–La Popa detrital zircon data contain a distinct ca. 132 Ma zircon component attributed to the western Mexican Alisitos and Guerrero arcs that is conspicuously absent in the Tornillo basin (Fig. 15). This indicates that, despite the apparent similarities in their detrital zircon age distributions, the Tornillo river drainage did not tap into western Mexico and was a distinct, separate system from the Parras–La Popa drainage basin. This indicates that the Chihuahua fold belt formed an effective drainage divide between the two drainage basins from the Late Cretaceous to Eocene.

After passing through the Tornillo basin, the Tornillo river presumably continued flowing east-southeast toward the northwestern Gulf of Mexico. The location where the Tornillo river intersected with the Gulf Coast remains uncertain, but one possible location includes deltaic–coastal marine depocenters of the paleo–Rio Grande embayment in south Texas, where the Tornillo river would have contributed to Wilcox Group sedimentation (Mackey et al., 2012; Blum et al., 2017). Alternatively, the Tornillo river may have fed into the northern part of the La Popa basin in northeastern Mexico (separate from the Parras–La Popa drainage in northern Mexico), in which case, sediment delivered from the Tornillo river would have been sequestered inboard of the Gulf Coast margin and would not have contributed to Wilcox Group sedimentation (Galloway et al., 2011). Further detrital zircon provenance analysis from Paleogene Gulf Coast strata in northeastern Mexico and south Texas would be able to test the possible connection between the Tornillo basin and the paleo–Rio Grande embayment and may have cascading implications for broader North American paleodrainage reconstructions to the northern Gulf Coast. Regardless, Tornillo Group provenance analysis and delineation of the Tornillo river through the U.S.-Mexico border region bring an underrepresented segment of Late Cretaceous–Eocene sediment routing from the North American hinterland to the Gulf of Mexico into sharper focus.

The Tornillo river was part of a protracted history of sediment-routing systems that flowed parallel to structural trends generated by the Mexican Border rift in the U.S.-Mexico border region. Following the cessation of rifting in the Late Jurassic, persistent east-directed sediment routing in southern Arizona and northern Sonora along the paleo-axis of the Mexican Border rift was separated by the Mogollon highlands from a northeast-directed routing system into the Sevier foreland (e.g., Lawton et al., 2003, 2014; Szwarc et al., 2015), until both systems were interrupted by marine transgressions during the Albian–Cenomanian (Lawton et al., 2020). Subsequent inversion and renewed Laramide-Mexican uplift of the Bisbee basin and Chihuahua trough during the Late Cretaceous–Eocene resulted in significant topographic relief in the Arizona-Sonora and New Mexico–Texas–Chihuahua border regions that followed the structural trend of the Mexican Border rift (e.g., Hennings, 1994; Haenggi, 2002; Clinkscales and Lawton, 2015, 2017; Chapman et al., 2020). The axial-trunk river of the Tornillo sediment-routing system likely occupied a narrow drainage corridor positioned along the northern flank of this belt of inversion structures, including Laramide basement–cored uplifts in the Bisbee basin and the Chihuahua fold belt (Fig. 14C). Longitudinal rivers within this inversion-flank drainage corridor were fed via transverse tributaries from Laramide-Mexican inversion structures along the entire length of the corridor to the Tornillo basin, highlighting the continued importance of inherited Mexican Border rift structures in Tornillo drainage. On the north side of the corridor, a drainage divide near the Mogollon Rim apparently separated the Tornillo river from north-flowing rivers such as the California river (Figs. 14C and 17). Although paleodrainage in southern Arizona was briefly north-directed in the Late Cretaceous (e.g., Pecha et al., 2018), the establishment of the Tornillo river along an inversion-flank corridor adjacent to inherited rift features suggests that Mexican Border rift architecture remained a first-order control on sediment routing in the border region from the Jurassic through at least the Eocene.

Despite being occupied by the Tornillo river, the inversion-flank drainage corridor was rendered topographically irregular by Laramide-Mexican uplifts along its reach, including Laramide uplifts in southern Arizona–New Mexico and monoclines within the Tornillo basin (Fig. 14C; e.g., Lehman et al., 2018; Favorito and Seedorff, 2021). Therefore, the Tornillo river and its tributaries must have navigated through complex, structurally controlled Laramide topography both within the drainage corridor and the adjacent inversion belt. The recent history of river drainages of the High Atlas of Morocco (>10⁵ km² in total) serves as a potential partial analog for the paleogeography of the Tornillo river. During early rift inversion in the High Atlas, river drainages were structurally controlled, with tributaries taking tortuous longitudinal and transverse paths within the orogen before making their way to narrow, axial systems on the flanks of the mountain belt (Babault et al., 2012). Even in the modern Atlas, in which transverse drainages are more direct and largely ignore structure within the High Atlas, rivers still flow to relatively narrow, predominantly axial systems along a highly fragmented inversion flank. The strong tectonic control on Miocene–recent drainage partitioning in the Andean broken foreland (~10⁵ km²) provides another analog for the way in which Laramide rivers in the Tornillo drainage likely responded to emerging, structurally partitioned topography by deflecting around emerging structures, flowing longitudinal to established ranges, and passing through structural gaps within the inversion belt (Seagren et al., 2022).

It has been hypothesized that the locations of large rivers are, to a first order, controlled by tectonic setting (Potter, 1978). Much of the upstream course of the modern Rio Grande closely follows the axial valley of an active rift system. Orogenic systems and their flexural foreland basins exert a first-order control on big river drainage patterns (e.g., Garcia-Castellanos, 2002). Modern megarivers like the Amazon, Mississippi, and Niger Rivers traverse entire cratons from distant orogenic highlands, but their downstream reaches are located in narrow remnant troughs of aulacogens and old fault/graben systems (e.g., Potter, 1978; Frostick and Reid, 1989). In this study, we posit that the Late Cretaceous–Eocene Tornillo river is an example of large paleoriver with a path that was controlled by its longitudinal position along the flank of a major tectonic inversion belt. We estimate that the Tornillo river was >1400 km in length and had a paleodrainage basin area of ~250,000 km² upstream from the Tornillo basin. This estimate does not include the full length or paleodrainage area of the Tornillo river, as the river likely continued to the Cretaceous–Paleogene shoreline in south Texas and/or northeastern Mexico, but it nevertheless suggests the Tornillo river was comparable in drainage basin size to modern big rivers (Fig. 18). Notably, the Tornillo drainage basin area was small relative to its length when compared to most modern big rivers. Narrow drainage basin aspect ratios are common for modern big rivers that exhibit a high degree of structural control on their drainage networks, such as the Nile, Mekong, and Salween Rivers (Fig. 18). Therefore, the narrow drainage basin reconstructed for the Tornillo river further reflects the primary control of Laramide uplift and Mexican Border rift inversion on Tornillo river drainage.

We demonstrated that inherited rift architecture exerted a first-order control on successor sediment-routing systems. The Tornillo river provides a case example demonstrating the first-order impact of inherited continental rift architecture on the establishment and long-term dispersal pathways of large drainage basins. Furthermore, the reconstructed route of the Tornillo river suggests that it did not simply flow axially along the flank of an inversion belt, but rather navigated a complex, structurally partitioned landscape. From its upstream reaches, the Tornillo river likely encountered a narrow and fragmented drainage corridor and must have skirted or transected the emerging topography of volcanic complexes, inverted basement-involved fault blocks, and growing folds. This complex drainage history illustrates the spatial and temporal heterogeneities in axial and transverse sediment routing that can occur in continental sediment-routing systems during rift inversion.

New detrital zircon U-Pb data from Campanian–Lutetian fluvial strata of the Tornillo basin in west Texas comprise dominant 130–87 Ma and 87–52 Ma age components representing the middle Cretaceous Cordilleran arc and Laramide magmatism, respectively, with subordinate Paleoproterozoic–Mesoproterozoic (Yavapai-Mazatzal and Granite-Rhyolite Provinces), Ectasian–Tonian, Tonian–Pennsylvanian, Permian–Triassic, and Jurassic–Early Cretaceous components. Detrital zircon MDAs bracket the stratigraphic location of the Cretaceous-Paleogene boundary to the uppermost Javelina Formation at Dawson Creek and to the lowermost Black Peaks Formation at Tornillo Flat and support placement of the Paleocene–Eocene thermal maximum horizon at or just below the south wall sandstone at Tornillo Flat. Detrital zircon MDAs also point to the presence of a basinwide ~3–4 m.y. early Danian hiatus within the lowermost Black Peaks Formation, possibly indicating a precursor phase to the monocline uplifts that disrupted the basin in the Eocene.

Depth-profiled core-rim age pairs comprise several prominent clusters, including Paleoproterozoic–Mesoproterozoic and Jurassic cores with Cretaceous–Paleogene rims and Ectasian–Tonian cores with Paleozoic rims. Paleoproterozoic–Mesoproterozoic and Jurassic cores with Cretaceous–Paleogene rims represent zircons from Yavapai-Mazatzal basement and Jurassic Cordilleran-Nazas arc magmatic rocks that were subsequently inherited by Laramide magmatic rocks of the porphyry copper province in southern Arizona, New Mexico, and northern Sonora. Zircons that have Ectasian–Tonian cores with Paleozoic rims were originally derived from the Appalachian orogen, and then were transported west to Mesozoic strata of the Colorado Plateau, or they were derived from the Coahuila terrane of north-central Mexico. In either case, zircons with Ectasian–Tonian cores and Paleozoic rims were shed into the Mexican Border rift and subsequently recycled out of Jurassic and Lower Cretaceous siliciclastic strata of the inverted Bisbee basin and Chihuahua trough.

These results demonstrate that strata of the Tornillo basin were predominantly sourced from first-cycle volcanic and plutonic complexes of the Laramide porphyry copper province in southern New Mexico, Arizona, and northern Sonora; first-cycle plutonic rocks of the Cordilleran arc in southeast California; Yavapai-Mazatzal basement; and recycled Mexican Border rift sedimentary rocks of Arizona, New Mexico, Sonora, and Chihuahua. This large drainage basin (>1400 km long, ~250,000 km²), with headwaters in southern California, Arizona, and northern Sonora, routed sediment to the Tornillo basin via an axial-trunk river that occupied a position on the northern flank of emerging Laramide topography, which included porphyry copper province volcanic and plutonic complexes and basement-cored uplifts of the inverted Bisbee basin. Upon entering west Texas, this river ran along the flank of the inverted Chihuahua fold belt, from which it received detrital input from recycled Mexican Border rift sedimentary rocks of the inverted Chihuahua fold belt via transverse tributaries. A lack of detrital zircons from western Mexico suggests that the emerging Laramide topography in Arizona, Sonora, and New Mexico and the Chihuahua fold belt formed a regional drainage divide along the southern margin of the Tornillo drainage basin. This drainage reconstruction indicates that the Tornillo river flowed adjacent to topography formed by the inverted early Mesozoic Mexican Border rift basins. The river was funneled along a narrow inversion-flank drainage corridor parallel to the paleorift axis, demonstrating that the inherited Mexican Border rift architecture exerted a first-order control on sediment routing to the Tornillo basin.

We thank the University of Texas–Austin Equinor Fellows Program and the sponsors of the University of Texas–Austin Bureau of Economic Geology Quantitative Clastics Laboratory (http://www.beg.utexas.edu/qcl) for financial support. Field work and laboratory analyses for this study were generously funded by a Statoil/Equinor grant to Daniel Stockli, the UTChron laboratory, and the Chevron (Gulf) Centennial Professor Endowment to Stockli. We would like to thank Don Corrick and the National Park Service for assistance in permitting and field access. Lisa Stockli provided training and assistance to Cullen Kortyna with laser ablation–inductively coupled plasma–mass spectrometry (LA-ICP-MS) data collection. Eugene Szymanski assisted in field work and sample collection. Cullen Kortyna also thanks Doug Barber, Megan Flansburg, Margo Odlum, Eirini Poulaki, Kelly Thomson, and the Stockli Research Group for training and assistance in mineral separation and detrital zircon U-Pb data reduction, visualization, and interpretation. Last, we thank Trystan Herriott and an anonymous reviewer for their constructive reviews, which greatly improved this manuscript.