The details and mechanisms for Neogene river reorganization in the U.S. Pacific Northwest and northern Rocky Mountains have been debated for over a century with key implications for how tectonic and volcanic systems modulate topographic development. To evaluate paleo-drainage networks, we produced an expansive data set and provenance analysis of detrital zircon U-Pb ages from Miocene to Pleistocene fluvial strata along proposed proto-Snake and Columbia River pathways. Statistical comparisons of Miocene-Pliocene detrital zircon spectra do not support previously hypothesized drainage routes of the Snake River. We use detrital zircon unmixing models to test prior Snake River routes against a newly hypothesized route, in which the Snake River circumnavigated the northern Rocky Mountains and entered the Columbia Basin from the northeast prior to incision of Hells Canyon. Our proposed ancestral Snake River route best matches detrital zircon age spectra throughout the region. Furthermore, this northerly Snake River route satisfies and provides context for shifts in the sedimentology and fish faunal assemblages of the western Snake River Plain and Columbia Basin through Miocene–Pliocene time. We posit that eastward migration of the Yellowstone Hotspot and its effect on thermally induced buoyancy and topographic uplift, coupled with volcanic densification of the eastern Snake River Plain lithosphere, are the primary mechanisms for drainage reorganization and that the eastern and western Snake River Plain were isolated from one another until the early Pliocene. Following this basin integration, the substantial increase in drainage area to the western Snake River Plain likely overtopped a bedrock threshold that previously contained Lake Idaho, which led to incision of Hells Canyon and establishment of the modern Snake and Columbia River drainage network.

We resolve the Neogene history of the Snake River, as modulated by the Yellowstone Hotspot, using detrital zircon U-Pb provenance analysis.

Mantle dynamics and their impact on the thermomechanical state of the overlying lithosphere are often invoked to explain topographic change (Stanley et al., 2013, 2015); however, few geologic records document the interplay of competing drivers of landscape evolution. Miocene and younger strata in the Columbia and Snake River basins preserve a rich record of landscape evolution and coincident faunal shifts along the Yellowstone Hotspot track and have spawned considerable research during the past century on regional landscape evolution (Lindgren, 1898; Livingston, 1928; Wheeler and Cook, 1954; Beranek et al., 2006; Stearley and Smith, 2016). Central to the history of the Snake River is Hells Canyon, through which the Snake River currently carves the deepest canyon in North America (Fig. 1), and its relation to draining Lake Idaho, a massive lake centered in the western Snake River Plain during late Miocene–Pliocene time (Wood and Clemens, 2002; Fig. 2).

On the basis of anomalous tributary orientations and a thick sequence of lake strata upstream of Hells Canyon (Fig. 1), early researchers hypothesized that the Snake River likely changed its course in recent geologic time, but these studies disagree upon the mechanism of landscape evolution and interpret vastly different pathways for the ancestral Snake River (Lindgren, 1898; Livingston, 1928; Wheeler and Cook, 1954). Proposed paths include a diversion around the Wallowa Mountains via the Baker and Grande Ronde Valleys (Livingston, 1928) and a route across eastern Oregon, Nevada, the Sierra Nevada, and into the Sacramento River drainage (Wheeler and Cook, 1954; Fig. 2). More recent studies of fish and rodent fossils from the western Snake River Plain and Columbia Basin suggest that Lake Idaho was isolated from the Columbia Basin until the late Pliocene, which lends support to the Sacramento drainage route of the Snake River (Repenning et al., 1995; Smith et al., 2000; Stearley and Smith, 2016). However, detrital zircon age distributions of Miocene–Pliocene Ringold Formation lacustrine deposits from the Columbia Basin suggest that the Snake River was part of the Columbia River drainage network over the same time period (Staisch et al., 2018).

In this work, we use a suite of statistical methods and modeling for detrital zircon provenance of modern and ancestral river sands collected throughout Oregon, Washington, Montana, and Idaho, USA. We synthesize our results with existing paleontological evidence for basin isolation and integration to propose a sequence of radically shifting Snake River pathways consistent with both the paleontological and sedimentological evidence, driven by topographic adjustments to dynamic volcanic and tectonic environments resulting from passage of the Yellowstone Hotspot and magmatic densification of eastern Snake River Plain continental crust.

The Columbia River originates in Columbia Lake in British Columbia, Canada, and flows for roughly 2000 km and drains over 660,000 km2 before it enters the Pacific Ocean near Astoria, Oregon. It is one of the larger river networks in North America and ranks high among global river networks in terms of mean flow, mean elevation, and slope (Vörösmarty et al., 2000). The Snake River, which originates in northwestern Wyoming, represents nearly 40% of the total Columbia River watershed. On its path toward the Columbia Basin, the Snake River carves Hells Canyon, which is the deepest canyon in North America (Fig. 1).

Tributaries immediately upstream of Hells Canyon, such as Indian, Pine, and Wildhorse Creeks, are classic “barbed tributaries” that flow dominantly southward, away from the northerly path of the Snake River along this reach (Fig. 1), and indicate that the original path of the trunk stream was opposite its current trajectory (Wheeler and Cook, 1954). Farther upstream, the Snake River flows through the broad and low-relief Snake River Plain (Fig. 2), a major topographic depression bounded in the north and east by the northern Rocky Mountains and in the south by the Basin and Range province (Fig. 2).

