Thermal budgets of magma storage constrained by diffusion chronometry: the Cerro Galán ignimbrite

The long-term thermochemical conditions at which large bodies of silicic magma are stored in the crust is integral to our understanding of the timing, frequency, and intensity of volcanic eruptions, and provides important context for volcano monitoring data. Despite this realization, however, individual magmatic systems also may have unique time-temperature paths, or thermal histories, that are the result of many complex and, sometimes, simultaneous/competing processes, ultimately leading to an incomplete understanding of their long-term thermal evolution. Of recent interest to the volcanology community is the length of time large volumes of eruptible and geophysically detectable magma exist within the crust prior to their eruption. Here we use a combination of diffusion chronometry, trace element, and thermodynamic modeling to quantify the long-term thermal budget of the 2.08 Ma, 630km3 Cerro Galán Ignimbrite (CGI) in NW Argentina, one of the largest explosive volcanic eruptions in the recent geologic record. We find that diffusion of both Mg and Sr in plagioclase indicate that erupted magmatic material only spent decades to centuries at or above temperatures (~750°C) required to produce and store significant volumes of eruptible magma. Calculated plagioclase equilibrium liquid compositions reveal an array that is controlled overall by fractionation of plagioclase + biotite + sanidine, although high-resolution trace element transects record a diversity of long-term storage conditions with some plagioclase recording periods of co-crystallization with biotite and sanidine, while others do not. Despite these chemical differences in long-term storage, we find diffusion models record a unimodal distribution which, when combined with prior work revealing zircon residence times of ~105 years, and calculated zircon saturation temperatures of 807 ± 8°C, provide compelling evidence that the CGI magmatic system spent most of its upper crustal residence in a largely uneruptible state and was ultimately remobilized shortly before eruption. INTRODUCTION Large silicic magma reservoirs, responsible for producing the biggest explosive volcanic eruptions in the geologic record (Mason et al., 2004; Wilson et al., 2021), exist as local thermal anomalies in the otherwise cold mid to upper continental crust (Huber et al., 2019; Turcotte and Schubert, 2002). Maintaining magma volumes that are sufficiently liquid-rich to erupt reflects a balance between conductive cooling, advective gain or loss of heat via eruption and addition recharge magma, and latent heat of crystallization (e.g., Degruyter and Huber, 2014). Because of this, rapid addition of magma is thought to be necessary to grow and sustain large volumes of magma in the crust (Annen et al., 2015; de Silva and Gregg, 2014; Gelman et al., 2013). Individual magma systems (and different regions with a single magma system) may also have unique time-temperature paths, or thermal histories, that are the amalgamation of many complex and, sometimes, simultaneous processes (e.g., recharge, eruption, second boiling, thermal buffering, magma ascent, etc.). Despite these complexities, quantifying the thermal state and evolution of a magma reservoir remains important as it has been shown to be responsible for controlling eruption timing, frequency, and dynamics (e.g., Degruyter and Huber, 2014) and thus an understanding of the long term thermal history of a crustal magma reservoir is important for understanding volcanic processes. The thermal histories of crustal magma systems, and the conditions that give rise to large bodies of eruptible magma, have been investigated for several decades using numerical modeling (e.g., Annen et al., 2006; de Silva and Gregg, 2014; Gelman et al., 2013; Jellinek and DePaolo, 2003; Karakas et al., 2017) and also more recently these issues have been investigated by petrological and geochemical techniques. The latter are based on a range of approaches, including single mineral geochronology (Andersen et al., 2017; Barboni et al., 2016; Klemetti et al., 2011; Szymanowski et al., 2017), thermobarometry and phase equilibria (e.g., Walker et al., 2013) or a combination of these approaches. In addition, in recent years, diffusion chronometry in a range of different mineral phases (e.g., Bradshaw, 2017; Cooper and Kent, 2014; Rubin et al., 2017; Shamloo and Till, 2019) has also emerged as a powerful way to constrain the thermal evolution of magmas. Although diffusion chronometry is more often used to quantify timescales of short duration magmatic processes – typically those associated with the buildup to eruption (e.g., Couperthwaite et al., 2020; Ruth et al., 2018; Shamloo and Till, 2019), mineral and element pairs with slower diffusivities at temperatures relevant to silicic magmatic systems can also be used to study longer term magmatic processes (e.g. Cooper and Kent, 2014; Bradshaw, 2017; Rubin et al., 2017). Based on these studies, two broad end-member models have emerged. In the first scenario crustal magma reservoirs spend the vast majority of their time in a state of near solidus temperatures, where they are not eruptible, and experience punctuated thermal events that generate and/or erupt eruptible volumes of magma (e.g. Cooper and Kent, 2014; Rubin et al., 2017; Szymanowski et al., 2017). This view is also supported by geophysical approaches, which rarely observe melt-dominated (greater than 50%) magma reservoirs (Lundstrom and Glazner, 2016), although some exceptions to this may exist (Laumonier et al., 2019). In contrast it has also been argued that magma reservoirs spend the vast majority of their time at temperatures that allow for a significant and eruptible melt fraction throughout much of their history to be present (e.g. Barboni et al., 2016; Kaiser et al., 2017; Tierney et al., 2016). These alternatives imply distinctly different behavior within the associated magmatic plumbing systems and are thus important to constrain. However, considering these scenarios as a simple binary is also likely to obscure important details about the complexity of magmatic systems, and larger magmatic reservoirs may exhibit both of these storage types within a single system (e.g., Andersen et al., 2017; Bradshaw, 2017; Mucek et al., 2021). Many thermal history investigations have also focused on systems that erupt relatively small volumes of magma (e.g., less than 101-2 km3; Barboni et al., 2016; Cooper and Kent, 2014; Rubin et al., 2017; Tierney et al., 2016; Till et al., 2015). Larger systems have also been studied (e.g., greater than 500 km3; Andersen et al., 2017; Bradshaw, 2017; Shamloo and Till, 2019; Szymanowski et al., 2017), although many of the available constraints utilize relatively small data sets, and thus are less suitable to see if there is complexity in the thermal histories experienced by different magmatic components within an individual magma reservoir. One way to investigate this issue further is to conduct more detailed studies of a single large eruption, and to combine diffusional constraints with other petrologic approaches to provide additional context within which to interpret results. In addition, use of multiple trace elements for diffusion modeling can also make diffusion constraints more robust (e.g., Chamberlain et al., 2014; Morgan and Blake, 2006; Shamloo and Till, 2019; Till et al., 2015). Here we utilize diffusion modeling of Sr and Mg in plagioclase and detailed documentation of core to rim changes in plagioclase chemistry in the 2.08 Ma, 630 km3 Cerro Galán Ignimbrite (Folkes et al., 2011c; Kay et al., 2011). We also combine these results with thermodynamic modeling (e.g., Rhyolite-MELTS; Gualda et al. 2012) and other petrological observations to constrain the long-term thermal evolution in this system. Geologic Background The Cerro Galán caldera, located in northwest Argentina on the eastern edge of the Puna plateau, is a part of the larger Central Volcanic Zone of the Andes (Figure 1). Between 5.6 and 2.08 Ma the Cerro Galán magmatic system produced >1200km3 of high-K, crystal rich (4050%), chemically homogenous (68-71wt% SiO2) rhyodacite (Folkes et al., 2011a, 2011c) that have been classified as a ‘monotonous intermediates’ (Hildreth, 1981). The largest and most recent eruption, the 2.08 Ma Cerro Galán Ignimbrite (CGI), produced 630km3 of ignimbrite that extends outward from the caldera 40km in all directions and up to 80km north of the caldera (Folkes et al., 2011c). The CGI is a massive, crystal-rich (40-50% crystals), and generally pumice poor ignimbrite with a mineral assemblage of plagioclase + sanidine + quartz + biotite + Fe-Ti oxides + apatite + zircon. The majority of pumice have between 69-71 wt% SiO2, however two types of pumice are evident (Wright et al., 2011): 1) the majority of pumice (~95%) are “white pumice” that contain 44 – 57% crystals, no microlites, > 76 wt.% SiO2 groundmass, and higher Ba concentrations; 2) a volumetrically small proportion of pumice (~5%) are “grey pumice” that contain 35 – 59% crystals, abundant microlites, groundmass SiO2 concentrations of 69 – 74 wt.%, and ~150ppm lower Ba concentrations at equivalent SiO2 concentrations as white pumice. Fe-Ti oxide geothermometry on white pumice record temperatures of 790-820°C (Wright et al., 2011), however Folkes et al., (2011) note that due to the high fO2 (+1-1.7 NNO) calculated for the CGI these temperature estimates may have uncertainties of 50-100 °C. Trace element ratios in the CGI, when compared to previous ignimbrites from the Galán magmatic system, indicate magmas that experience closed system evolution (Folkes et al., 2011b). Appearance of sanidine in the CGI, coupled with the disappearance of amphibole present in previous eruptions from the Cerro Galán system has been experimentally demonstrated to be the result of a shallowing of the magmatic system over time, and this is also supported by volatile contents in quartz-hosted melt inclusions (Grocke et al., 2017). Previous reconnaissance studies (e.g., Folkes et al., 2011b; Wright et al., 2011, Bradshaw, 2017) of the thermal evolution of the Cerro Galán Ignimbrite indicate that there are multiple populations of plagioclase with one population recording magma storage for long durations at high temperature (> 750 °C; type 1), while the other (and more common) records short time at high temperature (type 2). Additionally, these different plagioclase populations contain unique trace element signatures with type 1 plagioclase exhibiting a positive correlation between Ba and An, while type 2 exhibits a negative correlation. This has been interpreted as being the result of plagioclase recording crystallizing environments in which the magma system switches from conditions that inhibit sanidine crystallization (type 2) to those that promote it (type 1) as the magma system shallows with time (Bradshaw, 2017). It is also possible, however, that these trends may be heavily influenced by the presence of biotite (found up to 15 vol% vs. 1.5 vol% sanidine in the CGI) and reflect a magma reservoir with diverse crystallization conditions, further establishing the need for an in depth, multidisciplinary assessment of the time-temperature evolution of the Cerro Galán magmatic system.


