Volcano-pluton connections at the Lake City magmatic center (Colorado, USA)

Establishing the mechanisms and timescales associated with the generation of silicic magmas is fundamental to understanding the evolution and eruption of voluminous rhyolitic magmas, the generation of Earth's continental crust, and the extent to which these processes are, or are not, related. The idea that large volumes of eruptible rhyolitic magmas are intrinsically related to upper-crustal granitic and granodioritic plutons has long been debated, but many workers have converged on models that support a relationship (e.g., crystal mush, Bachmann and Bergantz, 2004; Hildreth, 2004; transcrustal magmatic system, Cashman et al., 2017). These models suggest that rhyolites are generated through low-pressure crystallization of, and segregation from, less evolved intermediate magmas. This process results in crystal-dominated residues (plutons) and ephemerally stored melt-dominated eruptible rhyolites (e.g., Gualda et al., 2012; Gualda and Ghiorso, 2013; Cooper and Kent, 2014). On the other end of the spectrum are models that suggest there is no connection (e.g., Glazner et al., 2004, 2008; Tappa et al., 2011; Mills and Coleman, 2013). In these models, intermediate-to-silicic plutons are built incrementally over 10⁵–10⁶ years during periods of low magma flux, while large-volume rhyolites are generated in the lower crust, transported, stored, and erupted within 10³–10⁴ years, with little to no differentiation in the upper crust. These two models are end members that imply that extrusive and intrusive silicic rocks in a volcanic-plutonic complex are either wholly related (former model) or discrete (latter model), but some systems, particularly large and long-lived ones, may represent a hybrid scenario, with periods of pluton growth that are linked to volcanic activity and others that are not.

Many studies addressing this debate have relied on comparisons of volcanic and plutonic systems from global data sets (e.g., Glazner et al., 2008, 2015; Keller et al., 2015) or on individual plutonic or volcanic deposits (e.g., Mills and Coleman, 2013; Deering et al., 2016; Streck, 2014; Schaen et al., 2021). Localities where seemingly coeval volcanic and plutonic sections are juxtaposed give a more direct glimpse into these relationships, but well-exposed and preserved examples are rare. This study focuses on the Lake City magmatic center (southwest Colorado), which is one of these exceptional places: The Lake City caldera was produced during the eruption of the compositionally zoned (high-silica rhyolite to trachyte) Sunshine Peak Tuff (SPT) and was subsequently intruded by intermediate composition magmas (syenite and monzonite) during post-caldera resurgence (Hon, 1987; Kennedy et al., 2012). Both the reversely stratified ignimbrite units and resurgent hypabyssal intrusions, which have been suggested to be a cumulate residue (Hon, 1987; Kennedy et al., 2016; Lubbers et al., 2020), are exposed within the caldera, providing an excellent opportunity to explore relationships between the extrusive and intrusive magmas.

In this work, we aim to build upon and complement several previous studies on field relations (Hon, 1987; Hon and Lipman, 1989; Kennedy et al., 2016), geochronology (Reynolds et al., 1986; Bove et al., 2001), whole-rock and phenocryst geochemistry (Kennedy et al., 2016; Lubbers et al., 2020), and phenocryst textures (Lubbers et al., 2020) within the Lake City magmatic system. We employ high-precision zircon geochronology and geochemistry of individual zircon crystals, determined using U-Pb chemical abrasion, isotope dilution, thermal ionization mass spectrometry, and trace-element analysis (CA-ID-TIMS-TEA, hereafter referred to as TIMS-TEA, Schoene et al., 2010), and zircon trace-element modeling to assess the volcanic-plutonic connection and the origin and longevity of eruptible high-silica rhyolite at the Lake City magmatic center.

The models described above provide hypotheses that can be further tested with zircon geochronology and geochemistry. Those in favor of volcanic-plutonic connections suggest there should be little difference in the ages of zircon grains crystallized in extrusive and intrusive units, zircon age dispersion of samples through an intrusive section should be similar, and the geochemistry of zircon grains in extrusive and intrusive units may be overlapping or complementary. Opposing models require no age or geochemical relationship between zircon grains in extrusive and intrusive units and suggest that different portions of intrusive units will show different age dispersion and, potentially, geochemistry. These expectations are equally relevant for the hybrid models.

At the heart of these arguments is whether the magmas that produced intrusive and extrusive deposits crystallized contemporaneously. The existence of an age relationship between them is a first-order test: if the units do not overlap in age, then there is no support for a direct relationship. From there, zircon geochemistry can provide additional detail on the mineralogical and chemical evolution of the system. Upper-crustal magma bodies are generally thought to exist on timescales of <500 k.y. (e.g., Petford et al., 2000; Michel et al., 2008; Schoene et al., 2012; Wotzlaw et al., 2013; Schaen et al., 2021); therefore, discerning such age relationships requires high-precision dating (ideally ≤~100 k.y. [2σ]). However, zircon crystals commonly display complex zoning patterns in cathodoluminescence (CL) images, which require interrogation through intragrain compositional analysis. Thus, the ideal methodological approach to assess these models with zircon would be a technique that provides high precision, in situ data.

Unfortunately, no such method exists. In situ geochronological methods (secondary ion mass spectrometry [SIMS] and laser ablation–inductively coupled mass spectrometry [LA-ICPMS]) cannot currently provide ages with the requisite precision (limited to ~2%–4%, or ~400–800 k.y. [1σ] on rocks of similar age and U content to those studied here; e.g., Lidzbarski, 2014), and spot sizes are rarely small enough to analyze single zones in a crystal. For young systems, TIMS-TEA can produce ages with such precision (0.001%–0.01% [2σ]), but the method requires full dissolution of zircon grains or grain fragments; so there is little to no spatial resolution achievable. All of these methods produce results that are necessarily weighted averages, with volumetrically larger zones imparting greater influence on the results; thus, results are generally weighted toward younger ages. Consider, for example, the idealized case of a crystal with zones that are concentric and of equal thickness—edge zones will be more voluminous and contribute more to the result than core zones. Of course, zircon zoning is rarely so simple; therefore, this relationship can be much more complex.

Despite these complications, the age precision achievable with TIMS-TEA makes it our preferred method for investigating the intrusive-extrusive relationship. Such precision makes assessment of dispersion within a sample non-trivial, and the averaging effect of the method likely results in underestimates of the duration of zircon crystallization. Nonetheless, as long as units have overlapping ages, there is first-order support for a relationship between them.

