Zircon and the role of magmatic petrogenesis in the formation of felsic-hosted volcanogenic massive sulfide (VMS) deposits: a case study from the mid-Paleozoic Yukon-Tanana terrane, northern Canadian Cordillera

Abstract

Magmatism is a critical component in sustaining hydrothermal convection and metal transport during the formation of volcanogenic massive sulfide (VMS) deposits. Previous studies of magmatic petrogenesis in VMS systems have demonstrated that ore-related volcanic rocks have distinct whole-rock geochemical and isotopic signatures (i.e., high HFSE, REE, Th, εHf-Nd, zircon saturation T) relative to barren volcanic rocks, which supports models of elevated crustal heat flow during periods of ore deposition; however, the petrologic characteristics and intrinsic parameters (e.g., T, fO₂) related to these magmatic events in VMS districts remain poorly understood. Arc–back-arc assemblages from the mid-Paleozoic Yukon-Tanana terrane are well-characterized and include the Finlayson Lake VMS district, which is host to several felsic-hosted deposits (e.g., Kudz Ze Kayah, GP4F, Wolverine) that were generated in a peri-Laurentian continental back-arc tectonic setting. In this study, zircon from back-arc and coeval arc rocks in the Yukon-Tanana terrane was used as a proxy for primary magma formation conditions that generated VMS-proximal and VMS-distal stratigraphy. Our results indicate that zircon grains in VMS-proximal environments have unique textural, geochemical, and isotopic characteristics (e.g., low-aspect ratios, greater abundance of zircon-phosphate intergrowths, Th/U > 1, Zr/Hf > 80, T_zrc > 780 °C, εHf_i > –7) that are clearly distinguished from zircon in VMS-distal rocks in both the back-arc and arc settings (Th/U < 1, Zr/Hf < 80, T_zrc < 780 °C, εHf_i < –7). These signatures correlate to VMS-proximal magmas that were hotter, less fractionated, and contained greater juvenile melt contributions compared to VMS-distal magmas and reflect a series of high-flux magmatic events that directly correspond to the early tectonic development of Yukon-Tanana terrane. Moreover, this study underscores the importance of mineral-scale petrology, geochemistry, and geochronology in defining the primary magmatic conditions that generated VMS-related felsic rocks and highlights the utility of zircon as a prospectivity tool in both grassroots and brownfields VMS exploration.

Introduction

The primary geochemical signatures of volcanic and plutonic rocks that host VMS deposits have been used for exploration and deposit discovery for nearly four decades (e.g., Lesher et al. 1986; Hart et al. 2004; Piercey 2011). Many of these studies have focused on understanding the petrological evolution and tectonic setting of VMS-associated rocks to examine the interrelationship between tectonics, heat flow, and VMS deposit formation. Coincident with this work has been research on subvolcanic intrusions, particularly in Precambrian terranes, where the close spatial, geochemical, and temporal association of subvolcanic intrusive rocks of various compositions with VMS-hosting sequences has led many authors to hypothesize that subvolcanic intrusions were the heat engines that drove hydrothermal circulation during VMS deposit genesis (Campbell et al. 1981; Franklin et al. 1981, 2005; Galley 1993; Hart et al. 2004). However, numerous studies have shown that many VMS districts contain no coeval intrusive suites, whereas in other districts, U–Pb geochronological work has shown these intrusions post-date mineralization and were not related to VMS formation (e.g., Galley 2003; Manor et al. 2022b). Recent work on incremental pluton assembly further illustrates that many large intrusions were unlikely to have formed in a single magma pulse (e.g., Glazner et al. 2004; Paterson et al. 2011; Menand et al. 2015), suggesting VMS models that assume a single magma intrusion driving hydrothermal circulation may require re-evaluation. Furthermore, some workers have shown that in many districts, high-temperature VMS-related felsic magmas associated with mineralization have evidence for juvenile input, likely from underplated basalt at the base of rifts even in areas underlain by continental crust, suggesting that this deeper-seated basalt underplating may be the main driver of VMS hydrothermal activity (Piercey et al. 2008; Piercey 2011; Denisová and Piercey 2022a; Manor et al. 2022b). Despite whole-rock lithogeochemical and isotopic (Hf–Nd) support for the underplating model, little work has been undertaken at the mineral scale to test this hypothesis. Moreover, physiochemical conditions of magma generation (e.g., fO₂-T), and potential metal abundance in the VMS-associated (and barren) melts are poorly understood.

Zircon (ZrSiO₄) is a robust, high-temperature accessory mineral that records the primary crystallization age and intrinsic petrologic parameters of magmas, such as co-crystallizing minerals, crystallization temperature, oxygen fugacity, sulfur speciation, and mantle–crust components (Belousova et al. 2002; Hawkesworth and Kemp 2006; Kemp et al. 2007; Ferry and Watson 2007; Dilles et al. 2015; Smythe and Brenan 2016; Loucks et al. 2020). The composition of zircon has shown to provide detailed insight into the formation of and exploration for mineral deposits (e.g., porphyry Cu; Ballard et al. 2002; Chelle-Michou et al. 2014; Dilles et al. 2015; Smythe and Brenan 2016; Rezeau et al. 2019; Kobylinski et al. 2020; Lee et al. 2020; Viala and Hattori 2021). The use of zircon with respect to VMS deposits, however, has been limited in its use to crystallization temperatures of felsic host rocks (i.e., Ti-in-zircon thermometer; Codeço et al. 2018) and as a tool to assess tectonic evolution of the host rocks (Rosa et al. 2009; Zhu et al. 2017). In this study, we report a comprehensive petrological, geochemical, and isotopic characterization of zircon to address the role of regional magmatism in the formation of felsic-hosted VMS systems. We utilize samples from back-arc and arc-affinity rocks throughout the Late Devonian to Early Mississippian Yukon-Tanana terrane, including the Finlayson Lake VMS district; this district is ideal to study for a number of reasons, including: 1) there is a well-documented stratigraphic and tectonic setting for mineralization with precise temporal controls; 2) there is a range of deposits with high Zn grades, large tonnages, and likelihood to yield abundant accessory minerals based on previous lithogeochemical studies (Piercey et al. 2001, 2003, 2008; Denisová and Piercey 2022a; Manor 2022; Manor et al. 2022b); 3) the Finlayson Lake district is unique in that it preserves both VMS-proximal back-arc stratigraphy and VMS-distal back-arc and arc assemblages; and 4) there are coeval VMS-barren rocks outside the district that are part of similar Yukon-Tanana terrane stratigraphic assemblages (e.g., Colpron et al. 2006a, b). Thus, this study utilizes zircon chemistry as a petrochemical tool to identify characteristics of contrasting basement compositions and petrochemical conditions of melt formation in both ore-hosting and barren settings and provides further insight into metal contents of VMS-related magmas based on inclusions in zircon that were trapped early in the petrogenetic history of these melts.

Geological setting and metallogeny of the Yukon-Tanana terrane

The Yukon-Tanana terrane comprises Devonian to Permian continental margin assemblage rocks of arc and back-arc rocks that formed diachronously along the ancestral western Laurentian margin (Fig. 1; Mortensen and Jilson 1985; Mortensen 1992; Piercey et al. 2004, 2006; Colpron et al. 2006a, b; Manor et al. 2022b) and were accreted onto the Laurentian continental margin in the Mesozoic (Nelson et al. 2006, 2013). The terrane has been divided into four tectonostratigraphic assemblages (Colpron et al. 2006a, b; Nelson et al. 2006; Piercey et al. 2006): (1) Snowcap assemblage (Neoproterozoic to Early Paleozoic), which is considered the basement to Yukon-Tanana terrane (Piercey and Colpron 2009; Kroeger et al. 2023); (2) the arc/back-arc related Finlayson assemblage (Late Devonian to Early Mississippian); (3) Klinkit assemblage (Middle Mississippian to Lower Permian); and 4) Klondike assemblage (Middle to Late Permian). This study focuses on rocks of the Finlayson assemblage, and we refer the reader to Colpron et al. (2006a) for a summary of the other assemblages.

Fig. 1

Terrane assemblage map for the northern Cordillera of British Columbia, Yukon, and Alaska showing the extent of the mid-Paleozoic Yukon-Tanana terrane and Late Devonian and Early Mississippian intrusive suites (e.g., Grass Lakes and Simpson Range plutonic suites, respectively; modified after Colpron and Nelson 2011; Yukon Geological Survey 2020, 2022). Blue boxes indicate sampling regions for this study and associated numbers after the region labels indicate sample locations (refer to Table 1 and Figs. 2 and 3). Red starts are locations of Zn-Pb VMS or SEDEX mineral districts (from Mortensen et al. 2006; Nelson et al. 2013; and references therein). YT Yukon Territories, NWT Northwest Territories, B.C. British Columbia, AK Alaska (USA), Dev. Devonian, Miss. Mississippian

The Finlayson assemblage contains Late Devonian to Early Mississippian rocks that represent the initial arc–back-arc development along western Laurentia (Piercey et al. 2004, 2006; Dusel-Bacon et al. 2006; Nelson et al. 2006; Colpron et al. 2006a, b; Murphy et al. 2006; Manor et al. 2022b). Finlayson assemblage rocks are predominantly mafic and felsic metavolcanic rocks, which are cross-cut by the coeval Late Devonian Grass Lakes plutonic suite (Manor et al. 2022b) and the Early Mississippian Simpson Range plutonic suite (Figs. 1, 2; ca. 357 to 345 Ma; Mortensen 1992; Grant 1997; Murphy et al. 2006; Piercey et al. 2006; Manor et al. 2022b). Granitoids of the Simpson Range plutonic suite intrude both the Snowcap and Finlayson assemblage rocks throughout the Yukon (Figs. 1, 2; Mortensen and Jilson 1985; Stevens et al. 1996; Colpron et al. 2006a, b; Murphy et al. 2006; Yukon Geological Survey 2022). Samples from this study include Finlayson assemblage rocks from four areas, including the Finlayson Lake VMS district, and the Glenlyon, McQuesten, and Stewart River map areas (Fig. 1).

Fig. 2

A Geologic map of the south-central Finlayson Lake district, Yukon-Tanana terrane (modified after Piercey et al. 2003, Yukon Geological Survey 2022, and Manor et al. 2022b). The Money Creek thrust fault crops out south of the extent of the figure. Symbols represent sample locations; numbers correspond to Table 1. Numbers indicate locations of VMS deposits in the region. VMS volcanogenic massive sulfide. B Schematic stratigraphic cross section of the Big Campbell and Fire Lake thrust sheets in the Finlayson Lake VMS district (modified after Piercey et al. 2008 and Manor et al. 2022b). Colored units indicate protolith rock types. Numbers correspond to VMS deposits in the district as in A. Sample symbols are the same as in A. Vertical orange and red lines to the right of the figure indicate U–Pb geochronology constraints on stratigraphy from Manor et al. (2022a, b). Zigzag line represents an angular unconformity between the Grass Lakes and Wolverine Lake groups (Murphy et al. 2002). GLPS Grass Lakes plutonic suite, KZK Kudz Ze Kayah, fm formation