While currently an integrated basin, the eastern and western Snake River Plains are physiographically and geologically distinct basins. The eastern Snake River Plain trends generally northeast and, unlike the western Snake River Plain, does not appear to be fault-bounded (Wood and Clemens, 2002), but it subsided because of magmatic densification of the lithosphere since ca. 10 Ma (McQuarrie and Rodgers, 1998) and Basin and Range extension (Parsons et al., 1998). The western Snake River Plain is an intracontinental rift basin that is bounded by normal faults active during the past 9.5 m.y. (Wood and Clemens, 2002). The Miocene onset of the northwest-trending graben system may be related to magmatism of the eastward-migrating Yellowstone Hotspot (Clemens, 1993) and cross-cuts the older (10.5–15.5 Ma) north-trending Oregon-Idaho Graben (Wood and Clemens, 2002). Neogene basin-filling sediments are up to 2–3 km thick, including much of the stratigraphic units that record the history of Lake Idaho.

Lake Idaho was a large lacustrine system seated in the western Snake River Plain from <10.1 Ma to 2.5 Ma. The stratigraphy of the lake system suggests that there were two major lake phases recorded by the 10.1–6.4 Ma Chalk Hills Formation and 4.3–2.5 Ma Glenns Ferry Formation (Wood and Clemens, 2002). The ~100-m-thick Chalk Hills sediments are laterally variable and composed of fluvial, lacustrine, and volcanic facies that include basal sand and gravel derived from the Idaho Batholith and local volcanic units, noncalcareous mudstones that become more prevalent upsection, tuffaceous muds, and other volcanic and volcaniclastic beds (Malde and Powers, 1962; Middleton et al., 1985; Wood and Clemens, 2002). The stratigraphy of the Chalk Hills suggests that the lake system was transgressive before ca. 7.5 Ma and recessed sometime after 6.5 Ma (Wood and Clemens, 2002). A gentle angular unconformity between the Chalk Hills and Glenns Ferry Formations suggests that the lake system receded between 6.4 Ma and 4.3 Ma, but the mechanism of recession is not well constrained. The lower Glenns Ferry Formation, locally called the Terteling Springs Formation, is a transgressive lake sequence that includes shoreline oolitic sands that may be indicative of an alkaline, closed lake basin (Swirydczuk et al., 1979; Wood and Clemens, 2002). The majority of the sequence is broadly siliciclastic lacustrine mudstone with minor carbonates and interbedded with volcanic ash and lava beds and is up to ~600 m thick (Middleton et al., 1985). The lacustrine sequence is abruptly overlain by deltaic sands that are known as the Pierce Gulch Sand member of the Glenns Ferry Formation, possibly owing to rapid lake level lowering and erosion of lake-margin sands (Wood and Clemens, 2002; Wood, 2004).

Coincident with much of the tenure of Lake Idaho in the western Snake River Plain is the ca. 9.5–3.4 Ma Ringold Formation, which overlies the Miocene Columbia River Basalt Group and records fluvial and lacustrine sedimentation within the Columbia Basin (Lindsey, 1996; Staisch et al., 2018). Best exposed along the White Bluffs in central Washington, the lower Ringold Formation includes the Taylor Flat and Wooded Island members, which are characterized by coarse channel deposits, fluvial sandstones, and overbank-paleosol deposits that are overlain by the Savage Island member, a series of laterally variable lacustrine and paleosol phases (Lindsey, 1996). Near the Saddle Mountains, U-Pb zircon dating of interbedded tephra constrain the coarse lower Ringold depositional ages to between 9.5 Ma and 6.8 Ma and the uppermost lacustrine phase to between 6.8 Ma and 3.4 Ma (Staisch et al., 2018).

Materials and Methods

Detrital zircon provenance analysis allows potential source terranes to be identified by investigating the age and relative abundance of U-Pb zircon ages. In this study, we present new U-Pb zircon age analyses of eight modern river systems to characterize the detrital zircon populations of known source areas (Fig. 2). We also collected 28 new samples from fluvial sandstones throughout Washington, Oregon, Idaho, and Montana to characterize the source of ancestral river systems in the U.S. Pacific Northwest (Fig. 2). We supplement our new modern and ancestral detrital zircon age spectra with existing data sets collected in Idaho and Montana (Geslin et al., 2002; Link et al., 2005; Beranek et al., 2006; Stroup et al., 2008; Staisch et al., 2018). These data sets include samples from 25 modern rivers and 13 samples from ancestral fluvial strata (Table S11).

