Large detrital zircon data set investigation and provenance mapping: Local versus regional and continental sediment sources before, during, and after Ancestral Rocky Mountain deformation

The exponential expansion of detrital zircon U-Pb data sets (Gehrels, 2014; Sundell et al., 2020) have spurred a revolution in sediment provenance analysis by creating new opportunities and challenges in data synthesis and integration (e.g., Pullen et al., 2014; Sharman and Johnstone, 2017; Chapman and Laskowski, 2019; Leary et al., 2020; Sundell et al., 2020). World-wide databases of detrital zircon ages (Puetz and Condie, 2019) enable comparison of measured samples to potential sources across the globe. However, this raises the challenge of identifying the most likely sources among an ever-expanding set of potential sources (e.g., Hervé et al., 2003; Morón et al., 2019; Clark et al., 2020; Lawton et al., 2021) particularly during periods of supercontinent assembly. The common prospect of multiple potential sediment sources containing the same age modes presents the challenge of quantifying the contributions from those sources particularly when using individual age groups to infer source contribution. Finally, the increasing number of detrital zircon samples requires integrating information about age distribution, source identification, and source contributions in a spatial context along with complimentary geological data. In consideration of these challenges, this contribution develops new data handling and interrogation methods to better synthesize, integrate, and interpret these data.

We present a novel, integrative method of visualizing the modeled data, in concert with other geological data, to reconstruct Carboniferous–Triassic changes in sediment routing patterns in the Laurentian interior. Our approach here employs two different methods of mixture modeling to (1) develop models of sediment sources’ age distributions, (2) identify the sources that best match those models, and (3) determine the proportion of detrital zircons from each source. We then contour the proportions with the aid of contextual geologic data, creating sediment provenance maps to identify the changes in sediment routing in and around the Ancestral Rocky Mountains (ARM). The utility of this approach lies in the leveraging of large data sets to analyze and interpret age distributions rather than solely focusing on individual age groups, and visualization of those data integrated with other contextual information (e.g., paleocurrent data, stratal thickness, etc.).

The upper Paleozoic–lower Mesozoic stratigraphic record in western Laurentia, which embodies all the challenges described above, is an apt interval to evaluate our modeling and visualization approach. It illustrates the interaction between continent-scale sediment routing systems and areas of exclusive local sourcing driven by intraplate deformation. Intraplate deformation during the late Paleozoic covered a broad region of western Laurentia (a.k.a. the ARM; Fig. 1), isolating some ARM basins from distally derived sediment, and creating potential impediments to transcontinental sediment transport (Leary et al., 2020; Lawton et al., 2021). Nevertheless, Appalachian detritus derived from the eastern margin of Laurentia persisted in basins to the west of the ARM during and after ARM deformation (Dickinson and Gehrels, 2008). This apparent persistence, in part, may be due to rerouting of sediment pathways around the growing ARM. Alternatively, although not mutually exclusive, Appalachian-like detrital zircon age distributions may persist through local recycling of older strata with that signature, even if direct connection to the Appalachian orogen had been cut off. These scenarios can be discriminated by documenting sediment routing corridors before, during, and after ARM deformation and basin isolation. To do so, we endeavor to better visualize and understand the pre-, syn-, and post-ARM tectonically induced isolation and characterize distal sediment pathways through large data set sediment source modeling and mapping of western Laurentia. Other basins where an interplay between locally and distally derived sediment occurs in time (e.g., the Amadeus Basin, Australia; Maidment et al., 2007) or space (e.g., Andean foreland of Argentina; Capaldi et al., 2017) will serve as testing grounds for the efficacy of this approach.

The late Paleozoic (Late Mississippian–Permian) saw the amalgamation of the most recent supercontinent, Pangea. During this period Laurentia was marked by tectonic activity along its eastern, southern, and western margins (modern geographic reference frame), as well as basement-involved deformation within (Fig. 1). Along the eastern margin, the Alleghenian Orogeny (Pennsylvanian–Permian) resulted from full continent-continent collision (Hatcher, 2010) creating a Himalayan-style mountain belt. The Ouachita-Marathon orogen resulted from a series of crustal blocks that collided at different angles and varying degrees of intensity with the uneven southern margin of Laurentia to form a mountain belt that stretched from modern-day Arkansas to west Texas, USA (Royden, 1993; Keller and Hatcher, 1999). The western and southwestern Laurentian margins underwent a complex combination of collision, back arc spreading, and sinistral translation, which has been further obscured by subsequent Mesozoic and Cenozoic tectonic activity (Trexler et al., 2004; Lawton et al., 2017; Sturmer et al., 2018).

Collision along the continental margins in the Pennsylvanian–early Permian was accompanied by the uplift of a series of basement-cored blocks and development of adjacent basins (i.e., the ARM) that stretched between present-day Oklahoma and Utah, USA (Fig. 1) (Kluth and Coney, 1981; Dickinson and Lawton, 2003; Soreghan et al., 2012; Blakey, 2019; Miall, 2019;). This swath of predominantly contractional deformation initiated primarily in the Early Pennsylvanian (Leary et al., 2017; Sweet et al., 2021), occurred in Proterozoic crust (Karlstrom and Humphreys, 1998; Leary et al., 2020), and reactivated many pre-existing Proterozoic and Cambrian rift structures (Marshak et al., 2000; Craddock et al., 2017). ARM uplifts were thrust onto adjacent crust (Frahme and Vaughn, 1983; Brewer et al., 1983; Hoy and Ridgeway, 2002) creating flexural basins that accommodated detritus shed from the bounding uplifts (Barbeau, 2003; Soreghan et al., 2012). The core of the ARM straddled the trace of the transcontinental arch (Carlson, 1999; Leary et al., 2020), which was a continental-scale, long-lived, northeast-southwest trending feature (Fig. 1) that controlled deposition in the midcontinent and western Laurentia throughout the early–middle Paleozoic (Sloss, 1988; Linde et al., 2014; Linde et al., 2017).

Upper Mississippian fluvial rocks of the Grand Canyon area (a.k.a. the Arizona shelf) record the arrival of a characteristic detrital zircon age distribution interpreted to mark the sediment routing connection between the Appalachian Mountains and the western shores of Laurentia (Gehrels et al., 2011; Gehrels and Pecha, 2014; Attia et al., 2018; Chapman and Laskowski, 2019). This age distribution is defined by a positively skewed “Grenville” mode (0.95–1.3 Ga, peak at ca. 1.1 Ga), and subordinate Paleozoic–late Neoproterozoic modes (Gehrels et al., 2011; Chapman and Laskowski, 2019; Kuiper and Hepburn, 2020). The arrival and subsequent persistence of these age modes suggest that during the rest of the Paleozoic, a relatively continuous supply of eastern-derived sediment was routed to the western Laurentian margin. Prior to the introduction of this eastern-sourced age distribution, detrital zircon ages in older units were predominantly late Paleoproterozoic and early Mesoproterozoic (Fig. 2), reflecting first-cycle or recycled western Laurentian crustal detritus (Gehrels et al., 2011; Chapman and Laskowski, 2019). Appalachian or Appalachian-like detrital zircon age distributions, potentially contributed by other sources (e.g., Ouachita-Marathon sourced sediment), persisted in the successive sedimentary record of western Laurentia throughout the late Paleozoic–early Mesozoic and beyond (Dickinson and Gehrels, 2008; Gehrels et al., 2011; Liu and Stockli, 2019; Gao et al., 2020; Thomas et al., 2019).

Research indicates that the transcontinental sediment transport system began to grow from east to west across Laurentia in the Mississippian and was aided by shallow marine processes (e.g., tides, longshore currents, storms; Chapman and Laskowski, 2019). Other potential sediment sources from along the Laurentian margin bearing varying degrees of similarity to the Appalachian detrital zircon age distribution include the Ouachita-Marathon front (Thomas et al., 2019), the Arctic Ellesemerian mountains (Beranek et al., 2010; Anfinson et al., 2012; Leary et al., 2020), and the remnant Antler uplifts in Idaho and Nevada, USA (Linde et al., 2016; Beranek et al., 2016). Despite the influx and mixing of sediment from around Laurentia, many ARM basins, especially in the core of the ARM (Fig. 1), contain large regions of dominantly locally derived sediment (Leary et al., 2020). The temporal overlap of areas dominated by local versus distal provenance raises questions of what controlled the geographic distribution of sediment source through time. Further, what was the extent of ARM basin sediment source isolation and how were distally derived sediments delivered to the western margin of Laurentia during this time?