The Finlayson Lake district of southeastern Yukon contains rocks of the Yukon-Tanana and Slide Mountain terranes that have been displaced ~ 430 km from their original position in the Eocene along the Tintina fault (Fig. 1; Tempelman-Kluit 1979; Mortensen and Jilson 1985; Gabrielse et al. 2006; Murphy et al. 2006). Yukon-Tanana terrane rocks comprise variably deformed and metamorphosed Finlayson assemblage arc and back-arc associated volcanic, plutonic, and sedimentary rocks that were deposited upon or intruded Neoproterozoic to Late Devonian basement (Colpron et al. 2006b; Murphy et al. 2006; Piercey et al. 2006; Piercey and Colpron 2009; Manor et al. 2022b). Yukon-Tanana terrane rocks within the Finlayson Lake district are thrust imbricated and occur within the Big Campbell, Fire Lake, Money Creek, and Cleaver Lake thrust sheets (Fig. 2; Murphy et al. 2006; Manor et al. 2022b). Rocks in the lowermost Big Campbell thrust sheet comprise the Upper Devonian North River formation, Grass Lakes group, which includes the Kudz Ze Kayah and Wind Lake formations, and Lower Mississippian Wolverine Lake group, which includes the Little Jimmy, Little Wolverine, Go Creek, and Jasper Creek formations (Manor et al. 2022b). The Late Devonian Grass Lakes plutonic suite and North Lakes intrusion both cut the Grass Lakes group and are unconformably overlain by rocks in the Lower Mississippian Wolverine Lake group. The Fire Lake thrust sheet comprises the Upper Devonian Fire Lake formation, which formed in an oceanic setting outboard of the Big Campbell thrust sheet; rocks of the Fire Lake formation were thrust onto or adjacent to the Grass Lakes group in the Late Devonian (pre-360 Ma; Manor et al. 2022b). The Money Creek thrust sheet sits in the hanging wall of the Money Creek thrust and includes rocks from the North River formation and overlying Upper Devonian to Lower Mississippian Waters Creek and Tuchitua River formations; these rocks are intruded by granitic rocks of the Early Mississippian Simpson Range plutonic suite (Mortensen 1992; Grant 1997; Murphy et al. 2006). The Cleaver Lake thrust sheet contains relatively undeformed and unmetamorphosed rocks of the Late Devonian Cleaver Lake formation, along with mafic and ultramafic rocks related to the Klatsa oceanic slab and Klatsa metamorphic complex; these rock units are intruded by granitic rocks of the Simpson Range plutonic suite (Tempelman-Kluit 1979; Piercey and Murphy 2000; Murphy et al. 2006). The Big Campbell, Money Creek, and Cleaver Lake thrust sheets are each bound in their footwall by Jura-Cretaceous thrust faults that are part of the larger Inconnu thrust system (Manor et al. 2022b).

The Glenlyon region of Yukon-Tanana terrane is situated southwest of the Tintina fault and is dominated by Snowcap assemblage rocks that are overlain by Upper Devonian to Pennsylvanian volcanic rocks of arc affinity and associated clastic rocks of the Finlayson and Klinkit assemblages (Figs. 1, 3; Colpron 1999, 2001; Colpron et al. 2006a; Colpron and Reinecke 2000; Piercey and Colpron 2009). In this area, Yukon-Tanana terrane is bounded to the southwest by arc rocks of the Quesnellia terrane along the Cretaceous to Early Tertiary Big Salmon fault (Colpron et al. 2003; Simard and Devine 2003). The Yukon-Tanana terrane rocks in the Glenlyon region are deformed and metamorphosed to greenschist facies, but primary stratigraphy is preserved (Colpron 1999). Rocks of the Snowcap assemblage basement are overlain by Early Mississippian clastic rocks of the Drury and Pelmac formations. A thin horizon of mafic and felsic volcanic rocks comprises the Lokken member (ca. 350 Ma; Colpron et al. 2006a), which overlies the arkose grit of the Drury formation (Gladwin et al. 2002; Colpron et al. 2006a; Kroeger et al. 2023). These rocks were then intruded by syenitic to gabbroic rocks of the Early Mississippian Ragged pluton (Colpron et al. 2006a). Early to Middle Mississippian volcanic rocks of the Little Kalzas formation and intrusive rocks of the Little Kalzas plutonic suite (i.e., Simpson Range plutonic suite, ca. 349 Ma; Colpron et al. 2006a) define a magmatic arc that is interpreted to be coeval with arc magmatism and back-arc extension in the Finlayson Lake district (349–343 Ma, i.e., Wolverine cycle; Colpron et al. 2006a; Murphy et al. 2006; Piercey et al. 2006).

Fig. 3

Geologic maps of the A Glenlyon region, B McQuesten map sheet (i.e., Reid Lakes), and C Stewart River region (modified after Yukon Geological Survey 2022). Symbols represent sample locations

Yukon-Tanana terrane rocks in the northeastern McQuesten map area are located on the southeastern side of the Tintina fault (Figs. 1, 3; Colpron and Ryan 2010). The area is divided into two fault blocks that are juxtaposed along the extensional Willow Lake fault. To the south-southwest, rocks consist of variably deformed and metamorphosed Snowcap assemblage, Early Mississippian granitoids of the Simpson Range plutonic suite, and younger volcanic rocks. To the north-northeast, rocks are generally well-preserved and relatively undeformed and unmetamorphosed and consist of granitic rocks of the Early Mississippian Reid Lakes batholith and coeval intermediate to felsic volcanic rocks of the Reid Lakes formation (Finlayson assemblage). Granitic rocks of the Reid Lakes batholith are commonly quartz-phyric, biotite granodiorite; however, alkali feldspar megacrystic granodiorite to quartz monzonite is also present in the northeastern part of the intrusion (Colpron and Ryan 2010).

The Stewart River region is located south of Dawson and southwest of the Tintina fault and is situated immediately southwest of the Finlayson Lake district after restoration of the Tintina fault (Figs. 1, 3). Yukon-Tanana terrane rocks in this region are intensely deformed and metamorphosed and are defined by a well-developed regional foliation, isoclinal folds, and localized mylonite zones (Ryan et al. 2003; Gordey and Ryan 2005). The south and central parts of the map area are dominated by high-grade amphibolite facies rocks of the Snowcap assemblage, which are overlain structurally and/or stratigraphically by Late Devonian mafic volcanic rocks (i.e., amphibolites) and minor marble horizons of the Finlayson assemblage. These rocks are subsequently intruded by alkali feldspar granodiorite and tonalite orthogneiss that are correlative to the Early Mississippian Simpson Range plutonic suite (Gordey and Ryan 2005; Piercey et al. 2006).

VMS deposits in Yukon-Tanana terrane are hosted predominately in rocks of the back-arc environment of the Finlayson Lake district (Peter et al. 2007). The Fire Lake formation hosts the 10 Mt mafic-siliciclastic Kona deposit (Sebert et al. 2004). The Grass Lakes group hosts multiple deposits of continental affinity at different stratigraphic levels, including the ~ 1.5 Mt felsic-siliciclastic GP4F and ~ 18.1 Mt bimodal-felsic Kudz Ze Kayah deposits (Piercey et al. 2001; Sebert et al. 2004; Peter et al. 2007; Denisová and Piercey 2022a, b). In the Wolverine Lake group, the Little Wolverine formation hosts all known deposits or mineral showings, including the high-grade, 6.2 Mt felsic-siliciclastic Wolverine deposit and Fisher, Sable, and Puck zones (Bradshaw et al. 2008). Arc rocks in the Finlayson Lake district, Glenlyon, McQueston, and Stewart River areas contain VMS-like occurrences but are not associated with significant tonnages of VMS ores (Grant 1997; Piercey et al. 2001, 2003, 2006; Murphy et al. 2006).

Sampling strategy and terminology

The sampling strategy for this study targeted felsic rocks and sequences that were generated in arc and back-arc terranes located in present-day Yukon; within these tectonic environments, magmatic rocks were sampled to include those that host VMS deposits (n = 16; e.g., Finlayson Lake district) and rocks barren of significant VMS mineralization (n = 8; e.g., arc sequences within Finlayson Lake district and the rest of Yukon-Tanana terrane; Table 1; Manor 2022; Manor et al. 2022a, b; this study). Our approach utilizes new stratigraphic and U–Pb geochronological constraints in the Finlayson Lake district (Manor et al. 2022a, b) and elsewhere in Yukon-Tanana terrane (this study; Table 1) to define magmatic episodes that are coeval and/or distinct from VMS formation. Mineralized, or “VMS-proximal,” stratigraphy consists of volcanic and subvolcanic rocks within ~ 50 m of VMS deposits, which were further refined using U–Pb dates that correspond to discrete mineralized intervals (ca. 363.3, 362.8, and 355.2–355.0 Ma; Manor et al. 2022a). Barren, or “VMS-distal,” stratigraphy has been established both within the back-arc and arc environments; this group contains volcanic and plutonic rocks that have U–Pb dates that lie outside of 2σ uncertainty of previously defined VMS horizons. The VMS-distal stratigraphy in the back-arc is reported below based on the lithological group (i.e., volcanic vs. plutonic) to understand variations in extrusive versus intrusive magmatic processes. Table 1 summarizes the sample details for the study, which includes VMS-proximal volcanic rocks (n = 6) and VMS-distal volcanic and plutonic rocks (n = 10) from the Finlayson Lake district; VMS-distal plutonic rocks from the Simpson Range plutonic suite across Yukon (n = 6); and VMS-distal volcanic rocks from the Finlayson Lake district and Glenlyon region (n = 2). It must be noted that zircon separated for geochronological work and presented in Manor et al. (2022a, b) were also used for trace element and isotopic analysis presented in this study, and the presentation of the data represents a summary of observations between these studies.

Table 1 List of Yukon-Tanana terrane samples utilized in this study

Methodology

Mineral separation, scanning electron microscope (SEM) imaging, and electron probe microanalysis (EPMA) were undertaken in the Hibernia Electron Beam Facility in the Core Research and Instrument Training (CREAIT) Network at Memorial University of Newfoundland. Zircon was annealed at 900 °C for 60 h to eliminate minor radiation damage, enhance cathodoluminescence (CL) emission (Nasdala et al. 2002), promote more reproducible interelement fractionation during laser ablation-inductively coupled plasma mass spectrometry (LA-ICP-MS; Allen and Campbell 2012), and prepare the crystals for subsequent chemical abrasion (Mattinson 2005). Backscatter and CL imaging was performed on a JEOL JSM 7100F scanning electron microscope (SEM) equipped with a Schottky-type field emission gun, Deben Centaurus CL detector, and a Thermo energy-dispersive spectrometer (EDS). Major and minor element oxides in zircon were measured by a JEOL JXA-8230 SuperProbe electron probe microanalyzer, equipped with five wavelength dispersive spectrometers (WDS), a Thermo EDS, and a xCLent IV CL spectrometer, and a W electron gun. Laser ablation (LA)-ICP-MS and chemical abrasion (CA)-ID-TIMS trace element and geochronology analyses were performed at the Isotope Geology Laboratory at Boise State University, Idaho. Laser ablation was performed on each zircon grain using a New Wave Research UP-213 Nd:YAG UV laser (213 nm) using fluence and pulse rates of ~ 5 J/cm² and 10 Hz, respectively, during a 45-s analysis (15-s gas blank, 30-s ablation) for either a 25-µm-diameter spot or 10 × 50 µm raster line. Natural zircon (Plesoviče, Zirconia, Seiland) and synthetic glass standards (NIST SRM-610 and -612) were monitored for instrument drift. The CA-ID-TIMS analyses follow protocols reported in Manor et al. (2022a, b; ESM 1). Zircon grains were plucked from epoxy mounts following LA-ICP-MS analysis and subjected to a modified version of the chemical abrasion method of Mattinson (2005). Hafnium isotopes in zircon were determined in the Micro Analysis Facility in the CREAIT network at Memorial University following LA-ICP-MS and CA-ID-TIMS methods at Boise State University. Isotopic concentrations were measured using a ThermoFinnigan Neptune multi collector (MC)-ICPMS paired to a GeoLas Pro 193 nm ArF excimer laser under operating conditions of 10 Hz frequency, 5 J/cm² fluence, and a 40-µm-diameter spot size. Natural zircon (Plešovice, R-33, and FC-1) and synthetic zircon standards (MUNZirc3, i.e., b142; and MUNZirc4, i.e., b144) were used throughout the run to monitor for instrument drift. The full record of analytical methodology employed during this study is reported in ESM 1 and 2.

