Provenance of Pennsylvanian–Permian sedimentary rocks associated with the Ancestral Rocky Mountains orogeny in southwestern Laurentia: Implications for continental-scale Laurentian sediment transport systems

Over the last several decades, increasing attention has been focused on the importance of treating sediment production, transport, and storage processes within an integrated source-to-sink framework (e.g., Allen, 2008; Tinker et al., 2008; Covault et al., 2011; Romans et al., 2016; Wang et al., 2019). Vital to this approach is the determination of both the spatial and temporal extent of the source-to-sink system, or alternately stated, it is imperative to understand the degree to which the source is separated from the sink, both by physical distance and by the time elapsed between erosion and deposition. This can be relatively straightforward for small, short-lived basins and/or those with limited catchment areas (e.g., Sømme et al., 2011; Cheng et al., 2016; Sickmann et al., 2016). However, there is increasing documentation of the importance of continent/supercontinent-scale sediment transport systems in the geologic record (Gehrels et al., 2011; Blum and Pecha, 2014; Benyon et al., 2014; Blum et al., 2017) and in modern systems (Garzanti et al., 2014). Of particular importance in mixed terrestrial/marine systems are the roles of littoral and aeolian processes (Garzanti et al., 2014, 2018), and “axial” sediment transport (perpendicular to the steepest surface path across a source-to-sink system) in controlling sediment transport patterns and sediment budgets (Garzanti et al., 2014; Romans, et al., 2016). Recent work has demonstrated that transport by longshore current can result in mixing of provenance signals at long spatial and temporal wavelengths (e.g., Liu et al., 2007; Garzanti et al., 2014; Sickmann et al., 2016), but these processes are not as well considered in deep-time provenance records.

The late Paleozoic Ancestral Rocky Mountains (ARM) system presents an exceptionally well-documented opportunity to explore these concepts at a continental scale because the ARM preserves mixed marine/terrestrial depositional systems that were active across much of southern Laurentia over a period of ∼50 m.y. (Fig. 1; McKee and McKee, 1967; Mallory, 1972; Rascoe and Baars, 1972; McKee and Crosby, 1975). Moreover, basement ages of rocks within the ARM system from which detritus would have been locally derived are well dated (e.g., Sims et al., 2005; Karlstrom et al., 2004; Whitmeyer and Karlstrom, 2007), allowing relatively high confidence in provenance interpretations. Here, we present new U-Pb detrital zircon data from 72 samples throughout the ARM system and integrate these data with published data from an additional 241 samples to track the evolution of sediment provenance and transport pathways and their relationship to basin formation and orogenesis.

The Ancestral Rocky Mountains system consists of a series of basement-cored uplifts bounded by reverse faults that were active during the Pennsylvanian through early Permian periods. Movement along these faults resulted in the formation of asymmetric flexural basins proximal to uplifts (e.g., Barbeau, 2003); however, ARM-coeval subsidence and sediment accumulation is not confined to discrete basins but extends across much of southwestern Laurentia. The Ancestral Rocky Mountains were first recognized in central Colorado and Utah (e.g., Ver Weibe, 1930), but the system also extends across most of southwestern Laurentia including southern New Mexico, Texas, and Oklahoma (Kluth and Coney, 1981; Fig. 1). The ARM system was situated at or near equatorial latitudes during Pennsylvanian and early Permian time, during which Laurentia was rotated approximately 45 degrees clockwise from its current orientation (Scotese and McKerrow, 1990; Domeier and Torsvik, 2014; Cao et al., 2017; note that the modern reference frame and cardinal directions are used throughout this study).

The ARM system was bounded on three sides by convergent tectonic plate boundaries (Leary et al., 2017; Fig. 1). To the south and east of the ARM system, convergence between Gondwanaland and Laurentia created a series of paired fold-and-thrust belts and foreland basins including the Appalachian orogen and basin, the Ouachita belt, and the Black Warrior, Arkoma, Fort Worth, Val Verde, and Marfa Basins, in which subsidence propagated to the southwest over time (Graham et al., 1975; Ross, 1986; Dickinson and Lawton, 2003). By Pennsylvanian time, final Pangean suturing was complete in the southern Appalachians, and contractile deformation and foreland basin formation migrated to the Ouachita–Marathon fold-thrust belt and adjacent Delaware and Midland Basins of west Texas (Whiting and Thomas, 1994; Thomas, 1988). To the west, convergence along westward-dipping subduction zones resulted in a series of arc collision and accretion events that affected the paleo-Nevada continental margin (Dickinson et al., 2000). These include the Devonian–Mississippian Antler event that emplaced the Roberts Mountain Allochthon and late Permian–Triassic Sonoma collision that emplaced the Golconda Allochthon, both of which contributed to flexural subsidence and foreland formation (Speed and Sleep, 1982; Burchfiel and Royden, 1991; Dickinson et al., 2000). Although tectonic activity along this margin may have been concentrated during accretion events, the entire margin was subject to ongoing convergent deformation events from at least Mississippian through Permian time (Erickson and Marsh, 1974; Trexler et al., 1991, 2004; Cashman et al., 2011). The southwestern margin of the ARM system is poorly exposed in present-day northern Mexico and southern Arizona and has been strongly overprinted by several major episodes of tectonism (e.g., McQuarrie and Wernicke, 2005; Clinkscales and Lawton, 2018). Paleotectonic reconstructions (Dickinson and Lawton, 2001), recent late Paleozoic global plate models (Domeier and Torsvik, 2014), and kinematic analysis of the ARM deformation suggest that this margin evolved from oblique convergence to subduction during late Paleozoic time (Leary et al., 2017; Lawton et al., 2017).

The tectonic driver(s) of Ancestral Rocky Mountain basement deformation remain a topic of debate, and proposed models include collision along the Ouachita–Marathon belt (Kluth and Coney, 1981), left-lateral shear along major transcontinental faults (Budnik, 1986), flat-slab subduction along the southwestern continental margin (Ye et al., 1996), reactivation of pre-existing rift structures (Marshak et al., 2000), wrenching associated with the zippering closure of the Iapetus Ocean (Dickinson and Lawton, 2003), dynamic uplift/subsidence due to stress-field interaction with underplated Proterozoic mafic rocks (Soreghan et al., 2012), and oblique convergence along the southwestern Laurentian margin (Leary et al., 2017; Lawton et al., 2017). A complete discussion of the strengths and limitations of these models is beyond the scope of this paper; however, we emphasize several salient aspects of this system that are relevant to the current study: (1) The ARM system developed within an area underlain by Proterozoic Yavapai and Mazatzal basement rocks that were accreted to the Archean Wyoming Craton between 1.7 Ga and 1.6 Ga (Karlstrom and Bowring, 1988; Bickford et al., 1989; Reed et al., 1993) and were later intruded by ca. 1.4 Ga and ca.1.1 Ga alkaline “A-type” granitic intrusions (Anderson and Bender, 1989; Bickford et al., 2015). (2) The core of the ARM system in Colorado straddles the trace of the Transcontinental Arch, an enigmatic line of mid-cratonic positive relief that seems to have reduced accommodation and acted as a barrier to sediment transport for most of the Phanerozoic (Weimer, 1978; Garfield et al., 1988; Carlson, 1999; Amato and Mack, 2012). (3) ARM uplifts and resulting flexural basins strike generally to the northwest-southeast (modern reference frame) (Fig. 1; Ye et al., 1996), with preserved kinematic data suggesting northeast-southwest–oriented shortening (Brewer et al., 1983; Granath, 1989; Shumaker, 1992; Hoy and Ridgway, 2002; Cather et al., 2006; Cox, 2009). Uplift and basin orientations in several cases follow the trends of Proterozoic rifts (Karlstrom and Humphreys, 1998; Marshak et al., 2000), and multiple tectonic models have proposed that associated faults controlled the localization of ARM deformation (Marshak and Paulsen, 1996; Marshak et al., 2000; Soreghan et al., 2012; Leary et al., 2017). (4) ARM deformation did not substantially affect Archean crust of the Wyoming Craton (Karlstrom and Humphreys, 1998), although whether that can be attributed to crustal rheology or distance from active continental margins remains unknown. (5) The ARM system was amagmatic, which prohibits precise radioisotopic geochronology; however, marine strata within the ARM system are well dated biostratigraphically and can be precisely correlated to strata in the Donets Basin and the Ural Mountains of eastern Europe (Schmitz and Davydov, 2012). Available geochronology suggests that ARM-related subsidence (and inferred uplift) began more or less synchronously across most of the region during the Early Pennsylvanian (Leary et al., 2017).

Sedimentary facies within the ARM basins evolved from regional exposure and terrestrial facies during the Late Mississippian (e.g., Huddle and Dobrovolny, 1945; Evans and Reed, 2007) into locally thick accumulations of marine carbonate, shale, and evaporite during the Pennsylvanian (McKee and Crosby, 1975), then transitioned to terrestrial deposition during the early Permian (McKee and McKee, 1967), approximately consistent with the end of ARM deformation.

Samples analyzed in this study were collected from ARM-related sedimentary rocks from the core of the ARM system–the Denver, Eagle, and Paradox Basins–as well as peripheral areas in New Mexico, Arizona, and southeast California. Regional sample locations and stratigraphy are described in Figures 1 and 2; detailed sample locations and detrital zircon data are presented in Figures 3–8 and Table 1. Below, we briefly describe the sampled strata of each area. We refer readers to GSA Data Repository File DR1¹ for precise locations, measured sections, and interpretations of lithofacies, age, and depositional environments. Although previously published data from Grand Canyon samples (Gehrels et al., 2011) are reviewed in the results, the stratigraphy and sampling within that study is not described here. Most locations sampled in this study do not require substantial palinspastic restoration to account for Sevier, Rio Grande Rift, and/or Basin and Range deformation (e.g., DeCelles, 2004; McQuarrie and Wernicke, 2005; McQuarrie and Oskin, 2010) (Fig. 1). Note, however, that locations affected by these episodes of deformation are not palinspastically restored in maps within this study.

This section consists of Mississippian through middle Permian rocks (Figs. 2 and 4, panel 1). The Mississippian Redwall Limestone makes up the base of the section and is overlain by the predominantly siliciclastic Supai Group, Hermit Formation, and Schnebly Hill Formation (Blakey, 1980, 1990; McKee, 1982). We collected samples from very fine-grained to fine-grained sandstone of the Manakacha, Esplanade, Hermit and Schnebly Hill Formations.

In this section, the base of the Pennsylvanian–Permian section was deposited above a disconformable karst surface developed on the Mississippian Redwall Limestone (Huddle and Dobrovolny, 1945; McKee and Gutschick, 1969) (Figs. 2 and 4, panel 3). A wide range of stratigraphic schemes have been proposed for this region, which is transitional between the siliciclastic Pennsylvanian–Permian section to the northwest and the carbonate dominated Pennsylvanian section to the southeast (e.g., Huddle and Dobrovolny, 1945; McKee and Gutschick, 1969; Brew, 1970; Ross, 1973; Blakey and Knepp, 1989) (Fig. 2). Here, we adopt the stratigraphy of Blakey and Knepp (1989), in which the Pennsylvanian Naco Group (equivalent to the Horquilla Limestone of Southeast Arizona; Ross, 1973; Brew and Beus, 1976; Blakey and Knepp, 1989) is overlain by the Permian Hermit Formation, Schnebly Hill Formation, and the Coconino Sandstone. We collected three samples from fine-grained and medium-grained sandstones and a granule conglomerate within the Middle Pennsylvanian Naco Group; a medium-grained sandstone from the lower Permian Hermit Formation; and two samples from fine and very fine-grained sandstone of the middle Permian Schnebly Hill Formation.

The Paleozoic succession in southwestern Arizona lithologically resembles the Paleozoic section preserved within the Grand Canyon and Mogollon Rim and has been inferred to be stratigraphically and lithologically equivalent (Richard et al., 1993; Blakey, 2008) (Figs. 2 and 4, panel 2). The Plomosa Mountains are currently ∼200 km southwest of the mostly undeformed western Mogollon Rim section (above), but palinspastic restoration of Basin and Range extension (not shown on figures within this study) locates the original depositional position of these strata within ∼100 km of the edge of the Colorado Plateau (McQuarrie and Wernicke, 2005).

We collected one sample each from medium-grained sandstone in the basal Supai Group, coarse-grained sandstone of the lower Permian Hermit Formation, and fine-grained sandstone of the Coconino Sandstone exposed in the southern Plomosa Mountains southeast of the town of Quartzsite. These rocks are exposed in a steeply dipping, fault-bounded panel within a zone that has been subjected to multiple episodes of intense compressional and extensional deformation (Miller, 1970; Richard, 1992, 1993; Richard et al., 1993; McQuarrie and Wernicke, 2005; McQuarrie and Oskin, 2010; Spencer et al., 2015). Most late Paleozoic sedimentary rocks in this region have been intensely deformed and metamorphosed. The section sampled for this study is a rare exception and is preserved largely unmetamorphosed.

