Determining the source of placer gold in the Anaconda metamorphic core complex supradetachment basin using detrital zircon U-Pb geochronology, western Montana, USA

Placer deposits are broadly defined as any deposits of sand, gravel, and other detritus that contain accumulations of economically valuable minerals (Yeend and Shawe, 1989). Gold placer deposits are notably important deposits that result from weathering and release of gold from a bedrock source, gold transportation, and mechanical concentration of gold in streams and river gravels (e.g., Boyle, 1979; Loen, 1986; McCulloch et al., 2003). It is estimated that approximately two-thirds of the total world gold supply—and roughly half of the gold mined in California, Alaska, Idaho, and Montana (western USA)—has been produced from placer deposits (Fig. 1; Boyle, 1979; Yeend and Shawe, 1989). These deposits therefore play an essential role in the discovery and production of gold worldwide.

Because placer deposits are formed by normal surficial processes (such as fluvial and glacial erosion) that act on gold-bearing bedrock of varying richness, their geographic distribution is wide and their sizes and concentrations vary greatly (Fig. 1). These deposits occur predominately in Cenozoic and Quaternary rocks due largely to the destruction of older placers by erosion (Edwards and Atkinson, 1986). Additionally, the preservation of placer deposits in young rock units may reflect the role that recent tectonic and climatic events play in their formation (e.g., Roy et al., 2018). For example, the relatively young depositional ages (Oligocene and younger; Elliott et al., 1992) of many gold placers in Montana can be attributed to the shedding of detrital gold from the recently exposed gold-bearing core of the northern Cordilleran hinterland (Fig. 1).

The concentration and size of gold placer deposits depend almost entirely on the supply of source materials (Yeend and Shawe, 1989), and conducting provenance studies on the gold-bearing sediments can provide insight into the original bedrock source and/or rock units that hosted vein or skarn gold. Placer gold source regions that are actively shedding gold are of particular interest because their identification could allow for extraction of gold directly from the lode source. Additionally, it has been established that an effective technique to discover new placer deposits is through the identification of potential source areas (Edwards and Atkinson, 1986).

Most studies interested in gold placer provenance have investigated gold nugget morphology, surface texture, and bulk chemistry to determine an approximate transport distance from the source (e.g., Loen, 1986, 1995). Knight et al. (1999) conducted a study of gold particle shape and rim characteristics to determine the distance of fluvial transport of placer deposits in the Klondike District of Canada, concluding that gold morphology “shows a smooth, well-defined relationship to distance of transport from the lode source” (p. 635). Other techniques that have been used to determine placer gold source include clast type analysis and measurement of heavy mineral concentrations from within placer deposits (e.g., Loen, 1994). Although these techniques may give adequate preliminary insight into placer source characteristics, the development of more quantitative techniques to investigate provenance—specifically detrital zircon (DZ) geochronology—provides the opportunity to conduct a more rigorous investigation of placer gold source.

Detrital zircon geochronology has become an essential tool in the study of sediment provenance because of the ubiquity of zircon in most depositional systems and the increasing ability to determine U-Pb ages with reasonable efficiency, accuracy, and precision (e.g., Gehrels, 2014). Formed primarily in felsic igneous rocks, zircon is a heavy and resistant mineral that does not commonly break down when weathered into sedimentary systems; as a result, zircons can be recycled multiple times and be most recently sourced from strata that do not represent their initial depositional unit. It is for this reason that many detrital zircons may give insight into an original source, but not necessarily a proximal one (e.g., Thomas, 2011; Schwartz et al., 2019). Conversely, the weathering of metamorphic and igneous rocks can provide first-cycle grains that can allow direct interpretation of provenance. Although DZ U-Pb geochronology is an established technique used to determine the provenance of sedimentary rocks and sediments, it is not routinely applied to gold placer deposits. Davis et al. (1994) first proposed that geochronologic analysis of detrital zircons within placer deposits provides U-Pb ages that can be correlated to ages of surrounding igneous and sedimentary rocks, giving insight into potential source areas of the gold-bearing sediments. Three studies in the Witwatersrand basin of South Africa used U-Pb geochronology of detrital zircons to gain insight into the controversial origin of gold deposited in quartzite reefs (Ruiz et al., 2006; Koglin et al., 2010; Zeh and Gerdes, 2012). The technique has also been applied to fingerprint to source of gem and gold placers in the Mamfe Basin of southwestern Cameroon (Kanouo et al., 2012, 2018). More recently, researchers have also begun pairing the U-Pb ages of detrital zircons with corresponding geochemical signatures (such as Lu-Hf isotopic data), which serve as separate provenance tracers for mixed sediments (Zeh and Gerdes, 2012; Kanouo et al., 2018).

The Pioneer District placer deposits were the site of the first gold discovery in Montana in 1852 (Pardee, 1951), and their hypothesized proximity to their lode source and complex surrounding structural framework make them an ideal candidate for a test of the utility of DZ geochronology in determining the sources of zircons in placer deposits. In this research, we evaluate the provenance of detrital zircon samples from the Pioneer District deposits to assess whether the gold might have been derived from the same sources, thereby providing a case study of the use of DZ geochronology applied to placers.

There is debate surrounding whether the Pioneer District placer gold originated from a vein or skarn lode source, and the source location remains unknown (Pardee, 1951; Loen, 1986; McCulloch et al., 2003). Extracted primarily from Pleistocene till and alluvium, it is possible that the placer gold was sourced from the nearby Late Cretaceous Royal stock, initially concentrated in Oligocene Cabbage Patch Formation and Miocene Squaw Gulch conglomerate beds, and reworked during Pleistocene glaciation. We test this hypothesis by comparing new detrital zircon U-Pb ages (n = 1058) from the Pioneer District placer deposits with the age spectra of plausible source units that may have been in contact with the Royal stock. We also present new detrital zircon Lu-Hf (n = 61) isotopic data, which serve as an independent provenance indicator, from four samples within the Pioneer District placer deposits. Additionally, we evaluate the DZ U-Pb data using an inverse Monte Carlo unmixing model that, by comparing the mixed placer age spectra to the DZ spectra of potential sources, calculates probable relative contributions from each input source (Sundell and Saylor, 2017).