Results

Zircon morphology and growth zoning

Back-arc volcanic rocks and subvolcanic intrusions

Zircon in felsic volcaniclastic and porphyritic intrusive rocks are typically elongated prisms with aspect ratios up to 5:1, although examples of more equant (1:1 to 2:1) grains are present (and dominant in Wind Lake formation sample 17MM-001). This distinction is most notable between VMS-proximal and VMS-distal rocks, where proximal rocks contain zircon that are generally more equant and stubby, whereas distal rocks yield zircon that is more elongate with the highest aspect ratios (Fig. 4). One exception is the feldspar porphyry sample 18MM-102, where zircon is elongate with aspect ratios from 2:1 to 4:1. Regardless of rock unit, zircon crystals exhibit a wide range of igneous growth zoning patterns that include oscillatory, sector, and patchy zones, or combinations of each, with distinct CL-dark to CL-bright variations (Fig. 4). Xenocrystic cores are present in most populations (ca. 1060–2500 Ma; see ESM 1) and are CL-bright, commonly fractured, and/or have resorbed margins; however, in some rare cases, “cryptic” inheritance is recognized by older ²⁰⁶Pb/²³⁸U dates in LA-ICP-MS results (ca. 366–397 Ma) where these textures are absent, which may represent antecrystic or xenocrystic populations. Oscillatory zonation is present as micron-scale growth bands that mantle the entire crystal, or as linear or tabular bands that grew parallel to the c-axis (some grains are broken fragments; Fig. 4).

Fig. 4

Representative cathodoluminescence images of zircon, arranged in groups relative to VMS deposits in both arc and back-arc settings. Scale bar in the middle corresponds to the entire figure. Sample numbers are located above grains, whereas the associated laser ablation spot (e.g., LA8) and CA-ID-TIMS grain number (e.g., z2) are located below grain images

Back-arc plutonic rocks

Zircon crystals in rocks from the Grass Lakes plutonic suite and North Lakes intrusion have variable morphologies. Aspect ratios in the zircon grains range from equant (~ 1:1) to elongate (~ 5:1) with most between 2:1 and 4:1; there is no significant difference in crystal shape between phaneritic and porphyritic lithologies. The grains have abundant inclusions and contain irregular growth zoning from oscillatory to patchy with rare sector zones (Fig. 4). Xenocrystic cores (372–428 Ma and 1090–2780 Ma) have bright CL signatures and/or significantly resorbed margins (g16, 17MM-005; g31, 17MM-061; ESM 2).

Arc plutonic and volcanic rocks

Zircon in arc rocks have similar grain morphologies, regardless of location. Crystals are predominantly equant and stubby with lower aspect ratios (< 3:1) compared to back-arc-derived zircon (Fig. 4). Samples 18MM-108 and 19MM-06 are the only volcanic rocks in the sample suite, and both contain zircon grains with slightly more elongate crystal shapes (up to 4:1) than the plutonic rocks. Well-developed oscillatory growth zoning is present in all grains, and they also commonly contain sector zones either throughout the growth of the entire crystal or exclusively in the cores of crystals (Fig. 4); examples of the latter are only observed in zircon from rocks in the Glenlyon region (i.e., 19MM-06 and 19MM-07). No obvious xenocrystic textures are present in any of the picked grains, which is confirmed by LA-ICP-MS results where older dates may be antecrystic populations (358–375 Ma; ESM 1).

Inclusions and petrologic observations

Metal-rich inclusions in zircon

Zircon grains typically contain abundant silicate (quartz, plagioclase, orthoclase, white mica, biotite/phlogopite, thorite/uranothorite), phosphate (apatite, monazite, xenotime), Fe-Ti-oxide and spinel (magnetite, hercynite, ilmenite, rutile), and rare metal and volatile mineral-rich inclusions. Over 20 metal and volatile mineral-rich inclusions were discovered and characterized by SEM–EDS from zircon mineral separates and thin sections (n = 2). These inclusions were observed in zircon exclusively in back-arc rocks, including four VMS-proximal volcanic rocks (17MM-074, 18MM-102, 18MM-133, and P00-WV-1C), one VMS-distal volcanic rock (18MM-114), and two VMS-distal plutonic rocks (17MM-061 and 18MM-105). Metal- and volatile-bearing minerals are present as sulfides, silicates, and oxides, as rare metal alloys, or contained within melt inclusions/glass (Fig. 5; ESM 2). The inclusions typically have a diameter < 10–25 microns, are fully enclosed by zircon, and are not associated with any fractures (i.e., were captured as inclusion in zircon during crystallization from the magma). Mineral inclusions reveal the presence of base metals (Pb, Zn, Cu, Fe, ± Ni) and volatile (S, Cl, F) elements present as either native metal, alloys, or mineral phases. Elemental abundances are variable within individual inclusions but can be divided into six classes of metal or volatile element-rich inclusions: (1) Fe sulfide (± Zn); (2) Zn sulfide (± Fe); (3) Pb-Fe Cu sulfide; (4) Zn-Fe alloy; (5) Pb-Fe (± Re) alloy; (6) native Cu or Pb; and (7) Cl- or F-bearing minerals (Fig. 5; ESM 2). These inclusions are typically observed as individual occurrences in zircon; however, one rare case (sample 18MM-133) highlights the presence of up to six chemically and texturally distinct inclusions (Fig. 5J). Here, native Pb metal, Pb-Fe ± Re alloy, Zn-bearing pyrite, Fe-bearing sphalerite, and xenotime are all observed in one rounded inclusion.

Fig. 5

Representative backscatter electron (BSE) images of observed metal and volatile-rich inclusions in zircon. A Zn-Fe alloy–white mica, and spinel inclusions; B Pb-Cu sulfide inclusion; C melt inclusion containing native Cu with minor alkali feldspar crystals; D galena and Cu-Fe-S inclusions in melt bubble(?); E bleb of pyrrhotite with a melt inclusion(?) at contact with zircon, and a smaller melt inclusion with native Cu and pyroxene; F pyrrhotite bleb; G chlorine-rich inclusion with Al-glass? (chloroaluminite?); H pyrrhotite and K + Na-chloride inclusions; I sphalerite inclusion (18MM-133); J BSE image of multi-phase inclusion in ~ 25-µm bubble, which includes native Pb metal, Pb-Fe (+ Re?) alloy, pyrite (+ Zn), Fe-bearing sphalerite, and xenotime (P00-WV-1C); K inclusion with intergrowths of Fe-oxide (magnetite-wüstite?) and hercynite; L ilmenite, monazite-pyroxene/mica intergrowth, and glass inclusions. Mineral abbreviations: Afs alkali feldspar, Ap apatite, Gal galena, Hc, hercynite, Ilm ilmenite, Mnz monazite, Po pyrrhotite, Px pyroxene, Py pyrite, Qz quartz, Rt rutile, Sp sphalerite, Spl spinel, Xtm xenotime, Zrn zircon

Accessory mineral associations in back-arc and arc rocks

Zircon in back-arc rocks is typically associated with numerous other primary accessory minerals, particularly HFSE-REE phosphates that appear to control the magmatic REE-HFSE budgets of Yukon-Tanana terrane magmas. Apatite, monazite, xenotime, and rutile are the most common minerals associated with zircon; allanite and titanite are minor; and uranothorite and thorite are present but rare (Fig. 6). Many accessory minerals occur together in clots or intergrowths. Clots typically include rutile as a nucleating core and at least three minerals intergrown with each other, whereas simpler intergrowths can occur between any of the listed minerals. In VMS-proximal rocks, zircon and phosphate minerals are intergrown as inclusions or clusters and imply co-crystallization, where zircon, monazite, and xenotime are more abundant compared to VMS-distal rocks (Manor et al. 2022b). In contrast, zircon and phosphate minerals in VMS-distal rocks commonly occur as discrete crystals, have lower modal abundances, and are rarely present as zircon-phosphate intergrowths (e.g., Fig. 6A, D). Where intergrowths are present, in both VMS-proximal and distal rocks, zircon is commonly intergrown with apatite, monazite, xenotime, and rutile, and can have late xenotime present as diffuse or sharp rims on zircon; the latter texture is interpreted to be secondary in nature, either late magmatic or hydrothermal. Apatite is typically euhedral to subhedral, contains well-developed growth zoning, and can contain rare thorite and allanite inclusions. Monazite is present as inclusions in both zircon and apatite and as discrete subhedral grains that can have patchy and irregular grain margins and can be overprinted or associated with allanite. Secondary monazite fills fractures and cracks in gangue minerals and commonly occurs as isolated grains away from other igneous accessory minerals. Xenotime is euhedral to anhedral and commonly associated with zircon, occurring as inclusions or ragged secondary rims; it also occurs as inclusions in apatite and monazite, and itself contains rare inclusions of thorite. Rutile is present in all rocks sampled and is typically euhedral to subhedral as free grains or cluster-forming grains with zircon and phosphates; rarely, it occurs as blebby or wormy inclusions in zircon. Complex rutile-quartz intergrowths rimmed by pyrite rims are common in zircon from the Wolverine Lake group (e.g., Fig. 6H). Allanite is most common in VMS-distal plutonic rocks as primary minerals with monazite rims, whereas back-arc volcanic rocks typically contain secondary anhedral patches and/or rims around monazite or pyrite.

Fig. 6

Backscatter electron images highlighting representative accessory mineral assemblages for the Kudz Ze Kayah formation (A–C), Grass Lakes plutonic suite (D), and Wolverine Lake group (E–I). Stars (*) after the sample number denote VMS-proximal sample. A Intergrowth of primary apatite, zircon, and monazite. B Discrete crystals of monazite and zircon with late magmatic or hydrothermal xenotime rims. C Cluster of rutile, zircon, and monazite; the zircon crystal in the lower left contains a sphalerite inclusion (see Fig. 5I). D Cluster of intergrown monazite and zircon, at a grain boundary with apatite. E Complex clusters of rutile, zircon, and monazite. F Multi-mineral intergrown between zircon, xenotime, monazite, apatite, and pyrite. G Rutile inclusion in zircon, at the boundary of pyrite. H Rutile grains that have quartz intergrowths or dissolution textures (?), intergrown zircon and monazite, and late (likely hydrothermal) cubic pyrite. I Discrete apatite, zircon, and monazite crystals. Mineral abbreviations are the same as Fig. 5