The section is fault-bounded and consists of ca. 375 m of strata mapped as Pennsylvanian Horquilla Limestone deposited conformably atop Mississippian Escrabrosa Limestone (Wrucke et al., 2004) (Figs. 2 and 4, panel 4). The only siliciclastic material identified in this section was collected from 5 m of horizontally laminated, wavy-bedded, and current-ripple, cross-stratified, fine-grained quartzose sandstone; one sample was collected from this interval. Data Repository File DR2 provides a detailed report on fususlinid dating of this section.

The interval sampled from this section spans Middle Pennsylvanian through middle Permian (late Wolfcampian) (Figs. 2 and 4, panel 5). We collected one sample from a green, non-calcareous siltstone within the Horquilla Limestone and two samples from the Earp Formation (Creasy, 1967). Earp Formation samples were collected from a chert granule conglomerate in the middle Earp Formation (Red Chert Conglomerate of Armin, 1987) and from a fine-grained sandstone of the upper Earp Formation.

Samples were collected from three localities along a roughly northeast–southwest transect across the central western Paradox Basin (Figs. 2 and 6, panels 1–3). These three locations span distal, medial, and proximal basin positions along this transect. In the distal basin position, we collected samples of quartzose sandstone from two localities within the Paradox Formation and three within the Honaker Trail Formation. These units consist of nested 4^th and 5^th order sequence stratigraphic cycles of fine-grained siliciclastic and carbonate facies (Goldhammer et al., 1991). From the medial basin section (Cane Creek), we collected one sample of quartzose sandstone from the upper Honaker Trail Formation. In the proximal basin near Gateway, Colorado, we sampled two sandstones of the Cutler Formation Undivided, deposited in buttress unconformity with crystalline rocks of the ancient Uncompahgre Uplift to the northeast (Moore et al., 2008). Overall, these samples span the early Desmoinesian to Wolfcampian time. Seven samples were also collected from the eastern proximal Paradox Basin between Durango, Colorado, and just north of Ouray, Colorado (Figs. 2 and 6, panel 4). Sampled strata along this transect include the Lower Pennsylvanian Molas Formation; the middle Pennsylvanian Hermosa Group, which includes the Pinkerton Trail, Paradox, and Honaker Trail Formations; the overlying Rico Formation; and the lower Permian Cutler Formation (Campbell, 1981).

Pennsylvanian–Permian strata within the Eagle Basin were sampled along a generally east–west transect consisting of four detailed measured sections across the basin (Figs. 2 and 5, panels 1–4). These strata span Early Pennsylvanian through Late Pennsylvanian time. Sampled formations include the Belden (1 sample), Minturn (6 samples), the Eagle Valley (4 samples), the Eagle Valley Evaporite (1 sample), and the Maroon (4 samples). Three samples spanning Late Mississippian to Early–Middle Pennsylvanian were collected from the Central Colorado trough, the southeast extension of the Eagle Basin, northeast of Salida, Colorado, from the Leadville, Kerber, and Sharpsdale Formations (Fig. 5, panels 5–6). One sample from the Lower Pennsylvanian Belden Formation was collected near Crested Butte, Colorado (Fig. 5, panel 5).

We collected six samples from the Lower Pennsylvanian Glen Eyrie Member and fluvial strata from the lower, middle, and upper Fountain Formation (Suttner et al., 1984; Sweet and Soreghan, 2010) (Figs. 2 and 7). Above this, we collected one sample from aeolian Permian rocks that were originally mapped as the Lyons Formation (Keller et al., 2005; Sweet and Soreghan, 2010) but recently reinterpreted as equivalent to the Ingleside Formation (Sweet et al., 2015). To the south, the section is in fault contact with the Pikes Peak Granite (Keller et al., 2005), a 1.09 Ga intrusive suite (Schärer and Allègre, 1982).

We collected samples from Pennsylvanian and Permian rocks exposed in two areas within the southern Inyo Mountains of eastern California (Figs. 2 and 8). The northeastern area, north of Lone Pine, consists of the Pennsylvanian Keeler Canyon Formation and the overlying Permian Lone Pine Formation (Fig. 8, panel 1). Five samples were collected from the Lone Pine Formation. This formation is within the Owens Valley Group (Merriam and Hall, 1957; Stone and Stevens, 1987) and is subdivided, in ascending order, into Members A/B (1 sample), C (3 samples), D (1 sample), and the Reward Conglomerate (Stone and Stevens, 1987; Stone et al., 2000).

The other sampled area is located northeast of the town of Keeler, California (Fig. 8, panel 2). The late Paleozoic section here consists of the Keeler Canyon Formation and the informally named sedimentary rocks of Santa Rosa Flat, consisting of members 1–12 (Stone et al., 2009). Samples from calcic sandstones were collected from one location within the lower Salt Tram member of the Keeler Canyon Formation, one location in the sedimentary rocks of Santa Rosa Flat unit 11 (Ps11), and two locations in the sedimentary rocks of Santa Rosa Flat unit 12a (Ps12a) (Stone et al., 2009).

In addition to the locations detailed above, several isolated samples were collected from the periphery of the ARM system. Strata sampled in this manner include the sandstones within the Quadrant Formation (Pennsylvanian) and the lower member of the Phosphoria Formation (Lower Permian) in southwestern Montana, the middle Permian Rain Valley Formation in southern Arizona, and the middle Permian Kaibab Formation in southeastern California (one sample each).

All zircon samples were mechanically crushed and separated through density and magnetic methods. All processing followed standard procedures for concentrating zircon from a detrital sample (e.g., Gehrels and Pecha, 2014).

All zircon samples collected from the Lone Pine Formation and sedimentary rocks of Santa Rosa Flat (Keeler and Lone Pine Basins) were analyzed at the University of California–Santa Barbara (UCSB) Laser Ablation Split Stream (LASS) facility using a Nu Plasma high resolution multi-collector–inductively coupled plasma mass spectrometer (HR MC–ICPMS), a Nu AttoM single collector ICPMS, and an Analyte 193 excimer ArF laser-ablation system equipped with a HeLex sample cell using a 24 μm beam. Analyses were normalized against the 91500 zircon standard as a reference material (1062 Ma; Wiedenbeck et al., 1995). Data were reduced using Iolite 2.31 software in Igor Pro 6.3. Error assessment follows Kylander-Clark et al. (2013). These data are a subset of data presented in Lodes (2019).

Samples Quad1 and Phos, collected from the Quadrant and Phosphoria Formations in southwest Montana, were analyzed at the University of Washington using an iCAP RQ Quadrupole ICP-MS coupled to an Analyte G2 excimer laser with a spot diameter of 25 μm. Data reduction was conducted with Iolite, using their Geochron Data Reduction Scheme to calculate U-Pb ages uncorrected for common lead (Paton et al., 2010) and with the Andersen Routine of the VizualAge Data Reduction Scheme (Chew et al., 2014) for U-Pb ages corrected for common lead; detailed methods are available in Licht et al. (2018). Both approaches yielded similar zircon ages, and ages from samples Quad1 are Phos uncorrected for common lead.

All other zircon data presented in this study were collected at the Arizona Laserchron Center at the University of Arizona using a Photon Machines Analyte G2 excimer laser equipped with a HelEx ablation cell using a spot diameter of 20 μm and an Element2 HR–ICPMS (e.g., Gehrels and Pecha, 2014). See Data Repository for description of analytical methods and procedures. Typical random uncertainties in laser-ablation U-Pb analyses of zircon are 1%–2% with systematic uncertainties typically ca. 1%. All detrital zircon data presented in this study are presented in File DR3. Grain size was not individually measured, but in most samples, grain sizes of analyzed grains ranged from silt- to sand-sized.

All zircon data presented in this study are plotted as probability density functions, kernel density estimates, and traditional histograms. Probability density function curves were generated using DensityPlotter (Vermeesch, 2012) and/or DZStats (Saylor and Sundell, 2016). Kernel Density Estimation plots were generated using the same programs; a bandwidth of 20 m.y. was used for all kernel density estimate (KDE) curves. Histogram bin width is 50 m.y. in all plots.

This study presents new data from analyses of 72 detrital zircon samples collected over an area of ∼600,000 km². These data are compared to data from an additional 241 published samples collected over nearly the entire conterminous United States. Due to the volume of data considered and its broad geographical extent, traditional probability density function (PDF) or KDE curves for the entire data set cannot be effectively visualized and compared simultaneously, although traditional comparison is useful within individual basins (e.g., Figs. 3–8).

A variety of statistical tools has been developed to compare large detrital zircon data sets (Sircombe and Hazelton, 2004; Vermeesch, 2013; Satkoski et al., 2013; Vermeesch and Garzanti, 2015; Vermeesch et al., 2016; Saylor and Sundell, 2016; Sundell and Saylor, 2017; Vermeesch, 2018; Saylor et al., 2019). For regional comparisons in this study, we employ two methods. The first is a straight statistical comparison using the tools of Saylor and Sundell (2016). The second is a “bottom up” inverse modeling approach first applied to detrital geochronology by Sharman and Johnstone (2017).

Saylor and Sundell (2016) reviewed and tested a variety of statistical comparison methodologies for reliability, minimization of artifacts, and accuracy. For samples with n ≥ 300, cross-correlation of probability density functions showed the best performance, and we use this method in the current study.

One limitation of probability density plot cross-correlation is that similarity values calculated using this method can be sensitive to sample size below n = 300 (Saylor and Sundell, 2016). Samples analyzed as part of this study are mostly n≅100 or n≅300, although carbonate samples collected from the distal Paradox Basin had poor zircon yield, and many of the published analyses to which new data are compared also have n<100. Despite this limitation, cross-correlation of probability density plots provides the most reliable method for assessing similarity of the samples considered here. Likeness values calculated from kernel density estimation (KDE) are less sensitive to sample size but have limited capacity to distinguish differences in age spectra (Saylor and Sundell, 2016). Comparisons based on the K-S test are favored by some authors due to higher statistical robustness (e.g., Vermeesch, 2012, 2018; Vermeesch et al., 2016); however, these tests are very sensitive to sample size (Saylor and Sundell, 2016). KDE likeness statistics and Kolmogorov-Smirnov (K-S) test D-values are presented in the data repository but are not used in our interpretation.

Similarity matrices in which each analyzed sample is statistically compared to each other sample provide a detailed and granular way to determine the absolute similarity of two samples (e.g., Cowgill et al., 2016). However, large data sets such as the one presented here limit the practical functionality of this approach. To address this, we adopt an approach in which regional statistical comparison relies on comparing each individual sample to a composite data set for discrete intervals of geologic time. Here, we use aggregated samples from what we term the “Arizona shelf” (Fig. 1) divided into the four time intervals as this composite data set: Morrowan–Atokan (323–312.6 Ma); Desmoinesian–Virgilian (312.6–299 Ma); Wolfcampian (299–282 Ma); and Leonardian–lower Guadalupian (282–ca. 260 Ma). Absolute ages are assigned based on biostratigraphic correlations to well-radioisotopically dated sections (Menning et al. 2006; Schmitz and Davydov, 2012), and the ages of most strata considered in this study can be constrained to the stage level based on fusulinid and conodont taxa. The “Arizona shelf” in this study refers to a broad area of largely undeformed Pennsylvanian–Permian strata exposed across northern and central Arizona (e.g., Mckee, 1982; Blakey, 1990). Detrital zircon data from four localities in this area are employed here: the Grand Canyon (Gehrels et al., 2011), Sycamore Canyon, west of Sedona, Arizona (this study), the Mogollon Rim northeast of Payson, Arizona (this study), and a section of Pennsylvanian–Permian rocks exposed in the Plomosa Mountains of western Arizona that palinspastically restores to the Sedona–Mogollon Rim region (this study; McQuarrie and Wernicke, 2005). These sections span an age range of Late Mississippian to early Guadalupian (ca. 325–260 Ma). The Arizona shelf was selected as the statistical comparison index for this study because this area has large data sets (both sample numbers and analyzed grains) from each time interval, shows minor upsection changes in age spectra (see Results), and is adjacent to the Ancestral Rocky Mountain system. Age data for each sample and amalgamated Arizona shelf data are presented in File DR4.

We use DZnmf (Saylor et al., 2019) to factorize both probability density functions (PDFs) and kernel density estimates (KDEs, 20 m.y. bandwidth) for the 108 samples and sample groups. We compiled this data set from 208 detrital zircon samples from published data and data presented herein (File DR5 provides sample group details). We combined small samples (i.e., those with low grain counts) with geologically and geographically appropriate neighbors into sample groups to better characterize sink distributions.

Inverse source modeling requires no a priori source information but rather deconvolves constituent source distributions from the sink sample set distributions. This characteristic is particularly advantageous in circumstances where source age distributions are poorly characterized in comparison to sink distributions (Sharman and Johnstone, 2017). The data set used here is an optimal candidate for NMF source characterization because of the large number of sink samples, high dissimilarity among sink sample distributions, and relatively large grain counts (n > 150) of most samples and sample groups (Saylor et al., 2019). Factorization allows calculation of reconstructed sink samples (WH) based on the factorized age distributions (W), which are common to all samples and the weighting functions (H), which are unique to each sample.