From a geochemical standpoint, zircon geochemistry determined with TIMS-TEA will also be weighted toward compositions from volumetrically larger zones and zones with relatively extreme concentrations of a given element. The CL zoning of many zircon crystals is complex, and this method cannot resolve the rich histories these textures imply. Nonetheless, zircon grains that grew under distinct magmatic conditions, from distinct magmatic reservoirs, and/or at distinct moments in time, should be geochemically distinguishable from each other, even if they have some amount of shared history. Furthermore, the TIMS-TEA method averages the same volume of material for both the age and geochemical analysis. The result is that broad geochemical relationships and time-composition relationships can be reasonably considered using this method. There is also opportunity for combining methods, because sectioned grains can be analyzed by in situ methods for geochemistry prior to dissolution. It is worth noting, however, that sectioning of grains may come at the cost of some age precision, because it will leave a smaller amount of Pb to measure, and the 3D aspect of zircon zoning, which is often in sectors, will not be integrated in the in situ results. Examples of this approach and comparisons between in situ and TIMS-TEA data are given by Samperton et al. (2015) and Eddy et al. (2022).

It is also worth considering the robustness of the zircon record for addressing this problem at all. For all mineral phases, the record maintained by a given phase is restricted to the period of time when it is saturated in a magma. As such, the time that a system existed in a molten state could be underestimated from zircon geochronology, and geochemical trends in zircon could be missing, if there is a substantial amount of time that a magma spends above its zircon saturation temperature. Quantifying such missing time can be challenging because (1) parental magma compositions from which zircon saturation temperatures can be calculated are often poorly known; (2) zircon saturation thermometry of more mafic parent compositions may be misapplied; (3) the extent of thermal fluctuations between super- and sub-zircon saturation incurred by a magma (e.g., number and thermal effect of recharge events) is challenging to constrain; (4) magma liquidus temperatures for comparison with zircon saturation temperatures are difficult to estimate because they are strongly affected by volatile contents, which themselves are often poorly known; and (5) it is not possible to know the amount of resorption (i.e., lost geochemical information) that has occurred. All of these issues can be problematic for both intrusive and extrusive magmas to varying extents, and in this work we do not attempt to determine the potential missing time. As such, our estimates of zircon age dispersion and magma residence time may be minimums.

On the other hand, residence times of eruptible magmas may be overestimated, and pseudo-geochemical trends may be resolved, if a zircon-saturated system spends a significant amount of time in a subsolidus state. This may be important for explaining some of the discrepancies that have been identified between magma residence times estimated from zircon geochronology and those from geospeedometry in other phases (e.g., Allan et al., 2013; Cooper and Kent, 2014; Till et al., 2019).

Finally, sampling density and sampling bias can have an effect on the results that is independent of the method employed. Plutonic complexes and large-volume ignimbrite deposits can be enormous, and low sampling density or limits to sampling the entirety of a deposit (e.g., poor exposure and erosion) could hide the existence of age and geochemical relationships and could result in an underestimate of the duration of magmatic activity at a magmatic center (e.g., Matzel et al., 2006).

The Lake City magmatic center (Fig. 1) is located in the southern Rocky Mountains (southwest Colorado, USA) and sits in the western side of the Oligocene San Juan volcanic field. The Lake City caldera, formed during the ca. 23 Ma SPT super-eruption, is the youngest of the 25 Tertiary calderas in the Southern Rocky Mountain Volcanic Field (SRMVF; Lipman, 1976; Steven and Lipman, 1975). Its structure is well exposed, and the SPT makes up the bulk of the intracaldera fill; little extracaldera ignimbrite exists (Hon, 1987; Kennedy et al., 2016). The composition of the magmas erupted at the Lake City magmatic center is more varied than at other magmatic centers in the SRMVF; this variety may be related to the onset of extension and opening of the Rio Grande Rift in the region (Lipman et al., 1978). Large megabreccia wedges separate the ignimbrite into three distinct eruptive units (all volumes in dense rock equivalent)—Lower SPT (LSPT; 200–250 km³); Middle SPT (MSPT; 20–40 km³); Upper SPT (USPT; 40–60 km³); each unit is composed of multiple welded flow units (Hon, 1987). Post-caldera lavas, dikes, and hypabyssal intrusions outcrop in the center and on the eastern flank of the caldera; intrusions are associated with resurgence and ring fault reactivation (Hon, 1987; Kennedy et al., 2012). There is no field evidence for significant time breaks between deposition of the ignimbrite units (Hon, 1987; Hon and Lipman, 1989), suggesting the SPT eruption was relatively continuous, and ⁴⁰Ar/³⁹Ar ages and paleomagnetic data suggest that the entire caldera cycle (LSPT eruption to intrusions) lasted 80–330 k.y. (Reynolds et al., 1986; Bove et al., 2001).

The SPT deposit is reversely stratified, varying systematically in bulk composition from high-silica rhyolite (LSPT) to trachyte (USPT) (Hon, 1987; Kennedy et al., 2016). The same phenocryst assemblage of alkali feldspar, quartz, plagioclase, biotite, and clinopyroxene is found among all of the eruptive units, but the abundance of each phase varies between them (Hon, 1987; Kennedy et al., 2016). The crystal contents within each unit are also variable, and crystals show evidence of late-stage resorption, particularly in the USPT (Kennedy et al., 2016; Lubbers et al., 2020). In detail, the LSPT is unzoned, with a homogeneous high-silica rhyolite bulk composition (76 wt% SiO₂) and up to 40% crystals; the crystal content of pumice clasts in this unit ranges from 16% to 40%. The MSPT is heterogeneous, with a low-silica rhyolite bulk composition (74 wt% SiO₂), up to 45% crystals, and two pumice types—rhyolite (LSPT composition) and trachyte (65 wt% SiO₂). The USPT is heterogeneous, with trachyte bulk composition (67–69 wt% SiO₂), up to 50% crystals, four pumice types—rhyolite (LSPT composition), trachyte (MSPT trachyte composition), “high rare-earth element” trachyte, trachyandesite, and trachytic and trachyandesitic enclaves (Kennedy et al., 2016; Lubbers et al., 2020). Accessory minerals in all units include zircon, sphene, apatite ± chevkinite (Kennedy et al., 2016). Quartz grains are rounded or embayed in all units and are rare in the USPT trachyte and trachyandesite pumice clasts (Kennedy et al., 2016; Lubbers et al., 2020). The Al contents of quartz grains in the LSPT (<200 ppm Al) are also much more constrained than those in other units (up to ~350 ppm Al), and some alkali feldspar and quartz grains in the USPT have Ba- and Ti-enriched rims, respectively (Lubbers et al., 2020). Glomerocrysts containing one or two feldspar phases are common in the MSPT and USPT (Lubbers et al., 2020).