Arc rocks contain primarily zircon and apatite, minor titanite, allanite, and rutile, and rare monazite, uranothorite and thorite; no xenotime was observed. Zircon and apatite are typically present as discrete inclusions in primary igneous quartz, clinopyroxene, and amphibole, or within quartz-muscovite layers in volcaniclastic rocks. Both zircon and apatite can exhibit resorbed grain margins that contain rare very fine-grained monazite, uranothorite, and thorite; the resorbed margins are interpreted to be late magmatic or hydrothermal in nature. Apatite contains inclusions of monazite, allanite, and zircon. Allanite is common, especially in sample 18MM-107, where allanite up to 0.5 mm was observed with Th-REE-rich (thorite?) exsolution minerals and a radial fractures filled with secondary Th-poor allanite.

Electron probe microanalysis (EPMA)

Major element oxide results for zircon in the entire sample suite (n = 2358 spots) are presented in ESM 2. The predominant cations Zr⁴⁺ and Si⁴⁺ in the zircon tetrahedra are measured as ZrO₂, which ranges from 64.0 to 68.6 wt% (median = 66.1 wt%) and corresponds to 0.965 to 1.013 Zr atoms per formula unit (apfu), and SiO₂, which ranges from 31.4 to 33.9 wt% (median = 32.6 wt%) and corresponds to 0.980 to 1.024 afpu. The HfO₂ compositions range between 0.43 to 1.38 wt%, which corresponds to Hf⁴⁺ from 0.004 to 0.012 apfu. A wide range of P₂O₅ and Y₂O₃ contents are observed, ranging from 0 to 0.82 wt% (0–0.021 apfu P⁵⁺) and 0 to 1.21 wt% (0–0.020 apfu Y³⁺), respectively.

Chemical abrasion ID-TIMS and LA-ICP-MS U–Pb zircon geochronology

All U–Pb age constraints utilized in this study are shown in Table 1. Here, we present geochronology results for four arc rocks in Yukon-Tanana terrane south-southwest of the Tintina fault, including the Glenlyon, McQuesten, and Stewart River map areas (Fig. 1; ESM 1 and 2). Two samples were dated in the Glenlyon region: (1) 19MM-06, a felsic crystal-lapilli tuff from the Lokken member (i.e., Finlayson assemblage); and (2) 19MM-07, a strongly foliated biotite granodiorite of the Little Kalzas suite (i.e., Simpson Range plutonic suite). In 19MM-06, seven concordant CA-ID-TIMS zircon fractions yielded a weighted mean ²⁰⁶Pb/²³⁸U date of 354.255 ± 0.039 Ma (2σ; MSWD = 0.66), which is older than TIMS date of 350.4 ± 1.3 Ma reported in Colpron et al. (2006a), contains no inherited fractions, and is thus interpreted as the true eruption age of this volcaniclastic rock. In sample 19MM-07, LA-ICP-MS zircon data for this sample give a weighted mean ²⁰⁶Pb/²³⁸U date of 354.2 ± 1.8 Ma (2σ; MSWD = 1.3; n = 23); this date is older than a previous TIMS date of 349.0 ± 0.7 Ma by Colpron et al. (2006a) and interpreted as the crystallization age of the rock with limited inherited components. In the McQuesten map area, zircon from a hornblende-biotite granodiorite (19MM-01) from the Reid Lakes batholith (i.e., Simpson Range plutonic suite) was dated by LA-ICP-MS and yielded a weighted mean ²⁰⁶Pb/²³⁸U date of 347.6 ± 1.3 Ma (2σ; MSWD = 1.3; n = 33). This date is younger than the previously reported date of 355.9 ± 0.8 Ma from the same outcrop (Mortensen 2009) but overlaps with the single analysis outside the weighted mean (346.5 Ma), suggesting multiple populations of zircon may be present in this rock. Thus, our younger date is interpreted to be the most appropriate estimate for the crystallization age of this granodiorite. One granodiorite (02RAYP30D1) in the Stewart River region was dated by CA-ID-TIMS methods, where seven concordant zircon fractions gave a weighted mean ²⁰⁶Pb/²³⁸U date of 354.204 ± 0.056 Ma (2σ; MSWD = 1.13); this age is interpreted as the crystallization age of the granodiorite.

Trace element concentrations and Hf isotopes in zircon

Trace elements and Hf isotopes in zircon were assessed based on concomitant in situ U–Pb dating results, where each grain/analysis was subsequently divided into groups of autocryst, Pb-loss, or xenocryst. The presentation of results only includes grains in the autocryst and Pb-loss categories to determine the relationship to primary magmatic conditions. The Pb-loss category is interpreted to be magmatic in origin where only the U-Th-Pb isotopic systematics were altered but not the trace element abundances. Four primary groups are recognized relative to VMS-proximity and tectonomagmatic setting: (1) VMS-proximal, back-arc volcanic rocks; (2) VMS-distal, back-arc volcanic rocks; (3) VMS-distal, back-arc plutonic rocks; and (4) VMS-distal, arc plutonic and volcanic rocks. Zircon in all groups yield a wide range of trace element concentrations, but the defined populations above have distinct compositions (e.g., Th, U, Ti, total REE; Table 2; ESM 1). Only trace element results that fall above the lower limit of detection (LOD = 3 × detection limit; e.g., La) are plotted on chondrite-normalized REE plots. All trace element and Hf isotope results are presented in Figs. 7, 8, 9, 10, 11, and 12; Table 2; and ESM 1 and 2.

Table 2 Key magmatic parameters determined from zircon chemistry

Fig. 7

Binary fractionation plots for LA-ICP-MS results in zircon (unless otherwise noted). A Th/U vs. Zr/Hf (EPMA). B Th/U vs. Hf/Y. C Ti vs. Yb. D Ti vs. Th/U; horizontal grey dashed lines indicate approximate Ti-in-zircon crystallization temperature at a respective Ti concentration. Vectors in A–C show relative fractionation trends for respective accessory minerals. The overlapping polygons represent 90% contour for VMS-proximal (solid red) and VMS-distal (dashed black) groups

Fig. 8

Fractionation indices vs. U–Pb date (Ma). A Th/U, B Zr/Hf, C Yb/Dy, D U/Yb, E Hf, in ppm. Each sample is plotted as a single box plot, where points indicate outliers. Box colors and associated outlier points represent VMS-proximal (red) and VMS-distal (black) samples. Black arrows on the right indicate trends in data toward increased fractionation in the melt. Orange arrows indicate progression of samples up-stratigraphy in the Kudz Ze Kayah formation (blue) and Wolverine Lake group (green). FC fractional crystallization

Fig. 9

Key trace element and Hf isotope results for Yukon-Tanana terrane zircon. A εHf_i vs Th/U. B εHf_i vs. T_zrc (°C). C Th/U vs T_zrc (°C). Results show an increase in Th/U and εHf_i with crystallization temperature (T_zrc) due to increasing juvenile/mafic input (see discussion for details). The overlapping polygons represent 90% contour for VMS-proximal (solid red) and VMS-distal (dashed black) groups

Fig. 10

Eu_N/Eu_N* vs Hf (ppm) for A back-arc zircon in the Finlayson Lake district and B arc zircon from across Yukon-Tanana terrane. Trends for mineralization arc magmas the boundary between ore-forming and barren magmas (Eu_N/Eu_N* = 0.4) from Dilles et al. (2015) are superimposed on our results for reference. All back-arc zircon plot well below the ore-forming–barren divide and is attributed to highly reduced and plagioclase-dominant parental melts. Symbols for A are based on location to VMS deposit as in Figs. 7 and 9, whereas symbols in panel B are separated by region in Yukon-Tanana terrane and observed in the legend

Fig. 11

Age-corrected Hf isotopic compositions of zircon from Yukon-Tanana terrane rocks. A εHf_i vs. LA-ICP-MS U–Pb date, where grey error bars incorporate the maximum extent of 2σ error on a single zircon grain. B εHf_i vs. weighted mean CA-ID-TIMS U–Pb date. Both zircon (small symbols) and whole-rock (large symbols) from Manor et al. (2022b) are shown for comparison. Box and whisker plots for zircon populations show the 1Q (25th percentile), median, and 3Q (75th percentile) for each sample; red boxes denote VMS-proximal samples, whereas black boxes are for VMS-distal samples. Depleted mantle (bold black line) and crustal evolution (grey dashed lines) curves in A and B were calculated using ¹⁷⁶Hf/¹⁷⁷Hf_DM = 0.283238 and ¹⁷⁶Lu/¹⁷⁷Hf_DM = 0.03976 (Vervoort and Blichert-Toft 1999), and ¹⁷⁶Hf/¹⁷⁷Hf_CHUR = 0.283238 and ¹⁷⁶Lu/¹⁷⁷Hf_CHUR = 0.03976 (Bouvier et al. 2008), respectively, using an upper crustal isotopic reservoir with a ¹⁷⁶Lu/¹⁷⁷Hf = 0.0072 (McLennan 2001). Symbol colors correspond to those in Fig. 7

Fig. 12

Modelled oxygen fugacity results for zircon from Yukon-Tanana terrane parental melts, using the methods of Smythe and Brenan (2016) and Loucks et al. (2020). A log fO₂ vs T_zrc (°C) for Smythe and Brenan (2016); B log fO₂ vs T_zrc (°C) for Loucks et al. (2020). Small transparent symbols are individual zircon analyses that have an average 2σ error of ~ 2 log units. Large symbols are the mean of each sample population with associated 2σ errors. The nickel-nickel-oxide (NNO), fayalite-magnetite-quartz (FMQ), and iron-wüstite (IW) buffers are shown by black dashed curves. C ΔFMQ vs T_zrc (°C) for Smythe and Brenan (2016); D ΔFMQ vs. T_zrc (°C) for Loucks et al. (2020). ΔFMQ represents the variation of the fO₂ from the FMQ buffer. Fields highlight the oxygen fugacity of sampling groups relative to VMS deposits, where the distal field here only accounts for felsic volcanic rocks. KZK Kudz Ze Kayah formation, WLG Wolverine Lake group

VMS-proximal volcanic rocks (back-arc)

VMS-proximal rocks occur at two primary time intervals (363.3 to 362.8 Ma and 355.3 to 355.0 Ma) and contain zircon with the highest Th (37–3807 ppm) and lowest U contents (36–1088 ppm), and highest median Th/U of all zircon in the sample suite (1.3; range = 0.13–5.3). This group also contains high Ti (4.6–43 ppm; median = 16.4 ppm) and total REE concentrations (562–6511 ppm; median = 1553 ppm). Chondrite-normalized LREE ratios are generally steep with (Sm/La)_N = 81–2200 (median = 178), and HREE-MREE abundances illustrated by straight to slightly concave-down slopes where (Yb/Gd)_N = 6–69 (median = 11). These zircon analyses have prominent negative Eu and positive Ce anomalies and Eu_N/Eu_N* and Ce_N/Ce_N* between 0.04–0.42 (median = 0.08) and 1.6–135 (median = 7.6), respectively. Age-corrected Hf isotopic compositions (εHf_i) range from − 18 to + 2.5, with the least negative εHf_i median value in all groups (εHf_i =  − 6.2) and T_DM = 0.9–1.6 Ga (median = 1.2 Ga).