The optimal number of factorized sources is determined in DZnmf by identification of a “breakpoint” defined by a decrease in the rate of decrease of the final residual (E in Equation 1) versus the number of sources modeled (i.e., the factorization rank, Saylor et al., 2019). In order to robustly identify the breakpoint, we calculated factorized sources and final residuals for factorization ranks 2–30 for both PDFs and KDEs.

Samples from Grand Canyon National Park have been previously presented and interpreted by Gehrels et al. (2011); however, because this group of samples represents the most stratigraphically complete section, we begin by reviewing these data (Fig. 3). Detrital zircon spectra from the Upper Mississippian Surprise Canyon Formation through the Leonardian Coconino Sandstone are composed of largely similar age spectra. All samples from this section have broad age spectra with individual grains ranging from 3.0 Ga to 2.5 Ga (Archean) to late Paleozoic (ca. 300 Ma). Two samples from the Upper Mississippian Surprise Canyon Formation are dominated by 1.10 Ga peaks with numerous small Archean–Paleozoic peaks (Fig. 3). Above this, spectra from the uppermost Surprise Canyon Formation sample as well as those up to the upper Hermit Formation are dominated by 1.80–1.60 Ga grains with subequal peaks at ca. 1.80 Ga and 1.65 Ga; these samples also contain prominent 1.20–1.10 Ga peaks. A sample from the upper Hermit Formation alters this pattern and includes abundant Paleozoic and late Neoproterozoic grains as well as a prominent population at 1.10 Ga, with only minor older populations. Samples from the Leonardian Coconino Sandstone change once again with prominent peaks at 1.1 Ga and minor peaks at 1.43 Ga and 1.65 Ga.

Samples collected for this study from correlative rocks along the western Mogollon Rim northeast of Sedona, Arizona, (Fig. 4, panel 1) span Early Pennsylvanian (Atokan) to middle Permian (Leonardian) time. These samples show a high degree of similarity with PDF cross-correlation values all greater than 0.65. All samples contain Neoarchean grains between 2.5 Ga and 2.9 Ga (Fig. 4). All samples from this location also contain abundant Paleoproterozoic grains with peaks at 1.79 Ga, 1.75 Ga, and 1.65 (most prominent). Mesoproterozoic age peaks include 1.56 Ga (sample SY-A), 1.50 Ga (samples SY-D, SY-C, SY-A), and 1.48 Ga (sample SY-B). All samples contain abundant grains between 0.95 Ga and 1.30 Ga, although the shape and center of this peak varies somewhat (Fig. 4). A sample from the Hermit Formation (sample SY-D; sensu Blakey, 1990) is the only sample in this section with a substantial number of 650–550 Ma grains. All samples in this section contain 430–295 Ma grains.

Pennsylvanian through lower Permian samples from the central Mogollon Rim (Fig. 4, panel 3; 5PS-42, 5PS-58, 5PS-82; 1PS-4) yielded similar age spectra to those from the western Mogollon Rim section (above). Middle Permian (Leonardian) samples from the Schnebly Hill Formation (3PS-100 and 4PS-186) are substantially different from samples from the lower portion of this section in that they contain prominent Neoproterozoic peaks with a 1.05 Ga peak, which is the largest peak in sample 4PS-186. The most striking difference between Schnebly Hill Formation samples and those below is the abundance, particularly in 3PS-100, of Neoproterozoic grains that make up age peaks at ca. 600 Ma. Both samples also contain prominent Paleozoic age peaks and 1–2 grains near the depositional age.

The final section included in the “Arizona shelf” grouping (see Statistical Methods) was measured and sampled in the Plomosa Mountains near Quartzsite, Arizona (Fig. 4, panel 2). These samples are similar to those from the Sycamore Canyon and Mogollon Rim sections, including age peaks and upsection trends (Fig. 4), and as such are not described in detail.

A single sample from the predominantly carbonate section in the Mescal Mountains 130 km east–southeast of Phoenix, Arizona, (Fig. 4, panel 4) contains a broad distribution of Archean to Paleozoic grains with the most prominent ages centered on 1.65 Ga. A single grain overlaps the depositional age (latest Virgilian) at 296 ± 4 Ma.

Late Pennsylvanian (Virgilian) to middle Permian (late Wolfcampian) samples from the Whetstone Mountains in southeast Arizona yielded similar age peaks in each sample (Fig. 4, panel 5) and include minor populations of Archean grains between 2.9 Ga and 2.6 Ga, peaks at ca. 1.80–1.72 Ga with a large peak at 1.78 Ga in sample 1WM-302, major peaks at 1.65 Ga in all samples, peaks at 1.50 Ga in all samples, and a minor peak at 1.40 Ga at the base of the section. All samples have peaks between 1.2 Ga and 1.0 Ga. All samples contain grains between 700 Ma and 400 Ma, with a common peak at 430 Ma.

Detrital zircon spectra from the Vail, Colorado, section are typified by unimodal age distributions (Fig. 5, panel 4). The lowest sample in the Vail section (V20) has scattered individual grains between 2.7 Ga and 1.10 Ga, a minor peak at ca. 1.82 Ga, and a major peak at 1.72–1.70 Ga. Above this, the abundance of non-ca. 1.7 Ga grains decreases, and the age peak at 1.7 Ga narrows substantially (Fig. 5). The only exception to this is sample V10, collected from the middle of the Minturn Formation (Fig. 5). This sample contains minor Archean peaks, major peaks at 1.78 Ga and 1.69 Ga, multiple minor peaks between 1.45 Ga and 1.00 Ga, and two minor peaks between 500 Ma and 300 Ma.

In the Red Canyon section (Fig. 4, panel 3), a sample from the Eagle Valley Evaporite (EC01) contains a major peak at 1.7 Ga, secondary peaks at 1.43 Ga and 1.11 Ga, and minor peaks at 2.67 Ga, 1.88 Ga, and 380 Ma (Fig. 5). Above this, the sample from the lower Eagle Valley Formation (RC01) contains a single, narrow dominant peak centered on 1.71 Ga with scattered individual grains between 2.8 Ga and 400 Ma. Sample RC19 from the upper–middle Eagle Valley Formation contains a broad range of ages that includes major peaks at 1.75 Ga, 1.57 Ga, 1.45 Ga, and ca. 570 Ma as well as minor peaks at 1.30 Ga, 1.08 Ga, 980 Ma, and 420 Ma. This sample is moderately similar (PDF cross-correlation of 0.48) to sample V10 from the Vail section, which was collected from nearly identical stratigraphic positions.

Most samples from the Glenwood Springs section (Fig. 5, panel 2) are dominated by two subequal peaks at 1.70 Ga and 1.42 Ga. In samples with these peaks, scattered Archean grains as well as small peaks at 1.89 Ga, 1.08 Ga, and scattered individual Neoproterozoic and early Paleozoic grains are present but decrease in abundance upsection. One exception to this pattern is sample GS09, which was collected from the upper Eagle Valley Formation. This sample contains a much broader range of ages than samples above and below, and major peaks are centered on 1.73 Ga, 1.65 Ga, 1.43 Ga, and 1.15 Ga, with smaller peaks at 1.29 Ga and 1.07 Ga; scattered individual Archean–early Paleozoic grains are also present. Stratigraphically, this sample is roughly correlative with samples V10 and RC19 described above.

The stratigraphically lowest sample from a section near Meeker, Colorado, (Fig. 5, panel 1; MC05) is dominated by two subequal age peaks at 1.71 Ga and 1.42 Ga; scattered Archean–early Paleozoic grains make up small peaks at 1.06 Ma and 510 Ma (Fig. 5). Above this, sample MC17 contains a broad distribution of grains with many individual peaks (Fig. 5). The most prominent of these are peaks at 1.80 Ga, 1.73 Ga, 1.65 Ga, and 1.16 Ga. Sample MC20 contains a similar overall distribution of age peaks except for the appearance of a major 1.45 Ga peak as well as a minor peak at 420 Ma.

A sample from the Upper Mississippian Leadville Limestone (Fig. 5, panel 6) contains major populations at ca. 1.79 Ga and 1.43 Ga with minor peaks at 1.68 Ga, 1.38 Ga, and 1.13 Ga. Above this, samples from the Morrowan Kerber and Atokan Sharpsdale Formations (De Voto et al., 1986) are dominated by a single peak at ca. 1.70 Ga. The Kerber Formation also contains scattered minor Archean to late Paleozoic peaks, and these decrease in prominence up section in the Sharpsdale Formation.

Samples from the Middle to Upper Pennsylvanian Paradox and Honaker Trail Formations within the distal Paradox Basin (Fig. 6, panel 3) mostly returned poor zircon yield, and the number of grains analyzed from each range between 36–103. These samples are qualitatively similar and contain sparse 2.5–3.0 Ga grains, abundant 1.80–1.40 Ga grains (except for sample HT02), abundant 1.30–0.90 Ga grains, and several small populations of 500–300 Ma grains. Sample HT02 deviates from this overall pattern because its largest population is composed of 500–400 Ma grains. Upper Honaker Trail samples have notable multiple small peaks at ca. 1.80 Ga, 1.72 Ga, 1.65 Ga, 1.50 Ga, 1.32 Ga, 1.20 Ga, 1.15 Ga, 1.05 Ga, and 0.98 Ga.

Sample KA03, collected from the Honaker Trail Formation in the western medial basin (Fig. 6, panel 2), contains major age peaks centered on 1.70 Ga and 1.42 Ga with minor peaks at 1.00–1.12 Ga and 450 Ma.

Proximal basin samples from the Cutler Formation (undifferentiated) near Gateway, Colorado, (Fig. 6, panel 1) yielded age spectra dominated by 1.40–1.45 Ga grains with a minor secondary population centered on 1.73 Ga, although this secondary population is nearly absent in sample TP01.

In the eastern Paradox Basin (Fig. 6, panel 4), sample ML01 from the Morrowan–Atokan Molas Formation contains a broad range of ages including scattered 1.8–3.0 Ga grains, several small peaks between 1.62 Ga and 1.30 Ga, a major peak at 1.05 Ga, and two minor peaks between 500 Ma and 400 Ma. Stratigraphically above this, age spectra change drastically to spectra dominated by a major peak at 1.43 Ga, a secondary but still major peak at 1.73 Ga, and a small peak at 520 Ma. A few samples (SM01, SP01, DMR01) contain 1–2 single grain ages between 1.20 Ga and 1.00 Ga. Sample HC01, collected from the Cutler Formation (undifferentiated) north of Durango, Colorado, departs markedly from other samples in this area. Although the same age peaks are present, this sample is dominated by a peak at 520 Ma, and the 1.73 Ga and 1.43 Ga peaks are much less prominent.

Samples from Denver Basin strata analyzed in this study can be broadly organized into three similar groups arranged stratigraphically (Fig. 7).

1. Samples from Atokan strata within the Glen Eyrie and Fountain Formations (samples MS01, MS10, MS11) are similar and have broad age spectra with most grain ages concentrated between 1.90 Ga and 1.00 Ga (Fig. 7). Within this range, major age peaks are centered on 1.75–1.70 Ga, 1.45–1.42 Ga, and 1.10–1.08 Ga, with a minor peak centered on 1.21 Ga. Sample MS01 from the Glenn Eyrie Formation contains three individual grains between 600 Ma and 350 Ma.

2. Stratigraphically above this, the second group of samples (MS06, MS05, MS04) is statistically similar with PDF cross-correlation coefficients > 0.81. These samples contain one major peak centered on 1.09 Ga and a minor peak at 1.43 Ga that vanishes entirely up-section such that sample MS04 contains only the 1.09 Ga peak (Fig. 7).

3. At the top of the section, a sample from the Wolfcampian Ingleside Formation (Sweet et al., 2015) contains a broad distribution of ages including minor peaks centered at ca. 2.7 Ga, 1.82 Ga, 1.62 Ga, and 600 Ma and major peaks centered at ca. 1.05 Ga and 430 Ma.

One sample from the Salt Tram Member of the Keeler Canyon Formation (Stevens et al., 2001) contains minor 2.8–2.5 Ga grains, multiple 1.9–1.3 Ga peaks, a major peak at 1.1–1.0 Ga, and a prominent 450–400 Ma peak (Fig. 8, panel 2).

Age spectra from samples collected in Wolfcampian strata of the Lone Pine Formation are qualitatively self-similar (Fig. 8, panel 1). Most of these samples contain scattered Archean grains, and minor age peaks occur between ca. 1.85 Ga, 1.65 Ga, 1.43 Ga, and 600–550 Ma. Major peaks are centered on ca. 1.1–1.05 Ga and 400 Ma. There is no systematic up-section change in these data.