The intracaldera resurgent intrusion is syenite, with similar bulk composition (~65–69 wt% SiO₂) and crystallinity to the USPT. It has a phenocryst assemblage of alkali feldspar plus minor quartz and crystal clots containing clinopyroxene, biotite, and Fe-Ti oxides (Kennedy et al., 2016). Feldspar crystals range from euhedral to anhedral; anhedral crystals commonly have reversely zoned rims (Kennedy et al., 2016; Lubbers et al., 2020). Textural signs of magma mixing in the syenite include trachytic and trachyandesitic enclaves with crenulated-to-diffuse edges and feldspar crystals that coexist between an enclave and its host (Kennedy et al., 2016). Intrusions on the eastern flank of the caldera are monzonite (~62–69 wt% SiO₂), with a phenocryst assemblage and enclaves similar to the syenite and SPT units, but they lack quartz. Dacite lavas (~63–66 wt% SiO₂) erupted on the eastern margin of the caldera after its collapse and bookend emplacement of the resurgent intrusions.

Early work on the Lake City deposits found significant similarities between the SPT volcanics and the syenite intrusion, leading to the suggestion that they were derived from a similar source. In contrast, the dacite lavas and monzonite intrusion on the eastern side of the caldera are distinct and were inferred to have come from a separate source (Hon, 1987; Hon and Lipman, 1989). This hypothesis was supported by evidence in geophysical studies for a large, buried intrusion, which was interpreted to be the crystallized remnant of the magma chamber that produced the dacite and monzonite magmas (Grauch, 1985, 1987; Hon, 1987).

Kennedy et al. (2016) further assessed the geochemistry of bulk tuff and individual pumice clasts from the SPT and post-caldera intrusions. While they found no systematic compositional variations within the individual ignimbrite units, varying proportions of mechanically mixed pumice types in the MSPT and USPT units led them to suggest that at least two magma types from distinct magmatic systems are represented in the Lake City deposits—an evolved batch A, represented by the rhyolite and trachyte ignimbrite units, a post-caldera rhyolite intrusion, and the resurgent syenite, and a more mafic batch B, represented by the trachyandesite pumices and enclaves, the post-caldera dacite lavas, and the resurgent monzonite. They also found that the USPT and syenite units are enriched in compatible elements relative to the MSPT and LSPT. Based on rhyolite-MELTS and trace-element modeling, they suggest that 35%–70% fractional crystallization of a trachyte composition could produce the MSPT and LSPT compositions.

Ultimately, Kennedy et al. (2016) proposed a “crystal mush”–like model for the Lake City magmatic center, in which: (1) several independent pods of (primary) trachyte magma initially existed and crystallized to generate a crystal-rich mush with interstitial high-silica rhyolite melt; (2) the high-silica rhyolite was segregated from one of these pods and erupted (LSPT), leaving behind both primary trachyte magma and a more mafic trachytic mush; (3) intrusion of trachyandesite magma into the system drove remobilization and eruption of some remaining rhyolite and portions of the mush (MSPT); (4) additional residual rhyolite, primary trachyte magma, trachyte mush, and trachyandesite magma were erupted (USPT); and (5) residual rhyolite (LSPT composition) and syenite (trachyte mush) intruded the caldera center, while dacite lavas erupted and monzonite magmas intruded the eastern side of the caldera. Lubbers et al. (2020) expanded this model using information from glomerocryst textures, crystal-size distributions, major-phase (quartz, feldspar, and biotite) geochemistry and geothermometry, and additional rhyolite-MELTS modeling. They also concluded that: (1) the rhyolitic LSPT and MSPT magmas crystallized from melts extracted from a less silicic (trachytic/syenitic) mush at 50–70 vol% crystals; (2) the USPT magma was formed by partially melting and erupting the feldspar-rich cumulate mush, the remainder of which is represented by the syenite; and (3) the melting event and triggering of eruption was caused by injection of more mafic magma into the mush. In addition, they suggested that the initial primary trachyte magma body formed and crystallized incrementally and was thermally buffered at 40–50 vol% crystals. This buffering maintained a crystallinity in the system amenable to extraction of interstitial rhyolitic melt for 10³–10⁴ years, allowing for the generation of a rhyolitic cap that underwent fractionation to establish the chemical zoning evident in the MSPT and LSPT deposits. Thus, they suggest the stratification in the bulk composition of the erupted deposits is due to in situ fractionation, and the system was composed of two magma batches (equivalent to batches A and B of Kennedy et al., 2016). Furthermore, they suggested that recharge events that occurred after crystal-melt separation produced localized areas of melting that remobilized part of the mush (USPT) and generated high-Ba rims on sanidine and high-Ti rims on quartz.

We collected bulk tuff samples from the Upper (SPU, 40B), Middle (SPM), and Lower (48C) SPT, as well as the resurgent syenite (SY-2A, 1C, and 10B) and monzonite (AGQM) units. All samples except SPU, SPM, and AGQM were also studied previously by Lubbers et al. (2020) and Kennedy et al. (2016). Analyzing zircon crystals from pumice clasts, rather than bulk tuff, is preferred, but the fact that we find few xenocrysts in our analyses (see Zircon U-Pb TIMS section) is reassuring and supports our approach of using bulk tuff.

Zircon crystals were separated from each sample using standard techniques. Samples SPU, SPM, and AGQM were first crushed using a BICO Chipmunk Crusher and a disk mill. The resulting material was sieved and hand-washed to remove clay and silt-sized particles. The remaining material was magnetically separated using a hand magnet and repeated runs on a Frantz magnetic separator. Heavy liquids (methylene iodide) were used to separate zircon crystals from light minerals and then hand-picked from the separate. All other samples were separated by hand-crushing, sieving, winnowing in water, and hand-picking. Bulk zircon separates were annealed in quartz crucibles at 900 °C for 60 h, and a subset of crystals was selected for further analysis; efforts were made to sample the range of crystal sizes and textures. Crystals were separated into rough size categories (large, medium, and small), mounted in EpoFix epoxy, ground to approximate centers, polished, and imaged by CL with a Gatan MiniCL detector attached to a FEI XL30 field emission gun–scanning electron microscope (FEG-SEM) at Princeton University.