ARM deformation ceased in the early Permian, but post-tectonic infilling of some ARM basins continued into the middle and late Permian (Dickinson and Lawton, 2003), and remnant uplifts affected sediment routing through this area into the Mesozoic (Dickinson et al., 2010). By Permian time, increased aridity resulted in pervasive eolian sediment transport across western Laurentia (Leary et al., 2020; Lawton et al., 2021), and during the late Permian, much of this region experienced nondeposition and/or erosion (McKee and Oriel, 1967). Relatively brief Early–Middle Triassic deposition in marine and nonmarine settings (e.g., Moenkopi Formation) occurred in the flexural basin of the Sonoma Orogen along western Laurentia (Lawton, 1994; Dickinson and Gehrels, 2008), and was followed by another period of nondeposition and/or erosion. Arc collision along Laurentia’s western margin in the Late Triassic created back-arc rift-related accommodation, preserved transcontinental fluvial deposits, and set the tectono-stratigraphic stage for the remainder of the Mesozoic (Dickinson and Gehrels, 2008; Dickinson, 2018). Plutonic rocks exposed in southern California (USA) record development of the southern magmatic arc from late Permian to Late Triassic along the truncated southwest margin of Laurentia (Riggs et al., 2013; Cecil et al., 2018). In the Late Triassic, arc material from this southern margin was transported northeast into a northwest-directed fluvial system (Dickinson and Gehrels, 2008; Riggs et al., 2016). The headwaters of this northwest-flowing fluvial trunk resided in the remnant Ouachita-Marathon mountains and their exhumed foreland basins, which delivered sediment from south Texas to western Arizona and Utah (Dickinson and Gehrels, 2008; Dickinson, 2018).

We contribute 20 U-Pb detrital zircon samples to the existing body of published ages. New samples were collected to fill out a paucity of data within the core basins of the ARM, which include the Paradox and Eagle basins, and the Central Colorado and Taos troughs. Descriptions of samples, U-Pb age data for new samples, and ages and uncertainties for compiled data are provided in the supplemental data files (Supplemental Materials 1, 3, and 4)¹.

New samples presented here were prepared at the University of Houston, Texas. Rock samples (~2 kg) were crushed and disc milled down to 400 µm disc spacing. Mineral separation was performed using standard separation techniques: density separation via a water table followed by immersion in methylene iodide (3.28 g/cm³), and magnetic separation on a Frantz Magnetic Separator. We used an Olympus SZX12 microscope to separate zircons from non-zircon grains after rinsing grains in nitric acid. Zircons were then non-discriminately mounted on double-sided tape.

U-Pb age of zircons were determined via ablation by Photon Machine Excite 193 nm ArF laser systems coupled to a Varian 810 quadrupole inductively coupled plasma–mass spectrometer (ICP-MS) at the University of Houston (Shaulis et al., 2010). Individual zircon grains were lased with 240–200 shots at 10 Hz repetition rate and 30–25 µm spot size. We corrected for inter- and intra-element fractionation using the Bohemian Massif potassic granulite Plešovice (337.13 ± 0.37 (2σ, two standard deviations) Ma) (Sláma et al., 2008) as the primary zircon standard. The accuracy of the fractionation correction was verified by evaluation of FC5z from the Duluth Complex (1099.1 ± 0.5 (2σ) Ma) (Paces and Miller, 1993) as the secondary standard. Analyses of every 10 unknowns were separated by an analysis of the primary standard. In turn, every third analysis of the primary standard was coupled with a secondary standard analysis.

Ages were determined with U-Pb Toolbox, a MATLAB-based application, which calculates integrated isotopic ratios from raw counts per second that are exported from Quantum software, corrects for machine bias and fractionation, and filters the databased on user-defined parameters (Sundell, 2017). Zircons were filtered using a 15% uncertainty cutoff, and 20% and −10% discordance cutoffs (²⁰⁶Pb/²³⁸U–²⁰⁷Pb/²³⁵U ratio, the latter determined via measured ²⁰⁷Pb/²⁰⁶Pb and the accepted ²³⁸U/²³⁵U, 137.82) applied with a 900 Ma discordance transition. Instead of a common Pb correction (Stacey and Kramers, 1975), we accept grains whose 2σ envelope is <15% discordant by comparison of the ²⁰⁶Pb/²³⁸U and ²⁰⁷Pb/²³⁵U ages for ages <600 Ma.

Sediment source modeling and comparative statistical metrics provide quantitative tools to distinguish distinct source signals and identify similarities amongst age distributions in the ever-increasing volume of provenance data (Amidon et al., 2005; Vermeesch, 2012, 2013; Saylor and Sundell, 2016; Sharman and Johnstone, 2017). The approach adopted here merges the strengths of bottom-up and top-down modeling (Fig. 3) to avoid lengthy and often cumbersome descriptions of individual age modes for large data sets. It does so by characterizing possible source distributions, and the proportions of sources assigned to individual samples. Both bottom-up and top-down models include goodness of fit metrics that assess how well generated model distributions match the sample distributions that they are trying to recreate. Top-down and bottom-up sediment unmixing models are defined by whether known sources are used as an input, and how the algorithm determines source attribution to sink samples (i.e., samples from basin(s), which are being analyzed) (Sharman and Johnstone, 2017).

How do we determine what the actual detrital zircon sources looked like, especially in ancient basins? Bottom-up modeling answers this question and does so with no known source input. Rather, this method aims to determine unknown sources that contributed to the collection of sink data distributions (Sharman and Johnstone, 2017). Bottom-up modeling used here employs a non-negative matrix factorization (NMF) algorithm to deconvolve two matrices that define (1) a series of synthetic sources (a.k.a. factorized sources) and (2) corresponding weighting functions both from an input matrix composed of sample’s age distributions from a basin or basins (i.e., sink data set) (Sharman and Johnstone, 2017; Saylor et al., 2019). The deconvolved matrices are iteratively updated with the goal to minimize the error (e.g., final residual) between model and sample age distributions for a given number of factorized sources (Saylor et al., 2019). To identify the optimal number of factorized sources, we considered three ranks (i.e., factorized source sets) centered on the rank above which there is little or no decrease in the final residual (Fig. 4A) (Saylor et al., 2019). To discriminate among these three factorized source sets, we use geologic information (e.g., knowledge of realistic sediment sources within the study area, previously proposed sediment sources) to determine which set of factorized sources most accurately depicts the dominant detrital zircon sources in the study area. Thereby using the bottom-up approach to identify the empirical detrital zircon sources to be selected for top-down mixing.

Top-down models mix proportions of detrital zircon sources from sample data (i.e., empirical sources) to create a mixture distribution that most closely matches that of a given sink sample (e.g., Amidon et al., 2005; Saylor et al., 2013; Kimbrough et al., 2015; Licht et al., 2016; Mason et al., 2017). Here, we mix empirical source age distributions via the Monte Carlo approach of DZmix (Sundell and Saylor, 2017) to create a model age distribution from a random mixture of these inputs, and then compare that modeled distribution to a detrital zircon sample age distribution from our sink data set. This process is repeated many times (e.g., 10,000) for each sample in the sink data set, and the best matches are retained. If the input sources closely reflect the actual detrital sources, then the proportions of input sources assigned to the best fit models for each detrital zircon sample faithfully approximate actual detrital source proportions.

In the first step, we used DZnmf (Saylor et al., 2019) to deconvolve (a.k.a. factorize) a set of 191 Devonian–Triassic Laurentian detrital zircon samples and sample groups (i.e., model inputs) into a set of potential detrital zircon sources. We use the optimal number of sources suggested by the DZnmf (Fig. 4A) to guide our selection of source sets within the context of our geologic understanding of sediment routing in western Laurentia during this time (the process described below). The data used in modeling span North America and were distilled from 329 individual samples by grouping geologically and geographically adjacent samples. This curation of the data better characterizes age distributions by combining samples, which thereby increases sample size (n). To address the potential mismatched lumping of samples into sample groups, we also performed this exercise without any lumping (i.e., running 329 individual samples through DZnmf; see Supplemental Material 8 for results). Minor differences are noted in the results below, but none that we found to affect the overall results or interpretations of this study, especially with respect to empirical source identification. All calculations in the following sections are based on kernel density estimates (KDEs) with a 20 m.y. bandwidth. We use KDEs rather than PDPs, which use analytical uncertainty for bandwidth, to mitigate mismatches of analytical precision associated with the range of analysis types (e.g., laser ablation ICP-MS versus secondary ion mass spectrometry) (Vermeesch, 2012).

In the second step, we search for and identify empirical sources that match the factorized sources produced in Step 1 with DZnmf. Empirical source matches were selected based on their similarity to factorized age distributions (Fig. 5; Table 1). Here, we used coefficient of determination (R²) from cross-correlation to determine best empirical source match. The searching and identification of potential sources is often iterative (Fig. 3) because initial sources identified may yield low similarity, sending the researcher back into the literature to find a better match. We discuss the uncertainty and relevance of our detrital zircon source selection in the context of possible late Paleozoic and Triassic sediment sources.