Guadalupian samples also contain broad distributions of ages (Fig. 8, panel 2). These samples contain small Archean peaks and major peaks at 1.85 Ga, 1.05 Ga, and 400 Ma. Additional peaks present in single samples include a major 950 Ma peak in sample 6–20–17–3 and minor peaks at 1.65 Ga and 1.5 Ga in sample 6–20–17–4. Grains with 500–400 Ma ages are more abundant in the upper two samples in this group.

Results of similarity calculations between each sample and amalgamated samples of the Arizona shelf are presented visually in Figure 9. The aggregated Arizona shelf data set is also compared to aggregated data sets from six potential source terranes, which are discussed in detail in the subsequent discussion section. It is important to note that the cross-correlation values shown in Figure 9 reflect only the similarity of each individual sample or source terrane to the aggregated Arizona shelf data set in a specific timestep; there is no measure of similarity of individual samples to other individual samples in Figure 9. Cross-correlation values can range from zero (completely different) to one (identical). Complete results are presented in Files DR6 and DR7.

Samples from Morrowan–Atokan (323–321.6 Ma) strata show substantial statistical separation between the Arizona shelf and other Laurentian strata (Fig. 9) with all Laurentian, non-Arizona shelf samples returning values of <0.5. All samples within the ARM region have cross-correlation values <0.5, and most samples have values near 0.2. Comparison of the Arizona shelf during this time step to potential source regions yielded values of 0.1–0.5, with the northeast Ellesmerian orogen having the highest value (0.45).

Calculations from Desmoinesian–Virgilian (312.6–299 Ma) samples show low cross-correlation values (<0.4) between proximal ARM basin strata within the Paradox, Eagle, and Denver Basins. Strata in distal portions of the Paradox and Eagle Basins are more similar to the Arizona shelf, and these samples and those from the Bighorn and Wood River Basins have cross-correlation values >0.6. Samples from the Appalachian Basin, American midcontinent, and Ouachita-Marathon Basins have low cross-correlation values (< 0.3 except for two samples). Aggregated potential source regions yielded cross-correlation values of 0.1–0.5 with the northeast Ellesmerian orogen having the highest value (0.48).

Samples from Wolfcampian (299–282 Ma) strata have the highest homogeneity within the Arizona shelf (cross-correlation values mostly >0.7), and samples collected from southern Arizona have similarly high cross-correlation values. Two samples from northern Utah and central Idaho have cross-correlation values >0.6. Like underlying Pennsylvanian strata, proximal Paradox samples have cross-correlation values <0.2 with most near zero. In contrast with the underlying Pennsylvanian strata, Denver Basin samples have values from 0.3 to 0.6. Samples from the Lone Pine Basin in eastern California have cross-correlation values of 0.2–0.4. Potential source region cross-correlation values are 0.2–0.6, and like underlying strata, the highest value (0.56) was calculated for northeast Ellesmerian orogen data.

Leonardian–early Guadalupian (ca. 282–260 Ma) strata show higher cross-correlation values than older strata (Fig. 9), and most samples from within the ARM region have cross-correlation values mostly >0.6. Proximal Paradox Basin samples are more similar to the Arizona shelf than are samples from older strata with cross-correlation values of 0.2–0.7. Samples from Permian Basin strata of this age have cross-correlation values of 0.1–0.5. Potential source terranes have values of 0.1–0.6, with the northeast Ellesmerian orogen having the highest value of 0.55.

Samples used in DZnmf are distributed across North America and divided into five age divisions: Mississippian, Early Pennsylvanian, Middle–Late Pennsylvanian, early Permian, and middle Permian. Factorization of PDFs and KDEs produced similar results but differed in their identification of the optimum number of sources. Factorization of PDFs and KDEs yielded either 8 or 10 as the optimal number of sources, respectively. However, the 8-source factorization for both PDF and KDE models, which are nearly identical (Fig. 10), is more reasonable when considering modeling artifacts and geologic context (see discussion of NMF caveats below); we therefore focus discussion on results produced by the PDF model. The mean of quantitative similarity metrics for the PDF optimal source number model indicate that 95% of the reconstructed samples yield close matches to their empirical counterparts (likeness > 0.7 (94%), cross-correlation > 0.7 (98%), V values < 0.15 (96%), D value < 0.08 (92%); see Files DR8–DR10). The median cross-correlation coefficient is 0.87 with a maximum of 1.00 and a minimum of 0.61. We projected the factorized weighting functions (H in Equation 1) onto east-northeast–west-southwest and northwest–southeast North America-scale cross sections for each of the five age divisions (Figs. 10 and 11).

The eight factorized sources produced by DZnmf exhibit geologically reasonable age distributions (Fig. 10) and are described below. Source A contains two dominant age modes in the early Paleozoic and late Neoproterozoic and a lesser age population from 0.9 Ga to 1.1 Ga. Source B contains an asymmetrically shaped age distribution from 0.9 Ga to 1.4 Ga with a mode at 1.1 Ga and minor early Paleozoic populations. Source C contains a 1.8 Ga mode with a broad age distribution and a lesser Archean age population of 2.5–2.7 Ga. Source D contains a dominant age mode at 1.65 Ga flanked by lower amplitude modes at 1.5 Ga and 1.75 Ga, a broad suite of ages from 1.1 Ga to 1.4 Ga, and minor early Paleozoic and Archean age modes. Source E contains a dominant 2.7 Ga age mode flanked by minor Archean age populations and several minor and discrete Proterozoic age populations centered around 1.1 Ga, 1.3 Ga, 1.5 Ga, 1.65 Ga, and 1.85 Ga. Source F exhibits a unimodal age distribution at 1.7 Ga. Source G contains a dominant 1.4 Ga mode and a lesser 1.7 Ga mode. Source H contains a unimodal age distribution at ca. 530 Ma.

Archean crust in Laurentia is confined to several discrete regions within the northern United States, central Rocky Mountains, and Canadian Shield (Fig. 12). These include the Superior, Slave, and Wyoming Cratons and the Hearne and Rae provinces, which are interpreted to be fragments of Archean cratons surrounded by Proterozoic mobile belts (Whitmeyer and Karlstrom, 2007; Condie et al., 2009). Zircon derived from these cratons yields broadly similar ages with a few key exceptions. Dating of zircons from the Superior Craton indicates episodic zircon production in this region from 3.5 Ga to 2.6 Ga (Bickford et al., 2006; Schmitz et al., 2006). This resembles zircon populations recovered from the Canadian Superior province (Percival et al., 1994). Zircon from the Superior Craton form an age peak at ca. 2.7 Ga as do grains from the Hearne Craton (Condie et al., 2009). In contrast, zircon ages from the Wyoming Craton are more widely distributed with the most prominent peak centered at ca. 2.65 Ga and a greater abundance of > 3.0 Ga grains than in Archean terranes farther north.

Grains between 2.5 Ga and 1.8 Ga could have originated from a variety of Laurentian sources. During the 2.4–2.0 Ga time period, a series of rifting events occurred and was followed by a series of collisions between 2.0 Ga and 1.80 Ga that assembled the Slave–Rae–Hearne and Superior Cratons (Sims et al., 1993; Whitmeyer and Karlstrom, 2007). The remnants of these collisions are preserved as juvenile arcs and accreted material separating Archean blocks (Fig. 12; Whitmeyer and Karlstrom, 2007). Roughly concurrent convergent orogenesis at 1.85–1.78 Ga and 1.88–1.83 Ga occurred in the Trans-Hudson and Penokean events, respectively, in northern Canada and the north central United States (Van Schmus, 1976; Schulz and Cannon, 2007; Whitmeyer and Karlstrom, 2007).

Proterozoic basement rocks of southwest Laurentia were accreted to the Wyoming and Superior Cratons, are well dated, and likely contributed sediment to many of the strata analyzed in this study. The Yavapai province was initially defined only for basement rocks of northern Arizona (Karlstrom and Bowring, 1988; Karlstrom et al., 1993), but it is now understood to extend from Arizona to the southern Great Lakes region (Fig. 12; Hoffman, 1989; Whitmeyer and Karlstrom, 2007; Gehrels and Pecha, 2014). Yavapai province basement generally yields zircon crystallization ages between 1.8 Ga and 1.7 Ga (Gehrels et al., 2011; Karlstrom et al., 1987; Condie and Knoper, 1986; Bickford et al., 1989, 2008). The ca. 1.7–1.6 Ga Mazatzal province extends from southern Arizona into the central and southern mid-continent (Whitmeyer and Karlstrom, 2007; Gehrels et al., 2011; Gehrels and Pecha, 2014). Though Mazatzal province ages are generally younger than Yavapai ages, the boundary between the two provinces is somewhat diffuse, reflecting orogenic processes that overprinted the region several times as terranes were accreted to Laurentia’s cratonic core (Karlstrom et al., 1987; Karlstrom and Bowring, 1988; Duebendorfer et al., 2015). In central and southern Arizona, Mazatzal province basement ages range between 1.76 Ga and 1.62 Ga (Karlstrom and Bowring, 1988; Karlstrom et al., 1993). Granites in central Colorado were emplaced at 1.7–1.66 Ga (Bickford et al., 1989), and northern New Mexico (mostly Yavapai) basement comprises rocks dated at 1.8–1.63 Ga (Maxon, 1976; Karlstrom et al., 2004). Southern New Mexico (Mazatzal) basement has been dated at 1.73–1.61 Ga (Stacey and Hedlund, 1983; Bowring et al., 1983; Roths, 1991).

Mesoproterozoic intrusions of the 1.5–1.3 Ga granite–rhyolite province in the U.S. midcontinent and coeval alkaline plutons in western Laurentia (Bickford et al., 2015) may also be important sources for zircons in ARM basins. These intrusive bodies gained the name “anorogenic” due to their general lack of metamorphic fabrics and timing relative to more deformed Yavapai–Mazatzal rocks; however, closer study has revealed that these rocks formed while orogenic processes were still ongoing (Bickford et al., 2015). Due to their alkaline composition, they are still referred to as “A-type” granites (Bickford et al., 2015). Mesoproterozoic plutonism initiated in eastern Laurentia around 1.5 Ga, and most ages in this area are 1.5–1.45 Ga with an additional peak in age distribution at ca. 1.37 Ga (Bickford et al., 2015). Mesoproterozoic magmatism appears to have younged to the southwest, and granite–rhyolite province rocks in the U.S. midcontinent have been dated at ca. 1.5–1.35 Ga with major peaks at 1.47 Ga and 1.36 Ga (Bickford et al., 2015). From the Rocky Mountains west (the area most proximal to ARM uplifts), Mesoproterozoic pluton ages are 1.48–1.34, with most ages at 1.45–1.41 Ga (Bickford et al., 2015).

The Grenville province consists of a band of ca. 1.3–0.9 Ga crystalline rocks stretching discontinuously from Newfoundland to Texas and into Mexico (Fig. 12; Rivers, 1997; Bickford et al., 2000; Heumann et al., 2006; Whitmeyer and Karlstrom, 2007). This province formed by accretion of arc and continental blocks to (modern) eastern Laurentia at ca. 1350–1220 Ma during the Elzevirian orogeny (Moore and Thompson, 1980), 1200–1140 Ma during the Shawinigan orogeny (Rivers, 1997; Corrigan and Breemen, 1997), and culminated with the interpreted collision of the Amazonia Craton during the 1090–1020 Ma Ottawan orogeny (Rivers, 1997; Hynes and Rivers, 2010). Deformation continued in the northern part of the orogen until ca. 980 Ma (Hynes and Rivers, 2010). In Texas and southeastern New Mexico, Grenville-age rocks with zircon U-Pb ages between 1.38 Ga and 1.07 Ga are exposed in the Llano Uplift, the Van Horn region, and the Franklin Mountains (Walker, 1992; Bickford et al., 2000).

Within the ARM system, the late Mesoproterozoic Pikes Peak batholith, now exposed west of Colorado Springs, may have shed sediment into adjacent basins. This voluminous batholith is considered a type example of A-type plutonism but has been precisely dated at ca. 1.09 Ga, several hundred million years after the most voluminous phase of A-type plutonism in the western United States (Schärer and Allègre, 1982; Smith et al., 1999; Guitreau et al., 2016). Recent analyses of zircons from ring dikes intruding roof pendants within this batholith yielded ages between 1115 ± 12 Ma and 1078 ± 11 Ma (Guitreau et al., 2016). A granitic pluton in the Little Hatchet Mountains of southern New Mexico has also been dated at 1.07 Ga (Amato and Mack, 2012), and plutons of the Aibo granite in northern Sonora have been dated at 1.11 Ga (Farmer et al., 2005). Diabase intrusions of similar age are widespread throughout the southwestern United States (Heaman and Grotzinger, 1992; Stewart et al., 2001), although here we assume that diabase would be substantially less zircon fertile than granitic intrusions (e.g., Moecher and Samson, 2006). Bimodal volcanic rocks erupted at ca. 1.1 Ga in the midcontinent rift system and now exposed northern Minnesota could also have been potential sources for detritus zircon of this age (Green and Fitz, 1993).

