Establishing temporal constraints on major volcanic eruptions is limited by the precision of existing geochronometers. Prior work on the La Garita caldera (Colorado, USA), created by the eruption of the Fish Canyon Tuff, failed to resolve temporal differences between pre-, syn-, and post-collapse eruptive units. Here, we report 40Ar/39Ar geochronologic data supporting a ∼100 ka eruptive history of the La Garita caldera, and resolving the timing of the pre-caldera Pagosa Peak Dacite, syn-caldera Fish Canyon Tuff, and post-caldera dacite of Nutras Creek. Minimizing uncertainty in neutron fluence by rotating samples during irradiation and employing Bayesian statistical interpretation of analytical data enables resolution of the ∼60 ka pre-caldera eruptive history and a hiatus of 0–20 ka prior to the eruption of post-caldera lavas. The improved precision demonstrated using these methods provides previously unresolvable temporal constraints on physical processes in the La Garita magmatic system and underscores the potential of unraveling other closely spaced events in geologic time.


To understand cause-and-effect relationships, geologic processes, eruptive histories, and evolutionary events in Earth history, we must resolve the timing of closely spaced events. Volcanic units, whether they are sampled in proximal locations or as distal tephra units in a sedimentary record, are often used to understand these events, and they are commonly indistinguishable from their nearest neighbors. One of the most studied eruptions on Earth is the Fish Canyon Tuff (Colorado, United States), but it has not previously been resolved from its pre- and post-caldera eruptions. Here, we temporally resolve pre-, syn-, and post-caldera events in the Southern Rocky Mountain volcanic field to track magmatic cycles and understand the physical processes active in magmatic systems.

The Oligocene La Garita caldera and associated Fish Canyon Tuff (ca. 28 Ma) in southern Colorado record one of the largest volcanic eruptions known on Earth. Erupted ignimbrite volumes are estimated at >5000 km3, and the eroded caldera (area >2600 km2) contains pre- and post-caldera deposits of Fish Canyon magma interpreted as nearly contemporaneous eruptive products, bracketing the La Garita structural collapse (Lipman et al., 1997; Lipman, 2006; Bachmann et al., 2007). The eruptive phases of the La Garita system correspond to three lithostratigraphic map units (Fig. 1) of Lipman (2006): (1) Pagosa Peak Dacite (PPD)—pre-collapse “blob-and-ash” and massive lava-like deposits preserved adjacent to the southern caldera margin; (2) Fish Canyon Tuff (FCT)—intracaldera and outflow ignimbrite; and (3) dacite of Nutras Creek (NCD)—massive to flow-layered lava and breccia preserved in a small area on the north flank of the La Garita resurgent block. Estimated eruptive volumes for the pre-, syn- and post-collapse deposits are 200 km3, >5000 km3, and <1 km3, respectively, (Lipman et al., 1997; Bachmann et al., 2000; Lipman, 2006). All three units are petrographically and chemically indistinguishable crystal-rich dacites.

Sanidine from the FCT eruption is a widely used neutron fluence monitor for 40Ar/39Ar geochronology, and the “absolute” age of Fish Canyon sanidine (FCs) from ignimbrite outflow deposits has been determined and revised numerous times (Cebula et al., 1986; Renne et al., 1998, 2010, 2011; Kuiper et al., 2008). Bachmann et al. (2007) determined statistically indistinguishable 40Ar/39Ar ages at the 2σ confidence interval for sanidine from PPD, FCT, and NCD (Fig. 2). Likewise, the >400 ka zircon crystallization history from Wotzlaw et al. (2013) tracks the thermal and chemical evolution of the La Garita magma chamber, but the youngest zircons from FCT and NCD are indistinguishable in age (PPD zircons were not analyzed) (Fig. 2).

In this study, we resolve the timing of pre-, syn- and post-caldera deposits of the La Garita system through 40Ar/39Ar geochronology by exploiting the rotation capabilities of the U.S. Geological Survey nuclear reactor. We utilize Bayesian analysis to measure the probability of age differences between samples with similar eruptive ages. For the first time, we can resolve an eruptive history that begins with the PPD ∼60 ka prior to the eruption of FCT and likely continued for as much as 20 ka longer with the NCD.