INTRODUCTION
Large silicic magma reservoirs, responsible for producing the biggest explosive volcanic eruptions in the geologic record (Mason et al., 2004;Wilson et al., 2021), exist as local thermal anomalies in the otherwise cold mid to upper continental crust (Huber et al., 2019;Turcotte and Schubert, 2002). Maintaining magma volumes that are sufficiently liquid-rich to erupt reflects a balance between conductive cooling, advective gain or loss of heat via eruption and addition recharge magma, and latent heat of crystallization (e.g., Degruyter and Huber, 2014). Because of this, rapid addition of magma is thought to be necessary to grow and sustain large volumes of magma in the crust (Annen et al., 2015;de Silva and Gregg, 2014;Gelman et al., 2013).
Individual magma systems (and different regions with a single magma system) may also have unique time-temperature paths, or thermal histories, that are the amalgamation of many complex and, sometimes, simultaneous processes (e.g., recharge, eruption, second boiling, thermal buffering, magma ascent, etc.). Despite these complexities, quantifying the thermal state and evolution of a magma reservoir remains important as it has been shown to be responsible for controlling eruption timing, frequency, and dynamics (e.g., Degruyter and Huber, 2014) and thus an understanding of the long term thermal history of a crustal magma reservoir is important for understanding volcanic processes.
The thermal histories of crustal magma systems, and the conditions that give rise to large bodies of eruptible magma, have been investigated for several decades using numerical modeling (e.g., Annen et al., 2006;de Silva and Gregg, 2014;Gelman et al., 2013;Jellinek and DePaolo, 2003;Karakas et al., 2017) and also more recently these issues have been investigated by petrological and geochemical techniques. The latter are based on a range of approaches, including single mineral geochronology (Andersen et al., 2017;Barboni et al., 2016;Klemetti et al., 2011;Szymanowski et al., 2017), thermobarometry and phase equilibria (e.g., Walker et al., 2013) or a combination of these approaches. In addition, in recent years, diffusion chronometry in a range of different mineral phases (e.g., Bradshaw, 2017;Cooper and Kent, 2014;Rubin et al., 2017;Shamloo and Till, 2019) has also emerged as a powerful way to constrain the thermal evolution of magmas. Although diffusion chronometry is more often used to quantify timescales of short duration magmatic processes -typically those associated with the buildup to eruption (e.g., Couperthwaite et al., 2020;Ruth et al., 2018;Shamloo and Till, 2019), mineral and element pairs with slower diffusivities at temperatures relevant to silicic magmatic systems can also be used to study longer term magmatic processes (e.g. Cooper and Kent, 2014;Bradshaw, 2017;Rubin et al., 2017).
Based on these studies, two broad end-member models have emerged. In the first scenario crustal magma reservoirs spend the vast majority of their time in a state of near solidus temperatures, where they are not eruptible, and experience punctuated thermal events that generate and/or erupt eruptible volumes of magma (e.g. Cooper and Kent, 2014;Rubin et al., 2017;Szymanowski et al., 2017). This view is also supported by geophysical approaches, which rarely observe melt-dominated (greater than 50%) magma reservoirs (Lundstrom and Glazner, 2016), although some exceptions to this may exist (Laumonier et al., 2019). In contrast it has also been argued that magma reservoirs spend the vast majority of their time at temperatures that allow for a significant and eruptible melt fraction throughout much of their history to be present (e.g. Barboni et al., 2016;Kaiser et al., 2017;Tierney et al., 2016).
These alternatives imply distinctly different behavior within the associated magmatic plumbing systems and are thus important to constrain. However, considering these scenarios as a simple binary is also likely to obscure important details about the complexity of magmatic systems, and larger magmatic reservoirs may exhibit both of these storage types within a single system (e.g., Andersen et al., 2017;Bradshaw, 2017;Mucek et al., 2021). Many thermal history investigations have also focused on systems that erupt relatively small volumes of magma (e.g., less than 10 1-2 km 3 ; Barboni et al., 2016;Cooper and Kent, 2014;Rubin et al., 2017;Tierney et al., 2016;Till et al., 2015). Larger systems have also been studied (e.g., greater than 500 km 3 ; Andersen et al., 2017;Bradshaw, 2017;Shamloo and Till, 2019;Szymanowski et al., 2017), although many of the available constraints utilize relatively small data sets, and thus are less suitable to see if there is complexity in the thermal histories experienced by different magmatic components within an individual magma reservoir.
One way to investigate this issue further is to conduct more detailed studies of a single large eruption, and to combine diffusional constraints with other petrologic approaches to provide additional context within which to interpret results. In addition, use of multiple trace elements for diffusion modeling can also make diffusion constraints more robust (e.g., Chamberlain et al., 2014;Morgan and Blake, 2006;Shamloo and Till, 2019;Till et al., 2015).
Here we utilize diffusion modeling of Sr and Mg in plagioclase and detailed documentation of core to rim changes in plagioclase chemistry in the 2.08 Ma, 630 km 3 Cerro Galán Ignimbrite (Folkes et al., 2011c;Kay et al., 2011). We also combine these results with thermodynamic modeling (e.g., Rhyolite-MELTS; Gualda et al. 2012) and other petrological observations to constrain the long-term thermal evolution in this system.