We interpret some of the empirical sources to directly reflect the location of detrital zircon sources. However, other factorized sources are non-unique, and share similarities to multiple empirical sources. Nevertheless, even in these cases, a broad distinction can be drawn between sediment derived from basement sources adjacent to the ARM basins and sediment transported from terranes on the margins of Laurentia. This level of detail works well for the types of sediment sources that existed during the study interval, and discrimination of ARM-driven sediment routing and basin isolation. In cases of applying this approach where a broad distinction is not sufficient in addressing provenance questions, further data would need to be brought to bear.

In the third step, we complete the sediment source modeling loop by top-down modeling the sink sample data set that was used in bottom-up modeling (i.e., Step 1). Top-down modeling also uses the empirical source set developed in Step 2. Utilizing both of these data sets, we model the proportions of empirical sources for each sink sample via the MATLAB-based DZmix algorithm (Sundell and Saylor, 2017). However, we employ an edited version of this code that allow analysis of all the detrital zircon samples at once, and which is included in Supplemental Material 9.

Parameters for the model were set at 10,000 runs and we retained the best 0.1% fits to create mean model distributions. We compare bottom-up and top-down model results to assess consistency across model predictions (Fig. 6; Supplemental Material 1). DZmix produces best model fits for R², Kuiper V statistic, and Kolmogorov-Smirnov D statistic. We use the R² model results and compare these results to Kuiper V statistic to check consistency. When interpreting model results, we focus on the major contributors to age distributions, and discount sources whose contributions are within uncertainty of zero.

Visualization of provenance evolution in both time and space is a powerful approach to understanding and interpreting source-to-sink data. In this study, sediment sources are broadly defined by the zircon age distributions that they contribute. This research focuses on the degree of ARM-driven provenance isolation within a background of continental-scale sediment mixing. We use the character of U-Pb zircon age distributions found in the stratigraphic record to discern these local and distal sediment source signals. Therefore, sediment sources identified in bottom-up modeling, and prescribed in top-down modeling are lumped into two categories: (1) Local contributors of sediment, which constitute many of the basement-cored ARM uplifts, and (2) contributors of sediment distal from the core of the ARM, which tend to reside along the Laurentian margin (Figs. 5 and 7). The recipients of sediments from these two sources (i.e., detrital zircon samples) are visualized primarily based on the percentage of source attribution. This approach to grouping samples is designed to fit our focus of ARM-driven provenance isolation, and is not inherent to the data or the approach. This focus also guides the types of auxiliary data used to aid in contouring provenance maps.

Local sources, as defined here, contribute sediments that are characterized by unimodal age distributions that reflect Laurentian basement exposed in ARM uplifts, or recycled from similarly sourced lower Paleozoic strata (Fig. 7). On occasion, further distinction in the local source attribution is made by considering the composition of the sandstone, and where available, published descriptions of sandstone compositions are used to interpret direct from basement versus recycled sediment source (e.g., arkosic versus quartz-rich). The distal sediment source group is predominantly attributed to empirical source matches from the Laurentian margin and associated accreted terranes (Fig. 7). Most of these sediment source signals are characterized by multi-modal age distributions, and although some exhibit simple age distributions, they are not commonly associated with basement exhumed by ARM deformation (e.g., Archean; Leary et al., 2020). This binary categorization is intended to collapse some of the ambiguity in identifying specific sediment sources. Whereas most sediment sources can be confidently attributed to either local basement or distal margin, the specific location of sources for some samples may remain somewhat ambiguous.

Sediment source proportions from DZmix results (Supplemental Material 6) were contoured with the aid of stratigraphic thickness maps, paleocurrent data, and paleogeographic information (Supplemental Material 1) to create sediment source maps (Fig. 8). Paleocurrent data (Supplemental Material 7) were used to help infer direction of sediment routing. These data are plotted on provenance maps and colored based on the general mechanism of sediment transport (i.e., eolian versus fluvial/alluvial or deltaic/marine). We contributed 18 paleocurrent measurement stations to 201 published paleocurrents. New paleocurrent stations each contain 5–25 individual measurements (n = 288), and are comprised of trough crossbedding, planar crossbedding, and imbricated clast measurements.

We use a large contour interval (20%), which is intended to compensate for the degree of variability within the sediment source episode and the uncertainty of top-down model results. The few instances where closely spaced samples exhibit variability greater than the contour interval commonly reflect temporal changes in sourcing that are beyond the resolution of the selected time interval (i.e., sediment source episode, which is defined below). For example, in Atokan-Wolfcampian Episode 3, one of the several samples in the middle Paradox Basin yields a 73% local (i.e., basement or recycled basement) source contribution, and is surrounded by samples with ≥80% local contribution (Fig. 8C). We italicize this value in the source map to indicate the discrepancy and contour with consideration of all data points in the area (i.e., place the 80% local sediment source contour adjacent this point in the data cluster). In an absolute sense, this difference is not significant. However, it is important to understand why these discrepancies exist. We divide the study interval (Carboniferous–Triassic) into five episodes (Fig. 2A) and create sediment provenance maps for each, the details of which are covered in the Discussion section of this paper. The boundaries between episodes, defined here by trends observed in the provenance record, coincide with North American stages. However, the sedimentation regimes that define these episodes reflect tectono-sedimentary processes that do not necessarily turn off and on at stage boundaries. Therefore, episode boundaries are somewhat imprecise (Fig. 2A). Many formations from which samples were collected have age designations that span stage and episode boundaries. In these cases, where no further stratigraphic information may be brought to bear, we either assume a midpoint age or assign the sample to the episode containing other nearby samples that exhibit similar characteristics.

As sample numbers increase, it is increasingly cumbersome to present, describe, and discuss age distributions for each sample. The bottom-up modeling described above addresses this and allows us to focus on salient endmember features of the 191 samples and sample groups included in this study. Therefore, in the following section we discuss source distributions and their attribution to detrital zircon samples from basins. We also present both top-down model results and KDEs for selected detrital zircon data (Fig. 2C) to provide a basis for comparison between the more granular approach of visualizing each sample’s age distribution and the higher-level data analysis of the modeling and mapping approach proposed here (Figs. 3 and 8).

Bottom-up modeling with DZnmf indicates an optimal number of sources is 10 (i.e., low rank minimum final residual or breakpoint; Saylor et al., 2019); however, factorization of 9 and 11 source models yield similarly low residuals and are considered here (Fig. 4A). Mean and medians of R² values for factorization models and sample distributions do increase from 9 to 11 (Fig. 4B), but factorization 11 identifies an additional Laurentian crustal age (i.e., 1.74 Ga, see Step 2 below for further description) that we interpret to be an important component of ARM basin sediment supply. We therefore use results from the 11-source factorization, which yield generally high correlation and low residuals (i.e., the difference between mixed models and corresponding sink samples) for factorized source mixed models and sample age distributions (Fig. 4). The 11-source factorization models compared to sink sample age distributions yield mean and median R² of 0.90 ± 0.08 (1σ, one standard deviation) (Fig. 4B) and 0.91, respectively, and Kuiper V statistics of 0.07 ± 0.04 (1σ) and 0.07, respectively.

To determine what, if any, effect our sample groupings have on the factorization of sources, we also used a bottom-up modeling approach on the 329 individual samples (Supplemental Material 8), and compared the results with the results from the bottom-up model using the 191 grouped samples (Supplemental Material 5). The most notable differences in the 329-sample model are the absence of a ca. 2.7 Ga age distribution, the presence of a 1.085 Ga unimodal source, and a corresponding 1.085 Ga nadir in an age distribution that is otherwise consistent with a “Grenville” source (Fig. 5B). The splitting of the “Grenville” age mode does not appear to lend any insight into better empirical matches. Furthermore, nadirs in age distributions that correspond to modes in other distributions are common artifacts in NMF detrital zircon modeling (Leary et al., 2020), and we therefore discount the 1.085 Ga nadir as an artifact. This reduces the differences between the models using the 329 individual samples and the 191 grouped samples from three to two: (1) elimination of the Archean source and (2) addition of the 1.085 Ga source. The Archean source is a relatively minor contributor in comparison to the other empirical sources. Its elimination would most likely force top-down models to attribute other empirical sources containing approximately similar age distributions while yielding poorer model fits (as observed in the 329 model: Supplemental Material 8). The 1.085 Ga source is indeed a minor, but unique contributor to (mostly) Pennsylvanian–Permian sink samples in western Laurentia that has a realistic empirical source match. We experimented with including this empirical source in top-down models and provenance maps, but ultimately decided against it. Details of these models are covered in the Discussion: Consideration of the Methods section as well as Supplemental Materials 1 and 6.