Bimodal rift-related igneous suites exposed in the Amarillo–Wichita uplift have been dated at 539–530 Ma (Thomas et al., 2012; Hanson et al., 2013). A detrital sample collected from the Post Oak Conglomerate adjacent to the Amarillo–Wichita uplift and interpreted to have been sourced exclusively from basement rocks therein yielded a unimodal age peak at 535 Ma with 95% of grains dated between 560 Ma and 515 Ma (Thomas et al., 2016).

Neoproterozoic and early Paleozoic alkalic plutonic suites in north–central Idaho have been dated at 665–650 Ma (Lund et al., 2010) and 500–485 Ma (Evans and Zartman, 1988; Lund et al., 2010), and similar age igneous rocks are present along most of the western miogeocline boundary (Lund et al., 2010). Recent investigation of mafic rocks in the southern Yukon Territory yielded zircon ages of 488–473 Ma (Campbell et al., 2019).

Zircons separated from a hornblende–biotite syenite in the Wet Mountains of south–central Colorado yielded a crystallization age of ca. 523.98 ± 0.19 Ma (Schoene and Bowring, 2006). Zircons from syenites in the Powderhorn area of southwest Colorado yielded ages of 525 Ma and 583 Ma (Jaffe et al., 1959, p. 127). Rocks within this suite have yielded K-Ar and Rb-Sr ages of 579–565 Ma (Olson et al., 1977). Zircons from alkalic rocks in the Florida Mountains of southwestern New Mexico have been dated between 523 Ma and 503 Ma (Loring et al., 1987; Evans and Clemons, 1988).

Several episodes of orogenesis within the Appalachian system produced zircon-rich plutons and metamorphic rocks (Becker et al., 2005, 2006; Bream et al., 2004; Park et al., 2010). Early stages of the Taconic orogeny produced volcanic rocks at 550–500 Ma (Drake et al., 1989), and kyanite-grade metamorphism and magmatism occurred at 470–440 Ma (Hatcher et al., 1989; Drake et al., 1989). The Devonian to Early Mississippian Acadian orogeny resulted in abundant 423–315 Ma plutonism (Osberg et al., 1989; Bradley et al., 2000). Finally, 330–265 Ma plutonism within the Appalachian belt is associated with the Alleghenian orogeny, which marks the final stage of collision with Africa and Laurentia (Hatcher et al., 1989). There is little evidence of similar igneous activity in the Ouachita–Marathon belt to the southwest.

The results of NMF must be interpreted within the geologic framework in which samples are found. Interpretation of samples outside of this context is likely to result in incorrect conclusions. Below we present two potential pitfalls inherent in NMF, including (1) factorization artifacts, and (2) incorrect attribution of sources. A common artifact of the NMF modeling of this data set is the presence of nadirs (low points) in more complex distributions that align with unimodes in simple factorized source distributions. We interpret this behavior as a function of the algorithmic nature of NMF where, once a unimodal source is identified and factorized out, inclusion of that mode in complex distributions is not necessary, and hence the more complex distributions are left with anomalous nadirs. The ca. 1.4 Ga nadir in Source B (Fig. 10) is an inconsequential example of this phenomenon, but the 1.4 Ga and 1.7 Ga nadirs in Source D may contribute to the difficulty in identifying its real source match. As an example of the second pitfall, we conclude that NMF sediment source characterization incorrectly attributes the broad Appalachian-like source (B) to Denver Basin samples (MS04, MS05, MS06, and perhaps MS11; Fig. 7). Denver Basin samples exhibit zircon age unimodes at ca. 1.1 Ga, but we interpret these ages to be locally sourced from the Pike’s Peak batholith in the Ancestral Front Range based on the leptokurtic (narrow) nature of the peak. In PDF and KDE NMF model runs with more sources (> 8), a leptokurtic ca. 1.1. Ga unimodal age distribution is identified, but this results in a nadir in a source distribution that otherwise closely matches Appalachian sources (B; Fig. 10 and see Files DR8–DR10). These issues may be common features of inverse modeled source scenarios that contain both simple and complex distributions and should be considered when interpreting factorized sources.

We suggest real source matches for NMF-generated sources with varying degrees of confidence. Source A’s age distribution is similar to the northern Appalachian detrital zircon source invoked for Pennsylvanian sandstone in the Illinois and Forest City Basins (Kissock et al., 2018). In this interpretation, exhumation of Pan-African metasedimentary terranes (presently in the subsurface) within the northern Appalachian Mountains, containing Neoproterozoic zircons (Zartman et al., 1988; Heatherington et al., 1999; Wortman et al., 2000; Hibbard et al., 2002; Pollock et al., 2007; Fyffe et al., 2009), mixed with early Paleozoic (490–270 Ma) and Grenville (1.3–0.9 Ga) detrital zircons while en route to the midcontinent (Kissock et al., 2018) and possibly North American western margin. Source B’s age distribution is similar to the commonly invoked Appalachian detrital zircon age signature observed in the central Appalachian Basin (Fig. 13) from Mississippian through early Permian (Eriksson et al., 2004; Thomas et al., 2004; Becker et al., 2005, 2006; Gehrels et al., 2011; Kissock et al., 2018; Thomas et al., 2017). Source C’s age distribution closely resembles early Paleozoic rocks along North America’s western margin exhumed during the Antler orogeny (Beranek et al., 2016). Source D is the most difficult factorized source to find a realistic match for and is discussed in greater detail below. Source E’s dominant Archean modal age is similar to the cratonic Archean basement terranes of North America (e.g., Superior and Hearne Cratons; Condie et al., 2009) but also closely resembles Cambrian sandstone detrital zircon distributions in the midcontinent (Konstantinou et al., 2014). Source F’s 1.7 Ga unimodal age distribution is a close match to North American basement rocks in the Yavapai–Mazatzal terrane (Whitmeyer and Karlstrom, 2007; Gehrels et al., 2011; Laskowski et al., 2013). Source G also closely resembles a North American basement terrane commonly referred to as A-Type (a.k.a. “Picuris”; Daniel et al., 2013) granitoids with a minor contribution from the associated Yavapai–Mazatzal basement terrane (Anderson and Morrison, 1992; Bickford and Anderson, 1993). The age distribution of Source H is consistent with Cambrian age igneous rocks of North America, which are primarily associated with the Oklahoma Aulacogen (Thomas et al., 2016) but also present as minor alkali intrusions throughout the western interior (Olson et al., 1977; Loring et al., 1987; McMillan and McLemore, 2004; Schoene and Bowring, 2006).

Most of the factorized sources produced by DZnmf exhibit geologically reasonable spatial distributions consistent with the source interpretations discussed above. The Appalachian-like synthetic sources (A, B) are more prevalent to the east but are present in western depocenters from Mississippian through middle Permian. This is consistent with the development of a North American transcontinental sediment dispersal system by Late Mississippian time (Gehrels et al., 2011; Chapman and Laskowski, 2019). Both Ely Basin samples (Tonka, Battle) exhibit factorized Source C dominance, suggesting a recycled Paleozoic contribution from the Mississippian through Middle–Late Pennsylvanian (Beranek et al., 2016). However, factorized Source C dominance in the Mississippian Wood River Basin yields to Appalachian-like sources (A, B) and Source D dominance in the Middle–Late Pennsylvanian and early Permian. These observations support models of continued deformation along the Antler Orogenic front during the Pennsylvanian (Trexler et al., 2004; Theodore et al., 2004; Sturmer et al., 2018) but also suggest partitioning of a once-connected foreland basin in the Pennsylvanian that was potentially associated with onset of Ancestral Rocky Mountain (ARM) deformation.

Source D is the most difficult factorized source to find a realistic source match for. However, its predominance in the western portion of North America and increase in prevalence in the Early Pennsylvanian (Figs. 10 and 11) suggest a western-derived source that experienced exhumation coeval with early stages of ARM deformation. Its abundance in Arizona shelf samples could suggest nearby ARM exhumation within the Zuni–Defiance uplift. Conversely, the widespread, north-to-south distribution of Source D suggests a broad western margin source that was perhaps associated with continued deformation along an expanded Antler orogenic front to the north and south of the Ely Basin.

Punctuated and highly localized dominance of ARM basins by North American basement sources within them (Sources E, F, G, and H) highlights both the heterogeneous nature of locally sourced basins and the timing and pattern of ARM deformation. The Central Colorado trough and Denver Basin were the first ARM basins to receive the majority of their sediment from local basement uplifts in the Early Pennsylvanian. At this time the Denver Basin received its greatest proportion of Y–M basement-like factorized sources (Sources F and G). This proportion may be even greater if, in consideration of the issue surrounding incorrect attribution of the leptokurtic 1.1 Ga age mode to the Appalachian-like Source B described above, we consider Source B to be a locally derived basement source. This interpretation is bolstered by the increase in Source B in the Denver Basin during the Middle–Late Pennsylvanian expansion of ARM deformation (i.e., upper and middle Fountain Formation, sensu Sweet and Soreghan (2010); Fig. 7). The Paradox Basin samples are not dominated by basement-like factorized sources until the Middle–Late Pennsylvanian. During this time the Paradox Basin samples are dominated by a ca. 1.4 Ga zircon population (Source G), whereas the Central Colorado trough–Eagle Basin samples have a dominant ca. 1.7 Ga age mode (Source F). By the early Permian, areas of the Paradox Basin received exclusively locally derived basement detritus. Similarly, basement-derived (ca. 530 Ma; Source H) detrital zircons exclusively sourced areas in the Anadarko Basin from its bounding uplift (i.e., the Amarillo–Wichita Uplift). This pattern of basement-derived sediment dominance indicates spatial variability in ARM basin-uplift development and localized sediment sourcing from particular basement blocks but does not support a clear directional trend in deformation as suggested by other authors (Dickinson and Lawton, 2003).

Sediment sources in the middle Permian are more spatially homogenous and exhibit less inter-basin variability than in the early Pennsylvanian to early Permian (Figs. 10 and 11). Middle Permian factorized source attribution indicates increased sourcing by Source A and a marked decrease of basement-like sources. We interpret this regional change in sediment sourcing to reflect cessation of ARM deformation and basinal accommodation, the progressive infill of formerly marine ARM basins by far sourced aeolianites and fluvial strata, and reintegration of regional drainage networks.

Paleogeographic reconstructions of the late Paleozoic suggest that zircon could have been sourced from substantially outside the conterminous United States (Gehrels et al., 2011). To shed light on the potential sources for non-locally derived sediment, we compare Arizona shelf data to zircon data collected from the northeast Ellesmerian orogen, northwest Ellesmerian clastic wedge, northern and central Appalachian orogen, Ouachita–Marathon fold thrust belt (Gondwanaland sources), and Sonoran miogeoclinal rocks (Fig. 13).

Data from the northeast Ellesmerian orogen consist of 27 Neoproterozoic to Devonian sandstone from Franklinian Basin strata (Anfinson et al., 2012). These strata were deposited in the foreland basin of the Devonian–early Carboniferous Ellesmerian collision, which shed large volumes of sediment to the south (Schack Pedersen, 1986; Embry, 1988; Patchett et al., 1999; Gottlieb et al., 2014). These samples contain a broad Archean peak centered at ca. 2.7 Ga, abundant 2.0–1.6 Ga grains with peaks at 1.99 Ga, 1.86 Ga, and 1.65 Ga, a broad range of Mesoproterozoic and early Neoproterozoic grains with peaks at 1.49 Ga, 1.35 Ga, and 1.08 Ga, and 800–400 Ma grains with a minor peak at 650 Ma and major peaks at 540 Ma and 420 Ma.

To the west, 13 samples from the Late Devonian to Early Mississippian Ellesmerian clastic wedge, referred to as “northwest Ellesmerian orogen” throughout this study (Beranek et al., 2010), include a broad range of Archean through Neoproterozoic grain ages with a large peak centered at ca. 1.9 Ga. These samples also contain a prominent peak at ca. 430 Ma. Sediment eroded from the northeast Ellesmerian orogen has been interpreted to have transported through this area (Garzione et al., 1997; Beranek et al., 2010; Gehrels and Pecha, 2014), but it is also possible that some northeast Ellesmerian orogen sediment bypassed this portion of the basin through a sediment conduit farther east (Anfinson et al., 2012).

Detrital zircon data from the northern Appalachians (Fig. 13) are compiled from samples collected from Newfoundland and the Maine Appalachian orogen (Cawood and Nemchin, 2001; Bradley and O’Sullivan, 2017). In contrast with central Appalachian strata, the most abundant grains are ca. 500–400 Ma and define a major peak at 430 Ma. Northern Appalachian samples also contain abundant 1.2 Ga to 950 Ma grains with minor abundances of Archean to early Mesoproterozoic grains.