VMS-distal volcanic rocks (back-arc)

VMS-distal volcanic rocks in the back-arc environment contain zircon with distinct trace element and isotopic signatures compared to VMS-proximal volcanic rocks. Thorium contents are typically similar to slightly lower (23–1481 ppm), whereas U contents overlap (58–1091 ppm) with the VMS-proximal zircon. They have lower Th/U ratios of 0.1 to 2.5 (median = 0.58) and slightly lower Ti (1.6–33 ppm; median = 13) and overlapping total REE contents (480–4387 ppm; median = 1410 ppm) with VMS-proximal zircon. Zircon in this group typically exhibit REE patterns similar to VMS-proximal zircon with steep chondrite-normalized (Sm/La)_N (43–2560; median = 147) and (Yb/Gd)_N patterns (6–41; median = 13) and similar Eu and Ce anomalies, negative anomalies where Eu_N/Eu_N* = 0.01–0.32 (median = 0.09) and Ce anomalies that are slightly less-positive with Ce_N/Ce_N* = 1.7–49 (median = 6.5). The εHf_i of distal zircon overlap with proximal zircon (εHf_i =  − 18 to − 1.1) but do not exceed chrondritic values and have a significantly more negative median εHf_i of − 12 compared to VMS-proximal zircon grains. Model ages are slightly older than VMS-proximal group with T_DM = 1.0–1.9 Ga (median = 1.3 Ga).

VMS-distal plutonic rocks (back-arc)

Plutonic rocks in the back-arc are all part of the VMS-distal group based on crystallization age and stratigraphic setting. In this group, zircon have Th and U contents that typically overlap with the VMS-proximal group (Th = 21–2723 ppm; U = 29–1468 ppm) and give similar Th/U with some of the highest observed ratios of all zircons in this study (0.07–6.0; median = 1.0). The Ti (3.1–58 ppm; median = 12 ppm) and total REE (536–3578 ppm; median = 1299 ppm) concentrations are lower than all the back-arc volcanic rocks. Chondrite-normalized LREE patterns are steep with a wide range of (Sm/La)_N values (9.5–1360; median = 196) and have overlapping MREE-HREE ratios to volcanic rocks with (Yb/Gd)_N ranging from 5.7 to 54 (median = 11). Negative Eu and positive Ce anomalies are present in similar magnitudes to the volcanic rocks (Eu_N/Eu_N* = 0–0.32, median = 0.07; Ce_N/Ce_N* = 2.4–102; median = 13). Age-corrected Hf isotopes encompass the entire range of values reported for volcanic rocks in both proximal and distal categories (εHf_i =  − 18 to − 3.0, median =  − 7.5; T_DM = 1.1–1.6 Ga, median = 1.2 Ga). Within this range, however, two groupings are observed: (1) GLPS phaneritic and porphyritic granitoids (17MM-061, 17MM-062, 18MM-105) and the North Lakes diorite (96DM-065) crystallized between 362.6 and 361.7 Ma and have εHf_i values between − 18 and − 3.3 with medians between − 8.6 and − 6.1; and (2) GLPS granitoid (17MM-005) crystallized ~ 800 kyr after the other granitoids at ca. 360.9 Ma, yielding a much more negative range of εHf_i (− 17 to − 9.6; median =  − 13).

VMS-distal plutonic and volcanic rocks (arc)

Arc rocks contain zircon that generally have much more restricted trace element and Hf isotopic compositions compared to back-arc-derived zircon. Thorium concentrations (42–2893 ppm; median = 208 ppm) overlap with those in the back-arc, whereas median U contents are higher despite a similar range (109–2132 ppm; median = 378 ppm). As a result, zircon Th/U show a tight range between 0.26 and 1.4 (median = 0.6). Titanium (0.18–47 ppm; median = 4.7 ppm) and total REE (289–4160 ppm; median = 871 ppm) abundances are relatively lower than back-arc zircon. Chondrite-normalized ratios show relatively flatter LREE patterns and steeper MREE-HREE profiles with (Sm/La)_N = 1–886 (median = 63) and (Yb/Gd)_N = 14–139 (median = 34), respectively (ESM 2). Similarly, the negative Eu anomaly is typically less pronounced (Eu_N/Eu_N* = 0.07–1.18, median = 0.24), whereas the positive Ce anomaly is greater than back-arc zircon patterns (Ce_N/Ce_N* = 2.1–159, median = 22).

Hafnium isotopic compositions have more restricted ranges in εHf_i compared to back-arc zircon results; however, samples still span multiple epsilon units. Age-corrected εHf_i ranges from − 20 to + 2 with a median of − 13 and has Hf-depleted mantle model ages between 0.8 and 1.9 Ga (Fig. 11); however, there are multiple, geographically distinct populations within this range. In the Glenlyon region, sample 19MM-06 gives εHf_i values of − 4.6 to + 0.7, whereas 19MM-07 has slightly more evolved εHf_i (− 7.1 to − 4.3). The Stewart River Simpson Range plutonic suite granodiorite (02RAYP30D1) yielded more chondritic εHf_i values than any other rocks in this study, where εHf_i ranges from − 1.5 to + 2. In the Reid Lakes batholith, sample 19MM-01 has a restricted range of evolved εHf_i (− 15 to − 12). Plutonic samples in the Finlayson Lake district (18MM-107, P99-24, P99-82) have the most evolved Hf isotopic signatures (εHf_i =  − 20 to − 7), whereas the quartz porphyritic rhyolite (18MM-108) contains two distinct populations of εHf_i at a crystallization age of ca. 363.4 Ma (εHf_i =  − 19 to − 16, and − 12 to − 9.8).

Discussion

Petrogenesis and chemical evolution of zircon in Yukon-Tanana terrane magmatic systems

The zircon crystal structure can incorporate major elements of Zr, Si, and O, minor Hf and Y (< 1 wt%), and trace REE, P, U, Th, and Ti (10s to 1000s ppm; Hoskin and Schaltegger 2003). Substitutions of these elements occur within the Zr⁴⁺ octahedral or the Si⁴⁺ tetrahedral site in the zircon lattice (e.g., Hanchar et al. 2001) and the substitutions are affected by the fractionation from co-precipitating phases (Speer 1982; Finch et al. 2001), crystallization parameters of the host melt (e.g., temperature, redox state; Rubatto and Hermann 2007; Claiborne et al. 2018), and the influence of crust-mantle mixing during melt generation (i.e., Hf isotopes; e.g., Belousova et al. 2006; Kemp et al. 2007; Storck et al. 2020). Co-existing accessory minerals with elevated REE concentrations (e.g., apatite, monazite, titanite, allanite) are common in felsic rocks and can directly compete for REE and HFSE with zircon during crystallization (Bea 1996; Schaltegger and Davies 2017). Manor et al. (2022b) showed that modal abundances of zircon, monazite, and xenotime directly correlate with whole-rock compositions in Finlayson Lake district back-arc rocks and thus control the trace element budgets of rocks in VMS environments. In this study, we provide textural context for these MLA-SEM results and show that zircon contains inclusions of, or are intergrown with, apatite, monazite, and xenotime, suggesting that co-precipitation and fractionation of these minerals occurred in the early stages of crystallization in Yukon-Tanana terrane back-arc magmas (Figs. 5, 6); therefore, we test this hypothesis below using fractionation indices to determine the significance of co-fractionating minerals on zircon chemistry.

Co-precipitating accessory minerals and impact on zircon chemistry

Incompatible element concentrations (e.g., Hf, Ti, Th, U) and ratios (e.g., Th/U, Zr/Hf, Yb/Dy, U/Yb) in zircon can monitor changes in melt chemistry due to element fractionation by co-precipitating minerals (Claiborne et al. 2006, 2010; Samperton et al. 2015; Schaltegger and Davies 2017; Large et al. 2020; Lee et al. 2020). There are distinct, time-resolved shifts (of median values) from high Th and Ti, low U and Hf zircon in VMS-proximal rocks to low Th and Ti, high U and Hf zircon in VMS-distal rocks, which reflect the progressive differentiation of back-arc melts with time (Figs. 7, 8; Table 2). Similarly, VMS-proximal zircon records the highest Th/U and Zr/Hf of all zircon in the study, which is interpreted to reflect early zircon crystallization from relatively primitive crustal melts that have not undergone significant differentiation (e.g., Claiborne et al. 2006). As melts evolve and crystallize zircon, the residual melt should become enriched in Th/U and depleted in Zr/Hf, a trend that is observed only in VMS-proximal and some back-arc plutonic rocks. The dominant trend for back-arc zircon illustrates a decrease in both Th/U and Zr/Hf, which is interpreted as zircon crystallization coeval with Th-bearing minerals (in this case likely monazite ± apatite), and which impacted the Th budget available to co-crystallizing zircon (Fig. 7A; Bolhar et al. 2008; Grimes et al. 2015; Schaltegger and Davies 2017). There is also a good correlation between decreasing Ti and Th/U in back-arc zircon, which is interpreted to suggest that the back-arc melts cooled with progressive melt evolution (Watson and Harrison 1983; Ferry and Watson 2007). In contrast, in arc zircon, there is a relatively constant Th/U (including the highest U contents relative to back-arc zircon) with decreasing Ti, which suggests that temperature played a minimal effect on zircon chemistry in Yukon-Tanana terrane arc magmas (Fig. 9C).

Ratios for REE (e.g., Yb/Dy) can also be used to determine if MREE to HREE-bearing minerals (e.g., apatite, xenotime, amphibole, titanite) influenced zircon compositions (e.g., Chelle-Michou et al. 2014; Schaltegger and Davies 2017). Back-arc zircon grains have consistently low Yb/Dy (< 5) and Hf/Y (< 10) and are negatively correlated with Ti and Yb, which indicate zircon crystallization in the absence of amphibole (Lee et al. 2020); the results also do not clearly indicate whether or not apatite or xenotime imparted HREE or MREE signatures on zircon chemistry. Conversely, high Hf/Y (> 10) and Yb/Dy ratios (median =  ~ 5–10) and positively correlated Ti vs. Yb in arc zircon indicate that amphibole was an early fractionating mineral that controlled the HREE/MREE composition of both melt and zircon compositions in Yukon-Tanana arc settings (Figs. 7, 8). Therefore, fractionation indices in zircon highlight three important points about Yukon-Tanana terrane magmas: (1) in the back-arc, VMS-proximal zircon crystallized from melts with low degrees of fractionation that evolved with early fractionation of zircon ± apatite and slightly later monazite ± xenotime assemblages, whereas VMS-distal zircon were generated from more fractionated melts that are dominated by co-fractionating zircon ± apatite ± monazite ± xenotime ± allanite ± titanite assemblages; (2) there is a strong temporal control on zircon chemistry in the back-arc environment where variable fluxing of less fractionated and more fractionated melts is interpreted to represent changes in the tectonomagmatic environment, perhaps due to a basaltic underplate beneath the rift system (e.g., Piercey 2011; Manor et al. 2022b); this is particularly evident in the Little Wolverine formation, where Th/U and Zr/Hf both increase up-stratigraphy on a < 1 Myr timescale; and (3) back-arc magmas crystallized zircon in the absence of amphibole, whereas arc zircon crystallized coeval to, or following, early amphibole fractionation, which has important implications for water content of the melts in each tectonic setting (e.g., Kobylinski et al. 2020; Lee et al. 2020).