We compare the 11 factorized sources from the bottom-up modeling approach to several empirical sources with the aim to produce an empirical source set. The following sections describe and compare the factorized sources and empirical sources (Table 1; Fig. 5). Our goal is to find empirical sources that closely match factorized ones, and then to use those empirical sources in top-down modeling (i.e., Step 3). The empirical source set that we propose here and use in subsequent steps exhibits very high-moderate correlation to their factorized counterparts with a mean R² of 0.84 ± 0.13 (1σ) and a median of 0.88. It is important to note that we do not intend for the exact locations of the samples used as empirical sources, which are matched to factorized sources, to be the only locations from where those sources are attributed. Rather, the empirical sources are representative of detrital zircon source regions to varying degrees of geographic specificity, which is discussed for each source below. For, example, factorized source 10 is matched to our 530 Ma basement empirical source, which is perhaps the most spatially limited source proposed here. The sample we use for that empirical source is from Cambrian igneous rocks exposed in Colorado, but is used in the proposed source set to represent detrital zircons from several Cambrian igneous rock exposures across the midcontinent and western Laurentia. Many of the factorized sources discussed here are like factorized sources presented in Leary et al. (2020), which we present a brief comparison to in the discussion below.

Triassic arc. Factorized source 9 exhibits a dominant unimodal age distribution at ca. 230 Ma and exhibits a high correlation (R²: 0.88) to the zircon age distribution from volcaniclastic gravels in Triassic sedimentary rocks of southwestern North America (Riggs et al., 2016) (Fig. 5A). We interpret attribution of this source to represent volcanic material shed from the Mesozoic Cordilleran magmatic arc that stretched along the southwestern margin of Laurentia (Fig. 7) (Riggs et al., 2012). There is a minor age population at ca. 1.46 Ga in the factorized source that is not in the empirical source and therefore not incorporated into top-down models. Both top-down and bottom-up models identify this source only in Triassic basin samples.

Central Appalachian. Factorized source 2 exhibits a positively skewed broad age mode (a.k.a. Grenville age mode; 0.95–1.3 Ga) centered at ca. 1.1 Ga and is highly correlated (R²: 0.76) to the central Appalachian empirical source (Fig. 5B). In most cases, we interpret attribution of this source to indicate detrital zircons originally shed from the growing Appalachian Mountain belt in the late Paleozoic. We also compared a few other possible empirical source matches to factorized source 2, but all other comparisons yielded poorer correlations (Table 1). These other samples are representative of accreted peri-Gondwanan terranes along southern Laurentia (i.e., Coahuila and Sabine terranes; Fig. 7). Despite their comparatively lower correlation to factorized source 2, we acknowledge the possibility that peri-Gondwana terranes may have contributed detrital zircons bearing a similar age distribution (Thomas et al., 2019, 2021), especially in strata deposited along the southern margin of Laurentia. Due to this source area ambiguity, we use regional paleocurrent data and additional age modes to help guide our interpretation (see discussion of factorized sources 3 and 7 below).

Maya Block. Factorized sources 3 and 7 predominantly exhibit early Paleozoic and late Neoproterozoic age distributions, respectively (Figs. 5C and 5D). We broadly interpret factorized sources 3 and 7 as representing detrital zircons shed from terranes that collided along the southern and/or eastern Laurentian margin. We select empirical source matches that yield high (R²: 0.87) and moderate (R²: 0.58) correlations to the upper Paleozoic Macal and Santa Rosa formations, respectively, which are exposed in present-day Central America (Fig. 7). These strata are interpreted as a flysch sequence deposited along the northern edge of Gondwana during a complex collision between other peri-Gondwanan terranes, Gondwana, and Laurentia during the late Paleozoic (Weber et al., 2006, 2008). This empirical source match is consistent with recent provenance interpretation of similar age distributions in the Marathon and Permian basins that used both detrital zircon U-Pb and Hf isotope compositions (εHf(t); Thomas et al., 2019). However, similar late Neoproterozoic and early Paleozoic detrital zircon age populations are also observed in mid-continent basins, and they are attributed to a northern Appalachian source (Kissock et al., 2018; Thomas et al., 2020; Leary et al., 2020). The primary cause of the conflicting source attribution for these detrital zircon ages is the similar Neoproterozoic–early Paleozoic geologic and magmatic histories of many peri-Gondwanan terranes (Murphy et al., 2004; Pollock et al., 2012). Despite provenance attribution ambiguity, these source ages may provide a tool in discriminating between eastern and southern Laurentian margin provenance in future research. Detailed investigation of these late Neoproterozoic and early Paleozoic detrital zircon age populations that employ auxiliary provenance proxies (e.g., detrital zircon Hf isotope compositions in Thomas et al., 2019) may be key in differentiating sediment sources along the Laurentian margin. In the absence of a comprehensive auxiliary data set such as detrital zircon εHf(t), we use sediment provenance maps to make these interpretations (Fig. 8).

Antler recycled. Factorized source 5 exhibits a multimodal age distribution with a broad prominent positively skewed age mode at 1.8 Ga spanning much of the Paleoproterozoic. It also contains a minor Archean age mode (Fig. 5E). Factorized source 5 exhibits a very high correlation (R²: 0.97) to early–middle Paleozoic sedimentary rocks exposed in thrust sheets of the Antler Orogen in southern Idaho (Beranek et al., 2016). These sedimentary rocks are interpreted to have a northerly source (e.g., Peace River Arch of northwestern Alberta and northeastern British Columbia, Canada; Gehrels and Pecha, 2014; Beranek et al., 2016). We interpret the presence of this source in upper Paleozoic rocks to represent recycled material shed from the intermittently collisional western margin of Laurentia (Trexler et al., 2004; Sturmer et al., 2018), particularly the northern segment.

Early Paleozoic recycled. Factorized source 11 exhibits a dominant Archean mode at ca. 2.71 Ga, contains minor Proterozoic age mode, and is highly correlated (R²: 0.88) to Ordovician sandstones deposited from Michigan to Arkansas (USA) (Thomas et al., 2016) (Fig. 5F). This source is interpreted to primarily represent recycling of early Paleozoic strata. However, it could also be derived from the cratonic shield to the north; however, it is mostly attributed to quartz-rich and/or feldspar-poor strata (e.g., Vanoss Conglomerate; Thomas et al., 2016; Bushberg Sandstone; Chapman and Laskowski, 2019), which supports a recycled interpretation.

530 Ma basement. Factorized source 10 is a ca. 530 Ma unimode that is highly correlated (R²: 0.92) to Early Cambrian age igneous rocks (Fig. 5G) intruded into Laurentian crust during late-stage Rodinian rifting (McMillan and McLemore, 2004). These igneous rocks are best known for their exposure in southern Oklahoma as part of the exhumed ARM uplift. Other sites of Cambrian intrusion host igneous rocks of similar age occur across Colorado and New Mexico (McMillan and McLemore, 2004). We present zircon data from a Cambrian age syenite exposed in the Wet Mountains of Colorado and use this source to represent similar age rock exposed in basement or possibly sedimentary rocks sourced by that same basement and recycled from exhumed ARM uplifts.

1.43 Ga basement. Factorized source 8 contains a dominant ca. 1.43 Ga age mode and a minor ca. 1.73 Ga age mode (Fig. 5H), and exhibits a high correlation (R²: 0.94) to granitic rock intruded into the Laurentian midcontinent (Bickford et al., 2015). These Mesoproterozoic igneous rocks are commonly referred to as anorogenic granites, but mounting evidence indicates that igneous emplacement during this time was associated with a protracted convergent margin along eastern Laurentia (Whitmeyer and Karlstrom, 2007; Bickford et al., 2015). This source is interpreted to reflect predominantly Mesoproterozoic granitoids that were intruded into Yavapai and Mazatzal basement and exhumed in western Laurentia during ARM deformation. Alternatively, attribution of this source could represent recycled detritus from earlier Paleozoic rocks that received sediment from these basement terranes.

1.66 Ga basement. Factorized source 1 exhibits a dominant ca. 1.66 Ga mode with a shoulder at ca. 1.53 Ga and a broad minor mode centered at ca. 1.33 Ga (Fig. 5I). The shoulder and minor mode are separated by a 1.43 Ga nadir that corresponds to factorized source 8 (1.43 Ga basement source described above; Fig. 5H). The Kingston Mining District gneiss exposed in central New Mexico (Amato and Becker, 2012) yields a high correlation to factorized source 6 (R²: 0.81). We interpret ca. 1.66 Ga source to reflect detritus sourced by the Mazatzal terrane exposed predominantly in New Mexico during the late Paleozoic by ARM uplifts. This empirical source does not exhibit the minor 1.2–1.5 Ga age modes, which may contribute to misfits in top-down modeling (see the Discussion section).