Geologic Background
The Cerro Galán caldera, located in northwest Argentina on the eastern edge of the Puna plateau, is a part of the larger Central Volcanic Zone of the Andes (Figure 1). Between 5.6 and 2.08 Ma the Cerro Galán magmatic system produced >1200km 3 of high-K, crystal rich (40-50%), chemically homogenous (68-71wt% SiO2) rhyodacite (Folkes et al., 2011a(Folkes et al., , 2011c) that have been classified as a 'monotonous intermediates' (Hildreth, 1981). The largest and most recent eruption, the 2.08 Ma Cerro Galán Ignimbrite (CGI), produced 630km 3 of ignimbrite that extends outward from the caldera 40km in all directions and up to 80km north of the caldera (Folkes et al., 2011c).
Fe-Ti oxide geothermometry on white pumice record temperatures of 790-820°C (Wright et al., 2011), however Folkes et al., (2011 note that due to the high fO2 (+1-1.7 NNO) calculated for the CGI these temperature estimates may have uncertainties of 50-100 °C. Trace element ratios in the CGI, when compared to previous ignimbrites from the Galán magmatic system, indicate magmas that experience closed system evolution (Folkes et al., 2011b). Appearance of sanidine in the CGI, coupled with the disappearance of amphibole present in previous eruptions from the Cerro Galán system has been experimentally demonstrated to be the result of a shallowing of the magmatic system over time, and this is also supported by volatile contents in quartz-hosted melt inclusions (Grocke et al., 2017). Previous reconnaissance studies (e.g., Folkes et al., 2011b;Wright et al., 2011, Bradshaw, 2017 of the thermal evolution of the Cerro Galán Ignimbrite indicate that there are multiple populations of plagioclase with one population recording magma storage for long durations at high temperature (> 750 °C; type 1), while the other (and more common) records short time at high temperature (type 2). Additionally, these different plagioclase populations contain unique trace element signatures with type 1 plagioclase exhibiting a positive correlation between Ba and An, while type 2 exhibits a negative correlation.
This has been interpreted as being the result of plagioclase recording crystallizing environments in which the magma system switches from conditions that inhibit sanidine crystallization (type 2) to those that promote it (type 1) as the magma system shallows with time (Bradshaw, 2017). It is also possible, however, that these trends may be heavily influenced by the presence of biotite (found up to 15 vol% vs. 1.5 vol% sanidine in the CGI) and reflect a magma reservoir with diverse crystallization conditions, further establishing the need for an in depth, multidisciplinary assessment of the time-temperature evolution of the Cerro Galán magmatic system.

Bulk Rock Geochemistry
Eight samples from representative locations around the Cerro Galán caldera were chosen for bulk rock geochemical analysis. Samples were analyzed at Washington State University for both major (SiO2, TiO2, Al2O3, FeO, MnO, MgO, CaO, Na2O, K2O, P2O5) and trace element (La, Ce, Pr, Nd, Sm, Eu, Gd, Tb, Dy, Ho, Er, Tm, Yb, L, Ba, Th, Nb, Y, Hf, Ta, U, Pb, Rb, Cs, Sr, Sc, Zr, V, Cr, Ni, Cu, Zn, Ga) compositions. Major elements, high field strength elements (HFSE), and their corresponding uncertainties were characterized using a ThermoARL AdvantXP XRF following the method of Johnson et al 1999. Rare earth elements (REE) and remaining trace elements along with their uncertainties were measured by inductively coupled plasma mass spectrometry (ICP-MS) using an Agilent 7700 Q-ICP-MS according to Knaack et al., (1994). Bulk rock major and trace element concentrations, along with their uncertainties can be found in online supplementary material.

Electron Probe Micro Analysis
Major element (Na, Si, Al, Fe, Ca, K, Mg, Ti) spot analyses and backscattered electron (BSE) images of plagioclase were conducted using Cameca SX100 electron probe microanalyzer (EPMA) at Oregon State University with a focused beam of 5µm, 15kV accelerating voltage, and 30nA current beam current. Calibration used standard reference material NMNH 115900.
Uncertainties were calculated by repeated measurement of standard reference materials. Average measured values, accepted values, and their uncertainties can be found in Table 1.

Laser Ablation Inductively Coupled Plasma Mass Spectrometry (LA-ICP-MS)
Trace element analyses of pumice glass and plagioclase were conducted using a Photon Machines Analyte G2193 ArF Excimer laser system connected to a ThermoFisher Scientific by 50 µm with 5 µm between spot centers, pulse rate of 50 Hz, and analysis time of 10 seconds per spot. This approach allows for high spatial resolution to be maintained in the direction of the transect, while still allowing for higher count rates to improve precision. Ablated material was carried to the mass spectrometer using an Aerosol Rapid Introduction System (ARIS; Teledyne Photon Machines Inc., Bozeman MT, USA) microcapillary tube to decrease washout time, allow for easier ablation peak identification, and help improve detection of low concentration elements.
All total, 81 plagioclase transects across 9 samples were gathered. Profile locations within individual plagioclase were chosen following the advice provided in Shea et al., (2015) to mitigate the influence of sectioning effects and merging diffusion fronts in our diffusion models.
Trace element analyses of microlite free pumice glass were conducted using 30 µm diameter circular spots, a pulse rate of 7 Hz, laser energy of 6.42 J·cm -2 , and analysis time of 30 seconds per spot. A total of 16 glass analyses were utilized in this study.
Elemental concentrations were calculated from analytes using the methodology of Kent and Ungerer, (2006) and Longerich et al., (1996) and the software LaserTRAM-DB (Lubbers et al., 2021). Anorthite contents of plagioclase were calculated using measured Ca/Si ratios similar to the method of (Kent et al., 2008). BCR-2G was used as the calibration standard and was analyzed every 5 profiles or 10 spot analyses (i.e., ~15 minutes of analysis time) along with ATHO-G to monitor for drift in the mass spectrometer. BCR-2G, ATHO-G, NIST-612, and BHVO-2G were run as standard blocks at the beginning, middle, and end of each experiment.
Their concentrations and uncertainties can be found in the Supplementary Data. Analyses suggest precision for trace elements measurements in all materials are < 5% except for heavy rare earths, Co, and V in glass analyses, which are < 10-15%.

Diffusion Chronometry
This study utilizes Mg and Sr diffusion in plagioclase to quantify the duration that a given crystal (or area of a crystal) has experienced a certain temperature. Trace element partitioning in plagioclase is dependent on An content (Bindeman et al., 1998;Dohmen and Blundy, 2014;Nielsen et al., 2017) and follows an Arrhenius relationship such that: 1.
( ! ) = "# + Where R is the gas constant, T is temperature in Kelvin, XAn is the molar fraction of anorthite, and both A and B are constants (Bindeman et al., 1998;Nielsen et al., 2017). Likewise, the rate at which self-diffusion of both Mg and Sr occurs within plagioclase is dependent on An content (Cherniak and Watson, 1994;Costa et al., 2003;Giletti and Casserly, 1994;LaTourrette and Wasserburg, 1998;Van Orman et al., 2014). Having both the partition and diffusion coefficient values dependent on An content then necessitates a diffusion equation that incorporates these observations. As such we utilize the solutions from Costa et al. (2003) to model how Mg and Sr profiles will change with time:  and Casserly (1994), respectively. Equilibrium profiles were calculated by using Equation 1 and the composition of the most rimward analysis in the transect used to calculate an equilibrium liquid composition. We then assume that this liquid composition is what the plagioclase in the measured profile is striving to equilibrate with. While these analyses are not always exactly at the rim of the grain, we find that calculated equilibrium liquids for each profile are broadly consistent with one another (i.e., 0.4 ± .1 wt% Mg; 157 ± 37 ppm Sr), regardless of whether or not they are at the core or rim of a given grain. Furthermore, as long as the observed trace element profile in the plagioclase is trying to equilibrate with a single liquid composition (real or hypothetical) this assumption is valid (Costa et al., 2003).

Initial model conditions
Boundary conditions used as starting profiles for diffusion models were chosen using An profile shapes as a guide. This is justified for two reasons: 1) NaSi -CaAl interdiffusion in plagioclase is extremely slow (e.g., Grove et al., 1984;Liu and Yund, 1992) such that it does not change appreciably on magmatic timescales at 750°C; 2) Anorthite profile shapes are approximately the inverse of the calculated equilibrium profiles as partition coefficient is inversely related to An (e.g., Bindeman et al., 1998;Nielsen et al., 2017). This then means that if a diffusion profile at t = >10 4 years was tracked back through time to t = 0 years, it produces an initial boundary condition that mimics the width and height of an An profile. We also show that using boundary conditions that deviate too far from the inverse of the equilibrium profile produces solutions to the diffusion equation that do not approach the observed trace element profiles. (Supplementary Figure 2). Applying this logic, we create either a step-function profile or profile that closely mimics the shape of the An profile, and focus on regions where there are large changes (steps) in An concentration. While not as numerically complex as other models for estimating initial boundary conditions that require assumptions about the initial magma composition and how it evolves with time (e.g., Costa et al., 2003;Druitt et al., 2012;Mutch et al., 2021), we find that these boundary conditions accurately approach calculated equilibrium conditions and are therefore suitable for the purposes of this study. Furthermore, after a number of attempts we concluded that forward models of both Sr and Mg evolution in the CGI magmatic system were unable to reproduce the observed relationships between Sr, Mg, and An and if used for calculating initial Sr and Mg profiles in plagioclase, would introduce considerable additional uncertainty into our diffusion models. The benefit of using step function-like initial profiles is then two-fold in this situation: 1) less uncertainty introduced by inaccurate trace element evolution models; 2) diffusion model results reflect diffusion times that are maxima as it assumes that compositional changes in plagioclase are instantaneous relative to crystal growth.