Central Appalachian orogen data consist of 22 samples from Precambrian sedimentary and metasedimentary rocks within the Blue Ridge and Inner Piedmont zones of the central Appalachians (Bream et al., 2004; Eriksson et al., 2004) and 14 samples from Cambrian to Devonian sedimentary rocks within the Appalachian foreland basin (Eriksson et al., 2004; Park et al., 2010). Collectively, Appalachian foreland basin samples contain abundant ca. 1.25 Ga to 950 Ma grains with peaks at 1.18 Ga and 1.02 Ga, moderate abundances of 1.55 Ga to 1.25 Ga grains, and minor abundances of 500–400 Ma grains (Fig. 13).

Ouachita–Marathon data consist of five samples from Ordovician and Silurian strata exposed in the Ouachita Mountains of Oklahoma and Arkansas (McGuire, 2017). Together, these samples contain populations at 2.6–2.5 Ga, 1.4–1.35 Ga, and 1.1–1.0 Ga (Fig. 13).

Data from four samples from Cambrian through Devonian Cordilleran passive margin rocks now exposed in Sonora, Mexico, (Gehrels and Pecha, 2014) consist of a major Archean peak centered at ca. 2.7 Ga, major peaks at 1.85 Ga, ca. 1.7 Ga, and 1.43 Ga, and a subordinate peak at 1.11 Ga (Fig. 13).

The interpretations of sediment provenance presented here begin with data from individual basins, starting in the Denver Basin and progressing counterclockwise. We then integrate these interpretations into a broader interpretation of late Paleozoic sediment routing.

Age peaks from a sample of the Glen Eyrie Member (Fig. 7; MS01) are interpreted to represent multiple sediment sources including 1.7–1.9 Ga Yavapai crust and 1.4–1.5 Ga A-type granite sources. Zircons with these ages have also been identified in the Cambrian Sawatch Formation, which was deposited disconformably above the Pikes Peak Granite (Siddoway and Gehrels, 2014), and we link grains with these ages in the Glen Eyrie Formation to erosional recycling of Paleozoic strata. Grains making up a minor peak at 1.08 Ga are interpreted to have been eroded from the upper portions of the 1.08–1.09 Ga Pikes Peak batholith just to the west of the basin. Scattered individual Paleozoic grains are interpreted to represent either minor inputs of distal Appalachian-derived grains or recycling of similar age grains from the Devonian Williams Canyon Formation within the Paleozoic section overlying the Pikes Peak granite (Siddoway and Gehrels, 2014). Upsection, the lower Fountain Formation (MS10) contains similar age peaks, but the 1.7 Ga peak is less prominent, and this sample does not contain any Paleozoic grains. Sample MS11 from the middle Fountain Formation sits directly above an intraformational unconformity (Sweet and Soreghan, 2010); an increase in prominence of the ca. 1.08 Ga peak is interpreted as increased contribution of Pikes Peak batholith sediment to the basin. This peak is increasingly dominant upsection in samples from the middle Fountain Formation (MS06, MS05, MS04). We interpret these data to indicate that material eroded from the Pikes Peak batholith was the dominant source of sediment for this portion of the Denver Basin during Desmoinesian–Virgilian time. The transition from bimodal (1.09 Ga and 1.43 Ga) to unimodal (1.09 Ga) spectra in the upper Fountain Formation coincides with a major shift in paleocurrent indicators that indicate north-northeast sediment transport in the lower and middle Fountain Formation to east-southeast sediment transport during upper Fountain Formation deposition (Sweet and Soreghan, 2010). In the upper portion of this section, the Wolfcampian Ingleside Formation yielded a broad distribution of zircon ages distinct from strata below (Fig. 7). These are interpreted to represent distal Laurentian sources such as the Canadian Shield, Grenville, and early Acadian orogeny or age equivalent, respectively. The most plausible sources for grains of these ages are the northeast and northwest Ellesmerian and northern Appalachian orogen, or some combination thereof. The Pikes Peak batholith was likely still a major contributor of sediment to the Ingleside Formation, but equal abundances of grains within the 1.00–1.05 Ga and 1.05–1.10 Ga age windows suggests the emergence of a new source from which 1.00–1.05 Ga grains were derived.

Overall, these data suggest an unroofing sequence beginning with erosion of lower Paleozoic strata (Siddoway and Gehrels, 2014) and culminating in the sole dominance of Pikes Peak granite detritus in the upper Fountain Formation. A similar unroofing trend has been interpreted by the presence or absence of quartz arenite clasts within the Fountain Formation (Sweet and Soreghan, 2010). These data suggest the onset of Ancestral Front Range exhumation by Atokan time. These data also show a rapid shift in provenance from proximally to distally derived sediment during the Wolfcampian (early Permian). This shift suggests that the basin’s sediment budget was overwhelmed by aeolian material at this time (Sweet and Soreghan, 2010; Sweet et al., 2015; Sweet, 2017), a shift that suggests exhumation may have slowed or ceased within the Ancestral Front Range. These data also point to the Denver Basin being sedimentologically isolated from other Ancestral Rocky Mountain basins until Wolfcampian time.

West of the Ancestral Front Range uplift, samples from the Central Colorado trough are interpreted to show initiation of ARM exhumation during Morrowan time (Fig. 5, panel 6). Detrital zircon spectra show a substantial change from Mississippian (Leadville Limestone) to Morrowan (Kerber Formation), and the dominance of ca. 1.7 Ga grains in the latter is interpreted to indicate exhumation of Yavapai basement rocks, although the presence of Archean to Paleozoic grains suggests a broad mix of distal Laurentian sources as well. Above this, 1.7 Ga grains dominate the spectra from the Atokan Sharpsdale Formation. These spectra are interpreted to reflect intensification of ARM exhumation adjacent to the Central Colorado trough during the Atokan stage.

Eagle Basin samples were collected along a transect from the proximal eastern basin (Vail, Colorado) to the distal northwestern basin (Meeker, Colorado) (Fig. 5). Unimodal 1.7 Ga age distributions from the Vail section are interpreted to represent sediment derived solely from Yavapai crust exposed by the paleo-Front Range uplift directly east of the basin; 1.7 Ga rocks are currently exposed just east of this location in the hanging wall of the Gore fault (Sims et al., 2001). The 1.7 Ga peak in sample V10 is interpreted to represent the same sediment source dominant in the other samples from this section; however, Archean, Neoproterozoic, and Paleozoic grains in this sample are interpreted to represent distally sourced grains transported into the proximal Eagle Basin through aeolian processes (e.g., Schenk, 1987) and/or during sea level highstands.

To the west, grain ages from the Eagle Valley Evaporite (Fig. 5, panel 3) are interpreted to represent sediment input from the same Yavapai basement source (1.7 Ga) as well as A-type granitic sources (1.4 Ga). This sample (EC01) contains additional peaks at 1.1 Ga and 400 Ma that are interpreted as distally sourced grains from either the northern Appalachian or northwestern Ellesmerian orogen or a combination thereof. Above this, the basal Eagle Valley Formation is interpreted to represent derivation from a single 1.7 Ga local source with scattered single grains representing distal sources. The upper sample in this section contains a broad distribution of grain ages and is interpreted to represent sediment with mixed provenance including locally derived sediment (ca. 1.7 Ga and 1.45 Ga) as well as northern Appalachian/Arctic derived grains. A prominent group of ages at ca. 550 Ma is interpreted to represent grains sourced from the northeast Ellesmerian orogen or the northern Appalachians (Fyffe et al., 2009) (Fig. 5).

The Glenwood Springs section (Fig. 5, panel 2) is adjacent to the ancestral Uncompahgre uplift and was located in the west–central Eagle Basin during Pennsylvanian to Permian time (Johnson, 1987). This proximity is reflected by a major change in detrital zircon spectra between this section and the Eagle and Vail sections to the east. In contrast with the eastern basin, samples in this section show two dominant age peaks at 1.7 Ga and 1.43 Ga, which are interpreted to represent sources in Yavapai and A-type granitic sources exposed by exhumation of the ancestral Uncompahgre uplift. Crystalline rocks of this age are exposed along the northern flank of the Paradox Basin (see below) and are interpreted in the subsurface south of the Eagle Basin (Sims et al., 2001). One sample from this section (GS09) contains a much broader distribution of ages, including a major peak at 1.65 Ga that is interpreted to represent Mazatzal derived sediment, and multiple small Proterozoic peaks are interpreted to represent Appalachian and/or Ellesmerian sourced sediment.

The basal sample from Meeker (MC05; Fig. 5, panel 1), the westernmost section in the Eagle Basin, contains two dominant age peaks at 1.7 Ga and 1.4 Ga that are interpreted to represent Yavapai and A-type granite sources within the ancestral Uncompahgre uplift to the south. Above sample MC05, two samples from the upper Minturn Formation have zircon spectra with broad age distributions. Major peaks at 1.65 Ga and 1.43 Ga are interpreted to represent input from Mazatzal and A-type granite sources, respectively; subordinate age peaks between 1.8 Ga and 1.73 Ga are interpreted to represent a mix of individual Yavapai sources. A prominent peak in sample MC17 at 1.15 Ga is interpreted as Grenville or Ellesmerian derived material, and small numbers of grains making peaks at 420 Ma and 580 Ma in MC20 are interpreted as Appalachian/northeast Ellesmerian, respectively. Overall, data from this section suggest that during Desmoinesian time, this portion of the basin changed from sedimentation dominated by ancestral Uncompahgre, Front Range, and Sawatch uplift-derived material to sedimentation dominated by distally sourced material from either the Appalachian orogen (Gehrels et al., 2011) or the northeast Ellesmerian orogen (Anfinson et al., 2012; Link et al., 2014). An alternative scenario could be that spectra from these samples represent recycling of grains from lower Paleozoic strata unroofed by ARM uplift. However, we prefer the former scenario because of the similarity in spectra among samples MC17 and MC20 and age-equivalent strata to the northwest (Fig. 9; Lawton et al., 2010; Link et al., 2014). In our preferred model, the local/distal sediment mixing interface shifted from west or northwest of this section during early Desmoinesian time to east or southeast of this section by late Desmoinesian time.

Two detrital zircon transects across the Paradox Basin shed light on the sediment shed southward from the ancestral Uncompahgre uplift. In the western proximal basin near Gateway, Colorado, (Fig. 6, panel 1) samples TP01 and JB01 show nearly unimodal zircon ages centered at 1.43 Ga that are interpreted to record sediment eroded from small catchments eroding A-type granites as well as minimal Yavapai rocks to produce minor populations at 1.73 Ga. Ca. 1.4 Ga granitic rocks and Paleoproterozoic metamorphic rocks are currently exposed to the northeast of this location (Sims et al., 2001), and the Cutler Formation in this location was likely sourced from similar rocks. These samples are interpreted as Wolfcampian in age (see File DR1 for detailed discussion of age interpretation), and uni- and bimodal age distributions in these samples suggest that the ancestral Uncompahgre uplift still dominated the sediment budget of the proximal Paradox Basin at this time. Although this does not necessarily require active uplift (e.g., Paola et al., 1992; Heller and Paola, 1992), it does suggest at least recent, if not active, rock uplift and exhumation of the Uncompahgre uplift during Wolfcampian time.

Ages from one sample from Upper Pennsylvanian (Missourian) strata (sample KA03) within the western medial Paradox Basin (Fig. 6, panel 2) are interpreted to represent a mix of proximal and distally sourced sediment (Fig. 6, location 2). Major peaks at 1.7 Ga and 1.4 Ga are interpreted to represent grains sourced from exhumation of Yavapai and A-type granite sources within the Uncompahgre uplift to the northeast (Sims et al., 2001). A minor, broad peak between 1.2 Ga and 1.0 Ga is interpreted to represent distal, Grenville orogen-derived grains as are grains making up a minor peak at 430 Ma. Scattered Archean grains are interpreted to have been eroded from the Canadian Shield or the Wyoming Craton.

Desmoinesian–Virgilian samples from the distal portion of the western Paradox Basin (Fig. 6, panel 3) contrast with samples from medial and proximal basin positions. Spectra from these samples contain broad distributions of Archean to Paleozoic ages. In these samples, Archean grains are interpreted as sourced from the Canadian Shield and/or Wyoming Craton, 1.8–1.7 Ga grains are interpreted as sourced from Yavapai rocks, and 1.7–1.6 Ga grains are interpreted to represent Mazatzal sources along the southern margin of the basin or northern portion of the Zuni/Defiance uplift (Wengerd and Matheny, 1958). The 1.3–0.95 Ga grains are interpreted as sourced from Grenville rocks, and 500–300 Ma grains are interpreted as either northern Appalachian or Ellesmerian in origin. Overall, these strata suggest diverse sources of sediment within this portion of the basin, and these spectra contrast strikingly with those interpreted as having a single, proximal source within the proximal basin.