Zircon ages and geochemistry were determined using chemical abrasion–isotope dilution–thermal ionization mass spectrometry–trace-element analysis (CA-ID-TIMS-TEA) at Princeton University. Here we summarize our approach to preparing and analyzing zircon solutions and our data reduction with this method, and additional details can be found in File 1 in the Supplemental Material¹. After CL imaging, zircon grains were removed from epoxy mounts and placed individually in Savillex microcapsules with a 29M HF + 30% HNO₃ solution to undergo chemical abrasion to leach domains affected by Pb loss and melt or mineral inclusions. After chemical abrasion, solutions were discarded, and zircon grains were rinsed, spiked with EARTHTIME ²³⁵U-²³³U-²⁰⁵Pb isotopic tracer (ET535; Condon et al., 2015), and dissolved in 29M HF. After dissolution, the solutions were dried down and dissolved in 6N HCl to convert the samples to chloride form. The final solution was dried down, brought back up in 3N HCl, and put through anion exchange column chemistry to distill U, Pb, and trace elements. Two solutions were captured during column chemistry: a U-Pb aliquot and a trace-element aliquot (i.e., the wash solutions from distilling U and Pb). We analyzed the U-Pb aliquots on the IsotopX Phoenix62 TIMS at Princeton University to obtain an age. Corrections were applied to the resulting data to account for mass-dependent fractionation of Pb and U, Pb and U deadtime on the Daly photomultiplier, interferences on Pb, and common Pb. Details of these corrections can be found in File 1. For each analysis, we made a correction for initial secular disequilibrium in the ²³⁸U-²⁰⁶Pb system due to the exclusion of Th during zircon crystallization (e.g., Schärer et al., 1984) using a zircon/melt partition coefficient of 0.29. This value was calculated using data from Claiborne et al. (2018) for a whole-rock composition comparable to the MSPT. Given that the MSPT is hypothesized to be a mixture of high-silica rhyolite and trachyte, we consider this to be a reasonable approximation of the parent melt. The effect of using this single value relative to other values (e.g., other single values based on other units or unit-specific values) is negligible, but a comparison can be seen in File 2. Final results from U-Pb TIMS can be found in File 3.

We analyzed the trace-element aliquots using solution ICPMS. Solutions where made by re-dissolving the dried down solutions from column chemistry in 3% HNO₃ + 0.2% HF + 1 ppb In. These solutions were then analyzed on a Thermo-Fisher iCAP quadrupole ICPMS at Princeton University. Measured elements include Y, Zr, Hf, rare-earth elements (REEs), and Th; In was measured as an internal standard. We used a dilution series of a synthetic zircon solution to generate calibration curves and both a homogeneous solution of Plešovice zircon (Sláma et al., 2008) and a solution with a known Zr/Hf ratio to assess reproducibility. We measured procedural blanks to monitor for laboratory trace-element contamination. Sets of four unknowns were bracketed by individual measurements of the Zr/Hf solution and Plešovice standard, and a new calibration curve was made for every 20 unknowns. We converted solution measurements to zircon concentrations by assuming that all of the measured trace elements substitute for Zr⁴⁺, such that Σ = Zr + Hf + Sc + Y + Nb + Ta + REE = 497,646 ppm. The results from the analyses are presented in File 3 (see ^{footnote 1}). Procedural blanks show no significant laboratory trace-element contamination for the elements of interest. Repeat analyses of the Zr/Hf ratio in our Zr/Hf standard solution give a weighted-mean value of 50.53 ± 0.29 (2σ, mean square of weighted deviates [MSWD] = 0.30), within ~1% of the reference value of 50. Repeat runs of the Plešovice zircon solution are highly reproducible for REE heavier than Nd (Fig. S2 in File 2) and provide confidence in our reproducibility of natural zircon trace-element concentrations.