While complications have been reported (Renne et al., 2012), sanidine typically records eruption ages of volcanic units, and it is present in most major units of the La Garita system. We performed 40Ar/39Ar geochronology on a total of nine samples from PPD, FCT, and NCD. Samples and neutron fluence monitors were irradiated in the same annular position of an aluminum irradiation disc (Fig. DR1 in the GSA Data Repository1); here, we used sanidine from the Fish Canyon Tuff (FCs-EK; Morgan et al., 2014) as a neutron fluence monitor, assuming an age of 28.201 Ma (Kuiper et al., 2008). Samples and monitors were rotated at 1 rpm throughout an 8 hr irradiation in the U.S. Geological Survey TRIGA reactor (Denver, Colorado) to minimize variations in the efficiency of reactor-produced argon generation. Rotation during irradiation aims to smooth neutron fluence gradients, allowing for the same J value (a function of neutron fluence; see Equations 1 and 2) to be applied to an annular ring of an irradiation disc (e.g. Schwarz and Trieloff, 2007). However, Rutte et al. (2015) showed that rotating the canister does not guarantee equal fluence. Measurable fluence gradients can arise if the sample packet and rotation axis are not aligned parallel to the nuclear fuel elements. Here, ensuring parallel rotation yields a neutron fluence gradient of <0.1% across an 18.5 mm disc, based on F-values (the mean ratio of radiogenic to reactor-produced argon, 40Ar*/39ArK; after Dalrymple et al., 1981) for neutron fluence monitors placed at multiple positions around the disc.

Samples and standards were analyzed on a Thermo Scientific ARGUS VI mass spectrometer by multi-collection on four Faraday cups and one ion counter (for 36Ar). Ages herein were calculated using an age of 28.201 Ma for the FCs neutron fluence monitor (Kuiper et al., 2008), decay constants (λ) of Min et al. (2000), and an atmospheric 40Ar/36Ar value from Lee et al. (2006). Analytical details, summary tables, and an argon isotopic data table can be found in the Data Repository (methods and Tables DR1–DR4) and in Morgan (2019).

Other studies have identified complex age spectra in the FCs monitor (Phillips and Matchan, 2013; Phillips et al., 2017; Jicha et al., 2016) that may indicate issues with pre-eruptive argon retention, recoil of 39Ar, unresolved isobars, and/or isotopic fractionation during incremental laser heating. Similar complexities have also been identified in step-heating experiments on many samples considered here, but are not believed to have affected interpretations. See the Data Repository and Figure DR2 for further discussion and age spectra.


The F-value was assumed to reflect the time of an eruptive event (te), recognizing that past work has argued that this ratio may be influenced by other factors (e.g., Andersen et al., 2017; Phillips and Matchan, 2013; see the Data Repository for discussion): 
where λ is the total decay constant for 40K and J is a measure of the reactor production of 39Ar from 39K through comparison with a neutron fluence monitor: 
where tm and Fm refer to the assumed age and measured F-value of the neutron fluence monitor, respectively.
Given our assumption that neutron fluence gradients are negligible (and thus J is constant) within an annular ring of an irradiation disc (an assumption we evaluate below), the expected difference in F-values between our unknowns and the FCs monitor representing the FCT is related to the ages of these samples as: 
Here, we refer to the age and F-value of unknowns as tu and Fu, respectively. Expanding out J based on Equation 2, and expressing the timing of an unknown event as the sum of the age of our monitor and an offset, tu = tm+ Δt, we obtain: 
where R denotes the ratio of unknown to standard ratios, Fu / Fm. Given constant J, differences in mean age (e.g., non-zero values of Δt) between the units in our sampled stratigraphy and the eruption of FCT will have values of R that are credibly different than 1. We can also rearrange this to directly compute the difference in age: 


If age differences are present, measurement variability from analytical uncertainty can result in overlap in the distributions of F-values between samples and monitors. Here, we are interested in whether there are credible differences between their mean ages, and hence mean F-values (Equation 5). Traditionally, weighted means have been used to represent the mean of a set of analyses, but they are sensitive to the inclusion or exclusion of outliers. Thus, we adopt a Bayesian alternative to the t-test modeled after Kruschke (2013), using a Markov chain Monte Carlo (MCMC) approach (Foreman-Mackey et al., 2013).