Analyses of Leonardian strata across the Paradox Basin (Lawton et al., 2015; Fig. 5) show a major change in zircon spectra at the Wolfcampian–Leonardian boundary. Above this boundary, rocks of the Castle Valley and White Rim Formations contain broadly distributed ages in stark contrast to the bimodal spectra from the Cutler Group below (Fig. 5, panel 1). This reflects the shift from fluvial and alluvial depositional systems sourced from the Uncompahgre uplift to an aeolian system that overwhelmed the sedimentary environments and carried sand into the basin via longshore current and aeolian processes (Lawton et al., 2015). A prominent ca. 1.43 Ga peak is retained through the lower Castle Valley Formation (samples 12CV50, 10CVR, 10CVG) in the proximal basin (Fig. 5, panel 1), suggesting that proximal Uncompahgre sediment was a major contributor to the material that filled this portion of the basin during Leonardian time. However, 1.43 Ga grains are nearly absent from the upper Castle Valley Formation (sample 10CVW), which suggests the shutoff of Uncompahgre exhumation by this time (Fig. 5, panel 1; Lawton et al., 2015).

A transect across the eastern proximal basin (Fig. 6, panel 4) sheds light on the along-strike variability (compared to the western proximal Paradox Basin) of sediment sources contributing to the proximal basin. A sample from the Atokan Molas Formation yielded a major peak between 1.1 Ga and 1.0 Ga, a minor abundance of grains at 500–400 Ma, and a broad range of Archean–Neoproterozoic grains (Fig. 6, panel 4). This spectrum is interpreted to represent input from the Appalachian orogen via recycling of eroding Mississippian strata (e.g., Gehrels et al., 2011). Above this, all samples contain two peaks at 1.73–1.71 Ga and ca. 1.43 Ga (Fig. 6). From Desmoinesian through Wolfcampian samples, there is also an increase in prominence of a peak at 500 Ma, and grains of this age make up greater than half of the grains dated within the youngest Cutler Formation sample. These grains are interpreted to be sourced from diabase and syenitic rocks of this age in southwestern Colorado (Olson et al., 1977). Together, these data suggest an unroofing trend in which Atokan strata were sourced from erosion of the lower Paleozoic section, and Desmoinesian–Wolfcampian strata were sourced from progressive erosion of basement rocks. Increasing depth of basement exhumation and/or drainage network evolution altered the abundance of detrital ages upsection. These data in addition to previous sedimentologic observations (Wengerd and Matheny, 1958) suggest that the eastern portion of the basin was filled by sediment sourced from the San Luis highlands to the east and represent a distinct sediment source from the western basin. The absence of 550–600 Ma grains in the rest of the basin suggest that proximal sediment shed into the eastern basin did not reach the western basin.

Samples from the Arizona shelf and southeastern Arizona basins have high degrees of similarity. In statistical comparison of 17 total samples, only two samples yield PDF cross-correlation values <0.5 (see File DR6). Major peaks common to most of these samples include a 1.8–1.7 Ga peak that is interpreted to represent material eroded from Yavapai-aged rocks exposed in Ancestral Rocky Mountain uplifts to the north and east (Figs. 4–6); a 1.7–1.6 Ga peak that is interpreted to represent erosion from Mazatzal rocks in northeastern Arizona and New Mexico (Karlstrom and Humphreys, 1998; Whitmeyer and Karlstrom, 2007); and a peak slightly younger than 1.5 Ga that is interpreted to represent erosion from A-type granites (Whitmeyer and Karlstrom, 2007). Alternatively, it is possible that these ages represent recycling of ca. 1.5 Ga grains from Proterozoic quartzites exposed along the Defiance uplift in northeastern Arizona (Doe et al., 2012, 2013). Broad peaks typically centered at 1.1–1.0 Ga and commonly containing a subordinate peak between 1.2 Ga and 1.1 Ga are interpreted as “Grenville” grains. However, Grenville crystalline rocks are widespread as are Mesoproterozoic to Neoproterozoic strata containing Grenville zircon (e.g., Rainbird et al., 2017), and all potential source areas considered above except for Sonoran miogeoclinal strata could have contributed grains to produce these peaks (Fig. 13). Grains of this age also could have been recycled from Mississippian strata in central and western Laurentia (Chapman and Laskowski, 2019). Also abundant in most Arizona samples are 500–300 Ma grains with well-defined peaks at ca. 430 Ma. These grains are typically interpreted as “Appalachian” orogen grains but could also have been sourced from the northwest or northeast Ellesmerian orogen (Beranek et al., 2010; Anfinson et al., 2012). Also present in many of these samples, especially from the Wolfcampian and Leonardian samples, are minor peaks/populations between 650 Ma and 500 Ma. In classical North American zircon provenance studies, grains of this age range are typically interpreted to be sourced from the pan-African/Brasiliano Craton (Gehrels et al., 2011; Gehrels and Pecha, 2014; Alsalem et al., 2018; Chapman and Laskowski, 2019). However, these peaks are also well-represented within northeast Ellesmerian orogen sources (Fig. 13; Anfinson et al., 2012). In either case, ca. 2.7 Ga grains, 1.2–1.0 Ga grains, and 650–300 Ma grains within Arizona samples were likely eroded and transported many hundreds to thousands of kilometers from their sources and represent distally sourced sediment. Recycling of these grains from strata immediately underlying the Arizona shelf is precluded by the near absence of these grain ages within the lower Paleozoic section in this region (Gehrels et al., 2011). Arizona shelf and southern Arizona spectra collectively represent mixing of distal sources with “locally” derived Yavapai, Mazatzal, and Granite–Rhyolite grains.

Unlike samples near major ARM uplifts, Arizona samples show little upsection change (Fig. 13). The only exception to this is that 650–500 Ma grains become substantially more abundant in the western Mogollon Rim (Fig. 4, panel 1) and the Plomosa Mountains section (Fig. 4, panel 2); this section structurally restores in close proximity to the Mogollon Rim (McQuarrie and Wernicke, 2005) beginning in Wolfcampian time. A similar increase in abundance occurs beginning in Leonardian time within the Eastern Mogollon Rim (Fig. 4, panel 3). Interestingly, grains of this age are present within the Pedregosa Basin by at least Virgilian time (Fig. 4, panel 5). The increasing abundance of these grains up section and their early appearance farther south could suggest a southern, African, or South American source for these grains. However, we consider it more likely that these grains were sourced from the northeast Ellesmerian orogen or northern Appalachians (Fyffe et al., 2009). A substantial portion of the Schnebly Hill Formation and nearly the entire thickness of the Coconino Sandstone is made up of aeolian facies (Blakey, 1990), and southward paleotransport is well documented in the Coconino Sandstone (Reiche, 1938; Odyke and Runcorn, 1960; Peterson, 1988). Additionally, consistent marine carbonate deposition within the Pedregosa Basin as well as marine conditions within the Permian Basin system to the east (Greenwood et al., 1977; Frenzel et al., 1988; Soreghan, 1994; Wright, 2011) suggest that it would have been unlikely for Ouachita–Marathon derived sediment to bypass those depocenters and be transported north.

Samples collected from the Keeler and Lone Pine Basins are largely similar (Fig. 8). The age peaks at ca. 1.85 Ga are uncommon in late Paleozoic Laurentian samples presented in this study, and these ages could represent either sediment originating (1) within Ellesmerian orogen, in which these ages are abundant (Fig. 13), (2) from accreted rocks to the northwest of these basins such as the Antler and/or Klamath terranes (Gehrels, 2000; Grove et al., 2008; LaMaskin, 2012; Beranek et al., 2016), (3) from recycling of nearby lower Paleozoic passive margin rocks (Chapman et al., 2015), or (4) from erosion of the 1.8 Ga Mohave province (Whitmeyer and Karlstrom, 2007). A low abundance of 1.6–1.8 Ga grains in most of these samples is interpreted to indicate that material eroded from Ancestral Rocky Mountain basement uplifts was not a major contributor to the sediment budgets of the Keeler and Lone Pine systems. Age peaks at 1.1–1.0 Ga could represent either “Grenville” ages from the northern Appalachian orogen (e.g., Gehrels et al., 2011) or similar ages from the northeast Ellesmerian orogen and/or northwest Ellesmerian orogen (Beranek et al., 2010; Anfinson et al., 2012). The same provenance possibilities apply to peaks at 500–400 Ma (Fig. 13).

We interpret the detritus that filled Pennsylvanian–Permian Ancestral Rocky Mountain basins to have been sourced from two end-member sources: far-traveled sediment derived from the Ellesmerian orogeny and/or the northern Appalachian orogen and proximal material eroded from exposed basement rocks and potentially local recycled Mississippian strata within the Ancestral Rocky Mountain system. The zircon budgets of proximal portions of the Denver, Eagle, and Paradox Basins were dominated completely by ARM sources, whereas zircon in the medial and distal portions of the Eagle and Paradox Basins appear to consist of a mix of ARM material and distally derived material. Of the basins studied here, zircon within the deep marine basins of southeast California appears to contain the least amount of ARM-derived material (Fig. 9). These interpretations follow Gehrels et al. (2011) and Gehrels and Pecha (2014), with the notable exception of a diminished importance of Appalachian provenance in favor of sources in northern Canada.

Samples from Morrowan–Atokan (323–312.6 Ma) rocks show early signs of ARM-related exhumation. The appearance of major ca. 1.7 Ga peaks in Morrowan samples from the Grand Canyon (Gehrels et al., 2011) and Central Colorado Trough (sample BC01) as well as the Atokan appearance of ca. 1.65 Ga peaks in the Grand Canyon, Mogollon Rim (sample SY-D), and southwestern Arizona (sample 1QZ-19) are interpreted to represent the early signals of Yavapai–Mazatzal exhumation related to ARM uplift surrounding these basins. However, these spectra contain a mix of other ages, not the uni- or bimodal spectra that typify Late Pennsylvanian samples in these sections. Interestingly, samples from the Upper Mississippian Surprise Canyon Formation in Grand Canyon strata contain peaks at 1.7 Ga and 1.43 Ga (Gehrels et al., 2011), and the Upper Mississippian Leadville Limestone in the Central Colorado Trough (Fig. 5, panel 6) has peaks at ca. 1.42 Ga and 1.7 Ga, allowing for the possibility that ARM exhumation initiated during the latest Mississippian. Conversely, Atokan samples from the Molas Formation within the proximal Paradox Basin do not contain significant numbers of 1.7 Ga or ca. 1.4 Ga grains (sample ML01; Evans and Soreghan, 2015; Nair et al., 2018), both of which become the dominant grain ages during the Middle–Late Pennsylvanian. This suggests that substantial exhumation of Uncompahgre basement rocks did not occur until Desmoinesian time, although unroofing of lower Paleozoic rocks deposited above basement could have occurred earlier.

At the continental scale, Morrowan–Atokan Arizona shelf samples from this time interval are statistically distinct from samples proximal to nascent ARM uplifts as well as from samples within the Ouachita–Marathon foreland (Buratowski, 2014), the Illinois and Forest City Basins (Kissock et al., 2018), and the Appalachian Basin (Becker et al., 2005; Park et al., 2010) (Fig. 9). Compared to aggregated data from potential source terranes, aggregated Arizona shelf data are most statistically similar to the northeast Ellesmerian orogen (cross-correlation coefficient of 0.45; Anfinson et al., 2012) and the Sonoran margin (cross-correlation coefficient of 0.42; Gehrels and Pecha, 2014). We did not calculate cross-correlation values between the Arizona and mixtures of these end-member sources, but it is possible that such mixtures could produce higher statistical matches. Sample density surrounding the Arizona shelf is too low to definitively determine source or transport pathway.

The Desmoinesian–Virgilian (312.6–299 Ma) interval contains the greatest number of samples from the ARM system (Fig. 9). Locations proximal to ARM basement uplifts from this time interval show nearly unimodal or bimodal age spectra (e.g., Figs. 5–7) that suggest large amounts of basement exhumation and large volumes of sediment input into proximal ARM basins. Samples collected from locations distal to major basement uplifts contain substantial numbers of “local” Yavapai–Mazatzal grains that are interpreted to have been eroded from ARM basement uplifts; however, these samples also contain broad age spectra that are interpreted to have been sourced from outside southwestern Laurentia. In the western portion of the ARM system, samples distal to ARM uplifts make up a roughly north–south zone in which individual samples are moderately similar to aggregated samples of the Arizona shelf (Fig. 9). Most samples within this zone have PDF cross-correlation values > 0.5 when compared to the Arizona shelf, in contrast with most proximal ARM basin samples (< 0.4). One exception to this is one siltstone sample from northern Arizona (DryCreek; Soreghan et al., 2007) in which only 54 grains were analyzed. The small number of grains likely prevents meaningful statistical comparison of this sample to larger-n samples (Saylor and Sundell, 2016). These samples are also distinct from a sample collected in central Nevada (Battle; Gehrels, 2000) that contains a major peak at 1.86 Ga and a subordinate peak at 2.66 Ga. Similar prominence of these populations is not found in any of the distal ARM basin samples, although grains of these ages are present.