On constraining lock-up temperatures
The rate at which diffusion occurs is strongly dependent on temperature. As such, temperature is a critically important input in all diffusion models. In the case of long-term magma storage, it is difficult to accurately assign a temperature, and storage conditions are unlikely to be isothermal in any case. We instead follow the approach of Cooper and Kent (2014) where we constrain the maximum time that a given crystal could spend at a given temperature.
We further use the relationship between viscosity and temperature (e.g., Marsh, 1981) to constrain the maximum time that magma could spend in a state where it could be considered to be rheologically mobile. Bradshaw (2017) used Rhyolite-MELTS (Gualda et al., 2012) modeling to model changes in melt and magma viscosity throughout the evolution of the CGI magmatic system using the relationships put forth by Giordano et al., (2008) and Scaillet et al., (1998), respectively. As with many silicic magmas (e.g., Fish Canyon Tuff; Caricchi and Blundy, 2015), the CGI magma experiences a drastic increase in magma viscosity as it reached crystallinities of 30-40% and temperatures of 750 ± 10°C, such that it reached values near or exceeding the proposed magma extrusion limit (e.g., 10 6 -10 8 Pa×s; Scaillet et al., 1998;Takeuchi, 2011); which would imply that this temperature and crystallinity approximates an inflection point whereby the magma rapidly changes from being rheologically mobile to immobile (e.g., uneruptible). While we note that highly crystalline magmas experience non-Newtonian behavior (e.g., Mader et al., 2013;Pistone et al., 2013;Caricchi et al., 2007), and these viscosities are approximations of complex multi-phase fluid dynamics (e.g., Bergantz et al., 2017) and are heavily influenced by the presence of an exsolved volatile phase (e.g., Okumura et al., 2019;Pistone et al., 2013), further evidence to suggest that 40-50% crystallinity may be at the upper limit of magma eruptible magma viscosities is the observation that crystallinities of erupted (and therefore eruptible) volcanic rocks are rarely above 45-50% (Marsh, 1981), and that the crystallinities of Cerro Galán rocks themselves are only as high as 40-50% crystals (Folkes et al., 2011b).
Due to the crystallinity of the Cerro Galán system changing rapidly with small changes in temperature once it reaches 30-40 % crystals (Bradshaw, 2017), we argue that modeling diffusion at 750 °C is sufficient to broadly capture the amount of time the majority of the magmatic system spent at or above temperatures required to produce the observed crystal fraction as well as crystal fractions ranging from 30 ->80% (Supplementary Figure 3). While there is no "one size fits all" crystal fraction that large silicic magma reservoirs reach lock up at (Bergantz et al., 2017), we emphasize that changing the temperature significantly from 750 °C in our diffusion models would either result in quantifying timescales of magma storage at conditions required to produce magmas that are too crystal poor (i.e., increasing the temperature of our diffusion model) or crystal rich (i.e., decreasing the temperature of our diffusion model) relative to the threshold between eruptible and uneruptible magmas.
Recently it has also been shown that magmas or regions of a magma reservoir may be stored for significant periods of time below the solidus and erupted without prior rejuvenation or remobilization (Mucek et al., 2021). Furthermore, there is also evidence for re-incorporation of subsolidus magma reservoir rinds into voluminous eruptive products (Andersen et al., 2017;Rivera et al., 2016). Whereas we do not dispute the observation that magmas near the margins of large reservoirs likely experience much cooler and, occasionally, subsolidus temperatures for long durations of time, we stress that these portions of the system are unlikely to be representative of the overall magma storage conditions that contribute to the bulk of the 630 km 3 of erupted material associated with the CGI, as they are volumetrically minor relative to total erupted volume. As we are primarily concerned with the duration at which the CGI magmatic system spent at temperatures where there were large volumes of eruptible magma, we set our diffusion model temperatures at 750°C for reasons listed above. We interpret our diffusion models to then quantify the maximum duration a given grain, or portion of that grain, resided in a magma that was hotter than 750°C, and thus was stored at conditions at which the surrounding magma was sufficiently crystal poor to be mobile and eruptible.