Individual crystals in a sample are characterized by an F-value (Fcrystal) and an analytical uncertainty (σcrystal). For each sample (∼30 crystals per sample), we characterize the posterior probability, p(Fmean,ω | Fcrystalcrystal), of the mean, Fmean, and an overdispersion term, ω, that accounts for variance not explained by analytical uncertainties (e.g., due to geologic phenomena). The posterior probability is proportional to the product of the prior probability, p(Fmean,ω), and the likelihood, p(Fcrystalcystal | Fmean,ω) of our observations. We assign a Gaussian prior probability, p(Fmean) = 8.944 ± 0.042 (1σ), defined by the mean and standard deviation of all Fcrystal values measured here to acknowledge previous work that found the geologic units studied here to be indistinguishable in age (Bachmann et al., 2007). We assign an uninformative prior expectation for the overdispersion term, p(ω) = 1, as we have no existing expectations for the amount of overdispersion present. We follow Vermeesch (2018) and assume that the likelihood of any single crystal measurement is given by a normal distribution whose variance is the sum of that from analytical uncertainty and overdispersion: forumla.

We run the MCMC model for 5000 iterations after a burn-in period of 1000 iterations with 200 “walkers” (each walker takes a guided step through the parameter space), resulting in 106 samples that characterize the posterior probability (we refer to the 106 numerical, randomized samples drawn for each geologic sample as MCMC samples). We summarize the credible intervals (CIs) for each sample’s value of Fmean and ω2 with the median and 95% interval of MCMC samples (Table DR1). Histograms of MCMC-sampled values of Fmean, which approximate the posterior probability, for each geologic sample are shown in Figures 3A and 3B.

We evaluate our assumption of constant J throughout the irradiation disk by comparing the posterior probabilities of the FCs monitors (FCs-EK separate; Fig. 3B). Substantial overlap of the CIs of the mean, Fmean, for all FCs monitors requires immeasurable differences in J (Equation 1). While there is also significant overlap in the CIs of samples SRM08 and SRM09 (also samples of the FCT; Fig. 3B), we preserve these as independent samples to document our ability to differentiate age differences of 0. Similarly, we do not group rock samples from the same geologic unit to evaluate the reproducibility of our results and allow for geologic variability. To incorporate observed variability and analytical uncertainty from all fluence monitors, we compute the posterior probability of Fm from a MCMC model where we pooled all single-crystal analyses for the fluence monitors (FCs-C1, FCs-C5, FCs-C9; see Fig. DR1) in evaluating the likelihood (FCs-EK pooled; Fig. 3B).

To determine the posterior probability of R − 1 (Equation 5), we difference every entry from the MCMC model of unknown sample means (Fu, sample SRM01-09) from the MCMC model of monitor means (Fm, FCs-EK pooled; e.g., Kruschke, 2013). The posterior probability of R − 1, specifically the probability that R − 1 values are >0 or <0, reflects our confidence in differences in the mean age of unknowns and the FCT. We also compute a modeled age difference (Equation 5; Table DR1) based on a mean assumed age of Fish Canyon sanidine (ts) and the total decay constant (λ) (Min et al., 2000; Kuiper et al., 2008). While this doesn’t account for uncertainty in tm or λ, reported uncertainties in λ and tm would propagate into Δt on the order of ∼0.1% of Δt (assuming uncertainties in λ and tm are uncorrelated).

For most samples analyzed here, the intersample variability in Fmean values is consistent with the measured analytical uncertainty (Table DR1). The additional variance (e.g., due to geologic dispersion), ω2, is small with respect to the expected variance from analytical uncertainty (on average, the variance from analytical uncertainty is ∼3× ω2). This is consistent with MSWD values close to 1 (Table DR2). However, some samples do show signs of overdispersion. For sample SRM02, the variance from analytical uncertainty and ω2 are comparable, and for samples SRM05 and SRM08 and the pooled FCs-EK, the median values of ω2 are about one-half the variance expected from analytical uncertainty.