The Pioneer District of southwestern Montana, straddling the northern Flint Creek Range and the Deer Lodge Valley, is a site of extremely rich placer deposits from which ∼8500 kg (∼300,000 oz) of gold was recovered (Fig. 2; Loen, 1986). Mined between 1870 and 1986, the gold was extracted primarily from Pliocene fossil alluvial placers and Pleistocene glacial till and alluvium (Pardee, 1951; Loen, 1986).

Located along the eastern edge of the Cordilleran hinterland of western Montana (Fig. 2), the Flint Creek Range represents the northern footwall of the Anaconda metamorphic core complex (AMCC) (O’Neill et al., 2004; Foster et al., 2010). The AMCC exhumes metamorphosed Cretaceous–Paleocene plutonic rocks and Mesoproterozoic–Mesozoic sedimentary rocks from depths of ∼12 km, based on thermobarometry data (Grice, 2006). The metamorphic-plutonic footwall of the AMCC is separated from largely unconsolidated, synextensional sedimentary rocks of the hanging wall by the low-angle Anaconda detachment fault (O’Neill et al., 2004; Foster et al., 2010).

The northern footwall of the AMCC, proximal to the Pioneer District placers, is composed of folded Mesoproterozoic–Mesozoic metasedimentary and sedimentary rocks that have been intruded by several granite and granodiorite plutons (Fig. 2; Emmons and Calkins, 1913; Grice, 2006; Portner et al., 2011). The dominant plutons include the 69–60 Ma Royal stock and the 65 Ma two-mica (muscovite and biotite) Mount Powell batholith (Fig. 2; Marvin et al., 1989; Grice, 2006), both of which are exposed near the Anaconda detachment. Today, the Royal stock is located at ∼2.5 km above sea level in the Flint Creek Range, ∼8 km south of the Pioneer District placer deposits (Fig. 2). Many of the strata now exposed in the hanging wall of the AMCC are interpreted to have been deposited in an Eocene–Oligocene supradetachment basin (the “Flint Creek Basin”) that developed contemporaneous with slip along the Anaconda detachment (Janecke et al., 2005; Stroup et al., 2008). The late Oligocene–early Miocene Cabbage Patch Formation and mid-Miocene Squaw Gulch conglomerate beds were deposited atop the hanging wall of the Anaconda detachment and have a combined thickness of ∼800 m, constituting most of the Flint Creek Basin fill (Fields et al., 1985; Loen, 1986; Stroup et al., 2008; Portner et al., 2011). It has been hypothesized that the Cabbage Patch Formation and Squaw Gulch beds were the original depositional units of the Pioneer District placer deposits prior to subsequent glacial transport and redeposition as fluvial placers in the Pleistocene (Loen, 1986).

No researchers have reevaluated the Pioneer District gold placers since the AMCC was recognized as a Cordilleran metamorphic core complex by O’Neill et al. (2004). The discovery that the Flint Creek Range represents the footwall of a core complex adds an interesting tectonic element to the provenance study and a potential mechanism for exhumation of the gold-bearing unit(s). Foster et al. (2010) used biotite and muscovite ⁴⁰Ar/³⁹Ar thermochronology to determine that displacement on the Anaconda detachment initiated ca. 53 Ma and lasted until at least 38 Ma. Zircon fission-track ages of ca. 27 Ma (Foster et al., 2010) and mapped field relationships (Howlett et al., 2019) suggest that slip could have occurred for much longer than previously thought, into the Oligocene. This raises the possibility that the Pioneer District gold placers, if initially concentrated in Oligocene–Miocene Flint Creek Basin conglomerates, were shed from the AMCC footwall during active extension.

Pardee (1951) first concluded that after uplift of the Flint Creek Range, ancestral rivers excavated gold from its lode source and deposited it in “river gravels”, later to be destroyed by Pleistocene glaciations and redeposited further downslope. Loen (1986, 1994) expanded on the work of Pardee (1951) and concluded that the lode source of the gold underwent initial erosion in the late Oligocene but was mostly concentrated in the Miocene Squaw Gulch beds. After Miocene deposition, placer gold was redistributed by Pliocene rivers and Pleistocene glaciation (Loen, 1986). Loen (1994) used gold morphology and gold composition to conclude that the placer gold was likely sourced from mineralized veins, skarn, or replacement deposits associated with the Royal stock of the northern Flint Creek Range. Although no bedrock evidence of the source has been discovered, the most detailed geologic map of the region, produced by Mutch and McGill (1962), presents a mapped contact between the Royal stock and Mesoproterozoic–Mesozoic supracrustal rocks, which supports the hypothesis that the northern Flint Creek Range is a plausible source region for the placer zircons (Loen, 1994). McCulloch et al. (2003) concluded that the Pioneer District placer gold originated from skarn deposits by analyzing gold fineness, which is a measure of the proportion of gold in a gold-silver alloy expressed in parts per thousand (e.g., Loen, 1986). Fineness is commonly used in the determination of placer gold source because different lode types (e.g., skarn versus vein) have characteristic finenesses determined by their formation conditions (McCulloch et al., 2003).