We modeled trace-element compositions in zircon using an approach that links thermodynamic modeling of magma composition, zircon and sphene saturation thermometry, and partition coefficients. The aim of this modeling was to explore if and how zircon trace-element compositions would differ for zircon that crystallizes from different parent melts. To do this, we applied Perple_X thermodynamic modeling and free energy minimization software (Connolly, 2009, 2005) to model major-element chemistry and major-phase equilibria of a cooling and differentiating magma using the internally consistent thermodynamic database of Holland and Powell (2011), solid solution models based on the thermodynamic data of Tajčmanová et al. (2009), Green et al. (2016), and Jennings and Holland (2015), and the melt model of Holland et al. (2018). After running major-element calculations with Perple_X, we assessed accessory phase stability using saturation models for sphene (Ayers et al., 2018), zircon (Boehnke et al., 2013), apatite (Harrison and Watson, 1984; Bea et al., 1992), monazite (Montel, 1993), and xenotime (Rusiecka and Baker, 2019), though the latter did not saturate. Finally, we partitioned trace elements among all major and accessory phases using partition coefficients compiled and averaged from the literature (including via the GERM database, https://earthref.org/GERM) in the registered Julia package StatGeochem.jl (doi.org/10.17605/osf.io/tjhmw), which was also used to automate our Perple_X calculations (code available in File 4, see ^{footnote 1}). In order to represent the uncertainty in sphene abundance caused by the possible re-equilibration of Ti from other solid phases back into the melt as a response to sphene saturation, we repeated these accessory mineral saturation and trace-element partitioning calculations with Ti re-equilibration proportions varying from 0% to 50%. We ran four sets of simulations, each with a different parental melt composition, at 175 MPa, which is the melt extraction pressure suggested by Kennedy et al. (2016). We used bulk composition data from Kennedy et al. (2016) for the parent compositions, including: (1) the syenite used by Kennedy et al. (2016) in their fractional crystallization models (6Tsyd); this composition is a compromise between the syenite composition with the greatest cumulate signature and the most pristine (least evidence of mixing or mingling) rhyolite composition, aimed at best approximating the parent melt composition before crystal-liquid separation (Kennedy et al., 2016); (2) a high-silica rhyolite from the LSPT (1Tspmps); (3) a trachyte composition (1Tspupa); and (4) a high-REE trachyte composition (4Tspup). We do not have data for Dy-Tm for many of the starting compositions, and Gd is missing from the rhyolite and high-REE trachyte composition; therefore, our results lack these elements.

Detailed results for each sample can be found in File 3 (see ^{footnote 1}). A total of 183 zircon crystals were analyzed by U-Pb TIMS-TEA, 111 of which resulted in ages with precision better than ~0.15 m.y. (2σ). Larger uncertainties are primarily due to common Pb contamination and/or low radiogenic Pb. We focus on the 111 ages with uncertainties ≤0.15 m.y., because those with larger uncertainties are uninformative for the questions of interest here (see Introduction). Ages range from 23.137 to 23.481 Ma (Fig. 2). A few older outliers (>ca. 24 Ma) exist; we assume these outliers are due to the inclusion of older cores. Eruption (volcanic units) or emplacement (resurgent units) ages were calculated using the stochastic sampling approach of Keller et al. (2018). The results are all substantially older than the original ⁴⁰Ar/³⁹Ar eruption age (22.93 ± 0.04 Ma [2σ], Bove et al., 2001), but several samples overlap with a recalculated ⁴⁰Ar/³⁹Ar age (23.076 ± 0.10 Ma [2σ]) using the MMhb-1 standard value of Renne et al. (1998) (Fig. 2). For samples that do not overlap with this age, it is possible that either we have not sampled the youngest zircon crystals in those samples, the magmas that produced them were not zircon saturated at the time of eruption, or there are issues with the ⁴⁰Ar/³⁹Ar age. We favor the latter explanation given the improvements in Ar geochronology that have occurred since the work of Bove et al. (2001) and because the MMhb standard has been shown to be affected by alteration and intergrowth with clinopyroxene (e.g., Schoene and Bowring, 2006). Furthermore, although there is evidence for late-stage heating of the system, the zircon crystals we dated were all euhedral and showed no textural evidence of late-stage resorption (Fig. S3 in File 2), making it unlikely that undersaturation is the explanation. Nonetheless, the discrepancy is unresolved.

Eruption ages of the SPT samples and emplacement ages of the syenite samples overlap within uncertainty. The emplacement age of the monzonite unit (AGQM) is significantly younger than the other samples, but it overlaps within uncertainty with one trachyte (SPU) and one syenite (1C) sample. Individual samples show dispersion ranging from ~60–150 k.y.

Trace-element contents of zircon grains determined by solution ICPMS can be found in File 3 (see ^{footnote 1}). We include trace-element data for all analyzed grains even if they did not result in a precise age. Outliers were identified using a criterion of 1.5× the interquartile range (IQR) on both the entire data set and within each sample. Those grains that fell outside 1.5× IQR of both tests were considered outliers. We took this approach because the entire dataset shows a robust population of grains—the high-Hf population discussed below—that are considered outliers in intrasample tests. This is because relatively few of these grains were analyzed in each sample, which skews the intrasample median to lower Hf values. Given that all samples contain such grains, even if they are a relatively small subset of each sample, and we find no evidence of analytical issues to explain them, we do not think removing these analyses is justified. This is further supported by the fact that high-Hf grains are not identified as outside 1.5 × IQR in the SPM sample, for which an especially large number of these grains was analyzed.

Trace-element contents of zircon grains from the SPT and the syenite units overlap and generally follow trends expected with fractionation of other minerals (e.g., feldspar—Eu/Eu*; sphene, apatite—Dy/Yb) with increasing Hf (Fig. 3). Considering Hf alone (Fig. 4), two populations of grains stand out—we deem these the “high-Hf” and “low-Hf” grains. We use a cut-off of 11,200 ppm Hf to distinguish these populations. The gap between these groups is particularly noticeable in results from the syenite sample, but all samples have a high-Hf population (Figs. 3 and 4). When results for individual trace elements, rather than ratios, are considered, the high-Hf grains fall off trends formed by other grains (Fig. 5). The Hf contents of zircon grains from the high-silica rhyolite do not extend to the lowest values seen in other samples.

The trace-element contents of zircon grains from the monzonite sample are distinct, in some cases plotting with the other samples (e.g., Dy/Yb, Fig. 3) and in others not (e.g., Eu/Eu*, Figs. 3 and 5). The Hf contents of the majority of monzonite grains fall within the Hf gap noted above, but there are low- and high-Hf grains whose trace-element compositions plot with the other samples.

While trends appear in bivariate plots of geochemical indices, there are no trends between geochemistry and age for most samples; only the monzonite sample shows any compelling geochemical relationships with age (e.g., Figs. 3B, 3D, and 6B). A significant proportion of dated grains have low Pb* contents (69 grains with Pb* ≤ 5 pg), and there are no relationships evident between Pb* and other trace elements (e.g., Fig. S4 of File 2, see ^{footnote 1}).

A comparison of results from our trace-element modeling and TEA data can be found in Figures 7 and 8 and in File 2. The goal of this modeling is to further assess if our zircon compositions reflect crystallization from the same or different parent melts. In each set of simulations, we consider the range of zircon compositions from the first zircon to crystallize from that parent melt to that crystallizing at a melt SiO₂ content of 76 wt%. The latter is the SiO₂ content of the high-silica rhyolite (LSPT) bulk composition, which we take to represent the high-silica rhyolite melt; however, given that portions of the LSPT have high crystal contents (40%) and glomerocrysts that may have been entrained from a mush (Lubbers et al., 2020), this could be an underestimate of the melt SiO₂. It is important to note that the middle REEs are strongly affected by sphene crystallization and are therefore tied to the sphene solubility model we use and the amount of Ti that is available for sphene saturation. Given this sensitivity, we focus primarily on the heavy REEs.

Overall, simulation results overlap with the TEA data (Figs. 7 and 8; File 2, see ^{footnote 1}). A notable exception is that the heavy REE contents of zircon from the high-silica rhyolite parent range to significantly higher values than the TEA data and encompass a smaller number of the TEA analyses than the results from other parent melt compositions. Those from the trachyte parent (Fig. S5 in File 2) also range to much more depleted middle REE compositions than the TEA, and those from the high-silica rhyolite parent (Fig. 8) show much deeper Eu anomalies than most of the TEA data. Again, however, the middle REEs are strongly affected by our sphene model.

Previous work on the Lake City magmatic center resulted in a model for the generation of the SPT (Hon, 1987; Hon and Lipman, 1989; Kennedy et al., 2016; Lubbers et al., 2020) that, at its core, follows the crystal mush model (Bachmann and Bergantz, 2004; Hildreth, 2004). The crystal mush model implies that the generation and eruption of high-silica rhyolites is coupled to the formation of silicic or intermediate cumulates. At the Lake City magmatic center, field relations, bulk-rock geochemistry, phenocryst contents, and mineral chemistry have been applied previously to suggest the SPT was a vertically stratified system in which (1) the rhyolites were extracted from a crystal mush, the residuum of which is represented by the syenite; (2) the trachyte is a partially remobilized and erupted portion of this mush; and (3) the post-caldera monzonite represents a distinct magma that existed concurrently with the magma that produced the SPT (Hon, 1987; Kennedy et al., 2016; Lubbers et al., 2020).

Our zircon ages are consistent with this model. Eruption and emplacement ages of the extrusive and intrusive samples overlap, and intra-sample age dispersions are similar, demonstrating that each unit crystallized coevally. It is possible that we are underestimating the duration over which each of the melts existed, because we do not have constraints on how much time is missing from the zircon record. However, the overlap in ages between samples indicates that the magmas that produced them coexisted for at least 100 k.y. Previous studies proposed that the sequence from the SPT eruption to the end of resurgence lasted ~80–300 k.y. (Reynolds et al., 1986; Bove et al., 2001). Our ages shift and shorten this range, suggesting only 60–220 k.y. passed between the eruption of the high-silica rhyolite and the emplacement of the monzonite (Fig. 1).

The age data require that the intrusive units represented by our samples were emplaced and crystallized in <200 k.y. This differs from deep plutons that were emplaced incrementally over millions of years (e.g., Mahan et al., 2003; Coleman et al., 2004; Matzel et al., 2006; Schaltegger et al., 2009) but is similar to estimates for some shallow intrusive systems (Michel et al., 2008; Schoene et al., 2012; Barboni et al., 2015). The hypabyssal nature of the Lake City intrusive units indicates that they are likely more representative of the latter systems; however, only limited portions of the intrusive units are exposed, so it is also possible that we are not accessing the full extent of these bodies. This is potentially supported by geophysical measurements (Grauch, 1985, 1987) that may indicate the monzonite is a small exposure of a much larger intrusive body (Hon, 1987). The intrusions also cut the volcanic units but overlap with them in age, suggesting the SPT eruption mobilized mushy parts of the system that are now represented by these intrusive units. Nonetheless, such observations do not inherently indicate an age relationship between exposed and unexposed portions of intrusions (e.g., Gaynor et al., 2018).

Zircon TEA results also generally support the existing model because trace-element concentrations in zircon grains from the SPT and syenite units overlap considerably. In particular, zircon compositions from the syenite, trachyte, and low-silica rhyolite are essentially identical. Zircon grains from the high-silica rhyolite also plot with results from these units, though they fall into a more restricted range (Figs. 2–5). We discuss this further in the Longevity of Eruptible High-Silica Rhyolite Magmas section.

Finally, results from the syenite- and trachyte-parent trace-element simulations encompass those from TEA, while those from our rhyolite-parent simulations are less consistent with the TEA data (Figs. 7 and 8). These results further suggest that the grains in the volcanic and syenite units grew from melts of evolving composition, rather than from one of these compositional end members.

Zircon grains from the monzonite sample stand out. They overlap in age with the other samples and show similar, if greater, dispersion (Fig. 2); however, their trace-element characteristics are distinct. In some cases, they plot with zircon grains from other units, but in others they fall far off trend or sit in the Hf gap that appears in other samples (e.g., Eu/Eu*, Dy/Yb, Hf; Figs. 3A, 3C, and 4). These results support the conclusions of past work that the magma that produced the monzonite unit coexisted with the magmas that produced the other units but was a physically distinct batch.

Previous studies have suggested that the Lake City magmatic system experienced a vibrant pre- and syn-eruptive history involving multiple magma batches and mixing (Hon, 1987; Kennedy et al., 2016; Lubbers et al., 2020), and the results of our TIMS-TEA analyses provide additional insights into the dynamic nature of the system. First, there is no relationship between zircon age and trace-element composition in most of the samples. This is notable given that there are trends in geochemistry that may be interpreted as evolutionary trends of the system (e.g., lower Eu/Eu* with higher Hf). Insufficient age precision would make it difficult to resolve differences between samples and could be the root of this result. Alternatively, or concurrently, the existence of trends in geochemistry but not age may reflect the complexity of the magma reservoir or its cooling and mixing history. For example, there is ample evidence that the geochemistry of the SPT deposit is related to pre- and/or syn-eruptive mixing of several discrete magma batches (Kennedy et al., 2016; Lubbers et al., 2020). If each of these batches became saturated in minerals that exert control on magma trace-element contents at different points in time, the result upon mixing would be geochemical trends in zircon that are not correlated with time.

To illustrate this, we can first consider the most simplified case of a single homogeneous, scale-independent magma reservoir that is undergoing radial cooling and zircon saturation (Fig. 9). In this scenario, every location in the magma body experiences the same geochemical evolutionary history but will be at a different point along that path at any given time. Even at the same location in the body, different grains will begin crystallizing at different points in time, and existing grains will continue to grow; thus, these grains will be displaced in time and geochemistry, but they will fall along the same geochemical path. If zircon crystals from different regions within this reservoir are then amalgamated and erupted together, trends may be evident in their geochemistry, but not in age.

A similar or even more complex result could be realized if we make this example more realistic, with greater complexity in cooling history, composition, crystallizing mineral assemblage, and mixing, as suggested for the SPT. In turn, the existence of distinct zircon compositions (e.g., Figs. 3A and 6) and trends in zircon geochemistry with time in the monzonite (e.g., Figs. 3B, 3D, 3F, and 6), albeit slight, is consistent with the hypothesis that it represents a separate, isolated batch that did not mix extensively with other batches in the system (Hon, 1987; Kennedy et al., 2016). The existence of high-Hf grains in this sample does, however, suggest there was at least some interaction of this magma with the rest of the system.

A second aspect of our results worth considering is the somewhat perplexing high-Hf zircon population that occurs in all samples (Figs. 4–6). Given that our analyses involve whole grains (i.e., they are weighted averages), the high-Hf signature implies that either some parts of these grains must have been significantly enriched in Hf, or they grew predominantly or entirely from a high-Hf melt. However, the ages and textures of these grains are unremarkable relative to the low-Hf grains, and this lack of correlation between Hf, age, and texture, or CL pattern (File 2, see ^{footnote 1}, and Fig. 3) allows us to confidently rule out several potential explanations for the high-Hf population.

First, although these grains exist in all samples, they are not ubiquitously the youngest grains (Fig. 6), and there is ample evidence for late mafic recharge of the system (Kennedy et al., 2016; Lubbers et al., 2020). As such, it is hard to imagine that the high-Hf grains grew from a late, highly evolved melt that existed everywhere in the system, including in the seemingly independent magma batch that produced the monzonite. Similarly, the grains are not uniformly the oldest grains, and we see no commonality in their CL patterns (though we acknowledge the limitations of 2D images in considering textures), suggesting they are not all xenocrysts with high-Hf cores. It is also unlikely that they were all inclusions armored in other phases, because their textural and age characteristics do not allow them to be clearly separated into a class of their own. Finally, the high-Hf grains did not likely grow from the REE- and Hf-enriched “high-REE” trachytic melt identified by Kennedy et al. (2016), because the grains are not enriched in other REEs, and they appear in all samples, not just the trachyte. Mechanical mixing during the eruption may have mixed populations from different magma batches together, but the lack of “high-REE” trachyte pumice clasts in other eruptive units and the geochemistry of the high-Hf zircon grains make it unlikely that they came from the “high-REE” trachytic melt.

Three additional possibilities remain open:

1. The high-Hf grains are mantled with high-Hf rims, and the rims dominate the averaging of the geochemical signature but not the age.

2. The high-Hf grains or zones crystallized just prior to a mixing event, such that the youngest zircon to crystallize was actually less evolved.

3. The high-Hf grains grew from a separate evolved melt that was extracted from the primary trachyte but existed independently from the LSPT.

In considering these possibilities, the REE contents of the high-Hf grains are noteworthy. For example, Dy and Yb contents of low-Hf grains increase with decreasing Hf, but the high-Hf grains fall distinctly off these arrays—relative to the grains with the lowest Hf contents, they have similarly low Dy and moderate-to-high Yb (Fig. 5). This suggests that the high-Hf grains are recording crystallization either from a different magma or after a change in how the REE budget of the magma was controlled.

The accessory mineral population of the SPT includes zircon, sphene, apatite, and chevkinite (Kennedy et al., 2016). As a middle REE, Dy partitions more strongly into sphene than zircon ( = 100s–1000s, = 10s in natural high-silica rhyolites; Colombini et al., 2011; Padilla and Gualda, 2016), while Yb is similarly compatible in both phases (K_d = 100s; Colombini et al., 2011; Padilla and Gualda, 2016). Thus, if the zircon population is small and sphene is absent, zircon crystallization alone will not substantially affect the REE budget of the magma, and zircon compositions will passively record fractionation of the system. After sphene saturation, sphene crystallization will drastically deplete the melt in middle REEs, and zircon compositions will reflect this change.

Given this logic, the high-Hf grains may be recording zircon crystallization in a melt that was also crystallizing sphene, while the low-Hf grains are recording zircon crystallization without sphene. The presence of sphene could also explain the lack of relationship between Hf and U in our results (Fig. S4 of File 2, see ^{footnote 1}). Zircon Hf and U contents are typically positively correlated (i.e., higher Hf is associated with higher U and therefore younger age; e.g., Claiborne et al., 2010a, 2010b; Lidzbarski, 2014), but crystallization of another U-partitioning phase, such as sphene, could cause them to become decoupled.

If sphene crystallization is the explanation for the geochemical signatures we see, the high-Hf grains could be mantled grains (scenario #1); however, the evidence for late recharge of the system still competes with this possibility, as does the fact that all grains, not just high-Hf grains, lack a U-Hf relationship. If instead the high-Hf grains have thin rims of less evolved melts due to late mixing (scenario #2), the geochemical signature would still need to be dominated by high-Hf zones and decoupled from U contents. Sphene crystallization could contribute to this, but so could the influx of a new magma with different U content. Texturally, some of the high-Hf grains have slightly bright rims in CL images; however, this is not a global observation. Thus, such rim-related arguments for the high-Hf signature seem challenging to support, though surface analyses of grains using in situ techniques could provide more insight.

A somewhat simpler explanation could be the existence of a separate batch of rhyolite in the system that has not been previously documented. In this case, the magma would either have been more evolved at the time of zircon saturation and/or sphene would have also saturated. For this explanation to be consistent with the ages of the high-Hf grains, this melt would have had to have been isolated earlier in the history of the system but was then mixed during the eruption.

Ultimately, more work will be needed to answer these open questions, and combining data from in situ methods with TIMS-TEA could help resolve them. This does not, however, impact our general conclusion: the age distribution and trace-element characteristics of the SPT and syenite zircon crystals are consistent with a model in which the rhyolite is an evolved melt separated from a trachyte/syenite mush. Our results support previous models that indicate the Lake City magmatic system was dynamic and complex, and they may point to the existence of yet another evolved melt being involved in the system.

The longevity of large volumes of eruptible (i.e., melt-rich), high-silica rhyolite magmas remains a source of significant debate, with some studies suggesting they have protracted lifetimes in the crust (10⁵–10⁶ yr; e.g., Jollands et al., 2020; Audétat et al., 2021) and others suggesting their existence is fleeting (10¹–10³ yr; e.g., Shamloo and Till, 2019). A large variety of mineral phases and methods, including both absolute and relative dating approaches, have been used to assess these timescales, but each has important limitations on telling time. For example, U-Pb zircon geochronology is limited to timescales ≥10³ yr, while diffusion geospeedometry is limited to ≤10³ yr. Similarly, zircon can saturate much earlier in the crystallization history than phases such as quartz, and it is less easily resorbed; thus, each mineral phase has the potential to record and integrate different aspects of a system's history (e.g., Cooper and Kent, 2014).

Our age results suggest zircon in the high-silica rhyolite magma (LSPT) crystallized over a minimum of 160 k.y. This is consistent with studies suggesting long (≥10⁵ yr) residence times for high-silica rhyolite magmas. As previously discussed, however, the zircon ages and trace-element characteristics in this sample overlap with those from the other SPT and syenite samples, and we interpret this to indicate that they crystallized as part of the same system and prior to extraction and eruption. This suggests that the zircon crystals from the LSPT are not recording merely the lifetime of the eruptible high-silica rhyolite magma, but rather a more extensive history associated with the development of the overarching magmatic system—from the formation of the source mush to the eventual extraction and eruption of the high-silica rhyolite. This interpretation is consistent with conclusions from other zircon geochronological data sets, as well as petrologic and thermal models, for other large-volume, caldera-forming, silicic eruptions (e.g., Bishop Tuff, Gualda et al., 2012; Yellowstone, Stelten et al., 2013; Rivera et al., 2016; Shamloo and Till, 2019; Oruanui, Chamberlain et al., 2014).

Interestingly, despite overlapping in age and composition with the various other units, the range of LSPT zircon compositions is more restricted. This may merely be an artifact of low analysis number (11 crystals). If so, it suggests that zircon grains in the rhyolite wholly or predominantly crystallized prior to extraction of the LSPT from a mush, and/or the time between melt extraction and syenite emplacement was too short for zircon with distinct trace-element characteristics to crystallize in any of the units. The age overlap between the units supports these possibilities, as do our trace-element modeling results, which suggest that zircon grains crystallizing from an independent rhyolitic melt could be compositionally distinct from those crystallizing from the syenitic melt (Figs. 7 and 8). This is also consistent with our hypothesis that the crystallization duration recorded by LSPT zircon reflects a longer system history than just the residence time of the eruptible high-silica rhyolite magma.

Lubbers et al. (2020) noted that Al contents of quartz crystals from the LPST unit also fall into a more constrained range than those from the other units, suggesting the restricted zircon compositions could be a real phenomenon. If so, we have identified three possible explanations:

1. The LSPT magma was separated from the mush prior to substantial chemical mixing with the mafic intrusions (Hon, 1987; Kennedy et al., 2016; Lubbers et al., 2020), such that LSPT zircon did not crystallize from the less evolved, mixed melt. This would be consistent with existing models.

2. The LSPT was extracted from the trachyte and began crystallizing zircon when the melt had evolved substantially. This would also be consistent with existing models.

3. The LSPT was a separate magma batch that was unrelated to the other magmas and may or may not have coexisted with them. This would not be consistent with existing models.

Our data do not particularly support any of these possibilities. If the less evolved grains grew from late mixed melts, we would expect them to be the youngest grains; however, the zircon ages of all units overlap, and the age dispersion of only the most precisely dated, lowest Hf grains (uncertainties of tens of k.y., Fig. 6B) is wide enough to confirm that they did not all crystallize late. Zircon REEs from all samples also overlap with modeled compositions from the syenite- and trachyte-parent simulations (Figs. 7 and 8). Assuming these models are reliable, this also suggests that none of the grains grew entirely from one of the compositional end members or from an independent batch of rhyolite. Finally, it is unlikely that these grains are a mixture of some that crystallized early, from less evolved melts, and late, from a mixed melt. The ages of the lowest Hf grains span the entire range of time represented in the samples (Fig. 6), and comparison of TEA results and our trace-element models suggest the zircon compositions are likely reflecting crystallization from a fractionating melt, rather than from mixed melts (Figs. 7 and 8). It is possible, however, that a combination of effects is obscuring a distinct signal from any single one, given that TEA data amalgamate trace elements from whole grains. For example, older grains may merely be mantled with late-stage growth.

Further work will be needed to better understand the restricted compositional range of the LSPT zircons. As it stands, these zircon compositions are congruent, rather than complementary, to those from the other units. This strongly suggests the LSPT grains crystallized predominantly, if not entirely, in a thermally oscillating system that was in a primarily mush-like state for most of its history before the rhyolitic melt was extracted. As such, we interpret the age dispersion we see as reflecting this longer history of system evolution, not just the final pre-eruptive residence of an eruptible high-silica rhyolite.

1. Zircon U-Pb TIMS-TEA results are consistent with previous models for the Lake City magmatic center. Zircon ages and compositions overlap, suggesting the magmas that produced the rhyolite, trachyte, and syenite rocks exposed in the Lake City caldera coexisted, were genetically related, and may have mixed prior to eruption. The monzonite unit represents a distinct batch that existed in the system concurrently with these magmas and was affected by mixing with some of them.

2. We find trends in zircon geochemistry that do not correlate with age. This may reflect a complex history of zircon and accessory mineral, namely sphene, saturation effects and pre- and syn-eruptive mixing of rhyolite and trachyte magmas. A high-Hf zircon population may be evidence of a separate batch of rhyolite that was extracted from trachyte relatively early on, was saturated in sphene, and was mixed with other magmas in the system either before or during eruption.

3. Our results suggest that mixing occurred prior to or throughout the eruption sequence, not just after the eruption of the LSPT, as suggested by previous studies (Kennedy et al., 2016; Lubbers et al., 2020). This is consistent with textural observations of resorbing quartz and feldspar in the rhyolite but inconsistent with the homogeneous composition of the LSPT.

4. The new ages we present shorten the lifetime of the Lake City magmatic system, suggesting the time between eruption of the LSPT and emplacement of the monzonite at 60–220 k.y. They also suggest the rhyolite crystallized over a period of ~160 k.y. This timescale likely reflects a protracted lifetime of magma generation in this system that extends beyond the final residence of the eruptible high-silica rhyolite magma.

This work was supported by a Harry Hess postdoctoral fellowship at Princeton University to ASP and National Science Foundation grant (EAR 1249821) to CDD. Pete Lipman is thanked for his guidance in the field and Jordan Lubbers for sharing samples and data. Helpful reviews from Sean Gaynor and Laura Waters and editorial handling by Jonathan Miller are also appreciated. Finally, the first author would like to thank Calvin Miller for his insatiable curiosity, his perpetual support, and his absolutely boundless enthusiasm for everything from magmas to baklava.