All new samples underwent standard zircon separation procedures that included sieving and various stages of density separation in heavy liquids. For all samples except the Chalk Hills Formation (P12–15), we dated ~120 zircons using U-Pb laser ablation-inductively coupled plasma-mass spectrometry (LA-ICP-MS) at the U.S. Geological Survey in Denver. Zircon was ablated with a Photon Machines Excite™ 193 nm ArF excimer laser that was coupled to a Nu Instruments AttoM high-resolution magnetic-sector ICP-MS. The laser spot sizes for zircon were ~25 µm. Raw data were reduced off-line using the Iolite™ 2.5 program (Paton et al., 2011). Ages were corrected by standard sample bracketing with the primary zircon reference material Temora2 (ca. 417 Ma; Black et al., 2004) and/or secondary reference materials FC-1 (ca. 1099 Ma; Paces and Miller, 1993), Plešovice (ca. 337 Ma; Sláma et al., 2008), and Fish Canyon tuff (ca. 28 Ma; Schmitz and Bowring, 2001). Reduced data were compiled using Isoplot 4.15 (Ludwig, 2012). 206Pb/238U ages are reported for igneous zircon samples less than ca. 1300 Ma, and 207Pb/206Pb ages are used for older ages following the recommendations of Gehrels (2014).

Chalk Hills Formation samples (P12–15) were analyzed separately. For these, we dated ~100 zircons using U-Pb LA-ICP-MS methods at the Arizona LaserChron Center (laserchron.org; accessed September 2021) on a Thermo Element2 single-collector LA-ICP-MS using a Photon Machines Analute-G2 ArF 193 nm excimer laser ablation system with HexEx sample cell and a Nu Plasma high-resolution multi-collector. The laser spot sizes for zircon were ~30 µm. Raw data were reduced using AgeCalc (Gehrels et al., 2008; Gehrels and Pecha, 2014), and ages were corrected using zircon reference materials SL (ca. 563 Ma; Gehrels et al., 2008), FC-1 (ca. 1099 Ma; Paces and Miller, 1993), and R33 (ca. 419 Ma; Black et al., 2004) using standard procedures (Gehrels, 2014; Pullen et al., 2018). Plots and age calculations were carried out using Isoplot 3.75 (Ludwig, 2012). 206Pb/238U ages are reported for igneous zircon samples less than ca. 900 Ma, and 207Pb/206Pb ages are used for older ages.

Results from Modern River Sands

Detrital U-Pb zircon age results from 33 modern river sediment samples collected throughout the U.S. Pacific Northwest and northern Rocky Mountains (Fig. 2; Table S1; see footnote 1 for all supplemental figures and tables) show that each fluvial system contains unique age populations that reflect the exposed geology of the contributing landscape upstream at the time of deposition (Figs. 3, S1, and S2). Detrital zircons from modern river sands thus provide characteristic fingerprints associated with specific river systems. Distinctive zircon age groups in modern river samples include the “Lemhi doublet” (Dumitru et al., 2016), which is characteristic of sediment derived from central Idaho, such as in the Salmon and Clearwater Rivers, with zircon ages of ca. 1380 Ma and between 1650 Ma and 1800 Ma (Figs. 3 and S2). Mesozoic age distributions from rivers that drain the Cascade Range overlap with distributions found in the western and eastern Snake River Plain but lack Paleoproterozoic and Archean age populations (Figs. 3 and S2). Importantly, the modern Snake River is the only major river that has a strong presence of 2500–2800 Ma Archean, 1000–1200 Ma Grenville, and 1600–1800 Ma Yavapai-Mazatzal zircon populations (Figs. 3 and S2). Similarly, tributaries to the upper Snake River in Wyoming contain Archean, Yavapai-Mazatzal, and Grenville zircon ages that are recycled through Neoproterozoic and Phanerozoic sandstones exposed in the Idaho-Wyoming thrust belt (Laskowski et al., 2013) (Fig. S2). Grenville ages are absent from streams that drain central Idaho, which otherwise contain Paleoproterozoic and Mesoproterozoic zircons recycled through the Belt Supergroup (Figs. 3 and S2). Streams that drain into the eastern Snake River Plain have a strong component of Grenville ages that the western Snake River Plain rivers lack as well as a notable population of Archean zircon ages (Fig. S2).

Results from Miocene–Pliocene Strata

The distinctively upper Snake River-sourced zircon age groups are present in late Miocene–Pliocene strata from the Ringold Formation and downstream in the Columbia River Gorge and thus provide strong evidence of a fluvial connection between the Columbia Basin and the upper Snake River throughout the late Cenozoic. Grenville and Archean age groups are absent from Miocene Lake Idaho strata from the western Snake River Plain, with the exception of the Glenns Ferry Formation, which is the youngest lacustrine strata related to Lake Idaho (Figs. 3 and S2). Older (>4.3 Ma) lake and fluvial strata in the western Snake River Plain contain locally sourced zircon age groups such as the 42–55 Ma Challis volcanics and 80–100 Ma Atlanta lobe of the Idaho Batholith (Link et al., 2005; Figs. 3 and S2).

Methods of Intercomparison

We use several standard statistical metrics for detrital zircon provenance analysis that are described below. Provenance analyses were performed on the entire U-Pb zircon spectra (0–4000 Ma) as well as amended spectra, in which we removed all U-Pb ages younger than 17 Ma. The latter test was done to ensure that geographically broad aeolian transport of zircons from Yellowstone Hotspot eruption did not cause misidentification of fluvial connection. We find that removal of ages younger than 17 Ma does not change the intercomparison or modeling results or relative goodness of fit among zircon spectra substantially. Rather, it aided our ability to distinguish among source terrane zircon populations given the abundance of data sets that contain zircons younger than 17 Ma. We therefore report metrics of similarity and provenance using the amended spectra.