1.70 Ga basement. Factorized source 6 is a ca. 1.70 Ga unimodal distribution that exhibits a very high correlation (R²: 0.97) to the Caylor Gulch granodiorite exposed in the Front Range of Colorado (Guitreau et al., 2016) (Fig. 5J). This granodiorite is part of a broad suite of granitoid intrusions (1.72–1.68 Ga) that were emplaced during peak deformation of the Yavapai Orogeny (Whitmeyer and Karlstrom, 2007). We use the granodiorite sample to represent the abundant ca. 1.70 Ga age zircon source of the Yavapai basement terrane exhumed throughout the western interior during ARM deformation (Fig. 7).

1.73 Ga basement. Factorized source 4 contains a dominant ca. 1.74 Ga mode and a lesser 1.46 Ga mode that are both positively skewed (Fig. 5K). This factorization exhibits a moderate correlation (R²: 0.65) to our empirical source match with a slightly younger age mode (ca. 1.73 Ga). We interpret this source to broadly represent Paleoproterozoic basement rocks of similar modal age in the Yavapai crustal terrane (Jessup et al., 2006; Whitmeyer and Karlstrom, 2007), which are commonly associated with minor intrusions of Mesoproterozoic igneous rocks. The limited multi-modal character of the factorized distribution, albeit not compensated for in our empirical sample, suggests common heterogeneity in the catchment exposing basement of this age, proximal downstream mixing, and/or recycling of local basement sourced sedimentary rock.

Top-down models generally yield greater misfit than bottom-up models when comparing mixed models and empirical samples (Fig. 6A). This is due to the imperfect matches of the empirical source set to the factorized sources and propagation of uncertainty in the match. The mean R² between top-down mixed distributions and empirical samples presented above is 0.69 ± 0.14 (1σ) and the median is 0.68. The mean Kuiper V value for these comparisons is 0.19 ± 0.08 (1σ) and the median is 0.16. Attribution of local sources between top-down R² and Kuiper V models are highly correlated (R² = 0.92; Fig. 6C). Comparison of R²s for samples between bottom-up (DZnmf) and top-down (DZmix) indicate a poor correlation (R² = 0.46). Despite this difference in the goodness of mixed model fits between top-down and bottom-up approaches, source attribution for individual samples is remarkably similar across both. This is illustrated by the high correlation of local source attribution of both top-down and bottom-up mixture models (R² = 0.85; Fig. 6B). In other words, the mixing weight applied by the bottom-up and top-down approaches is similar, but the goodness of fit between sample and mixed distributions is not. In the case of top-down modeling, these results indicate a geologically realistic empirical source set that works well for sink data attribution, but one that would likely benefit by improving its fit to the factorized sources.

Basins within the ARM corridor exhibit an increase in the relative contribution of modeled local source attribution during ARM deformation (Fig. 2). However, presenting the results at the scale of a stratigraphic column, as in Figure 2, does not illustrate sediment fairways through time. We therefore complement this approach with a synoptic map view of these basins in the context of regional provenance trends.

To address the challenge of data visualization, we create provenance-contour maps that incorporate paleogeography, stratal thickness, and paleocurrent data. DZmix-modeled local source percentages are displayed next to data points (Fig. 8), which are shaded based on their goodness of fit (R²). These data and this approach provide geographic context for the changes in detrital zircon source identified in the following description of sediment provenance episodes. Maps shown in Figure 8 discriminate between local and distal sources and include symbols on data points that distinguish threshold contribution percentages for several distal sources (e.g., a light blue square that indicates >20% Antler recycled attribution).

Combining bottom-up and top-down modeling provides an effective approach to distill large detrital zircon data sets into a realistic empirical source set and assign source proportions (Fig. 3). Instead of describing each sample as discretized age bins, we approach detrital zircon age distributions as continuous functions, similar to seismic waves, and leverage the size of the data set to our interpretive advantage. However, source signals, or factorized sources, are only useful if they provide information that can be interpreted within a geologic context. This is where the iterative approach to identifying empirical sources and then comparing them to the factorized sources is particularly instructive. Top-down modeling also helps to mitigate bias. In the case of DZmix, it does so by randomly mixing known source inputs (a.k.a. the empirical source set) to produce a mixture model distribution. This modeled distribution is then compared to the actual detrital zircon age distribution to determine the best model fits. A benefit of using this approach over age binning of data is that the model can generally discriminate between sources containing the same non-unique ages based on the shape of the distribution and/or accompanying age modes. Moreover, it considers the proportions of those age modes. Combining bottom-up and top-down modeling helps to decrease bias, lay out an approach to source analysis for large provenance data sets, and provide guidance in sediment source identification and attribution.

However, the combined bottom-up and top-down approach is not without limitations. Despite developing an empirical source set that is very highly to moderately correlated (see definitions in Table 1) to the respective factorized sources (Fig. 5), top-down models yield a worse fit to samples than do bottom-up models (Fig. 6A). The greater misfit between top-down models and samples in comparison to bottom-up models is most likely due to the NMF algorithm attributing minor age modes to samples where they are absent in the empirical source samples. For example, factorized source 9 includes a minor age mode at ca. 1.46 Ga that is absent from the Triassic arc empirical sample (see the section Step 2. Comparison of Factorizations to Empirical Sources above). Nevertheless, correlation of source attribution between DZnmf and DZmix is high (i.e., mean R² of source attribution is 0.87 and local source attribution R² of 0.94; Fig. 6A). These metrics indicate high fidelity between source attribution of the two mixture modeling methods despite the worse top-down model fits. Further, these results suggest that uncertainties associated with mixed model-to-sample distribution comparisons, as well as factorized-to-empirical source comparisons propagate from bottom-up modeling to factorized-empirical source comparison and then to top-down modeling. It follows that there is room for improvement in the empirical source set, as well as in the approach itself.

Bottom-up modeling, particularly with tools such as DZnmf, provides a clarifying approach to evaluate which detrital zircon sources play a critical role in the basin record, and which do not. For example, it may be troubling to some readers that we use only three source ages (1.66, 1.70, and 1.73 Ga; Figs. 5I–5K) to represent detrital zircons from the Yavapai and Mazatzal terranes, as they host a protracted assemblage of igneous ages that span 1.60–1.80 Ga (Whitmeyer and Karlstrom, 2007). However, these ages directly correspond to factorized age distributions constructed from an 11-source factorization (factorized sources 1, 6, and 4, respectively; Figs. 5I–5K) that exhibit high to very high correlation (R² of 0.90 ± 0.08 (1σ) mean and 0.91 median) between sample and factorized model age distributions. In other words, the well correlated bottom-up model indicates that ages outside of the ones characterized in factorized sources 1, 6, and 4, but included in the Yavapai and Mazatzal terranes, were not major contributors to the sedimentary record examined here.

In some cases, non-unique age distributions remain an issue when applying these methods and point to the need for integrative data interpretation. In the case of discerning between similar age distributions of basement versus recycled sources, use of sandstone petrography (i.e., quartz-rich versus feldspar-rich) is a straight-forward approach to aid in discriminating sources. For example, the Denver Basin shows a relatively low percent of local source relative to other ARM basins proximal to uplifts during Episode 3 (e.g., Central Colorado Trough (CCT), Taos Tough). Although these strata (i.e., the Fountain Formation) are arkosic (Suttner and Dutta, 1986), top-down modeling of the 11-source set yields surprisingly high distal source attribution. This is due to the absence of an age distribution correlative to the 1.085 Ga Pikes Peak Batholith, which was exhumed in the Ancestral Front Range, in the 11-source set. Previous research makes a compelling case for the source of Pennsylvanian–Permian strata in the Denver Basin (Fountain and Ingleside formations) to have been shed by the adjacent Ancestral Front Range (Sweet and Soreghan, 2010; Leary et al., 2020). In the absence of a 1.085 Ga source, DZmix selected the central Appalachian source (Fig. 5B) as the best match to Fountain and Ingleside formation samples. DZnmf factorizes out a 1.085 Ga age distribution in the 12-source factorization (Supplemental Material 5) and we experimented with re-mapping Episode 3 with DZmix results from an 11-source set plus a Pikes Peak source (Supplemental Materials 1 and 6). Improvement of the DZmix model fits for the Denver Basin samples were proportional to the percent of the 1.085 Ga source attributed. The most drastic change was observed in an arkosic Denver Basin sample (Fountain Formation 2 SG; Supplemental Material 2) where the model fit improved over 2σ (from 0.58 to 0.91), central Appalachian source decreased (from 0.61 to 0.03), and 1.085 Ga source attribution was 0.84 (Supplemental Material 6). The provenance map contoured with modeled results from this alternate source set, which contains the 1.085 Ga source (Supplemental Material 1) shows that the Denver Basin appears to have a small area that was dominated by local sediment. The 11-source set model map in Figure 8C does not capture this. However, we decided not to include the 1.085 Ga source in the source set presented here because it is a relatively minor source, and although its addition improves some of the Denver Basin samples, it appears to be erroneously attributed to other samples distant from the Pikes Peak batholith, albeit in low proportions (Supplemental Materials 1 and 6), thereby inflating the influence of local source attribution across a broad area. Therefore, we decided that on balance we could not justify the addition of this source to our source set given the scale of our provenance maps. Instead, we highlight this result as a potential shortfall of the approach we present here.