Back-arc zircon have the largest variation in Th/U ratios and show evidence for co-precipitation of zircon, monazite, xenotime, and apatite (Figs. 6, 7), which are present as inclusions in and/or intergrowths with back-arc zircon. However, the consistent REE patterns and ratios, and elevated Th/U in zircon suggest these minerals, despite their apparent coeval textures, have had a limited control on zircon chemistry in VMS-proximal rocks. This observation is somewhat contradictory, because co-fractionating monazite should preferentially incorporate Th during crystallization based on its higher partition coefficients for Th compared to zircon, resulting in melt depletion in Th relative to U, which would result in lower Th/U in zircon grains that crystallized after monazite (Stepanov et al. 2012; Duc-Tin and Keppler 2015), features not observed in back-arc zircon chemistry. Above, we show that high Th/U and low Zr/Hf ratios can be achieved with zircon-only fractionation (Fig. 7A); however, this could also be explained by local variations in the mafic versus felsic components in the melt (e.g., Wang et al. 2011; Kirkland et al. 2015). Previous research has demonstrated that (1) the abundance of Th is at least twice that of U in mafic compared to felsic melts (Kirkland et al. 2015) and (2) zircon from mafic to intermediate rocks, on average, has higher Th/U than zircon from felsic granitic rocks (0.93 vs. 0.59, respectively; Wang et al. 2011). Zircon has also been hypothesized to crystallize from “open” system processes that allow for an unlimited supply of trace elements to the melt (Rubatto 2002). Both increased mafic component and “open” system processes point to a crustal magmatic system that is dynamic and can accommodate fluxing of new, mafic melts and increased dissolution efficiency during crustal anatexis (e.g., Wang et al. 2018; Storck et al. 2020). On a local scale, zircon crystals from these types of environments are likely to rapidly crystallize under non-equilibrium (i.e., disequilibrium) conditions which can explain the chaotic growth zoning and high variability of Th/U in zircon from VMS-proximal, back-arc melts (Wang et al. 2011; Kirkland et al. 2015); in contrast, equilibrium conditions allow for slower, protracted crystallization of zircon defined by partition coefficients, which is the case for zircon with lower and more consistent Th/U that were formed in VMS-distal arc melts. Zircon in VMS-proximal also have higher aspect ratios compared to zircon in VMS-distal rocks, which we interpret to reflect differences in crystallization temperature and/or rate of crystal growth in the respective melts (e.g., Corfu et al. 2003). Arc zircon, however, show well-developed oscillatory (± sector) zoning in relatively equant crystals that have more restricted ranges of Th/U, indicative of slower cooling in more equilibrated melts. If temperature is a dependent variable that directly correlates with Th/U (and zircon morphology), as others have suggested, we would expect other temperature-dependent variables to corroborate the hypothesis. We therefore utilize the Ti-in-zircon thermometer (Ferry and Watson 2007) and in situ Hf isotopes in zircon, which directly give us information about the crystallization temperatures and relative abundances of juvenile (i.e., hotter) melts that contribute to the formation of zircon in the Yukon-Tanana magmatic rocks.

Ti-in-zircon thermometry

The Ti-in-zircon geothermometer has been utilized in numerous studies to establish criteria for magma fertility in genetic and exploration models of mineralized systems such as porphyry Cu (e.g., Lee et al. 2020; Viala and Hattori 2021), iron oxide Cu-Au (Wade et al. 2022), and intrusion-related Au systems (Mathieu et al. 2022); however, its use is limited in studies of VMS deposits (Codeço et al. 2018). The minimum crystallization temperatures (T_zrc) can be estimated using Ti concentrations in zircon obtained by LA-ICP-MS methods using the methods of Ferry and Watson (2007):

$$\mathrm{log }(\mathrm{ppm Ti}-\mathrm{in}-\mathrm{zircon}) = (5.711 \pm 0.072)-(4800 \pm 86) / {T}_{\mathrm{zrc}}(\mathrm{K}) -\mathrm{ log }{a}_{\mathrm{SiO}2} +\mathrm{ log }{a}_{\mathrm{TiO}2,}$$

where T_zrc is dependent on the measured Ti concentration in zircon (ppm Ti-in-zircon), and a_SiO2 and a_TiO2 in the melt, which are unity (= 1) in the presence of rutile and quartz (Watson and Harrison 2005; Watson et al. 2006; Ferry and Watson 2007). All rocks in this study contain co-existing zircon, rutile, and quartz, thus a_SiO2 and a_TiO2 are assumed to be unity for our calculations. Specifically, rutile and quartz are present as inclusions or intergrowths with zircon (e.g., Fig. 6), or where these are not present, they occur in the same field of view as zircon (i.e., < 100 µm away). In Yukon-Tanana terrane, the minimum crystallization temperatures for zircon in the back-arc (T_zrc = 634–912 °C) are distinctly higher than those in the arc environment (T_zrc = 471–838 °C; Table 2; Figs. 7, 9). Back-arc plutonic rocks have zircon with slightly lower T_zrc (648–912 °C; mean = 767 °C) relative to VMS-distal volcanic zircon (634–907 °C; mean = 778 °C), whereas VMS-proximal zircon have the highest values (679–900 °C; mean = 790 °C); however, all average zircon values have 2σ uncertainties of ± 13 to 30 °C. A broad, positive correlation is observed between T_zrc and Th/U in back-arc zircon only, which supports the hypothesis that zircon morphologies and Th/U ratios are temperature-dependent in the back-arc setting (Fig. 9C). The arc zircon, however, display relatively constant Th/U over a range of ~ 250 °C and are thus not temperature dependent, which is interpreted to indicate zircon crystallization largely under conditions that reflect equilibrium partitioning of Th and U between melt and zircon (Fig. 9C).

Hafnium isotopes in zircon

Zircon in VMS-proximal back-arc rocks have higher age-corrected Hf isotopic compositions (median εHf_i =  − 6.2) compared to VMS-distal volcanic rocks (median εHf_i =  − 12) and show a positive correlation with Th/U and T_zrc (Fig. 9). The εHf_i also correlates with zircon crystallization temperature (as above), where low εHf_i zircon grains crystallized from colder, more fractionated magmas (i.e., VMS-distal) and higher εHf_i zircon crystallized from hotter, less fractionated magmas (i.e., VMS-proximal). Increased εHf_i in VMS-proximal versus VMS-distal zircon could be a result of local dissolution of pre-existing, xenocrystic zircon from increased juvenile magma input that subsequently increased the Lu/Hf in the melt and crystallized new, autocrystic zircon with increased ¹⁷⁶Hf/¹⁷⁷Hf (i.e., increase εHf_i; Wang et al. 2018). The back-arc plutonic rocks generally show the same trends as the zircon from volcanic rocks, but also contain zircon that overlap with the VMS-proximal field (Fig. 9). We therefore interpret the extrusive rocks to record singular crystallization events, whereas there was protracted zircon growth and crystallization in the plutonic rocks. Arc samples, however, do not follow a temperature-dependent correlation. Zircon grains from Glenlyon and Stewart River areas have the highest εHf_i of our sample suite (~ chondrite), yet have Th/U similar to more evolved arc zircon. This observation confirms our hypothesis that arc zircon form from equilibrium element partitioning, whereas back-arc zircon form from temperature-controlled disequilibrium crystallization. These results may also point to significant differences in the composition of crustal substrates and basement domains that magmas ascended through in the arc environments of Yukon-Tanana terrane.

Both back-arc and arc magmatic zircon display ranges of εHf_i within a single sample that spans ~ 5–15, and ~ 3–7 epsilon units, respectively (Fig. 11). Where available, whole-rock Hf isotope results from the same rocks overlap within the range of zircon isotopic values, indicating that the whole-rock compositions incorporate a mixture of Hf signatures (Fig. 11B; Manor et al. 2022b), especially inherited zircon that evolve to more negative epsilon values over time (Fig. 11A). In arc magmatic systems, epsilon unit variations of ~ 4–6 units for individual samples have been attributed to magma mixing processes due to variable abundances of juvenile magmas and evolved crustal components in the magmas (Broderick et al. 2015; Attia et al. 2020). Our results suggest that similar mixing between mafic and crustal melts was involved in the formation of Yukon-Tanana terrane magmas, with more pronounced variation in the back-arc region. To determine the relative abundance of juvenile and crustal materials to produce the observed range of Hf isotopic signatures and fractionation indices (e.g., Th/U, Zr/Hf, U/Yb, Hf), a binary mixing model was employed similar to Manor et al (2022b). We assumed a basaltic underplate model (e.g., Huppert and Sparks 1988) using endmember compositions of N-MORB basalt (i.e., asthenospheric mantle) and evolved Laurentian crust that comprises the basement to Yukon-Tanana terrane (i.e., Snowcap assemblage; Piercey and Colpron 2009). These parameters were used in the mixing equation of Faure (1986):

$${\mathrm{R}}_{\mathrm{x}}^{\mathrm{m}} \, \mathrm{=} \, \left[{\mathrm{R}}_{\mathrm{x}}^{1}{{\mathrm{C}}}_{1}{\mathrm{X}}\mathrm{+}{\mathrm{R}}_{\mathrm{x}}^{2}{{\mathrm{C}}}_{2}\left({1}-{\mathrm{X}}\right)\right]\mathrm{/}\left[{\mathrm{C}}_{1}{\mathrm{X}}\mathrm{+}{\mathrm{C}}_{2}\left({1}-{\mathrm{X}}\right)\right],$$