Metrics of Similarity

There are various statistical descriptors of detrital zircon dataset similarity, many of which have been evaluated for their ability to test sample similarity (Saylor and Sundell, 2016). Using the DZstats program and guidance provided by Saylor and Sundell (2016), we used five metrics of sample comparison: cross correlation from the kernel density estimate (KDE), likeness and similarity of probability density plot (PDP), as well as Kolmogorov-Smirnov (K-S) and Kuiper p-values calculated by repeated subsampling and comparison with subsampling n = 25 and 10 trials. To account for the impact of data set size (N) on likeness value (Satkoski et al., 2013), we report likeness values scaled by the number of dated zircons in each sample, N.

For each statistical metric, we chose a threshold value that signifies whether the comparison is positive (similarities in data sets) or negative (dissimilarity in data sets). For the K-S and Kuiper p-values, the threshold value is p<0.05 which indicates that two samples are not drawn from the same parent population at >95% confidence (Kuiper, 1960; Stephens, 1970; Press et al., 2007). For the cross-correlation, likeness, and similarity metrics, the threshold value choice is more arbitrary. To pick a meaningful value for the remaining three statistical metrics, we evaluated the cross-correlation, likeness, and similarity values among the Ringold Formation strata, which should be similar given that they are from the same geologic unit. We took the lowest values of Ringold detrital zircon intercomparison for each metric and rounded it down to derive a cutoff. For cross-correlation, likeness, and similarity, these cutoffs are 0.35, 0.60, and 0.60, respectively. Values of intercomparison above these thresholds were considered a “good” match. Given that all five statistical metrics used have their own intrinsic strengths and weaknesses (Saylor and Sundell, 2016), we consider a good match to require that at least four values suggest similarity between two compared spectra. We suggest this to be a relatively conservative approach to assessing whether detrital zircon data sets from ancestral rivers were derived from similar source terranes. Results of intercomparison statistics are in Table S3.

Multi-Dimensional Scaling

To graphically show intercomparison values, we use DZmds (Saylor et al., 2018) for multi-dimensional scaling (MDS). MDS constructs plots based on intercomparison metrics, in which samples that cluster together are more similar, whereas samples that are spread out are less similar (Vermeesch, 2013). For each metric of statistical comparison, a sample dissimilarity matrix and sample coordinates for MDS plots are calculated. To assess model fit, the coordinates are in turn used to calculate distance matrices, which can be plotted with the dissimilarity matrices to produce a graphical representation of fit called a Shepard plot. Stress is calculated as a metric of scatter on the Shepard plot, and values of <0.2 are considered a good fit (Vermeesch, 2013). Shepard plots with values >0.2 indicate that the MDS plots do not accurately represent the dissimilarity of the statistical metric (Fig. S3).

We compare detrital zircon data sets from sandstone samples collected from the Columbia Basin, western Snake River Plain, Rocky Mountain, and Columbia River Gorge to discern similarities between strata. For the available metrics, we found that the most consistent results and acceptable stress values were achieved by using Kuiper V and K-S D metrics of similarity (Fig. 4).

Statistical Results and Implications of Basin Interconnectivity

Quantitative comparison of detrital zircon age distributions using multiple statistical metrics for intercomparison (Saylor and Sundell, 2016; Table S3) and multi-dimensional scaling (MDS) (Saylor et al., 2018; Fig. 4) show that the Ringold Formation zircon age spectra are similar to those of the modern Snake River and the Rocky Mountains Sixmile Creek Formation (Stroup et al., 2008) (Table S3). MDS plots show that the two northern Ringold Formation samples (P1 and P3) are most similar to the Sixmile Creek Formation and ~Miocene Monida Pass fluvial gravels (Sears et al., 2009). The southern Ringold Formation sample (P4), on the other hand, is most similar to strata from the Columbia River Gorge, the upper Glenns Ferry Formation, and the Clarkston Heights gravels (Fig. 4). The Chalk Hills Formation is dissimilar to most other populations, including the Glenns Ferry Formation, which suggests a change of source terrane in the western Snake River Plain after ca. 4.3 Ma. Baker Valley detrital zircon samples are dissimilar from nearly all other sampled strata and modern rivers and do not contain zircon age populations indicative of Idaho or Rocky Mountain sources (Fig. 4; Table S3). Together, our quantitative assessment of the zircon age distributions supports a consistent upper Snake River presence in the Columbia Basin, but not through the western Snake River Plain, in the late Miocene–Pliocene. Similarities between the Ringold Formation, Sixmile Creek Formation, and Monida Pass gravels indicate that the upper Snake River took a more northerly route through western Montana and around the northern Rocky Mountains.