Source mapping, as presented here, provides a new method to visualize provenance data. This approach is different from previous attempts to map sediment sources (e.g., Wissink et al., 2016; Blum et al., 2017) in that we are not mapping proportions of binned ages, but rather proportions of sources which are themselves characterized by continuous functions. This approach also creates a synoptic view of provenance data that incorporates other geologic information (e.g., stratigraphic thickness, paleocurrents). Integration of additional data and hand-contouring yield a final product with an element of interpretation. Although this inherently introduces bias, it also facilitates interpretation of a larger, more diverse set of data. We use top-down modeled detrital zircon source contributions as the primary input; however, it would be possible to substitute these source proportions with one of many other quantitative or semiquantitative provenance proxies (e.g., proportions derived from clast-counts, bottom-up model source results, a selected binned age-group).

When creating a source map, data selection and filtering are of top importance and need to be customized to the purpose of the map. In the case presented here, the purpose is to investigate the distribution and magnitude of ARM-driven provenance in western Laurentia. Therefore, we lumped the modeled sediment sources into local or distal, contoured local source contribution with consideration of auxiliary data, and included threshold symbols for the dominant distal sources (e.g., a pink triangle for samples with >20% Triassic arc source attribution; Fig. 8).

The question of which data set is better to create source maps with, top-down or bottom-up modeled source proportions, is not addressed here. While we do not contour bottom-up results, local source attribution is highly correlated (R² = 0.94) between DZmix and DZnmf models (Fig. 6B), indicating that source maps would be similar. However, the answer to which approach is better is likely contingent on the purpose of the map(s) one is contouring. The purpose of this paper is to evaluate the distribution of local sources around the ARM, and to introduce a new approach of handling, interpreting, and visualizing large provenance data sets. For the latter purpose, showing the loop of using basin samples to create bottom-up factorized sources that are then used to inform selection of real-rock sources, which, in turn, are used in top-down modeling, is valuable in demonstrating the range of possibilities in this approach. We also prefer the real-rock link between source and sink inherent in modeling with the top-down model results. However, bottom-up results exhibit less of a misfit than top-down results (mean DZnmf R² = 0.91 versus mean DZmix R² = 0.70; Supplemental Material 6), leaving us agnostic as to which is “better” to map with. Further, it may be reasonable to say that there is no purpose in compounding the misfit between factorized and real-rock sources by modeling with an inferior source set (i.e., empirical source set used in top-down modeling). We hope to see many applications of this new approach, and utilization of many different provenance data sets, and, in doing so, the applicability and best practices for it will be developed.

Many of the factorized sources presented in this paper are similar to factorized sources presented in Leary et al. (2020), and most of our interpretations are similar as well (Table 2). However, a few sources presented here that exhibit similarity to Leary et al. (2020) sources are interpreted differently. For example, factorized sources 3 and 7 are interpreted here as Maya Block (or similar peri-Gondwana terrane) sources (Fig. 5). However, combined, these sources are similar to “Source A” from Leary et al. (2020), which was interpreted as a northern Appalachian source. The difference in our results, and therefore interpretations, is likely due to the expanded detrital zircon data set used here that includes additional new and previously published samples (Supplemental Material 2).

Factorized source 1, interpreted as 1.66 Ga basement, is similar to “Source D” in Leary et al. (2020), but lacks the 1.75 Ga mode and replaces the 1.1–1.4 Ga distribution with a minor peak at 1.35 Ga. We attribute the minor but clarifying differences in age distribution observed in factorized source 1 versus “Source D” to the larger size of the data set analyzed here, similar to the case mentioned above. This may have allowed the NMF algorithm to better characterize this source and/or further deconvolve a mixed “Source D” into its constituents. To this point, the primary recipients of factorized source 1 are several New Mexico samples that were added to our data set, but were not in the Leary et al. (2020) data set. With stated uncertainty, Leary et al. (2020) attributed “Source D” to either exhumation of the Zuni Uplift or recycling of foreland basin strata along an expanded Antler orogenic front. However, we favor a Mazatzal basement source interpretation, which underlies the Zuni Uplift.

We divide the Late Devonian–Triassic study interval into five episodes based on modeled sediment provenance trends and comparison to previously published literature (e.g., Dickinson and Gehrels, 2008; Gehrels et al., 2011; Chapman and Laskowski, 2019; Leary et al., 2020). Characterization of sediment provenance episodes is focused on the core ARM basins (Fig. 1), is defined by changes in the distribution of local versus distal sourced sediment, and reflects major changes in sediment routing and provenance in western Laurentia. These episodes include (1) arch-controlled sediment dispersal (Late Devonian–Middle Mississippian), (2) transcontinental sediment integration (Late Mississippian–Early Pennsylvanian), (3) ARM basin isolation (Early Pennsylvanian–early Permian), (4) tectonic quiescence and eolian sediment mixing (early Permian–middle Permian), and (5) transcontinental fluvial reintegration (Triassic: predominantly Late) (Fig. 2).

The eight samples and sample groups that comprise the Late Devonian–Middle Mississippian Episode 1 exhibit different detrital zircon age distributions based on their geographic location relative to the transcontinental arch (Fig. 8A). This pattern of sediment dispersal is documented as far back as the Cambro-Ordovician (Amato and Mack, 2012), the extent of which is outside the scope of this paper. During this interval western Laurentia received a mix of local (i.e., Yavapai and Mazatzal crustal terrane) and distal (predominantly Antler recycled) sources. Age distributions to the east of the transcontinental arch (Figs. 2Civ and 2Cv; Supplemental Material 6), as well as those deposited on the Montana shelf are predominantly attributed to the central Appalachian source (Fig. 5B) with a significant contribution from the early Paleozoic recycled source (Fig. 5F). Samples that fall along the transcontinental arch (e.g., sandstones within the Leadville Limestone) are either dominated by local sources (e.g., incipient Taos Trough sample) or contain a mixture of local and Antler recycled sources. To the west of the arch, the White River and Ely basins (Fig. 1) appear to be dominated by the Antler recycled source (Fig. 5E). We favor a model of sediment compartmentalization across the transcontinental arch (Linde et al., 2014) with intermittent cross-arch connectivity during sea-level highs (Chapman and Laskowski, 2019) to explain provenance distribution observed in the data (Fig. 8A). We interpret the source of local sediment within the study area to be directly arch-derived or recycled from local sediment that itself was originally arch-derived.

Episode 2 marks a period of major sedimentation reorganization in western Laurentia during the Late Mississippian–Early Pennsylvanian (Chesterian–Morrowan) (Fig. 8B), but it is complicated by a protracted hiatus in the Late Mississippian (Chesterian). McKee and Crosby (1972, plate 2) show that most rock beneath Pennsylvanian stratigraphy in western North America is Middle Mississippian (Meramecian) and older, which highlights the paucity of Upper Mississippian (Chesterian) sedimentary rock throughout the study area. This subcrop pattern follows the trace of the transcontinental arch, indicating its persistence as a positive feature at least up to Early Pennsylvanian (Morrowan) (Blakey, 2009).

There is a rapid, but slightly diachronous change in provenance seen in Upper Mississippian–Lower Pennsylvanian strata from the Arizona shelf, Paradox Basin, and Montana shelf (e.g., the Molas and Surprise Canyon formations; Fig. 2). Upper Mississippian (Chesterian) fluvial-estuarine rocks (i.e., Surprise Canyon Formation) (Billingsley, 1999) record the introduction of Appalachian-sourced sediment to the western Laurentian margin (Gehrels et al., 2011) (Fig. 8A). However, this was followed within ~5 m.y. by an increase in local-source zircons in the Lower Pennsylvanian (Morrowan) marine strata (i.e., Watahomigi Formation) (Hodnett and Elliott, 2018) (Figs. 2 and 8B). With no further context, these data could suggest that this increase in local provenance was due to early exhumation of ARM uplifts. However, Lower Pennsylvanian (Morrowan) fluvial deposits on the Montana shelf (i.e., Amsden Formation) and incipient Paradox Basin (i.e., Molas Formation, which stretches into the Atokan; Evans and Reed, 2007) exhibit similar facies, and detrital zircon age distributions as the Upper Mississippian (Chesterian) fluvial strata on the Arizona shelf (Fig. 8B). Like these fluvial Arizona shelf strata, the Lower Pennsylvanian fluvial deposits on the Montana shelf and Paradox Basin were deposited on top of karsted Mississippian carbonates. This observation spurs the question: Are these diachronous deposits related, and if so, how?