where ${\mathrm{R}}_{\mathrm{x}}^{\mathrm{m}}$ is the ratio of the mixture; ${\mathrm{R}}_{\mathrm{x}}^{1}$ and ${\mathrm{R}}_{\mathrm{x}}^{2}$ are the isotopic compositions of the endmembers (i.e., N-MORB: εHf_{360 Ma} =  + 14.5, Vervoort and Blichert-Toft 1999; and age-corrected Snowcap assemblage: εHf_i ~  − 26.2, Piercey and Colpron 2009); ${\mathrm{C}}_{1}$ and ${\mathrm{C}}_{2}$ are the trace element concentrations of the endmembers (i.e., N-MORB: Hf = 2.05 ppm, Sun and McDonough 1989; Snowcap: Hf = 6.7 ppm, Piercey and Colpron 2009); and X is the mixing factor that defines the volume ratio of mantle to crustal components. It must be noted that these results represent Hf components only and do not reflect bulk assimilation percentages per se. Isotopic modelling results for back-arc zircon show a wide range of juvenile (46–89%) versus evolved (11–54%) components recorded by εHf_i. Within this range, juvenile, mantle-derived melts contribute a median of ~ 75–78% Hf to magmas that generate felsic VMS-proximal stratigraphy, compared to ~ 61–69% Hf in VMS-distal rocks. In the arc environments, zircon from Simpson Range plutonic suite granitoids in the Finlayson Lake district and Reid Lakes batholith give wide component ranges (35–76% MORB and 24–65% crust) compared to arc volcanic and Simpson Range plutonic suite plutonic rocks in the Glenlyon and Stewart River regions (75–88% MORB and 12–25% crust). There are two primary outcomes of these isotopic and modelling results: (1) elevated εHf_i and relative proportion of juvenile components correlate well with high Th/U and Zr/Hf, and low U/Yb and Hf concentrations, which suggest both temperature and mafic input contributed to the distinct geochemical and isotopic signatures of VMS-proximal zircon; and (2) Hf isotopic compositions for arc zircon are highly variable and distinguishable by region in Yukon-Tanana terrane.

Hafnium isotopic compositions for zircon in Stewart River and Glenlyon regions are significantly less evolved than other arc-affinity rocks in the rest of Yukon-Tanana terrane (Fig. 11). These rocks contained zircon with slightly evolved to chondritic εHf_i (median =  − 6 to 0) and have relatively young Hf model ages (T_DM = 0.8–1.0 Ga). In both areas, crustal contamination is considered a common process in the formation of felsic rocks due to the presence of both inherited zircon grains and Proterozoic to Archean depleted mantle model Nd ages (Colpron et al. 2006a; Ruks et al. 2006; Simard et al. 2007) and the proximity to type-locality rocks of the evolved Snowcap assemblage (e.g., Piercey and Colpron 2009). More evolved isotopic signatures might be expected as in the Finlayson Lake district if ascending melts interacted with highly evolved Yukon-Tanana terrane basement domains; however, felsic rocks with the least evolved εHf_i occur at ca. 354.2 Ma and become more evolved in other parts of Yukon-Tanana terrane with time (Fig. 11B). Therefore, the more chondritic Hf isotopes in both the Glenlyon and Stewart River regions are interpreted to reflect the following: (1) there was more mafic input into the zircon Hf budget in both regions; (2) the assimilant and crust was less evolved in these regions compared to the Finlayson Lake district; and/or (3) there was a heterogeneous distribution of basement domains along the Early Mississippian Laurentian margin with a range in juvenile to evolved isotopic compositions.

Oxidation state of Yukon-Tanana parental magmas

Zircon is an optimal candidate for assessing the oxidation state of the parental magmas due to its robust ability to limit elemental diffusion during post-magmatic hydrothermal or metamorphic processes (Hoskin and Black 2000; Belousova et al. 2002). Both Ce and Eu show unique chondrite-normalized REE patterns, often as positive and negative anomalies, respectively, and are proposed to reflect the oxidation state of the crystallizing magmas (e.g., Ballard et al. 2002; Burnham and Berry 2012), where higher anomalies have been attributed to more oxidizing conditions in magmas (e.g., Kobylinski et al. 2020; Lee et al. 2020). Europium anomalies have also been proposed as a proxy for sulfur speciation and degassing processes in mineralized arc magmas (Dilles et al. 2015); however, Eu anomalies may also record the effects of plagioclase fractionation in the melt, especially in the absence of significant water contents (i.e., > 3 wt%); thus, negative Eu anomalies in zircon may be accentuated when both processes contribute (Dilles et al. 2015).

Zircon from Yukon-Tanana terrane rocks have distinct Ce and Eu anomalies depending on the tectonomagmatic environment (Table 2; Fig. 10; ESM 2). Cerium anomalies are typically lower in the back-arc zircon samples compared to the arc, but there is no noticeable difference between VMS-proximal and VMS-distal zircon (Table 2; ESM 2). The higher Ce anomalies in the arc zircon are attributed to an increase in magmatic water content that both creates more oxidizing conditions and reduces the liquidus temperatures of silicate magmas (Holtz and Johannes 1994), which correlates with the lowest T_zrc (> 550 °C) and suggests more water-saturated conditions during zircon crystallization in the Yukon-Tanana terrane arc magmas. Conversely, the low water content of back-arc magmas indicates higher temperature melts that would have had potential to provide added heat to the upper crust if the permeability and structural characteristics of the upper crust permit, thus being potentially prospective for VMS deposit genesis (e.g., Hyndman et al. 2005; Piercey 2011). Furthermore, lower water content in back-arc magmas likely enabled plagioclase to co-crystallize with zircon at much higher temperatures, and thus more likely overlap with temperatures of zircon saturation and crystallization (~ 950–650 °C; Fig. 7). Therefore, both Ce and Eu anomalies in Yukon-Tanana terrane rocks can distinguish relative differences between arc and back-arc magmas but only Ce may be used valid indicator of relative oxidation state due to highly reducing conditions.

Oxygen fugacity (fO₂) can be derived from Ce concentrations in zircon based on empirical calibrations (Trail et al. 2011, 2012; Smythe and Brenan 2015, 2016; Loucks et al. 2020) but has only been applied to natural zircon samples from arc settings (e.g., Kobylinski et al. 2020; Viala and Hattori 2021). We apply the model of Smythe and Brenan (2016) to model the primary fO₂ of both arc and back-arc zircon (and thus their parental melts) in Yukon-Tanana terrane:

$$\mathrm{ln}\left[\frac{\textit{x}_{\mathrm{Ce}^{4+}}^\mathrm{melt}}{\textit{x}_{\mathrm{Ce}^{3+}}^\mathrm{melt}}\right]\mathrm{=}\frac14\mathrm{ln}\,\textit{f}{\mathrm{O}}_2\mathrm{+}\frac{13136\left(\pm591\right)}{{\mathrm{T}}_\mathrm{zrc}}\mathrm{-2.064}\left(\pm0.011\right)\frac{\mathrm{NBO}}{\mathrm{T}}\mathrm{-8.878}\left(\pm0.112\right)\cdot\textit{x}{\mathrm{H}}_2\mathrm{O-8.955}\left(\pm0.091\right)\mathrm{,}$$

where Ce⁴⁺/Ce³⁺_melt was calculated based on the method of Smythe and Brenan (2015) using the lattice strain model of Blundy and Wood (1994) and partition coefficients derived from Ce concentrations in both zircon and the melt (i.e., whole-rock lithogeochemistry; Manor et al. 2022b; ESM 1); T_zrc is Ti-in-zircon crystallization temperature in Kelvin; NBO/T is the ratio of non-bridging oxygens to tetrahedral coordinated cations in the hydrous silicate melt (calculated from whole-rock lithogeochemical molar ratios); and xH₂O is the mole fraction of water present in the melt. Whole-rock lithogeochemistry and Ti-in-zircon thermometry inherently introduce additional error to the fO₂ calculations because of assumptions made for melt compositions and Si-Ti activity coefficients, respectively; however, these are best estimates for samples in this study. Water contents were assumed to be 1.0 ± 0.5 wt% and 4.0 ± 0.5 wt% for back-arc and arc rocks, respectively, based on estimates from modern back-arc and arc settings globally (Zimmer et al. 2010); an increase of 1.0 wt% H₂O is expected to raise the median fayalite-magnetite-quartz buffer (∆FMQ) by ~ 0.8 log units (Frost 1991). All arc plutons in this study contain primary amphibole, which alone suggests arc magmas contain ~  > 3.5 wt% H₂O (Naney 1983; Dilles et al. 2015). Uncertainties associated with this modelling span ~ 2–3 log units because of the large range in H₂O contents assigned (e.g., Kobylinski et al. 2020).

The fO₂ results from Smythe and Brenan (2016) were compared to the oxybarometer of Loucks et al. (2020), which is based completely on zircon compositions and does not require the incorporation of other datasets or estimation of water contents:

$$\log\;f{\mathrm O}_{2\;\left(\mathrm{sample}\right)}-\log\;f{\mathrm O}_{2\;\left(\mathrm{FMQ}\right)}=3.99\;\left(\pm0.12\right)\,\cdot\log\;\left[\frac{\mathrm{Ce}}{\sqrt{\left(U_{\mathrm i}\times\mathrm{Ti}\right)^Z}}\right]+2.28\;\left(\pm0.10\right),$$

where U_i is the age-corrected initial U concentration (using U–Pb zircon dates reported in this study and Manor et al. 2022a, b), and the superscript z indicates zircon. Internal uncertainties for our dataset were calculated using compositions of the Plesoviče zircon standard (n = 302; ESM 1), which yield 2σ uncertainty of ± 0.96 log units.

Oxygen fugacity calculations reveal that back-arc zircon formed from melts in reduced conditions (Smythe and Brenan: ∆FMQ − 5.1 to − 1.4 log units; Loucks: ∆FMQ − 2.3 to + 0.26 log units) relative to arc zircon (Smythe and Brenan: ∆FMQ − 1.2 to + 1.0; Loucks: ∆FMQ − 0.43 to + 1.2 log units; Fig. 12). Mafic volcanic systems from the modern seafloor and subaerial volcanoes show significant variation in oxygen fugacity of primary magmas in arc (∆FMQ + 0.5 to 2.0), back-arc (∆FMQ − 0.7 to + 0.5), and mid-ocean ridge settings (i.e., MORB; ∆FMQ − 2 to 0; Wallace and Carmichael 1994; Bézos and Humler 2005; Kelley and Cottrell 2009; Rowe et al. 2009; Zimmer et al. 2010). Crustal assimilation of Snowcap assemblage and/or North American margin strata has been proposed as an important factor in generating trace element and isotopic signatures in felsic rocks in the Finlayson Lake district (Piercey et al. 2003; Manor et al. 2022b), and it also contains graphitic or carbonate-rich lithologies, such as minor marble and graphite-bearing quartzite (Ryan and Gordey 2002; Colpron et al. 2006a) that could have reduced the original fO₂ of incoming peraluminous melts throughout Yukon-Tanana terrane by multiple log units (e.g., Nixon 1998; Mollo et al. 2010; Tomkins et al. 2012; Bucholz et al. 2018). Coincidentally, the presence of native metals (e.g., Fe, Zn, Pb, Cu), various Fe-Ti oxides (e.g., ilmenite, magnetite-wüstite?), and pyrrhotite as inclusions in zircon also supports a localized, highly reduced magmatic environment (~ iron-wüstite buffer) in the Yukon-Tanana terrane back-arc melts (e.g., Toplis and Carroll 1995; Distler et al. 2004). We therefore propose that the observed fO₂ from both models applied above to Yukon-Tanana terrane zircon (arc and back-arc derived) corresponds to crustal assimilation of graphite and/or carbonate-bearing lithologies that produced localized reduction by ~ 1–3 log units of parental melts.

Evaluation of magma petrogenesis in the formation of felsic-hosted VMS deposits