Unmixing modeling allows for a more sophisticated and nuanced approach to determining source contribution and is performed by comparing synthetic mixtures of modern river detrital zircon age spectra (Sundell and Saylor, 2017), for which we know the upstream contributing landscape, with the ancestral river sands. To determine the relative contribution of river networks to a single sample, we used the DZmix software (Sundell and Saylor, 2017). Modern river detrital zircon data sets were used as potential sources of Miocene–Pliocene samples (Fig. 3). We acknowledge the likelihood that modern river networks are not entirely equivalent to those of the past; however, the benefit of modern rivers is that we know the upstream drainage area, whereas upstream drainage area in past river systems can only be inferred.

For each analyzed ancestral river sample, we ran the unmixing models in DZmix with 50,000 trials of different relative river contributions and retained the top 1% of models that fit the data best (Figs. 5 and S4). All ancestral river samples are treated as individual paleo-river age spectra except for several geologic units with individual samples of limited zircon fertility. These include the Sixmile Creek Formation from Fries Place (P34–P36), the upper Glenns Ferry Formation (P7, 8, and 11), the lower Glenns Ferry Formation (P9–10), and western and eastern Chalk Hills samples (P12–13 and P14–15, respectively). Goodness of fit (GOF) between modeled and measured detrital zircon data sets was assessed using three metrics: cross-correlation, Kuiper V-statistic, and K-S D-statistic. We report averaged GOF values, which range between 0 and 1, with 1 being a perfect fit (Figs. 5 and S4). Pie charts show the relative river contributions for each best-fit DZmix model (Figs. 5 and S4).

Testing Hypothetical River Networks

Our proposed ancestral Snake River route is radically different from prior hypotheses (Fig. 2). We test our new hypothesis against two alternative ancestral Snake River routes (Livingston, 1928; Wheeler and Cook, 1954) using the detrital zircon unmixing models. The three proposed models include: (1) all rivers as they are today, including the full Snake River, called the “All-Snake” model. This model was designed to test whether the configuration of the current major rivers could replicate past river zircon spectra; (2) all modern rivers as they are today except for the Snake River, called the “No-Snake” model, which is designed to test the importance of the Snake River detrital zircon spectra; (3) the “East-Snake” model, which includes only the upstream reaches of the Snake River system, such as the South Fork Snake River (samples M30–31), Henrys Fork (sample M24), the East Lost Trough (samples M17–20), and the Portneuf River (sample M14). This model is designed to test whether Lake Idaho, seated in the western Snake River Plain, could have been isolated from lacustrine Ringold Formation strata for some time, and whether other portions of the Snake River system may have been routed differently to circumvent direct linkage between these two paleo-lake systems.

Estimating the Lake Idaho Drainage Area over Time

In a similar unmixing model approach, we analyze the provenance of western Snake River Plain strata using smaller tributaries of the Snake River as zircon sources (Fig. S5). In the western Snake River Plain unmixing models, we are not testing specific hypotheses, but rather are estimating the relative contribution of various portions of the Snake River Plain landscape to Lake Idaho over geologic time. For small tributaries, some modern sources were combined for simplicity based on the similarity of detrital zircon spectra and geographic proximity (Fig. S1).

Detrital Zircon Unmixing Results

Most model results show that the “East-Snake” model best replicates Miocene–Pliocene detrital zircon age spectra collected in the northern Rockies, Columbia Basin, and Columbia River Gorge (Figs. 5 and S4; Table 1). The northernmost Ringold Formation and Monida Pass samples show the largest improvement in model fit with the “East-Snake” model, whereas samples taken downstream of the White Bluffs and throughout the Columbia River Gorge have similar GOF values for the “East-Snake” and “All-Snake” models. Mixing results show a marked increase in Salmon and Clearwater River contribution required to match detrital zircon spectra from the White Bluffs (P4), in central Washington, and farther downstream (Fig. 5; P28–P29). An increase in Columbia, Methow, and Okanogan River contribution is required to match the westernmost Ringold Formation zircon age spectra (P2; Fig. 5). Several samples are essentially indifferent to the mixing model (P33, P37, P40, and P41) in that the GOF values for all models are similar. This may suggest that the Snake River is not important for replicating their detrital zircon spectra and therefore is not a likely contributing source. Importantly, the Clarkston Heights gravels (P40) that are along the modern Snake River, just downstream of Hells Canyon (Figs. 2, 3, and 5), are among these “indifferent” samples.

The unmixing models for western Snake River Plain stratigraphic units show a clear change in the provenance over time. Miocene strata, such as the Chalk Hills Formation (P12–P15), contain zircons that are locally sourced, such as those from the Boise River, Reynolds River, and southern margin of the western Snake River Plain (Fig. S5). Detrital zircon age spectra from the <4.3 Ma Glenns Ferry Formation samples require a contribution from the eastern Snake River Plain source terrane, which includes the eastern Lost Trough, Portneuf River, and a slight contribution from the upper Snake River tributaries (Fig. S5). However, the Terteling Springs Formation, which is considered a transgressive sequence at the base of the Glenns Ferry Formation (Wood and Clemens, 2002), is locally sourced from the Boise and Reynolds Rivers. This contrast with the distal sources of the remainder of the Glenns Ferry Formation is possible due to the far western and slightly isolated location of the sample or because the source contribution from distal river systems occurred slightly after deposition of the Terteling Springs Formation.