Several studies interested in the sedimentation histories of extensional basins in western Montana conducted provenance analysis of the Paleogene–Neogene rocks exposed in the Flint Creek Basin (e.g., Stroup et al., 2008; Portner et al., 2011). These studies were focused specifically on Eocene–Miocene basin development and did not consider the history of gold placer deposition, but they provide valuable DZ and sedimentological data for the Cabbage Patch Formation and other relevant strata (such as the Barnes Creek beds; Portner et al., 2011). Detrital zircon U-Pb geochronology results from a Cabbage Patch Formation sandstone include data from abundant grains (n = 26) that are Late Cretaceous in age (ca. 75 Ma) (Stroup et al., 2008), and paleocurrent analysis of Oligocene two-mica sands of the Flint Creek Basin indicate northwest-directed paleoflows (Portner et al., 2011). Paleoflow indicators in the late Miocene Barnes Creek gravel beds also record sediment derivation from the south-southwest, suggesting that sediment dispersal patterns in the Flint Creek Basin remained broadly consistent during the Cenozoic (Portner et al., 2011). This interpretation is also supported by the orientations of Pleistocene lateral moraines emanating from the northern Flint Creek Range, which are oriented predominantly northeast-southwest (e.g., Loen, 1986). Generally speaking, results from DZ geochronology, sandstone petrography, and paleocurrent analysis in the Flint Creek Basin are consistent with sediment derivation from the northern AMCC footwall (Stroup et al., 2008; Portner et al., 2011).

To constrain zircon provenance and assess whether the gold was derived from the same sources, four detrital zircon geochronology samples (101517AL1, 101517AL2, 101517AL3, and 101517AL4) were collected from the reworked placer deposits of the Pioneer District. Analysis of published maps from the Pioneer District (Pardee, 1951; Loen, 1986) and satellite imagery allowed for previously mined placer deposits to be located and targeting for sampling. Bulk samples were collected from reworked and unconsolidated sediments and placer tailings with sediment size ranging from medium sand to cobbles. We assume that density separation techniques used to mine placer gold did not bias DZ yield because the equipment is designed in such a way to allow for the capture materials of ultrahigh density (i.e., gold with density of ∼19.30 g/cm³) while those of lesser density (such as zircon with density of 4.7 g/cm³) are washed through the apparatus (Silva, 1986; McCulloch et al., 2003). Samples were collected from three different modern drainages over a lateral swath of ∼6 km to ensure representative age spectra (Fig. 2).