This study shows zircon chemistry can constrain important aspects of magma petrogenesis and its relationships to felsic-hosted VMS deposits. Zircon records differences in oxidation state and sulfur speciation for magmas in the back-arc (FMQ − 5 to 0, low Ce_N/Ce_N*) and arc (FMQ − 1 to + 1, high Ce_N/Ce_N*) tectonic settings (Table 2; Fig. 12). Results from the Loucks et al. (2020) oxybarometer also show that VMS-proximal zircon may have crystallized from slightly more oxidized melts (median ∆FMQ − 0.65) relative to VMS-distal zircon (median ∆FMQ − 1.8; Fig. 12B, D). In the Yukon-Tanana terrane back-arc, the highly reducing conditions of the magmas likely inhibited significant magmatic degassing, stabilized S²⁻, and subsequently limited the fluid-related metal and sulfur (as SO₄²⁻) contributions to the VMS-related hydrothermal circulation cell (e.g., Dilles et al. 2015). Most of the sulfur present in the back-arc melts was likely segregated during the ascent of the high-temperature melts if they achieved sulfide saturation. If this segregation occurred at relatively shallow depths in the upper crust, sulfides crystallized from these magmas could have created metal-rich rocks that circulating hydrothermal fluids could have stripped metals from (e.g., Patten et al. 2016). Hydrothermal leaching of these pre-enriched sulfide accumulations would have generated metal-rich fluids that could then be added to the hydrothermal systems (± degassing fluids) that ultimately precipitated the sulfide-sulfosalt assemblages observed in felsic-hosted VMS deposits in the Finlayson Lake district (e.g., Peter et al. 2007; Bradshaw et al. 2008; Layton-Matthews et al. 2008). Therefore, sulfide saturation and metal pre-concentration (i.e., sulfide pre-collection) and later hydrothermal leaching are interpreted to be important in the formation of VMS deposits in the Finlayson Lake district and may be critical to the formation of VMS deposits globally.

Trace element characteristics in zircon also indicate that the relatively dry, back-arc parent magmas associated with VMS deposits contained plagioclase as a liquidus phase, but not hornblende, in shallow upper crustal settings (Fig. 7C). This petrological distinction has been attributed to source melting conditions and depth of melting within crustal settings, in particular related to FI–FIV rhyolite assemblages (Lesher et al. 1986; Hart et al. 2004). Within this classification scheme, FI magmas form at > 30 km depth with garnet as a common mineral phase during melting; FII rhyolite magmas are generated in the presence of amphibole-plagioclase-dominant melt sources at ~ 10–15 km depth; and FIII–FIV rhyolites are formed from plagioclase-dominant (and amphibole-garnet-free) melt sources shallower than 10 km depth. Our observations from zircon suggest that back-arc magmas should plot within the FIII rhyolite fields of Hart et al. (2004), with respect to mineral stabilities in P–T space; however, the whole-rock compositions plot between FI to FII fields (Manor et al. 2022b; Denisová and Piercey 2022a). This discrepancy again illustrates the pitfalls of utilizing the FI–FIV rhyolite classification scheme for rocks younger than Archean that are generated by the melting of continental crust (e.g., Piercey et al. 2008; Piercey 2011), despite suggestions to the contrary (Hart et al. 2004). The geochemical results do, however, overlap with other post-Archean rocks that were formed from dry melting conditions in intracontinental back-arc and extended continental margin tectonic settings (Fassbender et al. 2022). Therefore, while we cannot confirm the depth of crustal melting, we can infer that the magmas formed in relatively shallow upper crustal chambers based on the predominance of plagioclase fractionation in both whole-rock and zircon geochemical results (Figs. 7, 10; Manor et al. 2022b).

Melts that generated erupted or emplaced VMS-proximal facies were higher temperature, less fractionated, and contained on average ~ 6–17% more juvenile Hf component compared to those in VMS-distal stratigraphy (Fig. 13). The most important petrological characteristics that define VMS-proximal stratigraphy coincide with specific periods of magmatism during back-arc evolution and are summarized by the following: (1) high-temperature felsic magmatism and rapid zircon growth (e.g., increased zircon Th/U and T_zrc, and whole-rock Zr_sat); (2) mineral fractionation in crustal melts, such as high Th/U and Zr/Hf, and decreased U/Yb and Yb/Dy, which signifies variations in co-fractionating minerals by temperature and parental magma composition; (3) elevated juvenile melt input (e.g., increased zircon εHf_i and Th/U, and whole-rock Nb/Ta), which likely elevated the regional geothermal gradient, influenced the efficient assimilation and dissolution of pre-existing crustal-derived minerals, and led to increased HFSE-REE and Lu/Hf in melts; and (4) high modal abundance of HFSE-REE-bearing minerals in the whole rock, such as monazite, xenotime, and apatite, which directly correlates to HFSE-REE concentrations of the crustal melts. Manor et al. (2022b) interpreted the highest whole-rock Zr_sat, Nb/Ta, and Hf–Nd isotopic compositions for VMS-proximal stratigraphy, which is supported in this study by zircon results (e.g., Fig. 13), from increased contributions from basaltic underplating within a continental back-arc regime (Fig. 14). In the back-arc of the Finlayson Lake district, the Kudz Ze Kayah formation and Wolverine Lake group record discrete magmatic cycles observed on timescales of 2–3 Myr, separated by a period of deformation at ca. 360–358 Ma (Manor et al. 2022b); both magmatic cycles culminate with discrete horizons of VMS mineralization in the Kudz Ze Kayah and Wolverine VMS deposits and are more broadly interpreted to represent fine-scale magmatic episodicity during the magmatic evolution of Yukon-Tanana terrane (Figs. 13, 14; e.g., De Silva et al. 2015; Attia et al. 2020). We interpret these magmatic cycles in Yukon-Tanana stratigraphy as fine-scale magmatic flare-ups in the back-arc that directly correlate with the presence of a basaltic underplate. With time and evolution of the back-arc rift, the underplate and its magmatic derivatives elevated the heat flow and the geothermal gradients sufficiently to sustain hydrothermal circulation and VMS formation in the upper crust. Moreover, this integrated study utilizing zircon and whole-rock geochronology, geochemistry, and isotopic methods indicates that incremental additions of high-temperature, juvenile melts to extending back-arc crust creates an environment with high heat flow that is optimal for VMS genesis. Therefore, the intrusive timing and characteristics of magmas in the Finlayson Lake district are critical to delineating prospective horizons for VMS mineralization and may be a process that is critical to the formation of other large tonnage and high-grade, felsic-hosted VMS deposits globally.

Fig. 13

Summary of key chemostratigraphic variations in back-arc zircon and whole rocks of the Finlayson Lake district. A Column that illustrates the time-resolved stratigraphic relationships defined by recent high-precision CA-ID-TIMS U–Pb geochronology (Manor et al. 2022a, b). The locations of VMS deposits are highlighted in red text, and horizontal red bars encompass the 2σ error of each VMS mineralization horizon. Ages of intrusive rocks from the GLPS (red stars) and North Lakes intrusion (purple star) are shown for reference. B–H Box and whisker plots for each parameter are shown relative to the weighted mean U–Pb crystallization or eruption age of the rock; red boxes and red outlier points denote VMS-proximal samples. B εHf_i; C abundance of juvenile input (in %), calculated from two-component mixing equation (see discussion); D Th/U; E) T_zrc (°C), as in Ti-in-zircon crystallization temperature; F U/Yb; G) Zr_sat (°C), as in whole-rock Zr saturation temperatures with 2σ error bars; H whole-rock Nb/Ta, with reference to crustal (Nb/Ta = 12) and mantle (Nb/Ta = 17) values (Manor et al. 2022b)

Fig. 14

Conceptual tectonomagmatic diagrams illustrating the difference between A crustal-dominated magmatism related to VMS-distal rocks (i.e., pre- and post-VMS mineralization) is associated with relatively low heat flux, longer crustal residence times, and poor dissolution efficiency of country rocks that are unfavorable to initiation of hydrothermal circulation; and B mantle-driven flare-up events generate VMS-proximal rocks (i.e., syn-VMS mineralization) from high fluxes of mafic magmas that heat the crust and efficient dissolution of country rocks during assimilation, which results in elevated isotherms and sufficient heat to sustain hydrothermal circulation for VMS genesis. Red dashed lines with numbers represent schematic isotherms; isotherms with the same temperature have the same number (1, 2, or 3). Blue arrows in B show hypothetical meteoric water circulation paths that develop in response to heat from below. C Illustration of Hf isotopic data for VMS-distal (17MM-002) and VMS-proximal (18MM-133) samples in the Kudz Ze Kayah formation, highlighting the large variation in εHf_i for a single sample; this variation is present in all back-arc samples. Labels (εHf-1, εHf-a) correspond to random locations of Hf isotopic development in A and B. These magmas with different εHf_i eventually mix and are preserved in crystallized rocks. D Summary of key features in both high-aspect ratio pre-/post-VMS zircon and low-aspect ratios (i.e., equant) syn-VMS zircon crystals that are generated because of differing magmatic environments. See text for details

Conclusions

Arc and back-arc rocks in the Yukon-Tanana terrane record the onset of magma-rich rifting and significant VMS mineralization along the western Laurentian continental margin in the mid-Paleozoic. Our results provide a detailed and comprehensive study using zircon as a tool to decipher the role of magmatic petrogenesis in the formation of felsic stratigraphy that directly hosts VMS deposits in the Yukon-Tanana terrane. We have illustrated that zircon displays a variety of textures that reflect their tectonomagmatic setting and petrogenetic processes. Zircon geochemical and isotopic constraints further distinguish VMS-proximal from VMS-distal stratigraphy, and of most importance are crystallization temperatures (i.e., Ti-in-zircon thermometer), fractionation ratios (e.g., Th/U, U/Yb, Zr/Hf), and Hf isotopic compositions. Integrated with whole-rock geochemistry and isotopes, these results indicate that greater abundances of high-temperature, mantle-derived juvenile melts (~ 75–80%) were critical in generating hot upper crustal magma chambers, and ultimately VMS-proximal stratigraphy, whereas VMS-distal magmas contained much less juvenile input (< 70%). Efficient crustal melting and dissolution of HFSE-REE-bearing minerals from increased juvenile input distinguishes VMS-related from VMS-barren volcanic rocks. Oxygen fugacity modelling suggests that zircon grains in both arc and back-arc rocks crystallized from magmas that were reduced through assimilation of graphite and/or carbonate-bearing basement domains, such as the Snowcap assemblage, and favored sulfide fractionation and pre-collection of metal-bearing minerals/alloys as inclusions within zircon grains. As a result, the reducing conditions inhibited sulfur degassing, and there were likely minimal contributions of magmatic fluids to the hydrothermal circulation systems in VMS deposits in the Finlayson Lake district. Zircon chemistry shows secular evolution and records three distinct mantle-driven flare-up events in the back-arc region of Yukon-Tanana terrane on timescales between 1 and 3 Myr. Magmatic fluxing thus has important implications on the timing and localization of VMS deposits in Yukon-Tanana terrane, particularly in the Finlayson Lake district, and also places important constraints on the timing and nature of magmatism during a period of complex subduction initiation-related tectonics in the Late Devonian to Early Mississippian. Zircon is thus an important tool to deciphering the magmatic parameters in ancient arc and back-arc settings that can provide insights into regional metallogenic and tectonic processes related to VMS genesis.