Based on the geographic and temporal distributions of strata analyzed for detrital zircon provenance and paleontology, we propose that the upper Snake River flowed into the Columbia Basin via an alternative route to the north that circumvented the western Snake River Plain (Fig. 2). New detrital zircon provenance results show that the upper Snake River was a major contributing source of the Ringold Formation, particularly north of the White Bluffs (Figs. 2, 3, and 5), and that the modern route through Hells Canyon postdates the youngest Ringold Formation sample (<3.5 Ma; P1). The increased contribution from the Clearwater and Salmon Rivers at and downstream of the White Bluffs suggests that the Snake River confluence with the Salmon and Clearwater Rivers was located near the White Bluffs during late Miocene–early Pliocene deposition of the Ringold Formation. Together these data suggest that the headwaters of the late Miocene–early Pliocene Snake River entered the Columbia Basin from the northeast (Figs. 6A–6B). While these results contrast with previously proposed routes (Livingston, 1928; Wheeler and Cook, 1954; Reidel and Tolan, 2013; Stearley and Smith, 2016), prior hypotheses assume that the upper Snake River flowed into Lake Idaho. Isolation of Lake Idaho from the upper Snake River is consistent with detrital zircon populations from western Snake River Plain strata showing no evidence for throughgoing hydrographic connection until after 4.3 Ma (Figs. 3 and S5). A northerly Snake River route bypassing the western Snake River Plain also could explain upper Snake River zircons within the Ringold Formation during a time when western Snake River Plain fossil fish were isolated from the Columbia Basin (Smith et al., 2000; Stearley and Smith, 2016).

Detrital zircon results from Miocene fluvial gravels on Monida Pass over the Tendoy Range and the Sixmile Creek Formation in the Ruby Graben (P34–P36, P38; Fig. 2) support the hypothesized northerly paleo-Snake River pathway. Unmixing models show major upper Snake River zircon contribution in these strata (Fig. 5), and MDS plots and sample intercomparison statistics suggest similarity with the Ringold Formation and therefore a similar source terrane (Fig. 3; Table S3). Furthermore, the northeast-flowing 6.0 ± 0.1 Ma Timber Hill Basalt was erupted in the east Snake River Plain and is interbedded with the Sixmile Creek Formation in the Ruby Graben, which is likewise indicative of a northward throughgoing route across the Tendoy Range (Kreps et al., 1992; Sears et al., 2009). Distinctive clast lithology within the Monida Pass gravels may match sources south of the eastern Snake River Plain (Sears, 2014). While the Monida Pass gravels and Sixmile Creek Formation were suggested to have been part of the Miocene Bell River and sourced from the Colorado Plateau (Sears, 2013), a lack of Mesozoic (ca. 150–250 Ma) arc-derived zircons typical of sediment sourced from the Mogollan highlands (Laskowski et al., 2013) does not support this hypothesis (Fig. 3). Basalt sourced from the eastern Snake River Plain continued to flow northward until 5.5–5.0 Ma (Fig. 6C), subsequent volcanic flows were confined to the eastern Snake River Plain (Fritz and Sears, 1993), and Ruby Graben strata ceased to be derived from central Idaho (Sears et al., 2009). Since 4.2 Ma, the eastern Snake River Plain has subsided 745 m (McQuarrie and Rodgers, 1998), which is nearly identical to the amount of relief between the modern Snake River and Monida Pass (2100 m asl) that now separates the Snake River Plain from the eastern Missouri River Basin (Fig. 2). Our proposed route for the paleo-Snake River again crosses the modern continental divide back into the Columbia Basin over the lower (1800 m asl) Deer Lodge Pass. While the modern topography profile suggests ~250–300 m of relief between the modern valley floor and Deer Lodge Pass, several mechanisms may be responsible for post-early Miocene incision and relief generation in southwestern Montana. For example, the Ruby Graben accumulated up to 300 m of alluvial fan and fluvial strata before eruption of the 6 Ma Timber Hill basalt (Sears and Ryan, 2003). Other Miocene graben systems nearby, such as the Big Hole and Deer Lodge grabens, hosted deep lake systems that accumulated several kilometers of lacustrine sediment (McLeod, 1987). The aggradation and Miocene lake systems within the SW Montana grabens suggest the base level has changed since then. By 2 Ma, the Miocene graben system was crosscut by northwest-trending normal faults related to passage of the Yellowstone Hotspot, which segmented the earlier grabens and diverted drainage systems (Fritz and Sears, 1993). Since 2 Ma, considerable incision has lowered rivers in the area based on the modern stream elevations relative to the 2 Ma Huckleberry tuff (Sears and Ryan, 2003). Together, these data sets suggest that the Snake River was redirected from its route across the Tendoy Range between 5.0 Ma and 4.2 Ma and that the continental divide has shifted westward since then (Fig. 6C).