Detrital zircon samples were prepared and analyzed using protocols consistent with those of the Arizona LaserChron Center (Tucson, Arizona, USA; http://www.laserchron.org) (Gehrels et al., 2008). Zircons were separated from ∼2.0 L bulk samples by pulverization in a jaw crusher, sieving, magnetic separation, density separation, and hand-picking. Zircons were mounted in epoxy and polished to a depth of ∼30 µm, and backscattered-electron images were obtained using a scanning electron microscope (SEM) for targeting during analysis by laser ablation–inductively coupled plasma–mass spectrometry (LA-ICP-MS).

Zircon U-Pb ages were obtained for ∼300 zircon grains per sample using a Photon Machines Analyte G2 Excimer laser (30 µm beam diameter) attached to a Thermo Element2 HR single-collector ICP-MS. The laser ablation process excavates pits that are ∼15 µm in depth, and the ablated material is subsequently transported in helium to the plasma source of an Element2 ICP-MS. U, Th, and Pb isotopes were measured simultaneously using Faraday detectors with 3 × 10¹¹ Ω resistors for ²³⁸U, ²³²Th, and ²⁰⁸Pb-²⁰⁶Pb, and using discrete dynode ion counters for ²⁰⁴(Pb + Hg) and ²⁰²Hg. Each analysis consisted of one 15 s integration on peaks with the laser off, 15 one-second integrations with the laser firing, and a 30 s delay to purge for the next analysis. All new analytical data and a detailed list of concordance filters are reported in Item S1 in the Supplemental Material¹. Kernel density estimate (KDE) plots were generated using the Python-based detritalPy toolset (Sharman et al., 2018).

Following U-Pb isotopic analysis, Hf isotope measurements were made for 61 grains <250 Ma using a Photon Machines Analyte G2 Excimer laser (40 µm beam diameter) attached to a Nu Plasma multicollector ICP-MS at the Arizona LaserChron Center. Grains <250 Ma were chosen due to the relatively well-constrained geochemistry of igneous rocks related to the development of the North American Cordillera. An average of 15 analyses were conducted for each sample, with measurements made from the same sample spots as LA-ICPMS U-Pb analysis to ensure that Hf isotopic data were determined from the same domain as the U-Pb age. Fragments of zircon standards MT, FC, SL, 91500, TEM, PLES, and R33 were used for standard sample bracketing during Lu-Hf analyses. Hf analyses are reported alongside detrital zircon U-Pb data in Item S1 (see ^{footnote 1}).

The North American Cordilleran hinterland consists of a shortened, thickened, metamorphosed, and magmatically infiltrated package of sedimentary rocks that contains an extremely wide range of zircon ages (e.g., DeCelles, 2004; Dickinson and Gehrels, 2008; Laskowski et al., 2013). The resulting DZ age spectra obtained from syn- and post-orogenic sedimentary deposits are very complex. Despite the complexities and implications of zircon recycling, U-Pb age spectra are still a powerful tool for determining provenance, and new modeling techniques and a number of first-cycle source areas with unique ages allow for direct interpretation of provenance. Additionally, analysis of Hf isotopes in DZs of a known age has emerged as a powerful tool for provenance analysis because it serves as an independent source indicator when paired with detrital U-Pb geochronology (Goodge and Vervoort, 2006; Stroup et al., 2008).

A list of potential source rocks and their characteristic DZ ages for the Pioneer District placer deposits are listed in Table 1. Potential sources include: (1) major volcanic fields of the northern Cordillera (Dillon, Absaroka, Challis, Lowland Creek, and Elkhorn Mountains volcanics) with zircon ages ranging from ca. 17 Ma to ca. 83 Ma; (2) voluminous plutonic rocks of Montana and Idaho (Pioneer, Idaho, Philipsburg, and Boulder batholiths) with major ages ranging from 50 Ma to 98 Ma; (3) proximal plutonic rocks in footwall of the AMCC (Mount Powell pluton and Royal stock) with ages between 60 Ma and 69 Ma; and (4) sedimentary strata of the Montana–Idaho segment of the Sevier fold-thrust belt, including but not limited to the Belt Supergroup, Paleozoic passive-margin strata, and the Cretaceous Kootenai and Blackleaf Formations and Beaverhead Group with DZ ages ranging from ca. 3.6 Ga to ca. 90 Ma (see Table 1 for references).

An unmixing algorithm determined the mixing of source samples through inverse Monte Carlo modeling by comparing detrital samples to randomly generated combinations of source distributions (“DZMix”; Sundell and Saylor, 2017). As a forward-modeling technique, detrital and igneous zircon age data for potential source units were input into the model and randomly combined to create a model age spectrum in the form of a KDE (Fig. 3A). Subsequently, a cross-correlation coefficient was calculated for the model and mixed DZ sample pair (Fig. 3B). This process was repeated 10,000 times to increase the likelihood that all possible combinations of source proportions are tested. The mean relative contributions and 1σ uncertainty of each provenance group were reported from the top 1% of trials ranked by the cross-correlation coefficient (Fig. 3C). In addition to using the cross-correlation coefficient comparison, the random mixtures of potential sources were compared by calculating a Kuiper test V value and Kolmogorov-Smirnov test D value (Sundell and Saylor, 2017). This allowed for a side-by-side comparison of relative source contributions calculated using different statistical algorithms. After running the inverse Monte Carlo model, it was possible to run an iterative optimization that determined a single best fit between the model and the mixed sample.

Potential source units that are hypothesized to have contributed sediment to the placer deposits were identified through a combination of fieldwork around the margins of the Royal stock and analysis of previously published regional geologic maps (this study; Mutch and McGill, 1962; Loen, 1986; Lewis, 1998). Nine distinct units were identified, including (from oldest to youngest): the Mesoproterozoic Ravalli Group (mixed metasedimentary rocks), the Paleozoic Amsden (limestone), Quadrant (quartzite), and Phosphoria (shale) Formations, the Lower Cretaceous Kootenai Formation (silt and sandstone), the Colorado Group (Thermopolis shale), and Blackleaf (shale) Formation, the Late Cretaceous Pioneer batholith (Royal stock age-equivalent), and the Upper Cretaceous–Paleocene Beaverhead Group (conglomerate). Additionally, a Cordilleran passive-margin compilation of DZ data from Laskowski et al. (2013) was included as a model input due to the mapped presence of the Mississippian Madison Formation (Mutch and McGill, 1962) and the high likelihood that passive-margin strata were previously eroded from upland areas and deposited in the Flint Creek Basin (e.g., Portner et al., 2011). All units (excluding the Beaverhead Group) are currently exposed in the northern footwall of the AMCC; we therefore hypothesize that they were also present during development of the placer deposits. It is also likely that some of the placer sediments were derived from sedimentary overburden that was unroofed during development of the AMCC. The compilation of potential sources for the placer deposits was narrowed down further by hundreds of trial runs in DZMix with differing model inputs. By presenting the randomly generated KDEs atop the KDE of the mixed sink sample, the DZMix output allows the user to identify age ranges that were missed by each combination of input sources. A process of trial and error was used to create a randomly generated KDE that most closely approximated the mixed sample.

For the final model run, new DZ distributions for the Pioneer District placer deposits (samples 101517AL1, 101517AL2, 101517AL3, and 101517AL4; Fig. 2) were combined to create a single sink input (n = 1058). The combination of the sink samples into one large-n data set for modeling purposes is justified by the overall similarity in age spectra between samples. In addition to combining the sink samples, we compiled the 10 potential source units listed above into bins based their respective geologic ages (Proterozoic, Paleozoic, Mesozoic, and Paleocene). This was done primarily to simplify and refine the modeling outputs because inputting high numbers of individual sources (i.e., all 10) commonly led to similar, small contributions from each. Following the consolidation of data in this way, the combined sink sample was compared to random combinations of the 10 potential sources in their respective temporal bins.

This forward-modeling approach to unmixing detrital zircon samples is superior to a strictly qualitative approach to provenance analysis in that it provides approximate percent contributions from each input source. In a single model run, input sources that are not being incorporated by the model can be generally considered implausible.

All four DZ samples contain Cordilleran grains <250 Ma (the youngest of which are ca. 25 Ma) as well as abundant grains (n >200 in each sample) that are >250 Ma. U-Pb ages from all major North American crustal provinces are present in each sample, with age-probability peaks centered at ca. 1850, ca. 1750, ca. 1640, ca. 1490, ca. 1100, and ca. 430 Ma (Fig. 4). Each sample contains ∼20 grains that are >2500 Ma.

Sample 101517AL1 was collected from processed placer tailings on the west side of the Pioneer Bar, downstream from an inactive dredge. The sample is a poorly sorted, matrix-supported pebble to cobble conglomerate, with clasts of quartzite, sandstone, and slate, in order of decreasing abundance. Zircon U-Pb ages (n = 261) range from 24.9 ± 0.3 to 3293.2 ± 8.4 Ma, with 181 of those ages (69%) being >250 Ma (Fig. 4B1). Significant peaks for grains <250 Ma in this sample occur at ca. 165, ca. 75, and ca. 25 Ma. Lu-Hf isotopic ratios were determined for 22 grains ranging in age from 25 to 171 Ma (Fig. 5). Epsilon-Hf (ε_Hf) values decrease from −2.8 at 171 Ma to −19.3 at 67 Ma, then increase to as high as 13.2 at 26 Ma.

Sample 101517AL2 was collected from processed placer tailings downstream from an inactive dredge, and is a poorly sorted, matrix-supported pebble to cobble conglomerate. In order of decreasing abundance, clasts include granite, sandstone, siltstone, quartzite, and schist with bronze mica. U-Pb ages (n = 265) range from 24.1 ± 0.2 to 3635.5 ± 7.6 Ma, with 204 of those ages (77%) >250 Ma (Fig. 4B2). This sample contains the most variability in zircon ages <250 Ma, with significant peaks at ca. 165, ca. 115, ca. 100, ca. 90, ca. 69, and ca. 25 Ma (Fig. 4B2). This is the only sample collected from the placer deposits that contains a prominent age-probability peak at 90–100 Ma (Fig. 4B2). Lu-Hf isotopic ratios were analyzed for zircons (n = 13) ranging in age from 24.9 to 167.8 Ma (Fig. 5). Epsilon-Hf values vary from −9.6 calculated at 167.8 Ma to –15.3 at 77 Ma. Similar to those of sample 101517AL1, ε_Hf values become much more juvenile with corresponding U-Pb ages in the Paleogene, with an ε_Hf value of 14.0 calculated at 24.9 Ma.

Sample 101517AL3 was collected from sandy, reworked tailings and contains clasts of quartzite, sandstone, biotite schist, metacarbonate, paragneiss, and granite. U-Pb ages (n = 264) range from 22.2 ± 0.2 to 3750.3 ± 7.7 Ma, with 208 of those ages (79%) being >250 Ma (Fig. 4B3). The age spectrum for grains <250 Ma is similar to that of sample 101517AL1, with the exception that this sample contains a greater abundance of ca. 69 Ma grains (30 grains). Lu-Hf isotopic ratios were collected for 14 grains ranging in age from 22.2 to 168.3 Ma (Fig. 5). Epsilon-Hf values range from −3.3 at 168 Ma to −21.9 at 65 Ma, followed by an increase in ε_Hf values to as high as 9.6 at 27 Ma.

Sample 101517AL4 was sampled from a matrix-supported glacial moraine deposit with cobble- to boulder-sized granite and quartzite. The sample was collected ∼5 km downstream from the Pikes Peak Creek placer mine. U-Pb ages (n = 268) range from 24.5 ± 0.4 to 3472.5 ± 10.1 Ma, with 228 of those ages (85%) being >250 Ma. The only prominent peak determined for grains <250 Ma is at ca. 69 Ma, and the sample contains only one grain at ca. 25 Ma (Fig. 4B4). Epsilon-Hf values were collected for eight grains ranging in age from 22.2 to 168.3 Ma (Fig. 5) and range from −4.8 at 167.4 Ma to −23.4 at 66.9 Ma. No Cenozoic zircon grains in this sample were analyzed for Lu-Hf.

The inverse Monte Carlo modeling results for the Pioneer District placer deposits are shown in Figure 6. Figures 6A and 6B display the top 1% of model fits retained relative to the combined, mixed samples. With a KDE bandwidth of 20 m.y., the source inputs (see Methods section) yield a cross-correlation coefficient (R²) of 0.873 ± 0.014 between the model and mixed sample. Studies that consider an R² >0.70 to indicate a reasonable representation of all potential sources (e.g., Sundell and Saylor, 2017; Garber et al., 2020) give us confidence that we are not missing a prominent source in our modeling. Figures 6C and 6D show the relative zircon contributions from each source input calculated by the cross-correlation coefficient and Kuiper test, respectively. The Ravalli Group of the Belt Supergroup displays the smallest relative zircon contribution of ∼4%. The Paleozoic passive-margin and Mesozoic foreland-basin rocks have relative contributions of ∼16% and 25%, respectively. The Royal stock of the northern Flint Creek Range has a calculated contribution slightly smaller than that of the passive margin at ∼11%. Lastly, the Beaverhead Group has the largest relative contribution of ∼43%. The only notable difference between the two contribution plots is that the Kuiper test resulted in a slightly larger contribution from the Mesozoic inputs and less from the Royal stock (Figs. 6C and 6D).

All detrital zircon samples from the Pioneer District placer deposits display strikingly similar age spectra (Figs. 4A and 4B). U-Pb ages from all major North American crustal provinces are present in all four detrital zircon samples, with major peaks in age spectra at ca. 1850, ca. 1750, ca. 1640, ca. 1490, ca. 1100, and ca. 430 Ma. This diverse range of ages suggests abundant zircon recycling, and the ages present in these samples are very similar to those found in the synorogenic Beaverhead Group conglomerate (Fig. 7) (Laskowski et al., 2013; Schwartz and Graham, 2017).

The most prominent age spectra peaks for grains <250 Ma occur at ca. 25, 68–80, 90–120, and 160–180 Ma. The ca. 25 Ma grains present in all samples likely represent ash fall from the Dillon Volcanics of southwestern Montana, an interpretation that is supported by positive ε_Hf(t) values that are characteristic of relatively recent derivation from mantle sources (Fig. 5). For the three samples that were obtained from reworked tailings (101517AL1, 101517AL2, and 101517AL3), the presence of these youngest grains could represent a ca. 25 Ma maximum depositional age (MDA) for the placer deposits. This assumes that no zircon mixing took place during gold extraction from the initial strata (i.e., the DZ ages of the now-disturbed placer deposits are representative of their original ages). It is noteworthy that younger zircon grains from overlying and/or adjacent units could have been mixed into the sampled tailings, which would result in unreliable MDA estimates. However, the localized mining efforts, negligible effects of density separation techniques on DZs during gold extraction (e.g., Silva, 1986), and previous conclusions that placer gold deposition began in the Oligocene (e.g., Loen, 1986) give us confidence that the youngest grains help constrain the age of the deposits. The Late Cretaceous grains are interpreted as first-cycle grains shed from the Royal stock of the northern Anaconda metamorphic core complex. This interpretation is supported by northwest-directed paleoflow measurements and petrographic observations of high quartz and feldspar content in the Oligocene–late Miocene strata of the Flint Creek Basin (Stroup et al., 2008; Portner et al., 2011). The orientation of Pleistocene glacial deposits in the northern Flint Creek Range suggests that similar north-northwest sediment dispersal continued into the Quaternary (e.g., Loen, 1986). The presence of first-cycle grains shed from the Royal stock confirms that material was being eroded from the hypothesized gold-bearing unit during placer development. Assuming that the DZs present in the placer samples are at least partially representative of the units that sourced the gold, this finding supports the interpretation that the Royal stock could have been the original source of the Pioneer District placer gold. Grains dated to 90–120 Ma have many possible explanations, including first-cycle derivation from the nearby Idaho batholith (Gaschnig et al., 2011), shedding of recycled grains from the now-eroded Beaverhead Group (Janecke et al., 2000; Schwartz and Graham, 2017), or contribution from the Cretaceous Blackleaf Formation, which has been mapped ∼10 km north of the placer deposits (Brooks, 2002; Brooks and Sears, 2009) and may have been eroded from the paleo–source area. Grains with U-Pb ages between 150 and 250 Ma are interpreted to have been originally sourced from the Coast Mountains batholith or Sierra Nevada arc segments (Paterson et al., 2011) and are likely recycled most recently from Mesozoic sedimentary units.

Systematic trends that can be observed in Lu-Hf signatures are not well studied for the complex plutons and dikes that intruded the hinterland of western Montana. However, various isotopic systems have been used to investigate the evolution of the Idaho batholith and Challis intrusions of eastern Idaho, which lie ∼100 km west of the AMCC (Foster and Raza, 2002; Gaschnig et al., 2011). One of the isotopic trends observed in the Idaho batholith is a steady decrease in ε_Hf values from –8 at ca. 90 Ma to approximately –25 at ca. 50 Ma, which was interpreted to represent progressive crustal thickening and incorporation of previously existing crust into the arc (Gaschnig et al., 2011). Epsilon-Hf values in the subsequent Challis intrusives and equivalent volcanics (ca. 48 Ma; Gaschnig et al., 2011) display much more juvenile ε_Hf values, ranging from –28 to –3 with an average around –11 (Gaschnig et al., 2011). This “pull up” in ε_Hf space was interpreted by Gaschnig et al. (2011) to represent crustal thinning by extensional collapse of the Cordilleran arc crust and increased mantle input. The DZ Lu-Hf results reported here display a nearly identical trend that extends from ca. 25 Ma to ca. 175 Ma (Fig. 5).

Documented plutonic sources in the Montana–Idaho segment of the arc are almost exclusively <100 Ma (Gaschnig et al., 2010; Gaschnig et al. 2011), which suggests that the older grains with wide-ranging ε_Hf values are recycled. The relatively large spread of ε_Hf values from intermediate to very juvenile for grains at ca. 160 Ma is consistent with derivation from the Sierra Nevada and Coast Range segments of the Cordilleran magmatic arc (Paterson et al., 2011). The highly evolved ε_Hf values for Late Cretaceous–early Paleogene (75–65 Ma) zircons support the hypothesis that many of these grains are first-cycle grains sourced from the Royal stock, given that similar values have been seen for temporally overlapping, arc-associated plutons (e.g., Gaschnig et al., 2011). Oligocene grains in three of the four samples are extremely juvenile, which is consistent with the interpretation that they were sourced as ash fall from the Dillon Volcanics of southwestern Montana (Fritz et al., 2007).

Application of the unmixing model supports the interpretation that the Pioneer District placer gold was sourced from a gold-bearing contact between the Royal stock and overlying supracrustal rocks with complex DZ signatures (Fig. 6). Modeling results make evident the significant number of zircons that were likely recycled into the placer deposits from the Beaverhead Group (Fig. 6). Another possibility is that the Flint Creek Basin shared a similar source with the Beaverhead Group. However, we deem the former interpretation more plausible due to the time gap between Flint Creek Basin and Beaverhead deposition, given that the Flint Creek Basin sedimentation occurred largely during the mid-Oligocene to late Miocene (Loen, 1986; Portner et al., 2011), while main-phase deposition of the Beaverhead Group occurred during latest Cretaceous and Paleocene time (Haley and Perry, 1991; Schwartz and Graham, 2017). Furthermore, the Beaverhead Group is a synorogenic unit that was locally sourced from the frontal fold-thrust belt and various Laramide intraforeland uplifts, resulting in spatially variable provenance trends (e.g., DeCelles et al., 1991; Haley and Perry, 1991; Garber et al., 2020). The differing provenance dependent on geographic location for the Beaverhead Group further discounts the hypothesis that the Flint Creek Basin had a shared source. The absence of strata of the Beaverhead Group in the hypothesized source area itself (Fig. 2) is probably due to its erosion during the Paleocene (e.g., Houston and Dilles, 2013; Schwartz and Schwartz, 2013) and contemporaneous with exhumation of the AMCC in the Eocene and Oligocene (Foster et al., 2010; Howlett et al., 2019; Reynolds et al., 2019). Composite cumulative distribution and KDE plots produced by the placer samples and the Beaverhead Group reveal obvious similarities in their DZ signatures (Fig. 7), corroborating the modeling results above and allowing for a qualitative assessment of how they differ. Most striking is the relative abundance of grains <250 Ma in the placer samples that are not present in the Beaverhead Group (Fig. 7). This discrepancy offsets the cumulative distribution plots between 0 and 250 Ma, but otherwise they follow a very similar trend until 3 Ga, which supports the hypothesis that the Beaverhead Group was an important source for the placer deposits.

Modeling results do not require significant sediment contribution from the Belt Supergroup, a result that may suggest that the relatively small modern extent of the unit in the northern Flint Creek Range is similar to its ancient extent (Mutch and McGill, 1962; Lewis, 1998). Paleozoic and Mesozoic rocks are approximately equally represented in the hypothesized source region, which is reflected in the modeling-result percentages. The slightly larger percent contribution from Mesozoic units could be a consequence of differences in the availability of zircons in the source units, given that the Paleozoic passive margin units, which are dominated by carbonates, contain fewer zircon grains. Additionally, some Paleozoic grains were likely recycled into Mesozoic units, which the model does not have the ability to recognize. The 5%–10% contribution from the Late Cretaceous Pioneer batholith is consistent with first-cycle zircons being shed from the Royal stock.

Taken together, the DZ age spectra and modeling results suggest sediment derivation and recycling from the Mesoproterozoic Belt Supergroup, the Paleozoic passive-margin sequence, Mesozoic foreland-basin sedimentary rocks, Cretaceous plutons, and the synorogenic Beaverhead Group. These interpretations are supported by the mapping results presented in Mutch and McGill (1962) and by the presence of quartzite, sandstone, metacarbonate, and granite clasts within placer deposits (Loen, 1994; this study). Although more robust constraints are needed to make conclusive interpretations on MDA relative to the timing of AMCC exhumation, we propose that placer erosion began during late-stage slip along the Anaconda detachment (early Oligocene; ca. 25 Ma) and that initial deposition occurred in the Cabbage Patch and Squaw Gulch conglomerates of the Flint Creek supradetachment basin (Fig. 8).

The use of DZ geochronology to investigate the provenance of gold placers hinges on two main assumptions, the first being that gold and zircon behave similarly during transport and deposition due to their high densities. This assumption is supported by the observation that placer deposits are commonly enriched in dense minerals like zircon, monazite, and garnet (Reid and Frostick, 1985). However, it has been demonstrated that sediment samples with identical provenance can display significant intersample variability due to the preferential entrainment of larger and lower-density minerals (e.g., Malusà et al., 2016). Therefore, it is important to consider that the high density of gold (∼19.30 g/cm³) could lead to it falling out of entrainment in a sedimentary system before the source-representing zircon grains (density of 4.7 g/cm³) that it was eroded with. This potential limitation is linked to the concept that large, low-density grains fall out of entrainment at the same time as significantly smaller high-density grains (the “principle of hydraulic equivalence”; Rubey, 1933). This source of error may become especially important to consider if the placer deposits of interest were transported great distances, such as in the Shotover/Arrow-Kawarau-Clutha system in Otago, New Zealand, where gold particles have undergone as much as 180 km of fluvial transport from their vein lode source (Youngson and Craw, 1999). These relatively large transport distances provide a sedimentary system more time to hydraulically separate the gold and source-representative zircons that were eroded together, making it more unlikely that the latter could be used to gain reliable insight into the gold’s source (Slingerland and Smith, 1986; Malusà et al., 2016; Roy et al., 2018). It is noteworthy that complex combinations of many hydraulic sorting mechanisms likely control density sorting in sedimentary systems and that the “pure settling behaviors” described by the principle of hydraulic equivalence may not always apply to the evolution of placer deposits (Slingerland and Smith, 1986). Sorting by density of entrained grains is controlled by many factors, including, but not limited to, transport velocity, grain-to-grain interactions, and bed roughness (Slingerland and Smith, 1986); it is therefore plausible that, in some systems, deposition of minerals of substantially different density can overlap spatially and temporally.

Another consequence of long-distance transport that may complicate our approach is morphological changes to gold (e.g., flattening) that affect gold entrainment and ultimately placer evolution (Slingerland and Smith, 1986). In the Shotover system of New Zealand, it was determined that no significant morphological changes occurred to placer gold within the first 10–15 km of transport (Youngson and Craw, 1999), which suggests that more proximal placer deposits would be less likely to be complicated by changes in gold morphology. The Klondike District of the Yukon Territory (Canada) contains gold placer deposits that are closer to their hypothesized source than those in the Shotover system, with transport distances ranging from 0.5 km to 25 km (Knight et al., 1999)—distances that we propose make the DZ provenance approach outlined in this study more suitable.

Considering these settling behaviors, it is reasonable to conclude that the DZs and gold in the Pioneer District placers shared a source. Located on average 2.5–7.5 km from their hypothesized source, these placers have relatively short transport distances, which reduces the likelihood of hydraulic separation of zircon and gold. Short transport distances are supported by coarse and commonly angular placer gold with an average grain size >1.5 mm and some nuggets weighing as much as 840 g (27 oz) (Loen, 1994). Furthermore, our data are consistent with gold derivation from the peaks of the northern AMCC footwall and proximal deposition in the AMCC supradetachment basin (Figs. 2, 8). This topographic landscape provides two conditions that would effectively encourage deposition of placer gold with the zircons that represent the source: (1) a steep channel gradient in resistant, gold-bearing rock that facilitates relatively fast, turbulent stream flows with no significant density sorting (Slingerland and Smith, 1986), and (2) a relatively abrupt transition to a low-gradient basin that provides the hydraulic conditions necessary to concentrate gold and zircon (Roy et al., 2018).

The second common assumption made in almost all DZ studies is that zircon fertility in igneous sources and the relative abundance of zircon in sedimentary sources do not bias the source-region interpretation (e.g., Dickinson, 2008). Varying bedrock mineral fertility near a placer lode source could lead to misinterpretation of age spectra or a complete absence of crucially important age domains. Increasing the number of grains analyzed per sample (n) is one way to reduce the likelihood of missing a potentially important source unit (Pullen et al., 2014). The large-n data sets collected for this study, combined with the well-understood configuration and well-constrained ages of possible plutonic source units, make it unlikely that a significant bedrock source is being missed.

A source of error in the DZ unmixing model may arise due to spatial variability in the DZ age distribution of a given unit. Because none of the potential sources input into the unmixing model were collected directly from the hypothesized source area, we are assuming that samples obtained from other localities have DZ ages that are representative of local units, thus not drastically affecting the modeling results. It has been documented that regional DZ signatures for some of the units in the western U.S. are spatially variable (Schwartz et al., 2019); therefore, sampling potential source inputs from the hypothesized source area itself would improve confidence in the unmixing modeling results. However, even if this were done, the top-down modeling approach must still make a priori assumptions about the potential sources and assume that the sink sample is a pure mixture. Therefore, it may be beneficial to complement future studies of this kind with recently developed “bottom-up” unmixing models, which create synthetic sources from mixed samples and do not require a priori assumptions of source units (DZ non-negative matrix factorization, “NMF”; Saylor et al., 2019).

The wide range of potential source-rock ages in the northern Flint Creek Range (Mesoproterozoic–Oligocene) adds an additional element of difficulty for DZ provenance interpretations. Other regions may have far simpler hypothesized source regions, making the use of DZ spectra and modeling easier.

Overall, our DZ results and DZ unmixing models display that even with abundant zircon recycling, the techniques have potential to provide insight into the source of high-density mineral deposits such as gold placers. The complementary collection of trace-element geochemical data allows a more detailed characterization of source because it can be paired with detrital U-Pb ages and serve as an independent provenance indicator (Gehrels, 2014). Similarly, Lu-Hf isotopic ratios, when paired with U-Pb ages, can provide information into the magmatic and tectonic evolution of potential source areas. An increasing ability for unmixing models to account for zircon recycling and fertility will greatly improve the accuracy of results.

There are many mining districts around the world where the application of these techniques may be of great value. U-Pb analyses of detrital zircons in the Archean Witwatersrand Basin of South Africa have given important insight into the controversial origin of the largest gold placer deposits on Earth (e.g., Gehrels, 2014), but further U-Pb dating and geochemical analysis of detrital zircons have the potential to provide more robust constraints on source units and transport history. Combined DZ U-Pb geochronology and Lu-Hf analyses from the gold-bearing Central Rand Group of the Witwatersrand suggest that zircon was derived from a 3.06 Ga magmatic arc to the north, a 2.94 Ga arc to the west, and >3.28 Ga granitoids, possibly from the east (Koglin et al., 2010). Although the hypothesized source regions are clearly well constrained using these techniques, the application of new modeling techniques, such as DZMix (Sundell and Saylor, 2017), would add another quantitative constraint on provenance and allow for the approximation of each source’s relative contribution. Furthermore, for extremely old placer deposits such as the Witwatersrand whose potential sources may be entirely eroded, it could be illuminating to use the “DZNMF” software (Saylor et al., 2019), which has the potential to reveal previously unknown sources by creating them synthetically. In other placer deposits—such as those along the Clutha River in the Central Otago goldfield of New Zealand—the application of the techniques presented in this research could complement existing lithostratigraphic and geomorphic approaches for understanding placer development (e.g., Youngson and Craw, 1999; Craw et al., 2013). More specifically, lithology-based conclusions that the gold-bearing gravels in this region were derived from recycled sediments during uplift of schist basement (Craw et al., 2013) could be made more robust with the incorporation of DZ U-Pb ages, geochemical analysis, and unmixing modeling.

Future efforts applying these techniques are not limited to provenance studies of ancient and reworked placers, but may also be used to determine potential source areas for modern placer sediments. We propose that a systematic investigation of detrital zircon grains from modern river placers—such as those present in Bonanza Creek of the Klondike Region, Yukon Territory, and the Shotover River in the Otago goldfield of southern New Zealand—could ultimately assist in determining lode sources that are actively shedding gold.

In the Bonanza Creek placers, gold occurs in gravel lags along modern valley floors and is hypothesized to have been shed from basement schists that underwent slow exhumation in the Cenozoic (Lowey, 2006; Knight et al., 1999; MacKenzie et al., 2008). These interpretations could be strengthened or challenged by conducting U-Pb geochronology and isotopic analysis on DZs from the modern sediments. For an even more robust analysis that would elucidate the tectonic history of the source region, low-temperature thermochronology (such as zircon fission-track or zircon [U-Th]/He) could be conducted on detrital grains within the placers (e.g., Cerveny et al., 1988).

The lode source of placer gold in the Shotover River, New Zealand, is relatively well constrained to Miocene quartz veins in the upper catchment of the system (Youngson and Craw, 1999; Craw et al., 2006). Therefore, conducting DZ analyses on modern placers in this system would not only further quantify these constraints, but also serve as an excellent test case for the effects of long-distance transport (>100 km) on gold and zircon settling. Additionally, DZ results could be combined with recent numerical modeling results from the region that investigate the sensitivity of placer formation to fluvial processes within dynamic climatic and tectonic environments (Roy et al., 2018). The pairing of DZ techniques with such numerical models, which simulate the effect of topographic uplift rate and storm intensity on bedrock of varying strength, would enable a comprehensive provenance analysis that determines not only likely source units, but also their vulnerability to erosion and likelihood of the system to concentrate productive gold placers (Roy et al., 2018). Regardless of geographic location and/or methods used, identification of a weathering lode source would enable more localized and potentially less environmentally impactful mining operations.

Field observations, detrital zircon U-Pb geochronology, DZ unmixing modeling, and DZ Lu-Hf analysis suggest that gold from the Pioneer District placer deposits of southwestern Montana was derived from vein and skarn lode sources in the northern footwall of the AMCC. Our data offer the first DZ-based support for previous interpretations that the Late Cretaceous Royal stock precipitated gold along its contact with overlying Proterozoic–Paleocene supracrustal rock (Loen, 1994; McCulloch et al., 2003). Following primary erosion, transportation, and deposition in the Flint Creek supradetachment basin during the late Oligocene–early Miocene, placer gold and the source-representative zircons were largely reworked by Pleistocene glacial and fluvial systems. Zircon signatures allow us to identify 10 likely units that were present in the catchment that sourced the gold. We conclude that one of these source units, the Upper Cretaceous–Paleocene Beaverhead Group conglomerate, was completely eroded from atop the AMCC footwall into the adjacent Flint Creek supradetachment basin. These results serve as confirmation of the utility of using DZ geochronology and emerging DZ unmixing modeling techniques to trace the source of zircon in placer deposits, which in turn can provide insight into the source of detrital gold. The worldwide occurrence of gold placer deposits with unknown source areas provides abundant opportunity to apply these techniques.

We thank Devon Orme, Saré Campbell, and Colin Kennedy for field assistance, and the entire staff of the University of Arizona LaserChron Center (Tucson, Arizona, USA) (support from U.S. National Science Foundation grant EAR-1649254) for analytical assistance. Kurt Sundell is thanked for insightful conversations on the DZMix software. Theresa Schwartz, an anonymous reviewer, Associate Editor G. Lang Farmer, and Science Editor David Fastovsky are appreciated for detailed reviews that improved the quality of this publication. This research was funded by Montana State University and a Geological Society of America Graduate Student Research Grant.