Timescales for pluton growth, magma-chamber formation and super-eruptions

Generation of silicic magmas leads to emplacement of granite plutons, huge explosive volcanic eruptions and physical and chemical zoning of continental and arc crust1–7. Whereas timescales for silicic magma generation in the deep and middle crust are prolonged8, magma transfer into the upper crust followed by eruption is episodic and can be rapid9–12. Ages of inherited zircons and sanidines from four Miocene ignimbrites in the Central Andes indicate a gap of 4.6 Myr between initiation of pluton emplacement and onset of super-eruptions, with a 1-Myr cyclicity. We show that inherited zircons and sanidine crystals were stored at temperatures <470 °C in these plutons before incorporation in ignimbrite magmas. Our observations can be explained by silicic melt segregation in a middle-crustal hot zone with episodic melt ascent from an unstable layer at the top of the zone with a timescale governed by the rheology of the upper crust. After thermal incubation of growing plutons, large upper-crustal magma chambers can form in a few thousand years or less by dike transport from the hot-zone melt layer. Instability and disruption of earlier plutonic rock occurred in a few decades or less just before or during super-eruptions. Analysis of inherited zircons and sanidines from Miocene ignimbrites in the Central Andes shows that plutons were emplaced for up to 4 million years prior to onset of volcanism and that disruption of plutonic rock occurs a few decades or less just before or during super-eruptions.

Generation of silicic magmas leads to emplacement of granite plutons, huge explosive volcanic eruptions and physical and chemical zoning of continental and arc crust 1-7 . Whereas timescales for silicic magma generation in the deep and middle crust are prolonged 8 , magma transfer into the upper crust followed by eruption is episodic and can be rapid 9-12 . Ages of inherited zircons and sanidines from four Miocene ignimbrites in the Central Andes indicate a gap of 4.6 Myr between initiation of pluton emplacement and onset of super-eruptions, with a 1-Myr cyclicity. We show that inherited zircons and sanidine crystals were stored at temperatures <470 °C in these plutons before incorporation in ignimbrite magmas. Our observations can be explained by silicic melt segregation in a middle-crustal hot zone with episodic melt ascent from an unstable layer at the top of the zone with a timescale governed by the rheology of the upper crust. After thermal incubation of growing plutons, large upper-crustal magma chambers can form in a few thousand years or less by dike transport from the hot-zone melt layer. Instability and disruption of earlier plutonic rock occurred in a few decades or less just before or during super-eruptions.
Large-volume silicic ignimbrites are co-genetic with emplacement of large granitoid plutons in the upper crust 1-5 . Ignimbrites provide snapshots in the evolution of silicic magmatic systems. Information about pluton emplacement, magma-chamber formation and magma dynamics that leads to super-eruptions comes from geochronology, petrology and crystal residence time studies 6-12 , complemented by numerical modelling [13][14][15][16][17] . Upper-crustal silicic igneous systems are part of transcrustal magmatic systems in which differentiated (silicic) melts can be generated in the middle and lower crust by reactive flow in mushes created by long-term influx of basalt 8, [18][19][20][21] . Buoyancy instabilities 8, 22 drive silicic magmas to shallow crustal levels, sometimes resulting in episodic explosive volcanism.
To understand the onset of pluton emplacement relative to volcanism, the formation of large eruptible magma bodies and the dynamic processes associated with super-eruptions, we integrate geochronology, crystal diffusion modelling and transcrustal magma transport modelling. We combine new 40 Ar- 39 Ar ages on tens of individual sanidine fragments and laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) ages of zircon crystals from early Miocene rhyolitic ignimbrites in northern Chile. These data (Supplementary Tables 1 and 2) are combined with existing (reinterpreted) 206 Pb-238 U chemical abrasion-isotope dilution-thermal ionization mass spectrometry (CA-ID-TIMS) geochronology of the youngest zircon populations identified by LA-ICP-MS. To interpret these data, we develop an internally consistent conceptual model with estimates of fluxes and timescales to build and maintain an eruptible magma body within a magma reservoir with considerations of magma supply and making space for pluton growth by ductile deformation of the crust. Terminology for magmatic systems is summarized in Methods.