Collectively, these data suggest the presence of an active sediment routing system that ran north to south from southern Canada to southern Arizona (Fig. 9). Previous detrital zircon work on Upper Pennsylvanian strata in southwestern Laurentia have interpreted the sediment that filled distal portions of ARM basins as allochthonous (Gehrels et al., 2011; Link et al., 2014; Beranek et al., 2016; Chapman and Laskowski, 2019), but whether that material transited the northern United States midcontinent from the Appalachian orogen or was derived from Arctic terranes and/or the northern Appalachian orogen was not clearly resolved (Gehrels et al., 2011; Link et al., 2014). In addition to PDF cross-correlation values, a transect of Desmoinesian–Virgilian samples from the Illinois, Forest City, and Wood River Basins and northern Wyoming (Fig. 14) suggests that Desmoinesian sediment transport from the northern midcontinent to the Bighorn and Wood River Basins is unlikely. Spectra from the Illinois and Forest City Basins contain few grains older than 1.5 Ga, in contrast with major peaks between 1.8 Ga and 1.6 Ga within Bighorn and Wood River strata. These data suggest a major Desmoinesian sediment divide between the Forest City and Bighorn Basins (Fig. 14B). This interpretation is also supported by the relative similarity of a Late Pennsylvanian sample from southeastern British Columbia (Fig. 14A) that suggests that the north–south sediment transport pathway likely extended at least this far north. We interpret Arizona shelf spectra to consist of a combination of far-traveled grains and grains eroded from the Zuni–Defiance uplift in northeastern Arizona based on the abundance of Yavapai and Mazatzal grains as well as age peaks centered on ca. 1.5 Ga (Doe et al., 2013). We interpret the presence of sediment divide between the Arizona shelf and the Paradox Basin during this time interval (Fig. 9) possibly associated with the Zuni–Defiance uplift. This interpretation is based on a mismatch between the 1.42 Ga population dominant within the eastern proximal Paradox Basin at this time (and later during Wolfcampian time in the western Paradox Basin) (Fig. 6, panel 2 and 4) and the prominent 1.50 Ga peak present in Arizona shelf samples (Figs. 3, 4, and 13) at this time. This mismatch implies that “A-type granite” zircons shed from the Uncompahgre Uplift did not reach the Arizona shelf and that 1.50 Ga grains within the Arizona shelf were derived from another source. One possibility is recycling of Proterozoic quartzites from the Zuni–Defiance uplift (e.g., Doe et al., 2013).

Statistical comparison of Arizona shelf data from this time interval to potential distal source terranes yields the highest similarity values with the northeast Ellesmerian orogen and/or Sonoran sources, but we cannot rule out the possibility of northern Appalachian sediment being transported across the Canadian Shield before being routed southward into the ARM system via longshore and aeolian processes (Fig. 15).

During Wolfcampian time (299–282 Ma) (Fig. 9), uni- and bimodal detrital zircon spectra from proximal Paradox Basin deposits are interpreted to indicate that sediment was sourced exclusively from the ancestral Uncompahgre uplift to the northeast. These spectra contrast with spectra from coeval Arizona shelf strata to the southwest, which contain broad distributions of grain ages (Figs. 3–4 and 9), and we interpret that a sedimentary divide, most likely topographic, separated these areas (Fig. 9) based on the mismatch between nearly unimodal 1.42 Ga zircon populations in the proximal Paradox Basin and prominent 1.50 Ga populations in Arizona shelf samples (see above). Samples from the Arizona shelf and southern Arizona have high statistical similarities, suggesting that sediment within this region was well-mixed. Similar to the Desmoinesian–Virgilian interval (above), a conduit of moderate statistical similarity extends northward at least as far as southern Idaho (Fig. 9), and we interpret this as the route by which far-traveled sediment entered the Arizona shelf and southern Arizona. Strata sampled in the Denver Basin from this period have greater similarity to the Arizona shelf than underlying strata (Fig. 9). This shift coincides with the transition from fluvial to aeolian deposition as the Denver Basin was overwhelmed by southward progression of a major aeolian system and a coeval change from local-to continental-scale sediment sourcing. This shift also likely coincides with the cessation of exhumation of the ancestral Front Range uplift (Fig. 7). Arizona shelf samples from this interval are also statistically distinct from Lone Pine Basin samples (Fig. 8). We interpret the zircon within the Lone Pine Basin to have been sourced from Ellesmerian or northern Appalachian orogens but to have bypassed the ARM system such that Yavapai–Mazatzal and A-type granitic zircons were not substantially incorporated (Fig. 9). As in the two previous time intervals, aggregated Arizona shelf data are most statistically similar to aggregated data from the northeast Ellesmerian orogen, but we cannot rule out contribution from the northwest Ellesmerian and northern Appalachian orogens (Fig. 15).

Leonardian to early Guadalupian (ca. 282–260 Ma) sedimentation within the ARM system was largely dominated by aeolian systems (McKee and McKee, 1967; Blakey et al., 1988; Lawton et al., 2015) driven by south-directed sediment transport (Parrish and Peterson, 1988; Peterson, 1988). This is reflected in the detrital zircon data (Fig. 9) by a higher degree of statistical similarity across the entire ARM system from northern Wyoming through central Arizona (Fig. 9). We interpret these results (and deposits) as part of a large, well-mixed volume of aeolian sand that moved southward into the region as Pangea drifted northward away from the equator (Scotese et al., 1999; Domeier and Torsvik, 2014) and the climate aridified (Tabor and Poulsen, 2008). We also interpret the switch from uni- and bimodal detrital zircon spectra to broad spectra with diverse grain ages in the Paradox Basin (Figs. 8–9) to suggest the cessation of major exhumation on the ancestral Uncompahgre uplift during Leonardian time, consistent with field evidence for cessation of active uplift by this time (Soreghan et al., 2012).

The interpretations above demand one or more mechanism(s) by which sediment was transported ∼2000 km along the dominantly shallow marine western margin of Laurentia during Pennsylvanian and Permian time (McKee and Crosby, 1975). We argue that coupled longshore current/aeolian transport is the most likely process. Silt-sized grains also could have been transported in suspension and deposited as loess (e.g., Soreghan et al., 2014). Aeolian cross bedding from a Morrowan–Atokan sequence of mixed aeolian and shallow marine strata in southwest Alberta records southwest to south-southeast sediment transport (Stewart and Walker, 1980), demonstrating that southward sediment transport along the northwestern Laurentian margin was occurring by at least Early Pennsylvanian time. Aeolian transport and deposition also became the dominant sedimentary process in Wyoming beginning during the Late Pennsylvanian (Link et al., 2014), and this aeolian system buried the Arizona shelf by Leonardian time (Blakey and Knepp, 1989; Blakey, 1990). Paleowind directions reconstructed from Pennsylvanian–Permian aeolian strata are consistent with global climate models that depict trade winds capable of driving such longshore currents and south-directed aeolian transport (Tabor and Poulsen, 2008).

An informative modern analog for longshore current transport operating at a similar spatial scale to that envisioned for the western Laurentian margin can be found in southwestern Africa, where longshore current transports sediment from the Namibian/South African Orange River ∼1800 km north at least as far as the Angola coast (Garzanti et al., 2014, 2018). This littoral transport system, which has been active since at least Eocene time (Bluck et al., 2007), is driven by onshore (southwest) winds and has been interpreted to be the result of low sediment accommodation along the Namibian and Angolan coasts, aridity, and high sand and gravel sediment supply from the Orange River (Bluck et al., 2007). Provenance signals such as sandstone composition, heavy mineral assemblages, and detrital zircon spectra can be traced throughout the length of this system (Garzanti et al., 2012, 2014, 2018). The Namib, Skeleton Coast, Cunene, and Moçamedes coastal dune fields are discontinuously distributed along the Orange River littoral cell and are separated by up to 450 km; these coastal dune fields are each terminated by rivers to their north that flush sand back to the littoral cell (Garzanti et al., 2018). Most sediment transport in this system occurs in the subtidal to intertidal zones (Bluck et al., 2007). The abrupt end of this littoral transport system occurs where the continental shelf decreases in width to only a few kilometers and sand is flushed into the deep ocean through submarine canyons (Garzanti et al., 2018).

Although the Orange River littoral cell is not a perfect analog for the sediment routing system interpreted along the Pennsylvanian–Permian western margin of Laurentia, major features of that system suggest the two may be comparable. The Orange River littoral cell transports sand from 28.6°S northwest to 15.7°S. The presence of Morrowan–Atokan aeolian and shallow marine sandstones, deposited at ca. 17°N (Domeier and Torsvik, 2014) in southwestern Alberta, suggests that the interpreted sediment routing system operated at similar latitudes as the Orange River system. These deposits, the Tyrwhitt, Storelk, and Tobermory Formations, are sedimentologically similar to Eocene deposits 150 km north of the start of the Orange River littoral cell that are interpreted as early remnants of that routing system (Scott, 1964; Stewart and Walker, 1980; Bluck et al., 2007). Much of the Pennsylvanian sedimentary record has been eroded from the Canadian Laurentian margin (Richards et al., 1994; Henderson et al., 1994), and much of what is preserved consists of coastal dune to outer shelf deposits (Richards et al., 1994). This implies that the intertidal, subtidal, and coastal dune systems, the primary sediment conduit for the interpreted routing system, may have been removed by erosion. Abundant unconformities within these deposits (e.g., Scott, 1964; Stewart and Walker, 1980; Henderson et al., 1994) suggest relatively low sediment accommodation with high sediment bypass. Although the climate from southern Canada to the ARM province was not as arid as the Namib and Angolan Coasts (Tabor and Poulsen, 2008), the presence of large aeolian units such as the Storelk and Weber Formations in southern Canada, partly aeolian Quadrant Formation in southwestern Montana, and an unnamed aeolian sandstone in the Eagle Basin (Schenk, 1987), suggests at least periodic aridity in the northern portion of this system during the Early to Middle Pennsylvanian (Fryberger, 1979; Stewart and Walker, 1980; Schenk, 1987), and sparse data from most of the area in which this system is interpreted suggest arid to semi-arid paleoclimate (Tabor and Poulsen, 2008).

The data and interpretations presented here are consistent with the long-standing idea that the Transcontinental Arch, extending from northern Arizona and New Mexico to the Lake Superior region of Canada (Keith, 1928; Carlson, 1999; Amato and Mack, 2012), influenced sediment transport networks during the Pennsylvanian and Permian (Fig. 15). Detrital zircon data suggest that little if any sediment was routed directly from the central Appalachian orogen across the northern U.S. midcontinent into the ARM system (Figs. 9, 14, and 15) during Pennsylvanian to middle Permian time. However, the current study cannot constrain how far north this barrier extended, and our data allow the possibility that sediment shed from the northern Appalachian orogen was transported north across the Canadian Shield to join the southward longshore/aeolian routing system interpreted here. The data presented here also cannot rule out the possibility of a Mississippian–earliest Pennsylvanian sediment pathway between the central Appalachian orogen and southwest Laurentia (Gehrels et al., 2011; Chapman and Laskowski, 2019). Mississippian Surprise Canyon Formation (Fig. 3) and Early Pennsylvanian Molas Formation (Fig. 6, panel 4) samples appear to show central Appalachian affinity, which suggests that the Transcontinental Arch may not have played a substantial role in Mississippian sediment routing (e.g., Chapman and Laskowski, 2019).

1. Detrital zircon deposited within Pennsylvanian–Permian Ancestral Rocky Mountain basins is a mix of grains eroded from southwestern Laurentian basement rocks exposed in ARM uplifts and grains from a wide variety of Laurentian sources transported hundreds to thousands of kilometers into the ARM system.

2. Similarity metrics suggest a north-to-south (modern reference frame) sediment routing system that moved detritus from at least as far north as southern Canada into the ARM system and filled basins as far south as southern Arizona. This system would have been active from at least Desmoinesian to Leonardian time. We propose that sediment was moved via a ∼2000-km-long longshore and aeolian transport system analogous to the modern Namibian coast (e.g., Garzanti et al., 2018).

3. Distal basin detrital zircon spectra show the closest statistical similarity to Arctic Ellesemerian orogen sources, although distal basin spectra could also be the product of combined northern Appalachian and locally derived grains.

4. Detrital zircon spectra suggest initiation of ARM exhumation adjacent to the Denver and Paradox Basins by Atokan (Early Pennsylvanian) and Desmoinesian (Middle Pennsylvanian), respectively, as well as local exhumation of the ancestral Front Range and Uncompahgre uplifts by Desmoinesian time. Zircon ages from one sample collected from the Central Colorado trough allow for the possibility of ARM exhumation as early as Late Mississippian. Aeolian systems prograded into the Denver and Paradox Basins during the Wolfcampian (early Permian) and Leonardian (middle Permian), respectively, overwhelming local sediment sources and suggesting the cessation of ARM-related exhumation.

This research was funded primarily through the generous support of Pioneer Natural Resources. Research at the University of Houston was funded through NSF EAR-1742952. Special thanks go to the Arizona LaserChron staff for help with detrital zircon analyses. The Arizona LaserChron Center is supported in part by NSF grant EAR-1338583. We thank Dustin Sweet and Alan Chapman for constructive reviews and Lithosphere Science Editor Sarah Roeske for additional guidance that helped us to improve this manuscript.