We interpret this diachronous package of fluvial rock to record progressive eastward infilling of incised fluvial valleys and onlap that developed on the complex-karsted topography of Mississippian limestones (Evans and Reed, 2007). The transgression was driven by a second order eustatic sea level rise that flooded vast continental areas in the earliest Pennsylvanian (Sloss, 1963; Vail and Thompson, 1977; Haq and Schutter, 2008). Local detrital zircon sources (1.66 and 1.70 Ga; Figs. 5I and 5J) observed in the marine Watahomigi Formation may indeed be a signal of early ARM exhumation (i.e., either from exposed basement or recycled from early Paleozoic strata), but are not from the area around the northern Uncompahgre Uplift. This can be extrapolated from the predominant distal source of detrital zircons in the approximately coeval fluvial rocks (i.e., Molas Formation) (Fig. 8B). Therefore, while fluvial systems that supplied the incipient Paradox Basin may have also delivered that distal sediment to the Arizona shelf, they did not transport locally derived sediment to it, or at least they did not until tapping basement sources downstream of the northern Paradox Basin area. Furthermore, most Paradox Basin samples receiving sediment from the Uncompahgre Uplift during the subsequent episode contain a significant component of 1.43 Ga source, which the marine rocks of the Arizona shelf (i.e., Watahomigi Formation) do not contain (Fig. 2 and Supplemental Material 6). Instead, during this time, we attribute local detritus identified in Arizona shelf strata to be from basement exposed either along the transcontinental arch or areas of early ARM exhumation (Fig. 8B), as unroofing of some ARM uplifts is documented during this time (e.g., Musgrave, 2003). This may include the southern portion of the Uncompahgre Uplift (also termed San Luis Uplift sensu Blakey, 2009).

The Episode 2 source map (Fig. 8B) depicts sediment dispersal systems controlled by the transcontinental arch. Both during the Late Mississippian, when it served as a large area of non-deposition and karsting, and during the Morrowan, when sediment covered it as seas rose. It was also the approximate site of possible early ARM deformation (Leary et al., 2017; Sweet et al., 2021) and local sedimentation in a few basins (e.g., CCT; Musgrave, 2003). While it would be possible to route sediment around nascent ARM uplifts, the presence of a continental scale arch, at least until the Early Pennsylvanian (McKee and Crosby, 1972) would appear to be an impediment to central Appalachian sediment arriving on the western Laurentian margin in the Late Mississippian. Therefore, how would transcontinental sediment pass this barrier?

We offer several possible scenarios that describe how sediment could have been transported from the eastern to the western margin of Laurentia (modern reference) during this time. (1) A southwest-directed pathway that delivered Appalachian sediment from the northeast, which effectively routed detritus around the northern extent of the arch. (2) Pathways that intermittently transported sediment over the arch by ocean currents during higher frequency sea level rises (Chapman and Laskowski, 2019). (3) Pathways that delivered sediment were routed through paleovalleys that cut across the transcontinental arch, similar to, but with different orientation than Late Mississippian paleovalleys identified in western Kansas and eastern Colorado (Sonnenburg et al., 1990; Bowen and Weimer, 2003; Wang and Bidgoli, 2019). (4) Alternatively, Carlson (1999) suggests that the transcontinental arch was not a continuous feature, but rather a discontinuous platform consisting of a series of basement-cored features. In this scenario, sediment pathways could have utilized gaps in the arch that may have been enhanced, or even created by early subtle ARM deformation.

Integration of the continental-scale drainage system during Episode 2 was short-lived in the core of the ARM. Along the transcontinental arch, pockets of Lower Pennsylvanian (Morrowan) strata in the CCT, Taos Trough, and western Orogrande Basin received predominantly local sediment (Fig. 8B) contributed by 1.70 and 1.66 Ga sources. Unlike subsequent deposition in the Pennsylvanian, these initial sandstones are quartz-rich (DeVoto et al., 1971; Baltz and Meyers, 1999; Musgrave, 2003; Amato, 2019) and are interpreted to indicate recycling of lower Paleozoic strata from along exposed areas of the transcontinental arch and/or nascent ARM uplifts.

Episode 3 is defined by the dominance of local basement sediment deposition around core ARM uplifts during the Early Pennsylvanian–early Permian (Atokan–Wolfcampian). The Paradox Basin, Eagle Basin, CCT, and Taos Trough contain large areas of ≥80% local source attribution. The Denver Basin appears to be an outlier regarding locally sourced sediment. However, although it is not as sedimentologically isolated as other nearby ARM basins, its uncharacteristically low local attribution illustrated in Figure 8C is in part an artifact of a shortcoming in the 11 source-set used to model it. This source set lacks the 1.085 Ga granitoid source, which is adjacent to the strata bearing this age zircon, and instead misattributes the Grenville-dominant central Appalachian source. For further detail of this misattribution, and why we decided not to include the 1.085 Ga source in the empirical source set applied here (Fig. 5), see the Consideration of the Methods section in the Discussion section above. Evidence of local-derived sediment in the Denver Basin resides in the Lower–Middle Pennsylvanian coarse-grain strata along the eastern basin margin, which document unroofing and subsequent arkosic sedimentation during this episode (Sweet and Soreghan, 2010). For the first time, the 1.43 Ga source (Fig. 5H) serves as a major contributor to the local source profile of many samples. Areas of local source attribution outside of core ARM basins expand to the northwest, commonly exhibiting 20%–50% local source attribution. The Arizona shelf is a mixing pot of local and distal sediment. It received distal sediment from a corridor to the northeast that appears to have transported detritus between the Emery and Piute uplifts to the west and the Uncompahgre Uplift to the east (Figs. 1 and 8C).

The prominence of basement uplifts as dominant sediment sources in the core of the ARM was not ubiquitous across all ARM basins. The Central Basin Platform, which separates the Delaware and Midland basins (Fig. 1), was only modestly emergent despite causing >2 km of flexural subsidence in the adjacent Permian Basins (Yang and Dorobeck, 1995). Instead, it primarily served as a high-relief carbonate platform (Ewing, 1993) between the two deep Permian Basins (i.e., Delaware and Midland basins; Fig. 1). Therefore, it never became a major contributor of siliciclastic sediment during ARM deformation, and adjacent basins instead received siliciclastic sediment from distal sources. This example of an ARM basin where adjacent local (siliciclastic) sediment production was suppressed highlights the continued presence of distal detritus and suggests that core ARM basins were overfilled by locally derived sediment. Furthermore, source mapping of this episode (Fig. 8B) suggests that most of the overflow of local sediment from core ARM basins was predominantly transported west (e.g., Arizona shelf) rather than east (e.g., Permian Basin area).

During this episode distal sources continued to dominate sedimentation outside the core ARM. Transcontinental sediment delivery to areas around western Laurentia did not shut down during ARM deformation, but was likely routed around it. Recycled Antler (Fig. 5E) and Maya Block (Figs. 5C and 5D) source contributions were location dependent with the greatest proportions located along the western Laurentian margin and Midland Basin, respectively. In contrast, the central Appalachian source (Fig. 5B) contributed to samples across the study area. Northeast to southwest sediment pathways running through the Wyoming and Montana shelves, northern Utah, eastern Nevada, and southeastern California provided routes for abundant central Appalachian sediment. The close match between measured samples and mixture models that incorporate the central Appalachian source support the continued delivery of these zircons to western Laurentia during Episode 3. The recycled Antler source dominated samples from the Ely Basin and contributed notably to samples in southwest Montana. We interpret this to indicate the intermittent tectonic activity along the western Laurentian margin (Trexler et al., 2004; Sturmer et al., 2018). In west Texas, the Maya Block and central Appalachian sources dominated the Marathon and Midland basins. We acknowledge the possibility that the central Appalachian source attribution in these samples may instead represent a similar Grenville age population contributed by peri-Gondwanan terranes to the south (e.g., Coahuila terrane as suggested by Thomas et al., 2019). In fact, Wolfcampian samples tend to exhibit poorer model fits. However, regardless of eastern or southern orogenic source, provenance mapping in this area offers a second possible distal source pathway into western Laurentia (i.e., from the southwest), although data are scant in the southwest during this episode.