Model uncertainties
The best fit diffusion model to the observed trace element profile was assessed using a standard chi-squared test that assigns a 'goodness of fit' for each iteration of the model. The smallest chi-squared value for a given model is then deemed to be the best fit diffusion time based on the input parameters. Uncertainties for each diffusion model were evaluated using a Monte Carlo approach where 1000 random profiles were generated for each modeled profile.
Random profiles were based on the analytical uncertainty at each analyzed point (i.e., for each point in the profile a normally distributed random number was generated based on the observed mean analysis value and its 1 sigma uncertainty). Best fit diffusion times were then fit to each random profile in the Monte Carlo simulation keeping the initial boundary conditions and temperature fixed. As we are quantifying a duration above 750°C, for reasons listed above, and not a measured temperature with an uncertainty, we keep this temperature fixed as well. Overall uncertainties for a given diffusion model were calculated by taking the mean and standard deviation for the Monte Carlo simulation. Most of the distributions from Monte Carlo simulations were determined to be log-normal, so prior to taking the mean and standard deviation they were first transformed such that they fit a normal distribution. Calculated means and standard deviations from the transformed distribution were then back transformed so that they were in correct time units (e.g., years and not log[years]). All data handling, figure creation, and diffusion modeling was completed using the programming language Python (version > 3.6) via the IPython environment (Pérez and Granger, 2007) and relies heavily on the following packages: numpy (Harris et al., 2020); scipy (Virtanen et al., 2020); matplotlib (Hunter, 2007); pandas (McKinney, 2010); statsmodels (Seabold and Perktold, 2010); seaborn (Waskom, 2021).
For a more complete explanation of the diffusion modeling process and the intricacies of the model, we have created a Jupyter notebook that is available for download (https://github.com/jlubbersgeo/diffusion_chronometry).

Whole Rock
Bulk rock major and trace element data for our samples are consistent with previous work (e.g., Folkes et al., 2011b;Wright et al., 2011). Pumice found within the CGI is relatively homogenous in SiO2 and other major elements (Figure 2A), but has been interpreted to be the result of fractional crystallization (Folkes et al., 2011b). Wright and others (2011) note, however, that there is a less dominant pumice species found within the CGI that displays elevated Ba concentrations, lower Sr/Ba, and Eu/Eu * . Rare earth element diagrams for the CGI show shallow Eu anomalies and middle -heavy rare earth concentrations of 10x chondrite consistent with other large silicic eruptions that have been interpreted to be from relatively cold, wet, and oxidizing conditions (Bachmann and Bergantz, 2008;Deering et al., 2008Deering et al., , 2010. Although there are some more subtle differences in Eu/Eu*, overall REE trends are similar for both pumice types Figure   2B).

Plagioclase
Plagioclase from the CGI has An values that are normally distributed around An36±4.
There are, however, small numbers of analyses that display An>60 and An<25. High An values are exclusively found in plagioclase cores. Plagioclase trace element compositions from the CGI, like An, show mostly unimodal trace element distributions and plagioclase compositions largely overlap between white and gray pumice (Figure 3). However, when plagioclase analyses are observed incorporating spatial information, we find that there are small geochemical heterogeneities on the intracrystalline scale (see section below on Ba-An distributions).

Diffusion of Mg and Sr in plagioclase
As discussed above, both Sr and Mg in plagioclase appear to be out of diffusive equilibrium with their plagioclase host. This is shown by the observation that both Sr and Mg profiles are significantly different from calculated equilibrium profiles (e.g., Figure 4) Magnesium in plagioclase diffusion models were completed for 44 grains and best fits range in duration from <1.5 to 358 years ( Figure 6). While this is a large range, 95% of times are less than 108 years. Strontium in plagioclase diffusion models were completed for 34 grains and range in duration from <10 years to 969 years with 95% of times being less than 452 years.
Diffusion coefficient and analytical resolution dictate the lower limit of resolvable timescales achievable when modeling 1D diffusion (Bradshaw and Kent, 2017). This effectively puts a lower 'resolution' boundary on all diffusion models, regardless of what the 'best fit' model suggests. At An50 (CGI 3 sigma upper value) this is 1.5 years for Mg and 10 years for Sr. Two grains modeled for Mg diffusion and one grain modeled for Sr diffusion return times less than these respective limits, and we simply report these as < 1.5 and < 10 years.
Representative diffusion models for Mg and Sr in plagioclase can be found Figure 7 and Using these initial Sr and Mg distributions and at a temperature of 750°C, our diffusion models report decadal to centennial timescales for transects that are in both the cores and rims of grains.
For grains where multiple transects were measured throughout, core transects consistently record longer diffusion timescales than rim transects (Figure 9).
There is significant overlap in the ranges of apparent timescales recorded by Mg and Sr ( Figure 6A), however there are also some significant discrepancies between the Sr and Mg best fit diffusion times in individual crystals ( Figure 6B). We find that as Mg diffusion times increase, Sr diffusion times do not increase in a predictable way (e.g., Figure 6C), implying that this observation is likely not the result of either systematic differences in partition or diffusion coefficients. Furthermore, despite Mg diffusing roughly 5 times faster than Sr at the appropriate CGI An content and at 750°C, Mg profiles consistently match the profiles of An and other slow diffusing elements (e.g., Ba) profiles more closely than Sr profiles (Figure 7; Supplementary   Figure 4). This observation alone suggests that little diffusion has occurred.
One explanation for the differences between Sr and Mg diffusion models might be differences in the variation of Sr and Mg that occurs with changes in An content during plagioclase formation. Elements that have large discrepancies in diffusion rates should produce noticeably different diffusion widths after a sufficiently long period of time (e.g., Morgan and Blake, 2006). Conversely, when their diffusion widths and magnitudes appear similar, they will yield different diffusion times reflecting their individual diffusion rates, subsequently leading to the interpretation that diffusion has not progressed significantly from initial boundary conditions (Shamloo and Till, 2019;Till et al., 2015). While CaAl -NaSi diffuses so slowly (Cherniak, 2010 and references therein) that it is not useful in quantifying timescales of volcanic processes, we apply similar logic to Till et al. (2015) inasmuch as Mg profiles in the CGI plagioclase simply could not have existed for >10 3 years at 750°C while still resembling the observed An profile more than equilibrium profile shapes ( Figure 10A) and maintaining a positive correlation with An ( Figure 10B). As some Sr diffusion models record diffusion times up to 1500 years, we suggest that initial step functions based on the shape of the An profile may be an overestimation of boundary conditions for Sr profiles in some instances and, since the Sr diffusion coefficient is significantly slower than Mg, lead to overestimations in the time a given grain has experienced at a certain temperature. We hypothesize that this may be the result of Sr and Mg partitioning in plagioclase having different sensitivities to changes in An (i.e., the "A" parameter in Equation 2; Bindeman et al., 1998;Nielsen et al., 2017). Furthermore, the faster diffusion of Mg means that overestimation of the diffusion times are likely to be smaller with incorrect boundary conditions, and thus Mg diffusion models are likely to: 1) be less sensitive in inaccuracies of diffusion model boundary conditions; 2) more accurately represent the duration at which the CGI magma existed at or above 750°C.
As mentioned in the Methods section, we suggest that the erupted CGI magma would have uneruptible viscosity at temperatures below ~750°C. Thus, our best fit diffusion timescales at this temperature for both Mg and Sr in plagioclase are interpreted to represent the maximum total amount of time that an individual grain, or portion of a grain, could have spent at magma reservoir conditions sufficient for the magma to be mobile. In addition, each grain is recording its own individual thermal history, and as diffusion is continuous from the moment a given chemical potential (i.e., zone boundary) within the mineral forms, then the calculated timescale provides a "thermal budget" for all the processes that might impact the thermal state of the crystal. This includes initial crystal growth -which can occur over a wide range of temperatures in the Cerro Galán magmatic system (Bradshaw, 2017), including temperatures over 750°C, as well as long term storage, and transient reheating events associated with magma recharge and eventiual eruption (e.g., Rubin et al., 2017).
Despite this potential complexity, our results collectively indicate that the population of CGI plagioclase we have examined, starting from when they began crystallizing in the CGI magmatic system, only spent decades to centuries at temperatures required to produce large volumes of eruptible melt. Moreover, despite analysis of a much larger number of crystals than previous studies we see broadly unimodal timescale distributions ( Figure 6A) with little evidence for plagioclase populations that have experienced systematically different thermal histories. This may relate to the relatively high uncertainties in our diffusion chronology estimates, but also may reflect that diffusion of heat in silicic magmas is relatively rapid (e.g., Jaeger, 1964) and heat transport is more rapid than transport of mosty chemical constituents via diffusion. Our Sr and Mg diffusion in plagioclase diffusion results also align with other plagioclase diffusion studies that investigated storage conditions of caldera forming eruptions (Bradshaw, 2017;Druitt et al., 2012).