The 10.1–6.4 Ma Chalk Hills Formation records an early Lake Idaho history (Wood and Clemens, 2002) with sediment locally sourced from north and west of the lake system (Guzzo and Link, 2018; Fig. S5) and does not support the input of a large and distally sourced river, such as the Snake River, into the late Miocene–early Pliocene Lake Idaho. The similarity of Chalk Hills fossils with modern aquatic biota in the Sacramento River network has been used to support a throughgoing fluvial connection between the Lake Idaho and Sacramento or Klamath River drainages in late Miocene–early Pliocene time (Stearley and Smith, 2016). Together, these data indicate that Lake Idaho was isolated from the Snake River and Columbia Basin between 10.1 Ma and 6.4 Ma but that it may have had a southern outlet (Fig. 6A). Whether the southern outlet was sustained or ephemeral, and whether the outlet reached the Pacific Ocean, is not yet clear.

The stratigraphy, paleontology, and provenance of the 4.3–2.5 Ma Glenns Ferry Formation signal a major change in the Snake River drainage network and late Pliocene incision of its modern outlet, Hells Canyon. Fish and rodent fossils in the Glenns Ferry and upper Ringold Formations are nearly identical and indicate a newly established fluvial connection between the western Snake River Plain and Columbia Basin after 3.4–3.0 Ma (Smith et al., 2000). Additionally, the appearance of an anadromous salmon (Stearley and Smith, 2016) in Lake Idaho supports a connection between the western Snake River Plain and the Pacific Ocean being established by 3.3–2.5 Ma (Fig. 6D). Detrital zircon unmixing models of the western Snake River Plain strata clearly show an increase in distally sourced zircons from the eastern Snake River Plain in the Glenns Ferry Formation (Fig. S5). The upper Glenns Ferry Formation stratigraphy of progressive cut and fill terraces and the prograding deltaic sedimentation of the Pierce Gulch Sand into Lake Idaho, with a maximum depositional age of 2.5 Ma, indicate gradual base level lowering (Sadler and Link, 1996; Link et al., 2002; Wood and Clemens, 2002).

The timing and location of volcanic centers and major volcanic episodes of the Yellowstone Hotspot match remarkably well with detrital zircon provenance, sedimentological, and paleontological change through the Miocene and Pliocene, which suggests a causal link between hotspot-related topography and the fluvial drainage reorganization. The timing of major volcanism and structural disruption around ca. 16.6 Ma signals the emergence of the Yellowstone Hotspot plume through North American continental crust that may have produced the topographic lows such as the Ruby Graben (Sears and Thomas, 2007) through which we propose the paleo-Snake River eventually flowed. As the Yellowstone Hotspot tracked eastward through North America, the Picabo volcanic center of the Yellowstone Hotspot was centered in the central Snake River Plain between 10.4 Ma and 6.6 Ma (Fig. 2). Based on the similar depositional age of the 10.1–6.4 Ma Chalk Hills Formation in Lake Idaho (Wood and Clemens, 2002), we suggest that thermally induced topography from the Yellowstone Hotspot and eastern Snake River Plain subsidence due to densification of the lithosphere (McQuarrie and Rodgers, 1998) isolated the western Snake River Plain from the upper Snake River and eastern Snake River Plain (Figs. 6A–6B). An unconformity between the Chalk Hills and Glenns Ferry Formations between 6.4 Ma and 4.3 Ma (Middleton et al., 1985; Wood and Clemens, 2002) matches the 6.6–4.5 Ma age of the Heise volcanic center of the Yellowstone Hotspot (Watts et al., 2011; Fig. 2). This overlap suggests that the eastward migration of the Yellowstone Hotspot barred major fluvial input into the western Snake River Plain and caused Lake Idaho to evaporate, drain, or recede substantially from its former extent (Fig. 6B). Lake Idaho lacustrine sedimentation resumed during a lull in Yellowstone Hotspot volcanic activity between the 4.5 Ma Kilgore tuff (Watts et al., 2011) and 2.1 Ma Huckleberry Ridge tuff (Lanphere et al., 2002), which is marked by the 4.3–2.5 Ma Glenns Ferry Formation (Fritz and Sears, 1993). The appearance of east Snake River Plain-sourced zircons in the Glenns Ferry Formation suggests that Yellowstone Hotspot-related topography that previously barred integration of the western and eastern Snake River Plain subsided by 4.3 Ma (Fig. 6C).

Recent interpretation of aquatic biota from the 4.5–7.5 Ma Cache Valley member of the Salt Lake Formation and comparison with the Chalk Hills and Glenns Ferry biota calls into question the longevity of hotspot-related topography and its impact on basin isolation (McClellan and Smith, 2020). McClellan and Smith (2020) suggest that rapid collapse of the dynamic topography from the Heise volcanic center is evidence that the Yellowstone Hotspot-related topography is characteristically short lived at other volcanic centers, such as Picabo or Twin Falls, and therefore is not a factor in basin isolation. However, the Heise volcanic center may be unique by being at the locus of subsidence in the eastern Snake River Plain, where hotspot-related topography may diminish more rapidly due to a combination of thermal relaxation and densification of the lithosphere. The additional mechanism of eastern Snake River Plain subsidence is either absent or minor for Yellowstone volcanic centers farther to the west.