Article
Molinos age is a new 40 Ar- 39 Ar analysis of sanidine. The ages are reported at the 95% credible interval including systematic uncertainties.
The Oxaya ignimbrites are rhyolites with mineral assemblages of plagioclase, quartz, biotite, FeTi oxides, ±sanidine, ±amphibole 27,30 and accessory zircon. The Cardones ignimbrite 30 is crystal-rich (around 30-40%) and contains two pumice varieties, one with high sanidine content but no amphibole and temperatures estimated 30 in the range 770-670 °C, and the other with low or no sanidine, minor amphibole and temperatures estimated in the range 850-750 °C. Barometry, thermometry and rare-earth-element geochemistry indicate that these silicic magmas were generated by differentiation from wet basaltic to andesitic magmas with temperatures of 950-850 °C emplaced in the middle and lower crust at depths of approximately 15 km or more 30 . The silicic magmas were then emplaced at depths of 6.0-8.7 ± 2.0 km (ref. 30 ) before eruption. Isotopic data indicate assimilation of older crust 31 .

Geochronology of sanidines and inherited zircons
Single-fragment 40 Ar- 39 Ar total fusion ages of sanidine crystals and ages determined by LA-ICP-MS spot analyses of inherited zircons show highly dispersed age spectra for each ignimbrite, with single-crystal ages spread over millions of years before eruption (Fig. 2). Although widely documented in zircon [32][33][34] , age heterogeneity is being recognized in erupted sanidine 26,35 .
The sanidine and zircon data indicate a magmatic history extending over 7.7 Myr from the oldest inherited zircons (27. interpret zircon and sanidine ages >22.7 Ma as evidence of pluton growth without associated volcanism. Data for three of the ignimbrites suggest that zircon and sanidine crystallized continuously mostly in the 1-Myr interval before eruption, and sanidines retained Ar because of cold storage. Data for Molinos suggest either a hiatus in zircon crystallization 1 Myr before eruption or that the Molinos magma body did not incorporate zircons from this period. We interpret the inherited zircons and sanidines as the consequence of crystallization of silicic magmas to form granite plutons in the upper crust before the onset of ignimbrite volcanism. The repose intervals between ignimbrite eruptions (Extended Data Fig. 2) are: 0.80 +0.07/−0.07 Myr for Poconchile-Cardones, 1.02 +0.08/−0.07 Myr for Cardones-Molinos and 1.28 +0.07/−0.08 Myr for Molinos-Oxaya. The inherited sanidines are qualitatively consistent with the magmatic history implied by the inherited zircons but show a younger and narrower age spread. The percentage of sanidine crystals in each sample older than the eruption ages are: Poconchile (64%), Cardones (60%), Molinos (38%) and Oxaya (45%).

Interpretation of inherited sanidine data
The 40 Ar-39 Ar ages older than eruption ages may have been modified by mixing or diffusion. Sanidine crystals are commonly zoned 30 , so the Ar isotopic composition could be mixtures of old crystal cores and young (eruption age) parts of the crystal, resulting in an intermediate age.
Between crystallization and eruption, radiogenic 40 Ar production is countered by diffusion, which depends exponentially on temperature 38 . We modelled variation of apparent age in sanidine owing to diffusive loss of 40 Ar as a function of crystal size and temperature to constrain storage temperatures. Diffusive loss of 40 Ar from sanidine will create rim-to-core age gradients in individual crystals, with the oldest ages preserved in the cores. Our analyses are typically of 0.5-mm to 1-mm fragments derived from much larger crystals up to 11 mm wide, so we modelled the ages of crystal centres (Fig. 3).
Model results (Fig. 3a) indicate that ages up to several Myr older than the eruption age can be preserved in the cores of large sanidine crystals if the sanidines are stored at temperatures below about 470 °C year, with the oldest sanidines remaining below 400 °C. Previously, the concept of cold storage based on zircon geochemistry and geochronology inferred temperatures at or below the solidus of rhyolitic magmatic systems (around 700 °C) 12,25 . Cold storage can thus occur at temperatures well below the solidus.
Preservation of old sanidine ages also indicates short magma residence times because diffusion of 40 Ar is fast at magmatic temperatures. We calculated the apparent age of sanidine crystal cores as a function of grain size for temperatures of 700 °C and 770 °C (ref. 30 ), assuming no previous diffusive 40 Ar loss during cold storage. To preserve the range of observed sanidine ages, magma residence times must be in years to centuries (Fig. 3b). For example, for an original grain diameter of 11 mm, ages older than 26 Myr can be preserved in crystal cores for about 50 years at 700 °C or 7-8 years at 770° C. These residence times agree with estimates for incorporation of inherited crystals into a magma body before eruption 8-10,12,24 . Some diffusive loss of 40 Ar during cold storage and subsequent residence in the magma probably explains a narrower range of sanidine ages than zircon ages.
Sanidine crystals may experience diffusive loss during welding and cooling of an ignimbrite. Diffusive losses at a temperature of 600 °C for 650 years at 200 m depth were estimated for the Cardones ignimbrite 39 . Post-eruption losses of 40 Ar, even for crystals that give ages >26 Myr, are minor (Fig. 3c).

Implications for plutons and volcanism
Our conceptual model of magma generation, magma transport, magma-chamber formation and onset of super-eruptions (Fig. 4)   Ar ages of sanidines. The data are ordered by age from left (youngest) to right (oldest). 2-sigma uncertainty (analytical precision) intervals are shown for individual crystal (sanidine and zircon) ages, whereas Bayesian eruption ages are 95% confidence intervals (full external precision).
Article by incremental intrusion of mafic magma 21 . A layer of buoyant silicic melt accumulates at the top of the hot zone. Rayleigh-Taylor instabilities develop 22 that eventually trigger rapid ascent of large volumes of silicic magma into the upper crust, initially forming granite plutons and leading to thermal conditions for magma-chamber formation. The Oxaya magmatic system lasted at least 7.7 Myr, a time comparable with the main plutonic episodes and development of large caldera systems 1-7,17,23 . Zircon ages are robust to thermal disturbance, so the absence of age clusters is consistent with continuous zircon crystallization related to upper-crustal magma emplacement. However, the low precision of LA-ICP-MS cannot preclude episodic magmatism at timescales <0.5 Ma (ref. 25 ). The magmatic history involved an initial stage of about 4.6 Myr of pluton growth in the upper crust, followed by roughly 3.1 Myr of episodic ignimbrite eruptions.
We interpret the first stage as the incubation period predicted from thermal models of incremental pluton growth and supported by geological and petrological evidence [13][14][15]17,23 . The incubation period is defined as the time for the upper-crustal magma reservoir to generate a region with a melt fraction exceeding 0.4 (ref. 14 ). During this stage, silicic magma is transferred from the middle crust (depths ≥ 15 km) to the upper crust (depths of about 5 to 8 km) 30 to form granite plutons in which older zircons and sanidines crystallized. Using equation (7) from Annen et al. 15 and a 1D approximation, we estimate an upper limit for magma reservoir growth of 3.5 mm year −1 over 4.6 Myr. We foresee pluton growth by displacement of hot ductile crust downwards and sideways 17 (Fig. 4). In incremental growth models, each increment of magma cools and crystallizes quickly while most of the growing pluton remains at much lower temperature. Models 13,40 of pluton growth indicate that temperatures below 450 °C can be sustained across the pluton for growth rates <1 cm year −1 , explaining preservation of old sanidine ages. The incubation stage involved only pluton emplacement, a conclusion supported by the lack of silicic volcanic clasts and detrital zircons with early Miocene ages in underlying Azapa sediments 37 .
Magma chambers supplied by magma ascent from the hot zone can assemble in the upper crust at the end of the incubation period (Fig. 4a). We attribute growth of large chambers as the consequence of highly episodic growth 40,41 . We interpret inherited sanidine and zircon as being incorporated from previously emplaced granite plutons at a late stage of chamber growth. These plutonic rocks, at temperatures well below the solidus, are disrupted and mixed into silicic magma chambers just decades before and during a super-eruption (Fig. 4b).
We apply a model 22 of episodic magma ascent owing to laterally confined Rayleigh-Taylor instabilities in slowly developing melt layers at the top of the middle-crustal hot zone beneath a subsolidus upper crust with a viscosity μ c (see Methods). On the basis of this model, we propose that repose periods between the main ignimbrite-forming eruptions are dominantly controlled by upper-crust rheology. Of interest are the conditions required for the timescale for the instability to grow into a diapir (τ d ) to match the approximately 1-Myr intervals between the main Oxaya ignimbrites (equation (3) in Methods). For a melt layer width (diameter) of 40 km and a density contrast between the buoyant melt and the overlying crust of Δρ = 300 kg m −3 , we calculate that τ d = 1 Myr requires a crustal viscosity of μ c = 5 × 10 19 Pa s, which is consistent with experimental data 42 at 500-700 °C. Large repose intervals that allow generation of magma volumes for super-eruptions can be explained by a strong crust and silicic magma accumulation in the middle crust.
We compared diking and diapirism 43,44 as mechanisms of magma transport (see Methods) and conclude that dike transport enables rapid formation of upper-crustal magma chambers implied by the geochronological and petrological observations. Conditions for dike formation are inferred to develop at the top of the incipient diapir resulting from Rayleigh-Taylor instability growth over a time comparable with τ d . We develop and justify (see Methods) a simple exchange flow model 45 of crust with viscosity μ c slowly subsiding over a large area and magma with viscosity μ m ascending through a narrow cylindrical conduit as an approximation of dike transport. Calculations of conduit radii, fluxes, magma chamber volume and assembly time are presented in Table 1 for Δρ = 300 kg m −3 , μ m = 10 5 Pa s and μ c = 10 19 Pa s, considering transfer of magma through the conduit from a layer initially 1 km thick of radius R. Three radius values for the collapsing cylindrical crust region are chosen at 15, 20 and 25 km to span the scale of large plutons and super-eruption caldera systems. Fluxes are large, explaining how  40 Ar- 39 Ar ages in the cores of sanidine crystals. We show model ages for the cores of crystals, in which the core is defined as the volume at the centre of the original crystal with a diameter 50% of the original crystal. We show 40 Ar- 39 Ar ages of crystal cores for two reasons. First, we know that the 40 Ar- 39 Ar measurements were made on sanidine crystal fragments and that the grains were fragmented during mineral separation. Second, diffusive loss of 40 Ar will cause rim-to-core age gradients; therefore it is the crystal cores that will provide the oldest ages in our observed 40  large-volume magma chambers can be assembled over periods of a few thousand years. Thermal models of episodic pluton growth and magma-chamber formation [13][14][15]17,40,41 are consistent with the exchange flow model. Episodic magma ascent at rates much higher than time-averaged rates are needed to form large upper-crustal magma chambers 40 and are illustrated by parametric models for the Jurassic Yerington batholith in Nevada 41 , which includes the Luhr Hill pluton with similar geochemical and mineralogical characteristics to the Oxaya ignimbrites. Intrusion rates greater than about 10 cm year −1 were required to form a large shallow magma chamber. For a cylindrical magma chamber with radius 20 km, this translates into 0.12 km 3 year −1 or more, comparable with the calculated exchange flow fluxes (Table 1) and about two orders of magnitude greater than the 0.001 km 3 year −1 needed to generate about 1,000 km 3 of magma in 1 Myr. Here we explicitly invoke slow extraction of silicic melt from mush in the middle crust. Episodic instability of an accumulating silicic melt layer beneath a high-viscosity crust leads to rapid magma transfer into the upper crust (Fig. 4).
This study combining geochronology on both zircon and sanidine, crystal diffusion modelling and magma transport modelling highlights the potential of a multifaceted approach to understanding the evolution and dynamics of large magmatic systems. Previous studies of zircon age distributions led to the concept of long periods of cold storage and rapid mobilization of silicic magma into eruptible magma bodies 24 . However, in our case, study preservation of old sanidines rules out cold storage within a non-eruptible mush system at or above the solidus 25 . Rather, the observations indicate catastrophic assimilation and mixing of subsolidus plutonic rock, which we propose occurred just before and during eruption (Fig. 4b).
We applied a Rayleigh-Taylor instability model to silicic melt accumulations in the middle crust to explain episodic ignimbrite volcanism in which observed repose times on the order of 1 Myr are controlled by upper-crustal rheology. A new model of exchange flow involving magma ascent along dikes and crustal subsidence can explain rapid upper-magma-chamber assembly sourced from middle-crustal mush zones, leading to super-eruptions.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-022-04921-9.   Results are shown for conduit radius, magma flux, magma chamber volume and magma chamber assembly times for a melt layer 1 km thick with crustal viscosity μ c = 10 19 Pa s. Calculations for three different radius values of the magma system (see Fig. 4) are presented.

Methods
Magma system terminology A transcrustal igneous system spans the mantle to the surface and includes magma chambers, igneous mush and fully solidified cognate igneous rocks below the solidus, as well as host rocks incorporated into the system by intrusive mechanisms. Mush is defined as a mixture of melt and crystals forming an interconnected framework in any proportions 8 , whereas magma is sufficiently melt-rich that any crystals are suspended. Transition between mush and magma occurs over crystal contents of typically 50-60%. A magma reservoir is that part of the system containing melt. Mechanisms of forming large bodies of magma (that is, magma chambers) 8 include: incremental intrusion at flux rates sufficient to sustain high enough temperatures for magma formation 14 , segregation of evolved melts at the top of or within mushes 19,21,46 , amalgamation of smaller-scale melt-rich layers within mushes to form larger-volume magma bodies 24 , reheating of mush or fully solidified igneous rocks 47 and fluxing of higher-temperature fluids transferred from hotter mafic magmas 48,49 .
Super-eruptions are defined as eruptions of magnitude 8 or greater 50 .
Geochronology 40 Ar- 39 Ar ages were obtained from sanidine phenocrysts. Pyroclasts were crushed, sieved and washed repeatedly in deionized water, before magnetically separated to isolate feldspar phenocrysts. The largest inclusion-free crystals were selected. Sanidine compositions within a single pumice are very homogenous except for some Ba zonation 30 . Sanidine phenocrysts were leached in an ultrasonic bath in 5% HF for 5 min to remove adhering groundmass glass, before being rinsed three times in deionized water in an ultrasonic bath. Dried sanidines were passed through a magnetic separator at low speed and low angle of tilt, to remove crystals with mineral or melt inclusions. Samples were then hand-picked under a binocular microscope to eliminate those with inclusions and any visibly altered crystals. Sanidine crystals were collected from individual pumice clasts and fiamme from the Cardones and Molinos ignimbrites to avoid contamination with accidental crystals picked up during eruption. However, for the fine-grained Poconchile and Oxaya ignimbrites, we isolated sanidines from bulk rock ignimbrite samples with very low lithic contents to minimize contamination. Many sanidines from the Oxaya Formation are rich in Ba or have Ba-rich growth zones, so we are confident that the sanidines are cognate with the magmatic system, an interpretation confirmed a posteriori. Our samples are typically around 2 mm in dimension but may be fragments of larger crystals, as phenocrysts up to 11 mm are common 27,30 . Laser fusion 40 Ar- 39 Ar geochronology was carried out on 50 individual crystals from each ignimbrite. Sample details and results are given in Extended Data Table 1. There is no link between age and composition with respect to K, Ca and Cl.
Pristine crystals were parcelled into Cu packets, or Al discs, stacked in glass vials and sealed in a large glass vial for irradiation. International standard Fish Canyon sanidine (FCs; with an age of 28.294 ± 0.0036 Myr) were used as fluence monitors for J-determination and packaged throughout the stack at known spacing (geometry) between samples. Samples and standards were irradiated at the Cd-lined (CLICIT) facility of the Oregon State University TRIGA reactor for 15.02 h. Single irradiated crystals (n = 30 per sample) were fused with a CO 2 laser and isotope data were collected using a MAP 215-50 noble gas mass spectrometer 51 .
Samples were analysed in a single batch; backgrounds and mass discrimination measurements (through automated analysis of several air pipettes) specific to each batch were used to correct the data. Air pipettes were run (on average) after every four analyses. Backgrounds subtracted from ion beam measurements were arithmetic averages and standard deviations. Mass discrimination was based on a power law relationship 52 using the isotopic composition of atmospheric Ar (ref. 53 ) that has been independently confirmed 54 . Corrections for radioactive decay of 39 Ar and 37 Ar used published decay constants 55,56 . Ingrowth of 36 Ar from decay of 36 Cl was corrected using the 36 Cl/ 38 Cl production ratio and methods of Renne et al. 57 and was negligible.
The samples were analysed by total fusion and step heating with a CO 2 laser. The mass spectrometer is equipped with a Nier-type ion source. The MAP 215-50 data were collected using an analogue electron multiplier detector. Mass spectrometry used peak hopping by magnetic field switching for ten cycles.
For age comparisons, contributions from sources of systematic uncertainty (that is, uncertainties in 40 Ar/ 40 K of the standard and 40 K decay constants) are neglected and only analytical uncertainties (referred to as 'analytical precision') in isotope measurements of samples and standards are included 62 . In this study, analytical uncertainties include contributions from uncertainties in the interference corrections that have variable effects owing to slight variations in sample composition.
Information on samples and geochronological data are summarized in Extended Data Table 1. The individual geochronological data of zircons and sanidines are provided in Supplementary Tables 1 and 2.
Zircon ages of three of the Oxaya ignimbrites have been previously presented 27 . The previous study determined eruption ages by single-crystal zircon U-Pb CA-ID-TIMS 206 Pb-238 U analyses of high-aspect-ratio zircons lacking any complex crystal shapes and evidence of older cores. In this study, zircons with complex resorbed cores that were previously excluded from the U-Pb CA-ID-TIMS 206 Pb-238 U analyses were analysed at the Geochronology and Tracers Facility, British Geological Survey using a Nu Instruments, Nu Plasma HR, multi-collector inductively coupled mass spectrometer. The Nu Plasma HR was operated in static mode, with simultaneous measurement of isotopes on either a Faraday detector or an ETP secondary electron multiplier (see Extended Data Table 2). Zircon age data are summarized in the supplementary material (Supplementary Table 2). Ages determined by LA-ICP-MS from inherited zircons (Fig. 2) show highly dispersed age spectra (mean squared weighted deviation > 7 for samples with more than 50 analyses).
Laser sampling used a New Wave Research 193-nm laser ablation system, incorporating an in-house-designed, low-volume sample cell with an ablation volume of about 3-4 cm 3 , which, when combined with approximately 1 m of tubing to the plasma torch, leads to a signal washout time of about 1 s. Ablation parameters were: 35 μm static spot, run at a repetition rate of 10 Hz, with a fluence of around 2.2 J cm −2 . Samples were ablated for 30 s with a 15-s washout/laser warm-up period between each analyses.
Data were acquired using the time-resolved analysis function of the Nu Plasma HR's software and processed using Iolite, a software package specifically designed to handle the large volumes of data produced by LA-ICP-MS. Iolite corrects data using the 'standard-sample-bracketing' technique, which applies a normalization factor (measured/known) to the data for the 207 Pb-206 Pb and 206 Pb-238 U of a primary zircon reference material (91500; 1,062 ± 0.4 Myr), analysed at regular intervals during each session. Two other zircon reference materials (GJ-1 and Mud Tank, 602 ± 1 Myr and 732 ± 5 Myr, respectively) were analysed during each session, to check for accuracy and precision. Propagated uncertainties were produced by Iolite and reflect the quadratic combination of the internal uncertainty (the reproducibility of the measured ratios) with the external uncertainty (the reproducibility of the bracketing reference material). Components relating to the systematic uncertainty of the method (age uncertainty of primary reference material, decay constant uncertainties and long-term variance of secondary reference material) are quadratically added, post Iolite.
Dispersed 40 Ar- 39 Ar and the previously published CA-ID-TIMS U-Pb age distributions (Fig. 2) preclude calculation of a weighted mean, leading us to adopt a Bayesian approach to eruption age estimation based on the algorithm of Keller et al. 29 . Bayesian eruption age estimation requires a previous estimate of the relative age distribution of crystallization (zircon) or apparent closure (sanidine) ages before eruption, which was estimated by bootstrapping 45 .
Incorporating all available 40 Ar- 39 Ar age distributions that feature well-resolved pre-eruptive heterogeneity, bootstrapping by kernel density estimation shows (Extended Data Fig. 1) a consistent, exponential form of the relative closure age distribution. This exponential form suggests an underlying survivorship process (for example, potentially consistent with geologic processes ranging from partial degassing of sanidine antecrysts to pre-eruptive Ar accumulated in a cold-storage regime). For excess Ar, the observed continuum of ages would not be expected. Using this bootstrapped age distribution, the resulting eruption age estimates based on 40 Ar- 39 Ar sanidine ages for the Cardones are indistinguishable within uncertainty from those based on U-Pb CA-ID-TIMS zircon crystallization ages, whereas the 40 Ar-39 Ar sanidine ages for Poconchile and Oxaya ignimbrites are just beyond uncertainty of each other. To account for this, we calculate an integrated 40 Ar-39 Ar and U-Pb age that accounts for the scatter.
We recalculated the previous CA-ID-TIMS zircon ages 27 using the Bayesian method, Incorporating constraints from both sanidine and zircon eruption age estimates, we estimated repose intervals between eruptions (Extended Data Fig. 2) using the superposition algorithm of Keller et al. 29 .
We first estimated empirically the form of the relative closure distribution, analogous to the relative crystallization distribution of Keller et al. 45 using a method equivalent to the 'bootstrapping' approach 63 . The results (Extended Data Fig. 1) showed a characteristic form of the closure distribution featuring a nearly exponential decrease in probability density with increasing time before eruption. The consistency and reproducibility of this form, to the first order, between all available well-resolved single-crystal volcanic sanidine 40 Ar-39 Ar age distributions (both from these Andean ignimbrites and the Mesa Falls Tuff 64 ) suggest that this exponential form may be underlain by a consistent physical process. A survivorship process wherein, for example, each sanidine has some finite probability of being reset by reheating in any given pre-eruptive time interval provides one simple mechanism for producing such a trend.
We applied the Markov chain Monte Carlo eruption age estimation algorithm in the Chron.jl software package 29 to each ignimbrite, using a half-normal relative crystallization distribution for all CA-ID-TIMS zircon ages, and our previously determined exponential relative closure distribution for all sanidine 40 Ar- 39 Ar ages (Extended Data Fig. 1). Systematic uncertainties were propagated using the 'optimization intercalibration' the constants of Renne et al. 60 for 40 Ar- 39 Ar ages, and the decay constants of Jaffey et al. 64 , along with the effective systematic uncertainty of the EARTHTIME Tracer 65,66 for U-Pb TIMS ages. Finally, to estimate repose intervals between each ignimbrite (Extended Data Fig. 2), we used Chron.jl to run a second 'stratigraphic' Markov chain Monte Carlo model, combining both the new eruption age estimates and the relative age constraints provided by the stratigraphy. Supplementary Table 3 shows all model outputs.

Diffusion modelling
Argon diffusion calculations were carried out using analytical solutions for simultaneous production and diffusion 67,68 . These solutions, which involve two infinite series, typically converge with fewer than 20 partial sums. We use measured argon diffusion kinetics for FCs 69 . We assume that all sanidine crystals form at 26.5 Ma, reside at a constant temperature until 21.8 Ma (the approximate eruption age of Cardones) and experience no 40 Ar diffusive loss after eruption. To estimate magma residence times, we assume that no argon diffusion occurs during cold storage (that is, that 40 Ar concentration profiles in sanidine crystals were uniform at the beginning of magma residence). Because some previous diffusive rounding of the 40 Ar concentration profiles probably occurred during cold storage, our estimates of magma residence times should be considered minima.
Uncertainties in cold storage temperatures owing to uncertainties in argon diffusion kinetics are fairly invariant and range from ±5 to ±6 °C (1σ), with the largest uncertainties corresponding to small grain sizes and low degrees of fractional argon loss (that is, older 40 Ar-39 Ar ages). Because magma residence times range over a few orders of magnitude, absolute uncertainties in magma residence times owing to uncertainties in argon diffusion kinetics also range over a few orders of magnitude. Generally, magma residence time uncertainty increases with increasing and degree of fractional argon loss (that is, younger 40 Ar-39 Ar ages), increasing grain size and decreasing magma residence temperature. Relative uncertainties in magma residence times, on the other hand, are essentially invariant with grain size or degree of fractional loss, and are about 7% for residence temperatures of 700 °C and about 3% for residence temperatures of 770 °C. For example, for an 11-mm-diameter sanidine grain that experienced 96% fractional loss (that is, has an 40 Ar-39 Ar age of 22 Ma), the residence time at 700 °C is 356 ± 25 years, whereas the residence time at 770 °C is 57 ± 2 years.

Magma transport modelling
We apply here an experimentally verified model for development of buoyancy-induced instability for a growing melt layer beneath a layer of much greater viscosity 22 . Here a silicic melt layer extracted from an underlying mush accumulates beneath the hot upper crust (Fig. 4). The fastest-growing wavelength of the Rayleigh-Taylor instability, λ, in the case of an unconfined layer is given by: in which ḣ is the growth rate of the unstable layer, g is gravity, Δρ is the density difference between melt and overlying crust, μ m is the silicic melt viscosity and μ c is the viscosity of the overlying hot just subsolidus crust. Representative values are Δρ = 300 kg m −3 , μ m = 10 5 Pa s and μ c = 10 19 Pa s (ref. 44 ). We take ḣ values of 1 and 5 mm year −1 based on models of reactive flow related to basalt underplating 21 , resulting in λ values of 600 and 1,400 km. Although approximate, these calculations show that the fastest-growing wavelength is much larger than the width of zones of magma generation beneath batholiths, taken here to be typically in the range 30-50 km. Thus we have applied the theory for confined instability growth 22 for μ m ≪ μ c to calculate a characteristic timescale, τ, for instability: c in which D is the width (diameter) of the layer. Observations of 16 experiments described in Seropian et al. 22 showed that the time, τ d , it takes an instability to transform into a detached diapir was about four times greater than τ, leading to: χ d which we applied in the main text.
One possibility is that the Rayleigh-Taylor instability grows to form a diapir that traverses the intervening plutonic crust. However, magma transport by diapirism is too slow to explain the rapid assembly of magma chambers before ignimbrite eruptions: we estimate using equation (8) in Burov et al. 44 that a 1,000-km 3 diapir with D p = 300 kg m −3 takes about 10 5 years to rise 10 km in a crust with an effective viscosity of 10 19 Pa s. This simplified calculation does not consider heat loss from the diapir, which locally reduces the viscosity of the surroundings 43 , which-in turn-enables faster ascent and could help assimilate older plutonic material into the diapir as it ascends. This mechanism could explain the range of zircon ages but is not consistent with abundant old sanidines. Thus dike transport provides an attractive mechanism to enable fast magma chamber assembly. In our conceptual model (Fig. 4), a conduit (dike or cylinder) is formed that allows an exchange flow 45 between the middle-crustal melt layer and an upper-crustal region in which a magma chamber forms. Here we predict that upward flow of magma along the conduit is balanced by the downward subsidence of the crust. We are interested in the case in which the cross-sectional area of the magma conduit is much smaller than the area of crust flowing downward (A m ≪ A c ) and the magma is much less viscous than the crust (μ m ≪ μ c ). In this scenario, the average speed of the crust downward (U c ) is lower than the average magma flow speed up the conduit, U c ≈ (A m /A c )U m , by many orders of magnitude. To develop a very simple model, we represent the subsiding crust as a large cylinder of radius R and the conduit as a small cylinder of radius r; note that for a dike with a length 1,000 times its width, its width is approximately a quarter of the radius r of a cylindrical conduit that would accommodate the same flux. Owing to the low crust velocity and the large viscosity contrast, the upward flow of magma is well approximated by flow through a cylinder with solid walls (Poiseuille flow): We have applied these equations to make the calculations (Table 1). The difference between an exchange flow along a cylinder and a dike is a matter of geometry, with viscous friction being a factor of a few greater in a dike with the same cross-sectional area as a cylinder. The length of the dike is another factor in governing friction and different choices could be made, but would have a minor effect on calculated magma fluxes. Thus the essential elements of exchange flow are captured by a cylindrical conduit. Even for a cylindrical geometry, our approximate calculations are intended only to illustrate the feasibility of the crust subsiding slowly over a large area, allowing an exchange flow with relatively fast ascent of magma from the mid to the upper crust.

Data availability
All data supporting the findings of this study are available within the paper and its Supplementary Information files. All isotopic and related geochemical data were placed in EarthChem: https://earthchem.org; https://doi.org/10.26022/IEDA/112268.

Code availability
Spreadsheets for carrying out the argon diffusion calculations can be found at: https://github.com/Thermochronology-At-Purdue/ Oxaya2021.