Plagioclase as a recorder of magma reservoir conditions
Plagioclase is a solid solution found in virtually all igneous rocks. It displays a myriad of growth, dissolution, and zoning textures, is capable of preserving long crystallization and mixing histories (e.g., Cooper et al., 2001;Davidson et al., 2001;Singer et al., 1995), and crystallizes over a wide P-T-X-fO2 range (e.g., Bowen, 1915Bowen, , 1913Drake, 1976), making it a key mineral to study for observing long-term magma reservoir processes. It also has many useful geochemical applications such as U-series dating (e.g. Cooper and Reid 2003), diffusion modeling (e.g. Costa et al. 2003;Druitt et al. 2012), and thermometry (e.g., Putirka 2008 and references therein). In the following section we explore how plagioclase can be used to infer long-term magma reservoir thermochemical histories using a combination of trace element partitioning and diffusion chronometry.
Plagioclase crystallizes over a wide range of temperatures in the CGI magma system (Bradshaw, 2017), suggesting that it should reflect crystallization from a melt that is becoming more progressively depleted with elements such as Sr and Mg that are compatible in plagioclase (e.g., Bindeman et al., 1998;Nielsen et al., 2017) or compatible in phases that are cocrystallizing with plagioclase. Using partition coefficients from Nielsen et al., (2017) and our measured plagioclase compositions we have calculated melt compositions in equilibrium with observed plagioclase compositions and compared these to the erupted pumice glass compositions ( Figure 11). Except for the lowest An compositions, the erupted glass composition for the CGI and the liquid compositions in equilibrium with CGI plagioclase do not overlap, and the vast majority of plagioclase equilibrium melt compositions are less evolved than the observed pumice glasses. While a decrease in temperature to near solidus (e.g., ~650-680) temperatures shifts plagioclase equilibrium melt compositions to lower Sr and Mg abundances, even at these low temperatures most calculated equilibrium compositions are still less evolved than the erupted glass compositions. Assuming the Kd values used are broadly correct this suggests that most plagioclase at Cerro Galán crystallized at a relatively early stage of magmatic evolution and in equilibrium with melts that are more enriched in both Sr and Mg (which is also consistent with the observed Sr-An and Mg-An relations we observe). Although there are uncertainties associated with any partitioning model for any phase, we find most calculated plagioclase equilibrium liquids and erupted glass compositions do not overlap within 2-sigma uncertainty.
These observations are also true when using partition coefficients from Bindeman et al., (1998), suggesting that calculated partition coefficients, themselves, are not a significant source of error.
Rather, we propose that the wide array of melt compositions in equilibrium with plagioclase compositions are the result of plagioclase recording magma reservoir conditions that do not reflect those immediately prior to eruption, but longer-term magma storage that is evolving due to fractionation of predominantly plagioclase, sanidine and biotite. Figure 11 also illustrates that nearly all CGI plagioclase formed in equilibrium with a liquid that was much more enriched in Mg. As biotite is the only phase within the CGI in which Mg is a major component (Mg# 58-66; Folkes et al., 2011b), we take this observation to imply that the majority of CGI plagioclase crystallized from a magma that was either lacking biotite or had crystallized less biotite than what is observed in erupted pumice.
Trace element contents also provide constraints for the thermochemical conditions that CGI plagioclase are recording. Barium is less compatible in plagioclase compared to lower temperature phases in the CGI (i.e., biotite, sanidine). Barium also diffuses sufficiently slowly (Cherniak, 2002) that measured Ba contents will not be significantly modified by diffusion.
Looking at the Ba-An relationship in a single grain thus allows us to determine whether or not biotite or sanidine were present as that portion of the plagioclase was growing, as a decrease in Ba with mineral growth implies a co-crystallization with a Ba compatible phase. Observed Ba-An relationships in plagioclase reveal a diversity of behavior and can broadly be grouped into six categories (fractionation with Ba compatible phases, fractionation without Ba compatible phases, recharge with Ba compatible phases present, recharge without Ba compatible phases present, little to no fractionation, and some combination of these trends; Figure 12) revealing that individual plagioclase crystals are not only recording different crystallization environments, but also have the potential to record multiple crystallization/dissolution events prior to their evacuation from the reservoir(s) in which they grew. This is also in contrast to the bulk rock major and trace element chemistry which suggest a more chemically homogeneous system (this study; Folkes et al., 2011b). To determine if there are statistically significant differences recorded in diffusion timescales between these groups, we performed a two-sample Kolmogorov-Smirnov (KS) test which tests the null hypothesis if two samples are drawn from the same distribution. In this case our samples were plagioclase that show signs of cocrystallization with biotite/sanidine and plagioclase that do not. We find that despite plagioclase recording very different chemical environments, there is no statistically significant difference in diffusion timescales between groups (i.e., KS statistic below critical value and p-value >> 0.05).
These results hold for both Mg and Sr diffusion results. Collectively these data are interpreted to show that although the magma system from which the CGI was erupted contained significant reservoir scale chemical heterogeneities through time, and that individual plagioclase experienced and recorded quite different petrologic environments, there are no substantive differences recorded in thermal history between these crystals.

Constraints from zircon
To further help quantify the long-term thermal state of the Cerro Galán magmatic system, we also compare zircon residence times (Folkes et al., 2011a), zircon saturation temperatures (Boehnke et al., 2013;Watson and Harrison, 1983), and plagioclase diffusive equilibration times. Folkes et al., (2011a) found that many zircons from the CGI have cores that crystallized at > 10 5 but < 10 6 years prior to eruption. This implies that parts of the CGI magmatic system must have been below zircon saturation temperatures for at least that duration of time. To calculate zircon saturation temperatures throughout the evolution of the CGI magmatic system we utilized Rhyolite MELTS for Excel (Gualda et al., 2012;Gualda and Ghiorso, 2015) as it quantifies many extensive properties of a magmatic system while it crystallizes from fully liquid to solid.
Starting compositions for MELTS models are average glass compositions of individual samples and can be found in the Supplementary Data. Log fO2 is constrained between +1 -1.7 NNO (Folkes et al., 2011b) and pressures range from 100 -200 MPa (Grocke et al., 2017). For each down temperature step in a given MELTS model a M parameter, [(Na + K + 2Ca)/(Al*Si)], was calculated for the remaining liquid composition. Zircon saturation temperatures were then calculated using (Boehnke et al., 2013) and we find that across all MELTS models zircon saturation occurs at 807 ± 7°C (Figure 13). Calculating zircon saturation temperatures from measured glass data yields near identical zircon saturation temperatures of 807 ± 5°C. While generally regarded as a robust and faithful recorder of magma reservoir processes in silicic systems, Bindeman and Melnik (2016) show that when exposed to temperatures above the zircon saturation temperature, even larger zircons (e.g., 200µm diameter) will dissolve in 10 3 -10 4 years. Therefore, for the CGI zircons to record such long residence times, while also remaining hundreds of microns in length, they could not have been above zircon saturation temperatures for a significant period. This also puts a limit on the thermal state of the CGI magmatic system, and when combined with plagioclase diffusion timescales, provides compelling evidence that prior to eruption the CGI magmatic system waxs predominantly stored at sufficiently low temperatures such that stored magmas were unlikely to be at an eruptible viscosity.

Periodic recharge
Plagioclase from the CGI exhibit mean An values of An35, however display frequent An>55 zones throughout many analyzed grains (Figure 4). An increase in An of this magnitude may due to a multitude of processes such as temperature or pressure increase Streck (2008) or rapid ascent of water saturated magmas Blundy et al., (2006), however, we propose that this increase in An is due to recharge of a less evolved magma, as increases in An are typically not found at the rims of a given grain. The high An zones are also associated with significantly higher concentrations of Mg (Figure 4) which may also be further support for influence from a less evolved magma. These observations also effectively rule out rapid ascent of water saturated magmas as a likely cause of high An zones, as rapid ascent driven changes in An would also likely be one of the last things a plagioclase experiences prior to eruption and be found at the rims of grains or as infilled sieved textures, both of which we do not see.
Our diffusion chronometry of Mg also puts thermal and temporal limits on these recharge events, as our data provide an estimate of the total time that a crystal could have been at temperatures in excess of 750°C. Thus, diffusion data define a thermal "budget" within which the analyzed plagioclase must form, be stored, and then erupted. Given the relatively short durations we estimate for residence at temperatures >= 750°C, it implies that although recharge events clearly occur, they may have had relatively little thermal input to the system overall, and are mainly reflected in changes in composition, not temperature.
Our data also show that grains that experience recharge occur throughout the fractionation history of the CGI system (i.e., some profiles have Ba-An relationships indicating co crystallization with biotite/sanidine, where others have Ba-An relationships indicating crystallization without biotite/sanidine). Rhyolite-MELTS models predict that biotite and sanidine begin crystallizing relatively late in CGI crystallization sequence (i.e., ~30°C lower than plagioclase). We use this to delineate between plagioclase that record relatively late stage (postbiotite/sanidine saturation), and longer-term magma reservoir processes that occur prior to biotite/sanidine saturation. We take the ubiquity of high-An zones in both these groups of plagioclase to imply that the Cerro Galán magmatic system has been periodically recharged with a magma that produces plagioclase of at least An60 and potentially as high as An75 composition.
This idea is also supported by: 1) the presence of a volumetrically minor, less evolved, pumice population (i.e., gray pumice) that Wright et al., (2011) used to suggest that there was a second, less evolved melt reservoir deeper in the crust that briefly interacted with CGI magmas prior to eruption; 2) volatile zoning in apatite indicating recharge shortly before eruption by a volatilerich magma into a reservoir that contained local variations in composition (Boyce and Hervig, 2008).