The five PPD samples analyzed indicate a protracted, pulsed eruptive history of the PPD. These five samples are separated geographically (Fig. 1). Two samples (SRM03 and SRM04) are from the southwestern area of mapped PPD (here named PPD south, or PPDS), and three samples (SRM05, SRM06, and SRM07) are from the northeastern area of mapped PPD (here named PPD north, or PPDN). The CI bounds on differences in the mean age for these samples range from 13 ka younger to 78 ka older than FCT (Fig. 3). The two PPDS samples are not distinctly older than FCT; the three PPDN samples define an interval from 15 to 78 ka prior to FCT. The PPDS samples are most likely younger than the PPDN samples. Prior studies inferred distinct eruptive sources for these areas, based on maximum thickness as distinct lobes, but were unable to discern relative ages from stratigraphic relations (Bachmann et al., 2000; Lipman, 2006).

The two FCT samples investigated here are from the lowermost and uppermost subunits of the FCT (respectively, units Tfc1 and Tfc6 of Lipman [2006]). The lowermost unit has long been used as a neutron fluence monitor for 40Ar/39Ar geochronology (Renne et al., 1998), but a recent sampling now in use as a standard (Morgan et al., 2014) was sourced from an intermediate subunit (unit Tfc5 of Lipman [2006]) and was shown to be indistinguishable from the lowermost unit in age. The CIs for all FCT samples and neutron fluence monitors overlap (Fig. 3). The eruption of NCD appears to have closely followed the caldera-forming FCT eruption, with both samples indistinguishable from FCT. The most likely age of sample SRM01 is younger than that of sample SRM02, but the range of credible ages for these samples overlaps. Because NCD overlies FCT, we conclude that NCD was concurrent with or postdated FCT by as much as ∼20 ka (Fig. 3).

These new ages are consistent with the mapped geologic field relations of Lipman (2006) and with existing geochronology (Fig. 2). For example, (1) PPD everywhere underlies FCT along the southern caldera margin; (2) only the southern PPD is associated with faults interpreted as syn-collapse features; and (3) the single exposure of NCD is stratigraphically above FCT on the La Garita resurgent block of intracaldera FCT. However, the extended eruption time of the PPD and the migration of PPD eruption activity from north to south was not previously interpreted.

The time scales presented here for pre- and post-caldera formation in the La Garita system are broadly similar to those of relatively young, albeit smaller caldera systems, which display prolonged but episodic activity with periods of dormancy that alternate with periods of eruption (Druitt et al., 2002; Hildreth, 2004). As with the La Garita system, for example, the Long Valley magmatic system (California, USA) had extensive precursor activity for >1 Ma prior to the caldera-forming eruption of the Bishop Tuff at 0.76 Ma (e.g., Andersen et al., 2017; Mark et al., 2017). The last eruption of Glass Mountain rhyolite (ca. 0.79 Ma) erupted ∼30 ka before the Bishop Tuff (Metz and Mahood, 1985; Simon and Reid, 2005; Hildreth, 2004), while in the La Garita system, the eruption of the PPD began ∼60 ka prior to the eruption of FCT.


The combination of advances in sample rotation in the reactor and the application of Bayesian analysis allowed us to temporally discriminate between the PPD, FCT, and NCD, which were previously indistinguishable by either U-Pb in zircon or 40Ar/39Ar analyses of sanidine. For the ca. 28 Ma samples analyzed here, intervals as small as ∼20 ka can be differentiated. Given that the PPD appears to have erupted intermittently over a protracted time span, the role of the PPD in triggering the FCT eruption (Bachmann et al., 2007) should be revisited.

Although discernible time differences such as these would vary in magnitude through Earth’s history and with sample and data quality, the advances presented here may allow for determinations of small Δt values for a range of Earth science questions—including rates of evolutionary change, tectonic uplift, and other volcanic eruptive activity.


We thank Pete Lipman for considering the geologic implications of our data and a helpful review of the manuscript, Brenhin Keller for advice on Bayesian modeling, Daniel Rutte for advice with sample rotation during irradiation, and Brett Davidheiser-Kroll for assistance with fieldwork. We thank Tom Sheldrake and two anonymous reviewers for their constructive and thoughtful reviews. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

1GSA Data Repository item 2019173, details of analytical procedures and incremental heating experiments, Figures DR1–DR2, and Tables DR1–DR4, is available online at http://www.geosociety.org/datarepository/2019/, or on request from editing@geosociety.org.
Gold Open Access: This paper is published under the terms of the CC-BY license.