Episode 4 is defined by an increase in eolian deposition (Blakey, 2008; Lawton et al., 2015) and a decrease in local sources throughout most of western Laurentia during the early–middle Permian (Leonardian–Guadalupian) (Fig. 8D). These changes are associated with waning ARM tectonics (Blakey, 2008) and an increasing aridity (Crowell, 1999; Preto et al., 2010; Lawton et al., 2021). This combination of factors facilitated increased mixing of abundant distal sources and limited local sources (Lawton et al., 2015, 2021), and homogenized the provenance signal across western Laurentia. The northern Paradox Basin appears to be the only exception as the Leonardian upper Cutler Formation (Trudgill et al., 2011) exhibits high local source attribution. This may be due to late-stage ARM exhumation of the northern Uncompahgre Uplift (Blakey, 2009). Across the rest of western Laurentia, samples and sample groups indicate <35% local attribution. Higher local sediment attribution occurs in a curved, mostly longitudinal band, and appears to come from the Uncompahgre Uplift. It is important to note that the data used to draw the contours that illustrate the longitudinal band are sparse, and one could imagine filling in the empty space on the map (Fig. 8D) with a different shape than is shown. The 1.43 Ga source (Fig. 5H) is dominant in the proximal northern Paradox Basin. However, the 1.66 Ga (Fig. 5I) source constitutes most of the local source attribution along the curved band. This may reflect a basement-derived detrital zircon age from the Uncompahgre Uplift that is not captured by regional basement maps (Whitmeyer and Karlstrom, 2007) (Fig. 7), sediment from another exposed remnant ARM uplift along that band (e.g., Zuni Uplift: Fig. 1), and/or recycling of recently deposited sediment (i.e., deposits from Episode 2).

Central Appalachian (Fig. 5B) and Maya Block sources (Figs. 5C and 5D) constitute the predominant distal sources, whereas the early Paleozoic recycled source (Fig. 5F) is present in a few samples, most of which are on the Wyoming and Montana shelves. Central Appalachian sourcing was the dominant distal source for a large area extending westward from ~105°W (present-day coordinates). The dominant distal source in the southeastern corner of the study area was the Maya Block, which is consistent with the location of the Midland, Delaware, and Marathon basins adjacent to peri-Gondwanan terranes. However, attribution of this source farther north is at odds with most fluvial and eolian paleocurrent data, which indicate southward sediment transport (Fig. 8D). We interpret this to indicate that either paleocurrent measurements presented here do not fully capture sediment transport direction during this episode, or that Maya Block attribution is serving as a surrogate for similar detrital zircon age distributions in the Appalachian Mountains (e.g., Kissock et al., 2018).

Episode 5 is separated from the early–middle Permian eolian sediment mixing by unconformities in Utah, Colorado, New Mexico, and Arizona (McKee and Oriel, 1967), and is marked by a return to fluvial-dominated deposition. This sediment transport system delivered detritus across Laurentia (Lawton, 1994; Dickinson and Gehrels, 2008), connecting the relict collisional margin of the Ouachita-Marathon Mountains to the Nevada shelf (Dickinson and Gehrels, 2008). These Triassic depositional systems span from the Sonoma Orogeny in Nevada to the initial development of the Cordilleran arc (Fig. 2) (Lawton, 1994), thereby marking another major change in tectonic regime and sediment transport in western Laurentia. Although this episode is broadly distinguished as Triassic, the data displayed in the source map (Fig. 8E) are mostly Late Triassic, and are therefore representative of that time. Due to a paucity of data from Lower–Middle Triassic Moenkopi Formation and equivalent strata, we combined these data with the abundant Late Triassic data.

The source map for Episode 5 (Fig. 8E) illustrates fluvial systems that connected the relict collisional margin of the Ouachita-Marathon Mountains to the Nevada continental shelf (Dickinson and Gehrels, 2008). Modeled samples exhibiting high distal source attribution stretch from southeast to northwest across western Laurentia (Fig. 8E). The Triassic arc source (Fig. 5A), for the first time, was a notable contributor to some detrital zircon samples along the southwestern Laurentian margin. The digitate pattern of Triassic arc attribution exhibited by samples along the southwestern highlands highlights the variability of provenance in these transverse fluvial drainages interacting with the trunk system that flowed to the northwest (Fig. 8E). Local sources along the southwestern margin of the study area are predominantly attributed to 1.43 Ga basement. The central Appalachian (Fig. 5B) and Maya Block (Figs. 5C and 5D) were the dominant distal sources, whereas the Antler recycled (Fig. 5E) and early Paleozoic recycled (Fig. 5F) sources were relatively minor contributors. Areas with high proportions of local sources appear to extend in tongues, east-to-west and southeast-to-northwest from basement sources that are inferred to have been exposed at this time (Dickinson, 2018). For example, the tongue which extends from southwestern Oklahoma is dominated by the 530 Ma source (Fig. 5G) and is most likely from the Amarillo Wichita Uplift (AWU) (Riggs et al., 1996; Dickinson and Gehrels, 2008).

While we favor the AWU as the most likely source of Cambrian detrital zircons in northern Texas and eastern New Mexico, we propose a possible alternative source for the northeastern Utah sample. It is possible that this sample was instead sourced by either Cambrian age igneous rocks still exposed in parts of the largely buried Ancestral Front Range, or by Pennsylvanian–Permian CCT strata that were themselves sourced from those same Cambrian rocks. Recent provenance work reveals that Pennsylvanian–Permian strata in this area were, at times, exclusively sourced by proximal Cambrian rocks (Smith et al., 2023). The exclusivity of these Cambrian detrital zircons in some CCT strata demonstrate potential significance of Cambrian age igneous rocks and recycled sedimentary rocks west of the AWU as sediment contributors. The source map for Episode 5 (Fig. 8E) illustrates these two hypotheses by including a red swath across an area of no data to depict a possible sediment pathway from the AWU (as proposed by Dickinson and Gehrels, 2008). Alternatively, the map also includes a few areas of the remnant Ancestral Front Range uplift to indicate possible sourcing from a more proximal region.

We present two novel elements for interpreting large detrital zircon geochronology data sets. (1) Combination of bottom-up and top-down source modeling, and (2) sediment-source mapping. The first approach leverages large data sets and source modeling to provide guidance in developing a vetted empirical source set and determining source attribution for large detrital zircon data sets (Fig. 3), rather than focusing on lengthy descriptions of binned age groups. While the data set used here is large, when filtered for time intervals of interest and plotted on maps, large areas contain no data, whereas others are densely populated. Source mapping is an integrative and interpretive approach that addresses this data distribution heterogeneity, aiding in geographic visualization of sediment dispersal systems, and demonstrating the utility of this approach across areas of sparse data. As with other geological applications, where data are limited, creation of contour plots allow extrapolation in a way that honors available data.

We applied these techniques to a late Paleozoic–early Mesozoic Laurentian data set and show that 191 samples and sample groups (created from 329 samples) can be factorized into 11 sources for which real examples (i.e., empirical sources) can be found. Contour mapping of top-down mixture model results and synthesis of previously published data reveal that western Laurentia experienced five episodes of sedimentation defined by changes in provenance and sediment routing particularly around the core of the ARM (Fig. 1). These include: (1) arch-controlled sediment dispersal (Late Devonian–Middle Mississippian), (2) transcontinental sediment integration (Late Mississippian–Early Pennsylvanian), (3) ARM basin isolation (Early Pennsylvanian–early Permian), (4) tectonic quiescence and eolian sediment mixing (early Permian–middle Permian), and (5) transcontinental fluvial reintegration (Triassic: predominantly Late). Source maps uniquely reveal the degree and extent of ARM-associated tectonically driven source isolation. Furthermore, areas of intraplate deformation diverted transcontinental sediment dispersal systems around affected areas but did not shut those dispersal systems off. Sediment source modeling and provenance mapping indicate that influx of central Appalachian sediment from the northeast remained relatively consistent in the late Paleozoic. Beginning in Episode 3, Maya Block sources began to dominate provenance signatures inboard of the Marathon fold and thrust belt. The Triassic ushered in a new episode of transcontinental sediment dispersal, this time from the south and southeast, which included magmatic arc material.

Funding for field and lab work, as well as publication costs, was provided by National Science Foundation (EAR-1824557), and graduate students grants from American Association of Petroleum Geologists, Geological Society of America, Sigma Xi, the University of Houston (UH), and University of British Columbia. Many thanks to field assistants Lokin Casturi, Crystal Saadeh, Kenny Lambert, and Noah Karsky. Special thanks to Minako Righter and Yongjun Gao for assisting with the U-Pb analyses and maintaining the UH laser-ablation–inductively coupled plasma–mass spectrometry lab facilities, to Tim Lawton, Bill Thomas, Scott Ryan, and two unnamed reviewers for their insightful and constructive reviews on this and/or earlier versions of this manuscript, to Tim Kusky, associate editor, for comments and guidance, and to Don Rasmussen for sharing his decades of experience in the Paradox Basin of western Colorado and Utah (USA). Thanks to James Hagadorn and Christopher Holm-Denoma for assistance navigating the U.S. Geological Survey data repository.