Long term thermochemical conditions
The results from our plagioclase trace element and diffusion modeling show that the plagioclase measured in this study are recording a range of long-term and late stage magma reservoir processes within the Cerro Galán magmatic system, and that overall, eruptible magmas are only present for periods of decades to centuries or less, regardless of whether or not plagioclase show signs of long-term residence (i.e., crystallization without biotite and sanidine) or later-stage crystallization (i.e., co-crystallizing with biotite and sanidine). This is important as it implies that short diffusion times are not simply reflective of young crystals and that the longterm thermal history of the CGI magmatic system was one was dominated by storage temperatures below 750°C. Below, we address how these timescales fit with previous literature on magma remobilization timescales and what it means for the long-term thermochemical state of the magma system.
Previously, the Cerro Galán magmatic system has been classified as homogenous on the basis of bulk rock and major phase chemistry (Folkes et al., 2011b). Whereas homogeneity of a given system is ultimately defined by the scale of investigation (e.g., Fish Canyon Tuff, Bachmann et al., 2002), the CGI is homogenous inasmuch as it does not display bulk rock major element zoning characteristic of many large silicic eruptions (e.g., Carpenter Ridge Tuff, Bachmann et al., 2014;Youngest Toba Tuff, Chesner and Rose, 1991;Kneeling Nun Tuff, Szymanowski et al., 2019). The CGI, however, also has been interpreted as the result of assimilation-fractionation of a 50-50 mix of local metamorphic basement and mantle derived basalts (Folkes et al., 2011b;Kay et al., 2011), suggesting that there may have been large scale thermochemical zoning related to progressive differentiation at some point in its history (Bachmann and Bergantz, 2008). Huber et al., (2012) find that these reservoir scale heterogeneities may be removed in crystal rich, eutectoid systems with high degrees of exsolved volatiles on the order of years or less, which is significantly shorter than the time required to reactivate a given mush body via repeated injection of hotter magmas. While longer than homogenization times, other numerical models suggest that reactivation times for a hydrous, crystal rich, magma systems may be on the order of hundreds to thousands of years (Hartung et al., 2019). These results align with our diffusion modeling results that show CGI plagioclase spent only decades to centuries at or above temperatures required for production of eruptible magma. Based on heterogeneous crystal cargo that display frequent zoning (Figure 4, Figure   7, Figure 8), evidence for long term storage with repeated recharge of less evolved, volatile rich material (Boyce and Hervig, 2008;Wright et al., 2011;Figure 13), and short times at remobilizing temperatures ( Figure 6, Figure 10, Figure 11), we argue it is highly plausible then, that the Cerro Galán magma system was in a chemically heterogeneous but uneruptible state for a significant portion of its history and was only homogenized and remobilized shortly before eruption.