We can reconcile these data sets by comparing the fossil and zircon data sets and paying close attention to the time transgressive Yellowstone Hotspot locations relative to Cache Valley (Fig. 2). Several fish species in Cache Valley have morphological differences to those found in the Chalk Hills and are therefore divergent (McClellan and Smith, 2020), which suggests a period of isolation between the Chalk Hills and Cache Valley biota. This time period coincides with the Yellowstone Hotspot Picabo volcanic center (10.4–6.6 Ma) that is located to the northwest of Cache Valley. The Glenns Ferry fossil assemblages more closely resemble Cache Valley species. This is similar to our detrital zircon provenance results, which show an increased contribution from eastern sources, such as the Portneuf River, located east of the Picabo center but west of the Heise volcanic center (6.6–4.5 Ma). Cache Valley is close to the Portneuf River (Fig. 2), and thus its connection to Lake Idaho during Glenns Ferry time (<4.3 Ma) is consistent with provenance analysis.

The large increase in drainage area into the western Snake River Plain after 4.3 Ma likely increased the volume and elevation of Lake Idaho and led to its eventual demise by overflow and incision of its outlet. We suggest that the coincident 3.3–2.5 Ma stratigraphic change in Glenns Ferry Formation sedimentation (Wood and Clemens, 2002) is evidence that the overfilled Lake Idaho overtopped its bedrock threshold, or sill, and initiated the incision of Hells Canyon (Fig. 6D).

New detrital U-Pb zircon data and analyses clarify how and when the western Snake River Plain was incorporated into the greater Snake and Columbia River drainage system and constrain the timing of Hells Canyon incision between 3.4 Ma and 2.0 Ma. Previous work on the fossil aquatic biota assemblages from western Snake River Plain and Columbia Basin lacustrine strata provided important groundwork for identifying changes in major river systems during the late Cenozoic; however, detrital U-Pb zircon analyses and modeling uniquely identify past fluvial interconnectivity and pinpoint specific river connections and pathways over time.

We propose a northerly route of the late Miocene–early Pliocene Snake River that satisfies provenance, paleontological, and stratigraphic data sets (Fig. 2). Synthesis of these diverse data sets regionally allows us to correlate major events in the fluvial record, including incision of Hells Canyon, with the volcanic history of the Yellowstone Hotspot. We identify thermally induced topography of the Yellowstone Hotspot plume as a mechanism for eastward drainage divide migration, which—with magmatic densification of the eastern Snake River Plain lithosphere—was responsible for drainage network reorganization and long-term landscape evolution of the U.S. Pacific Northwest and Rocky Mountains.

Critical input and discussions with Brittany Guzzo, Spencer Wood, Kurt Sundell, Nick Zentner, Matthew Morriss, Alison Duvall, Brian Yanites, and Gene Humphreys greatly improved our approach, analysis, and interpretations. Permission to sample on private land was graciously granted by Carl Hurlburt. We thank Nora Nieminski, Spencer Wood and John Shervais for reviewing this work and providing critical feedback for the improvement of the manuscript. Funding was provided by the U.S. Geological Survey (USGS), College of Idaho Museum, and Idaho State University College of Science and Engineering. Author contributions are as follows. L.M. Staisch: conceptualization, methodology, software, validation, formal analysis, investigation, data curation, writing of the original draft, visualization, project administration, and funding acquisition. J.E. O’Connor: conceptualization, investigation, writing—review and editing, project administration, and funding acquisition. C.M. Cannon: investigation, writing—review and editing. C. Holm-Denoma: methodology, resources, writing—review and editing. J. Lasher: investigation, resources. P.K. Link: Chalk Hills investigation, review and editing, funding acquisition. J.A. Alexander: investigation, resources. Supplementary information, figures, and analytical and statistical results are available at the USGS ScienceBase repository https://doi.org/10.5066/P9T2K62G. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

1Supplemental Material. Figure S1: Detrital zircon age spectra from modern rivers; Figure S2: Detrital zircon age spectra from fluvial and lacustrine sandstones; Figure S3: Shepard plots from Multi-Dimensional scaling (MDS) analysis comparing distance and disparity for four metrics of detrital zircon similarity; Figure S4: DZmix results for three hypothesized river networks; Figure S5: SRP sample location map and detrital unmixing results; Table S1: Modern and ancestral river detrital zircon sample locations, ages, and references; Table S2: U-Pb zircon age results for new modern and ancestral river sands; Table S3: Intercomparison results between modern and ancestral river sediments; Table S4: Best-fit DZmix results estimating the relative contribution of hypothesized sources to measured detrital zircon age spectra of ancestral river sands; Table S5: Best-fit DZMix results that estimate the relative contribution of Snake River Plain tributaries to Miocene-Pliocene Lake Idaho strata. Please visit https://doi.org/10.1130/GSAB.S.16654777 to access the supplemental material, and contact editing@geosociety.org with any questions.
Science Editor: Rob Strachan
Associate Editor: Eric Roberts
Gold Open Access: This paper is published under the terms of the CC-BY license.