CONCLUSION
We applied diffusion modeling of Mg and Sr in plagioclase from the 2.08 Ma, 630 km 3 , Cerro Galán Ignimbrite (CGI) to investigate the long-term thermochemical magma storage conditions of large silicic magma reservoirs. Our results indicate that, although bulk-rock geochemistry implies a magmatic system that is chemically homogeneous at the time of eruption, plagioclase are recording diverse crystallization environments (i.e., late stage co-crystallization with Ba-compatible phases such as biotite and sanidine as well as earlier crystallization without Ba-compatible phases) throughout their crystallization history that largely do not reflect those at the time of eruption. Furthermore, despite this evidence for long-term storage within the reservoir, CGI plagioclase are only recording timescales of decades to centuries at temperatures required to maintain eruptible volumes of magma. We interpret these results to reflect that the CGI magmatic system existed in a predominantly immobile, but chemically heterogeneous state for most of its residence within the crust. These results help to further constrain the long-term thermochemical storage conditions of magma systems capable of producing catastrophic caldera forming eruptions. Figure 2: A) Total alkalis vs silica (TAS) diagram for CGI pumice compared to other CAVZ volcanics from the last 10Ma. Note the overlapping between the two pumice types as well as the unimodal distribution of data in both silica and total alkalis. B) REE diagram for CGI pumice (gray and white) compared to other CAVZ volcanics from the last 10 Ma. It can be seen here that CGI pumice have REE trends that are more consistent with cold, wet, oxidizing than hot, dry, reducing rhyolites from other large volcanic eruptions (Deering et al., 2010). Figure 3: Panel of selected trace elements in CGI plagioclase filtered by what pumice type they are from (i.e., gray or white). Similar to bulk rock pumice data, plagioclase from the CGI overlap in composition for all trace elements measured and show unimodal distributions. Note, however, the elevated Mg concentrations found at An60 and above. Figure 4: Plagioclase from the CGI exhibiting high-An cores. High-An cores in CGI plagioclase also consistently contain significantly higher Mg concentrations (e.g., > 200 µg/g). Some high-An core boundaries have Mg profile that also mimic the shape of the An profile well suggesting little to no diffusive equilibration (A), while others show significant deviation from the shape of the An profile (B). . Gray area is the range that would be predicted using Cs = Kd*Cl where Kd is calculated using Equation 1 and observed median glass composition. Based on the observed global positive correlations between Sr and Mg vs. An, CGI plagioclase are not in diffusive equilibration. Note, that even though diffusively equilibrated plagioclase do not fall within the gray region in (B, D), their slope still broadly matches that of the gray region. Gray region location on the charts is ultimately determined by an assumed liquid composition. B and D then suggest that CGI plagioclase have formed in a liquid that is elevated in Sr and Mg relative to erupted glass compositions.     Different colored lines (A) and symbols with regression lines (B) indicate different durations of diffusion ranging from initial boundary conditions to equilibrium and illustrate that as Mg profiles near equilibrium (e.g., long durations of time at high temperature), they begin to have a negative correlation with An. As negative correlations are never observed between Mg and An in CGI plagioclase, we take this to mean that 1) initial Mg profiles based off An profiles accurately represent the initial Mg concentration at the time of plagioclase formation (i.e., initial profiles in diffusion models) 2) CGI plagioclase thermal histories are more accurately recorded by Mg diffusion than Sr diffusion 3) CGI plagioclase have not spent long durations of time at or above 750 °C. Figure 11: Comparing melt compositions in equilibrium with CGI plagioclase (colormapped circles) to observed glass compositions (red circles) for Ba, Sr, and Mg. Equilibrium melt compositions are calculated using Equation 1 for the plagioclase partition coefficient. We can see that the observed glass is much more depleted in Sr (from plag fractionation), Mg (from biotite fractionation), and Ba (from biotite, sanidine, and, to a lesser degree plagioclase, fractionation). The elevated Mg equilibrium liquid comps suggest that much of the plagioclase that we analyzed crystallized (at least partially) before biotite fractionation. As biotite in the CGI are typically greater than a couple of mm and euhedral we interpret this as evidence for plagioclase experiencing extended periods of time in the magma reservoir. We can see in this wide array of trends that plagioclase is crystallizing in multiple environments; some with Ba compatible phases (e.g., biotite and sanidine), and some without. Combined with Figure 11, this suggests that while most of the plagioclase equilibrium liquid compositions at the time of transect formation are less evolved (potentially pre biotite), individual grains are experiencing heterogeneous chemical environments where biotite/sanidine are coming in and out of the liquidus. Figure 13: Results from MELTS and zircon saturation modeling. For each down temperature step, a given MELTS model's remaining liquid is either zircon undersaturated ( !"# < 0) or saturated ( !"# > 0), where !"# is the difference between the system temperature and the calculated zircon saturation temperature. MELTS models suggest zircon saturation temperature of the CGI is 807 ± 8 and for our observed glass data this is 807 ± 5. This suggests two things: 1) MELTS models are accurately (within error) portraying zircon saturation temperatures for the CGI magmatic system and 2) zircon saturation in the CGI magmatic system is relatively low temperature.