Long-term temperature cycling in a shallow magma reservoir: insights from sanidine megacrysts at Taápaca volcano, Central Andes


 Hybrid dacite magmas from Taápaca volcano in the Central Andean volcanic zone (18ᵒ, N. Chile) contain sanidine crystals of unusual size (1 to 12 cm) and abundant mafic enclaves of variable composition throughout the entire eruptive history (1.5 Ma to recent) of the volcano. They are rich in mineral inclusions and strongly zoned in Ba with distinct growth bands separated by resorption interfaces. Resorption is followed by a sudden increase in Ba with compositional contrasts up to 2.3 wt% BaO. We argue that resorption and the sharp jumps in Ba-concentration reflect distinct heating and melting events, suggesting that different growth zones formed at different temperatures. Amphibole-plagioclase thermo-barometry based on mineral inclusions gives variable temperatures of ∼720 – 820 ᵒC at shallow pressures (0.1 – 0.3 GPa) for individual growth zones. Using these temperatures for diffusion modelling, Ba-profiles from x-ray scanning profiles and grey scale gradients based on accumulated BSE images across these interfaces allow to estimate crystal residence and reactivation times prior to eruption. This temperature control allowed the application of a “non-isothermal” diffusion algorithm to obtain diffusion times for individual diffusive boundaries that range from 0.4 to 490 ky and add up to total residence times of 9 to 499 ky for different crystals from different stages of eruption. A combination of temperatures, pressure, diffusion times and R-melts modeling of the parent rhyodacite suggests storage conditions for the Taápaca reservoir at near eutectic composition at shallow depth (4 – 10 km). Temperatures never fell below the magma solidus but frequently cycled between 720 °C and 820 °C, i.e. between eruptible and non-eruptible state with crystallinity circling around ∼40 – 50 vol%, for tens to hundreds of thousands of years. We define this as “Long-term Transitional Temperature Cycling” or LTTC storage. Frequent recharge events of basaltic andesite magma, as represented by abundant mafic enclaves, orchestrated the temperature cycling, resulted in multiple heating events that caused frequent resorptions and interrupted crystal growth, and kept the reservoir thermally ‘alive’. Recharge events became more frequent only ∼3 – 11 ky before the eventual eruption that carried a particular set of sanidine megacrysts to the surface. Thus, after many earlier recharge events that did not result in eruption, a final event involved mixing at a critical recharge rate to mobilize, entrain, and erupt a particular set of megacrysts from the resident rhyodacite in a hybrid dacite host. This process, happening not more than a few centuries before an eruption, has been repeated at similar timescales at different stratigraphic stages throughout the 1.5 My long history of Taápaca volcano. The observed mineral zonation patterns and size of sanidine crystals from the resident magma reservoir below Taápaca volcano are identical to those observed in the megacrysts from granite intrusions that also show typical age ranges of zircon crystallization which are comparable to the residence times extracted here from Ba zonation. Taápaca sanidines thus may represent an erupted equivalent and provide “smoking gun” evidence of temperature cycling during the formation of such K-feldspar megacrysts in granites.


INTRODUCTION
Constraining rates of magmatic processes, durations and conditions of crystal and magma storage are key insights into the pre-eruptive history of magma systems. Such studies therefore serve as essential tools in understanding potentially dangerous and active volcanoes, and may help in eruption forecasting, hazard assessment and risk management. Some of the most hazardous volcanoes are silicic volcanoes, which have relatively long eruption histories with low eruption frequencies that can erupt with no or little anticipation (e.g. 2008 rhyolitic eruption of Chaitén volcano, Chile; Alfano et al., 2011;Saubin et al., 2016). Such volcanoes and in particular, large silicic caldera-forming magma systems, are apparently fed from reservoirs with extended magmatic histories where magmas of variable crystallinity can be stored for thousands to hundreds of thousands of years (e.g. Reid et al., 1999;2017;Brown & Fletcher, 1999;Charlier & Zellmer, 2000;Simon et al., 2008;Wotzlaw et al., 2013;Szymanowski et al., 2017;Keller et al., 2018;Schaltegger et al., 2019;Karakas et al., 2019). Even longer times of magma accumulation, residence, and crystallization of a few Ma were documented by high-resolution zircon dating from intrusive rocks (e.g. Coleman et al., 2004). Yet the time to (re-)activate and assemble large volumes of magma may be relatively short (Bachmann & Bergantz, 2008;Druitt et al., 2012;Cooper & Kent, 2014;Rubin et al., 2017) compared to the magnitudes of storage http://www.petrology.oupjournals.org/ times. This has led to new concepts for the evolution of silicic magma systems, where magmas reside as non-eruptible crystal mushes for longer times and may become eruptible relatively fast via recharge-induced heating and mush reactivation (Bachmann & Bergantz, 2008;Reid, 2008;Annen et al., 2008;Huber et al., 2012;Ellis et al., 2014;Cooper & Kent, 2014).
During storage, evolution and activation of magma reservoirs, the associated temperature, pressure and compositional changes are imprinted into the composition of growth zones in the minerals that crystallize from the magmas (Costa & Chakraborty, 2004;Morgan & Blake, 2006;Chakraborty, 2008;Ginibre et al., 2007;Costa & Morgan, 2010). In order to constrain the rates of processes controlling the evolution and eruption of magmas with their crystals, two principal methods have been developed in the past two decades. One method (1) is absolute dating of crystals by short-lived isotopes of U-decay-series (e.g. Condomines et al., 1988;Hawkesworth et al., 2000;Hawkesworth et al., 2004;Schmitt, 2011) that provides insights into their time of formation relative to eruption times and thus gives the minimum times of magma storage. For example, applications of direct U-Th-Ra mineral dating (Condomines et al., 1988) e.g. at Santorini, Greece  and Kilauea Volcanoes, Hawaii (Cooper et al., 2001) has constrained histories of crystallization, magma transport and residence spanning over up to hundreds of ky. A second method (2) is the estimation of pre-eruptive residence times of crystals using diffusion chronometry of zoned crystals, which has been applied to several minerals e.g. plagioclase (Zellmer et al., 1999;Costa et al., 2003), pyroxene (Morgan et al., 2004;Chamberlain et al., 2014;Petrone et al., 2016), olivine (Costa & Dungan, 2005;Costa & Chakraborty, 2004;Sundermeyer et al., 2020), quartz (Chamberlain et al., 2014) and sanidine (Chamberlain et al., 2014;Iovine et al., 2017;Rout & Wörner, 2020).
Both methods suffer from limitations. Cooper & Kent (2014) reviewed results from these methods and compared diffusion ages to U-series ages for a large range from small to large, and mafic to silicic magma systems. Particularly for long-active silicic systems, it appears http://www.petrology.oupjournals.org/ that durations of diffusion are usually significantly shorter than absolute pre-eruptive crystal ages recorded by U-series dating. This discordance indicates that the crystals seem to have spent a significant portion of their "life" at low temperatures where diffusion is sluggish. Thus, estimating magma/crystal residence times from diffusion modeling critically depends on the thermal history of a crystal, which can be constrained using various thermometers (e.g. Putirka, 2008;2016). Only if the thermal history suggests that the crystals spent most of their storage time well above the diffusion threshold, and this temperature is well-constrained, then the diffusion times will correspond to the actual crystal residence times (e.g. Morgan & Blake, 2006;Petrone et al., 2018).
We will argue that the dacites are products of mixing/mingling between a sanidinemegacryst bearing rhyodacite and mafic andesite recharge magma. With detailed mineralogical, geochemical and isotopic investigation, we will identify temperature to be the dominant factor controlling the size and zonation of the megacrysts. As the main motivation behind this study, we apply core-to-rim non-isothermal diffusion modeling (Petrone et al., 2016(Petrone et al., , 2018Di Stefano et al., 2020;Rout et al., 2020) to extract diffusion and storage times with the help of intracrystal thermometry. This approach gives more reliable diffusion time estimates than using a constant temperature value for the entire growth history. These results provide evidence for a sustained crystal residence conditions with cyclic high-T excursions from below to above the The 'modern' cycle of active continental margin magmatism in the Central Andes started around 180 Ma. Jurassic intrusive and extrusive arc rocks are exposed in the Coastal Cordillera. Starting from Jurassic times, the focus of volcanism migrated for about 150 km from the coastal Cordillera in the West towards the East and widened across the western Altiplano after the onset of crustal thickening in Miocene times (e.g. Scheuber et al., 1994;Wörner et al., 2018). Reasons for arc migration are debated and could be due to massive tectonic erosion of the forearc and changing subduction angles through time. Presently, the subducting slab descends below the CVZ at a relatively shallow angle of ~25 to 30ᵒ (Cahill & Isacks, 1992) at a convergence rate of 75 -80 mm/a (Samoza, 1998). The volcanic front of the CVZ lies approximately 120 -150 km above the subducting slab, and magmas from the mantle wedge must traverse remarkably thick continental crust (>70 km) (James, 1971;Zandt et al., 1994;Allmendinger et al., 1997;Scheuber & Giese, 1999;Yuan et al., 2002;Tassara et al., 2006).
The crust below the Altiplano consists of the Palaeoproterozoic (2.0 -1.8 Ga) Arequipa terrane represented by metamorphic and igneous rocks exposed in southern Peru and in a series of inliers in northern Chile and western Bolivia (Wörner et al., 2000a, Loewy et al., 2004McLeod et al., 2013). Proterozoic amphibolites and gneisses of the metamorphic Belen Complex are exposed in northernmost Chile along the western slope of the Andes just 35 km to the S of Taápaca volcano, and granulites and charnockites of Cerro Uyarani occur on the western Altiplano, as reported by Wörner et al. (2000a) and are potential equivalents to the basement rocks underlying the volcanoes in northern Chile.

History and compositions of lavas erupted at Taápaca volcano:
The eruptive history of the TVC has been described by Kohlbach & Lohnert (1999) who identified three stages of evolution based on contact relationships, mineral mode, flow morphology and vent location. These stages have been refined further and presented in more http://www.petrology.oupjournals.org/ detail by , who recognized four evolutionary stages based on geochronological and morphological criteria on the volcanic complex, events of gravitational sector collapse and subsequent migration of the main vents (Fig. 1). In this study, we use the stratigraphy and eruption stages defined by  and complemented by additional 40 Ar/ 39 Ar sanidine dating .
The volcanic edifice of Taápaca is constructed on a deeply eroded topography on the western margin of the Western Cordillera and consists of an arcuate cluster of domes of which the highest reaches 5,850 m asl at the present summit. The main edifice volume is estimated to be 35 km³; the eruptive products of Taápaca cover an area of 250 km² .
The initial eruptive stage (Stage I) formed a gently outward-dipping stratovolacano consisting of moderately porphyritic (< 15% crystals) two-pyroxene andesite lava flows (~60 wt% SiO 2 ) containing small amounts of sanidine and hornblende. Stage I is estimated to be older than 1.5 Ma. Stage II forms the main volume of the TVC and ranges in age from 1.5 to 0.5 Ma. The oldest 40 Ar/ 39 Ar ages from the earliest stage II is 1.46 ±0.07 Ma , from plagioclase) and 1.27±0.004 Ma (Wörner et al., 2004, from sanidine). It consists of viscous dacite lava flows, which formed a steeply-dipping stratovolcano. A major collapse event is documented within stage II by a voluminous debris avalanche. Stage III (0.5 -0.47 Ma) consists of small volume lava domes and block-and-ash flow deposits concentrated mainly in the central part of the dome complex. The eruptive products differ from stage II only by higher amounts of sanidine megacrysts and mafic enclaves. Partial collapse of the southern part of the ancestral stage II edifice and stage III domes mark the start of the youngest and http://www.petrology.oupjournals.org/ morphologically most complex stage IV (0.45 Ma -recent). The dacites generated during Late-Pleistocene -Holocene eruptive activity of the TVC are petrographically undistinguishable from those of units II and III. Stage IV forms the present central edifice of TVC, characterized by extrusion of voluminous domes and associated block-and-ash flows, blast deposits, tephra fallout, and pyroclastic flows and lahars.
Numerous debris avalanche deposits resulted from frequent edifice collapse events at TVC.  recognized two types (triggers) of debris avalanches at Taápaca: (1) debris avalanches that are a consequence of gravitational instability of the edifice due to extensive hydrothermal alteration, and (2) by intrusion of a cryptodome causing deformation.
The latter collapse triggering mechanism is evident from blast deposits associated with the debris avalanche, which results from a rapid decompression after abrupt mass unloading.
Despite such repeated mass unloading events, a notable change in the composition and mineralogy of erupting products after edifice collapses has not been observed.
All dacites from stage II to IV contain mafic enclaves and sanidine megacrysts next to phenocrysts of plagioclase, amphibole, biotite, titanite, Fe-Ti-oxide and accessories (apatite, zircon, xenocrystic clinopyroxene). Quartz and traces of small, anhedral clinopyroxene also occur as xenocrysts. Abundance and size of both, mafic enclaves (1 -6 vol% and 2 -15 cm in size) and sanidine megacrysts (1 -7 vol% and 1 -12 cm in size) increase from stage II through IV of the eruptive history of Taápaca volcano .

Samples and petrography
In this study, the sanidine megacrysts from dacite host rocks from stage II (1.5 -0.5 Ma) to stage IV (0.45 Ma -recent) were examined. These porphyritic dacites (Fig. 3a) are dominated by up to 30 vol% phenocrysts of mostly plagioclase, which are typically up to 5 mm in size.
Sanidine megacrysts are highly variable in size but generally larger than 1 cm. The goundmass of fresh dacites range from black in colour, dense, partly glassy to light grey and friable textures due to incipient degassing and associated crystallization of microphenocrysts. The groundmass primarily constitutes of plagioclase and amphiboles. Oxidized samples show transitions from light grey to pink and brownish-red colors due to syn-eruptive oxidation of dome rocks.
Plagioclase occurs as crystallites, microcrysts and phenocrysts, with phenocrysts dominating the population. Common textures in plagioclases include sieve textures with oscillatory zoning and multiple resorption interfaces. The phenocrysts are generally unaltered and rarely contain mineral inclusions like euhedral titanite, euhedral to subhedral amphibole and sporadically http://www.petrology.oupjournals.org/ zircons. Amphibole crystals (Fig. 4) are euhedral to subhedral, ranging from crystalite to phenocryst sizes, with some crystals showing prominent zoning after resorption. Most amphibole phenocrysts occur as completely unaltered and only rarly show breakdown rims to oxide and pyroxene. Biotite (up to 3 mm) occurs as fragmented euhedral to subhedral phenocrysts. They are infrequently altered and rarely contain inclusions of apatite and zircon.
Oxide and clinopyroxene crystals are generally euhedral to subhedral microcrysts. Quartz crystals are up to 2.5 mm in size, rounded, fractured, and strongly resorbed.
Sanidine megacrysts (Fig. 3b -j) range in shapes from euhedral with sharp crystal faces to rounded, and range from 10 to 120 mm in size in the longest dimension. The megacrysts contain large amounts of euhedral to subhedral inclusions ( Fig. 4a) of plagioclase, amphibole, titanite, Fe-Ti oxide, glass, sporadically biotite and apatite, rarely quartz and minute crystals of zircon. These inclusions are often aligned parallel to the crystal rims tracing the crystal shapes and distinct growth zones (up to 5 mm wide), which are already evident in the thin sections ( Fig. 3j). Glass inclusions appear in partially filled channels and holes that irregularly cross the sanidine growth zones. Plagioclase inclusions are up to 4 mm in size, with larger inclusions occuring within the crystal cores. Amphibole and titanite inclusions are up to 1.5 mm and 3 mm in size, respectively. Both amphibole and titanite also dominate the mineral population within the glass inclusions (up to a few mm in size). The size, shape, composition and zonation patterns (for more detail on zonation based on elemental mapping, see below) of these megacrysts are indistinguishable from K-feldspar megacrysts that are frequently observed in granitic rocks (e.g. Vernon, 1986;Moore & Sisson, 2008;Slaby et al., 2008;Johnson & Glazner, 2010).
The mafic enclaves (2 -15 cm in size) are rich in acicular amphibole and contain subordinate plagioclase, both in a poorly versicular glassy groundmass. Besides plagioclase and amphibole microcrysts, the enclaves also contain relatively large and strongly resorbed http://www.petrology.oupjournals.org/ xenocrysts of quartz and plagioclase as well as amphibole and biotite with breakdown textures along their rims.

Sample preparation
24 euhedral sanidine megacrysts from eruption stages II to IV of the Taápaca volcanic complex were studied. The crystals were impregnated with a silicate-binding agent and then heated in the furnace at 90 ᵒC for ~12 hours. This procedure was repeated five times to saturate cracks and harden the crystals to prepare them for further steps of slicing by diamond saw. The crystals were cut at 0.5 mm from the crystal's center perpendicular to the (010) crystallographic plane and perpendicular to the cleavage on this plane. The larger half of the crystal was ground planar using different grain sizes of SiC, lapped and polished. Polished surfaces were glued on to a glass slide and the crystals were cut carefully by a low speed saw to obtain a 400 μm thick slab remaining on the slide. The top side was again ground, lapped and polished. To avoid contamination during electron microprobe work, the samples were cleaned and heated at 75 ᵒC prior to carbon coating to ensure electrical conductivity at the surfaces during measurements.
For isotopic analyses, we prepared a second slab from one of the largest crystals, ca 1 mm thick and embedded and fixed it between two glass slides. This provided a position just off-centre and opposing the thin section and allowed to precisely locate the second slab with respect to the Ba growth zones mapped by electron microprobe. The slab was cut into a rod vertical to the growth zones and that was cut into small cubes (1 mm 3 ) from core to rim. By superimposing the location of these cubes with respect to the Ba-intensity map analysed by EMP on the directly opposing thin section plane, all samples could be directly located with respect to high-Ba and http://www.petrology.oupjournals.org/ low-Ba growth zones. The sanidine cubelets were liberated from the glass by dissolving the resin and cleaned from inclusions by handpicking.

Electron microprobe and LA-ICPMS
Microprobe measurements were done with a JEOL JXA-8900R microprobe at GZG, Universität Göttingen. Quantitative analyses used an accelerating voltage of 15 KV with a 15 nA beam current and a 5 -15 μm beam diameter. Counting times for all elements (except Ba and Sr) were 15 sec on peak and 5 sec on the background. Ba and Sr in feldspars were measured at a maximum of 120 sec and 60 sec on peak and background, respectively. Calibration standards for the feldspars were: sanidine for K, albite for Na, anorthite for Si, Ca and Al, celsian for Ba, SrTiO 3 for Sr, and hematite for Fe. Calibration standards for the amphiboles were: sanidine for K, albite for Na, olivine for Si and Mg, anorthite for Al, wollastonite for Ca, hematite for F, TiO 2 for Ti, Rhodonite for Mn, Cr 2 O 3 for Cr, and NiO for Ni. The relative standard deviation for major oxides was below 2% and the calculated absolute errors for minor oxides were below 0.03 wt%.
Accumulated back scattered electron (BSE) images were prepared by superimpositions of ten individual BSE image accumulations acquired in COMPO mode with 20 kV accelerating voltage and 20 nA beam current and in slow scanning mode with an acquisition time of 120 sec per accumulation. Detailed element mappings concentrated on inclusion-poor regions before identified by the COMPO-images. During each mapping, five elements, Ba, Ca, Fe, Na and K were measured simultaneously. Additionally, a BSE-image of the same area was also recorded.
The mapping sections/areas were between 0.12 cm × 0.5 cm and 0.5 cm × 0.5 cm, with multiple http://www.petrology.oupjournals.org/ (up to 3) sections per crystal. A stage scan was used with an accelerating voltage of 20 kV, beam current of 40 -50 nA and a beam diameter of 10 μm. The grey-scale values in the BSE images show a linear correlation with the Ba intensities in the Ba-maps (see below) and thus can be used as a proxy for Ba-content.
In addition, x-ray line scans for Sr, Ca and Ba were also acquired in COMPO-mode stage scan conditions to resolve the compositional gradients across the zone boundaries in 16 crystals. 50 accumulations per scan were acquired perpendicular to the zone boundaries at 10 kV accelerating voltage and 30 nA beam current, with an acquisition time of 2 sec per pixel.
LA-ICPMS analyses were also done at GZG, Universität Göttingen with a Perkin Elmer SCIEX Elan DRC II instrument coupled to a Nd-doped YAG UV laser (266 nm). Trace-element line profiles were selected according to previously mapped elemental compositions and done parallel to electron microprobe traverses for selected samples. Beam diameters were 80 to120 μm, and total ablation times along the profile ran for up to 1,200 sec per laser scan at a laser frequency of 10 Hz and dwell times of 20 -50 ms. For further and more precise location of profile and spot measurements, the thin sections were repolished after laser ablation, again carbon coated and BSE maps were repeated to precisely locate LA-ICPMS measured positions with respect to the distinct Ba growth zones. Elemental contents are recorded as counts per second (cps) for the line scans with relative standard deviation <8%. Certified "NIST 610" glass was used as the reference material. Electron microprobe and LA-ICPMS data are given in electronic appendix 1, 2 and 3.

Isotope measurements
Oxygen isotope compositions of dacites and sanidine separates were measured by Laser fluorination IRM gas chromatography-mass spectrometry at the Department of Isotope http://www.petrology.oupjournals.org/ Geology at GZG, Universität Göttingen. Laser fluorination uses an Nd-doped UV laser for insitu analysis and a CO 2 laser for bulk materials. The oxygen isotopic composition was measured with an irmGCMS Finigan MAT262 gas mass spectrometer. The instrument and analytical techniques have been described by Wiechert & Hoefs (1995), Wiechert et al. (2002) andHagen (2002). The variations in the oxygen isotopic ratios are expressed as the δ notation and relative to VSMOW. The used standards, garnet UWG-2 at the CO 2 laser line and MORB glass at the UV laser-assisted in-situ measurements, were defined as δ 18 O SMOW = 5.7. The precision of the measurement was determined to be 0.14 -0.17% (Wiechert et al., 2002). With the UV laser insitu measurements, the oxygen isotopic composition of the megacryst sample TAP 28/1 was measured at six compositional zones with measurements at three to five spots per zone. Laser settings were as follow: laser energy: 260 mJ, voltage: 30 kV, frequency: 50 Hz, measuring time: 60 sec for standard (3,000 laser-pulses) and 70 sec for the sample (3,500 laser pulses).
The intra-crystal oxygen isotopic composition of four other megacrysts (18 K1, 11 K1, 18 B3 and 34 K2) was measured with a CO 2 laser line using 1 mg of material that was separated from specific individual zones (from inclusion-free homogeneous sites) in crystal-core, -mantle and -rim zones. The choice of zones ascertained to include at least one Ba-poor and one Ba-rich zone per crystal. For the bulk rock analyses of the host dacites, about 1 mg rock powder (size <65 microns) of each sample was used. The samples were chosen so as to represent all stages of eruption.
Isotope ratios of 87 Sr/ 86 Sr were meaured on sanidine cubelets separated across the coreto-rim profile of the sanidine megacrsyt TAP 28/1. Analyses were done by TIMS on a Finnigan MAT262-RPQ + multicollector mass specrtometer at Department of Isotope Geology, GZG, Universität Göttingen. 3 to 4 mm large pieces were cut parallel to the (001) surface and perpendicular to the seven prominent growth zones of the crystal. Each of these small pieces were crushed and carefully hand-picked to select pure sanidine without mineral inclusions.
http://www.petrology.oupjournals.org/ They were leached in diluted HCl, rinsed and cleaned in milipure water using an ultrasonic bath. 100 mg of sample powder was dissolved in 6 ml HF: 32 HNO 3 (1:1) and heated for 16 hours at 200 ᵒC within Savillex beakers. The solution was completely evaporated at 140 ᵒC followed by repeated (twice) dissolution and evaporation in 4 ml 6N HCl. For the last evaporation, it was re-dissolved in 2.5 ml 2.6 N HCl, stored in PE vials and centrifuged. For separation, the sample solution was rinsed with 2.6 N HCl through columns of ion exchange resin BIORAD AG 50W-X8 Resin, 200 -400 mesh. The strontium-rich elution fraction was then collected and evaporated. For measurement, Sr was dissolved in 0.5 N H 3 PO 4 and mounted on Re-double filaments (~1 μg). Measured Sr isotope ratios were corrected for mass fractionation to 87 Sr/ 86 Sr = 0.1194 and normalized to values for the standard NBS987 (0.710245). Standard values over the period of the study were 0.71025 ±3 (14 analyses).
External 2σ errors on the 87 Sr/ 86 Sr ratio are estimated at <0.004%. Total procedural Sr-blanks are 0.4 ng and thus, neglegible for our measurements with sample content averaging at ~900 ppm Sr. The Sr-isotope data is provided in electronic appendix 4.

Mineral Chemistry and isotopic composition
Amphibole compositions in the dacites are bimodal (Fig. 5a) with low Al-Ti Mg-hornblende and high Al-Ti Mg-hastingsite. Mg-hornblende forms up to 2 mm large crystals in the dacite and is the exclusive type of amphibole found as inclusions in sanidine megacrysts. By contrast, Mg-hastingsites are common in the dacite as small (<0.3 mm), poorly zoned crystals but also as overgrowth after resporption on the Mg-hornblendes in the dacite (Fig. 4). Mg-hastingsite also is the exclusive type of amphiboles in mafic enclaves. Like the amphiboles, the plagioclase http://www.petrology.oupjournals.org/ compositions (Ab 45 -Ab 80 ) in the dacite are also bimodal ( Fig. 5b) with low-Fe (<2,000 ppm) and high-Fe (2,000 -4,500 ppm) concentration. Low-Fe plagioclase occurs as individual crystals in the groundmass of the dacite and as cores of larger phenocrysts. Inclusions of plagioclase in sanidine megacrysts have a range in An similar to that of the high-Fe plagioclase in the dacite host, but they are, however, exclusively of low-Fe composition. By contrast, high-Fe plagioclase occurs as microcrysts in the dacite groundmass and as phenocryst rims, and is the only type of plagioclase in the mafic enclaves.
Compositions of sanidine megacrysts vary from Or 65 to Or 80 with An-contents ranging between 0.3 and 1 mol%. Spectacular compositional zonation in Ba (0.2 -2.5 wt% BaO) shows multiple oscillatory growth bands from core to rim ( the sanidine data but also reach to slightly higher values (7.1 -8.6 ±0.2‰). 87 Sr/ 86 Sr ratios are also found to be uniform at ~0.707 throughout the analyzed crystal 28/1 (Fig. 9, electronic appendix 4).
http://www.petrology.oupjournals.org/ Our objective is the documentation and interpretation of the megacrysts' compositional zonation as a basis for Ba-diffusion modelling with the ultimate goal of unraveling the conditions and duration of magmatic residence of the sanidine megacrysts. However, before we describe the zonation in detail, we first need to clarify the origin of the host dacite and its genetic relation to the sanidine megacrysts, then constrain the melt composition from which the sanidines formed, and finally define their P-T conditions of their crystallization.

Evidence for magma mixing and endmember compositions of Taápaca dacites
There are several lines of evidence that argue for magma mixing as a principal process in forming sanidine-megacrystic dacite magmas at Taápaca: (1) Mafic enclaves occur in abundance and are found in different stages of disintegration.
We consider the mafic enclaves as evidence for the ongoing process of hybridization and argue that the observed enclaves represent mafic recharge that has not yet been completely homogenized. This interpretation is in accord with findings by Humphreys et al. (2009) and Ruprecht et al. (2020). Their studies suggest that (i) mingling and mixing are two closely related and not contrasting processes, in which (ii) enclaves are continuously being eroded and immersed to homogeneity within their host and (iii) such erosion and disintegration of, for example, a 5 cm-large enclave would take between a few days and a couple of decades. This indicates that the observed mafic enclaves are related to the (latest stage of the) process that actually formed the hybrid dacites.
(2) We observe bimodal compositions of amphiboles in the forms of Mg-hornblende and Mg-hastingsite in the dacite. Mg-hastingsites are observed together with high-Fe http://www.petrology.oupjournals.org/ plagioclase in mafic enclaves whereas Mg-hornblende are the exclusive type of hornblende associated with low-Fe-plagioclase as mineral inclusions in sanidine megacrysts (Fig. 5). Both types of amphibole and both types of plagioclase populations are observed in the dacite groundmass. Additionally, Mg-hastingsite and high-Fe-plagioclase are also frequently observed as the post-resorption overgrowths in the Mg-hornblende (Fig. 4b) and low-Feplagioclase populations, respectively, in the dacite.
(3) Bulk rock dacite compositions form clusters in variation diagrams rather than crystal fractional trends (Fig. 2). At first glance, this does not seem to support a magma mixing relation but certainly excludes a formation by fractional crystallization. However, Blum-Oeste & Wörner (2016) observed that individual mafic enclaves, that are highly varied mostly in trace element composition form tie-lines with its specific host dacite that project to a common point of intersection. This intersection identifies a silicic mixing endmember as a rhyodacite (see Additional evidence is provided by the composition of plagioclase crystals: Plagioclase inclusions in the sanidine are exclusively low in Fe (< 2,000 ppm), as opposed to the bimodal low-and high-Fe plagioclases compositions in the dacite. This suggests that the parent melt of the sanidine would not have more than 2.9 wt% of FeO (calculated with a plagioclase-melt partition coefficient of 0.09; Ginibre et al., 2002). This value agrees with the FeO-content of 0.7 -2.7 wt% in the rhyodacite endmember estimated by Blum-Oeste & Wörner (2016) and the rhyodacites observed at Parinacota, and clearly contrasts with the higher FeO contents (3.0-5.2 wt%) observed in the host dacite.
We conclude that the sanidine megacrysts did not form from their host dacite but crystallized from a more silicic rhyodacite melt. Such melt, however, was never actually erupted from Taápaca volcano. This is highly relevant for the correct interpretation of the sanidine zonation patterns (Zellmer & Clavero, 2006;Higgins, 2011).

Conditions of crystallization
To constrain crystallization temperature, we applied thermometric formulations to mineral inclusions in the sanidine megacrysts, as well as to the phenocrysts in the dacites and the mafic http://www.petrology.oupjournals.org/ enclaves. The Holland & Blundy (1994) amphibole-plagioclase thermometer used here (Table   2) is based on the edenite-tremolite and edenite-richterite equilibria between amphibole and coexisting plagioclase in silica-saturated rocks. It obtains temperature values with uncertainties of ±28 ᵒC (1σ) and was applied to Mg-hornblende and low-Fe-plagioclase pairs from the inclusions in the sanidine megacrysts (Fig. 4a). In the host dacite, the same thermometer was applied to the Mg-hornblende and Mg-hastingsite populations (Fig. 4b) while pairing them with the dacite-hosted low-Fe and high-Fe plagioclases, respectively. We also applied the thermometer to the Mg-hastingsite and high-Fe plagioclase populations in the mafic enclaves as well. The compositions used in the thermometer were systematically measured at the rims of the phenocrysts and the mineral inclusions in the megacrysts.
Because of the large size of the sanidine crystals and the abundance of amphibole and plagioclase inclusion in them, it was possible to determine crystallization temperatures for individual growth zones with respect to the documented Ba-zonation. These temperatures give a relatively large range between 720 and 820 ᵒC that is clearly outside the error of the method (Table 2, Fig. 11). Temperatures for Mg-hornblende-plagioclase pairs from the host dacites also lie in a similar range between 740 and 830 ᵒC. By contrast, the Holland & Blundy (1994) thermometer gave significantly higher temperatures for Mg-hastingsites from the dacite and for the mafic enclaves (~860 -1050 ᵒC and ~950 -1100 ᵒC, respectively). These results document a bimodal distribution of amphibole crystallization temperatures for the compositionally different phenocryst populations (Fig. 11), representing both, a more mafic and a more evolved magma component for the hybrid dacite.
The thermometric data was also verified by using the amphibole-only thermometer after Ridolfi & Renzulli (2012) applied to the Mg-hornblendes and Mg-hastingsites. This thermometer yields temperatures within ±30 ᵒC from the amphibole-plagioclase thermometer, which is similar to the uncertainty range of the thermometers.
http://www.petrology.oupjournals.org/ equilibrium with alkali feldspar, biotite, Fe-Ti-oxides, titanite, apatite and quartz. In the case of Taápaca sanidine megacrysts, plagioclase and amphibole crystals are included in the sanidine with the other required minerals, but rarely quartz. In the dacite, quartz is also rare with only a few strongly resorbed grains observed in thin sections. R-MELTS modeling (see below) and experimental evidence (R. Botcharnikov pers. comm.) indicates that quartz should also be a phase in equilibrium with this set of minerals in the rhyodacite host of the sanidine. However, quartz inclusions, just like the phenocrysts in the dacite, likely dissolved upon mixing and heating prior to eruption. Therefore, the required assemblage for applying the plagioclase-   Fig. 11). Mg-hastingsites in the mafic enclaves also indicate pressures between 0.3 ±0.30 and 0.7 ±0.08 GPa. We appreciate that the application of the Ridolfi & Renzulli (2012) barometer can be problematic, as has been shown in many studies (see Erdmann et al., 2014;Putirka, 2016). In fact, the direct comparison between Mg-hastingsite composition and their crystallization conditions based on an experimental study of a Taápaca dacite composition suggests that, at high pressures (>0.3 GPa), the Ridolfi & Renzulli (2012) barometer underestimates the true crystallization pressure by ~0.1 -0.2 GPa (R. Botcharnikov pers. comm.). Since the range of calculated pressures for the Mg-hastingsites is already higher compared to the pressure derived from the Mg-hornblendes, this underestimation would indicate that the Mg-hastingsites may indeed have formed at greater depth from more mafic magmas compared to the Mg-hornblende that crystallized from the rhyodacite. We note, however, that a precise knowledge of the depth of crystallization and residence (or the lack thereof) will not affect the conclusions from our diffusion modeling and sanidine residence time estimates.
To confirm our results, we additionally used other barometric formulations based on Al in hornblende (after Johnson & Rutherford, 1989;Schmidt, 1992;Anderson & Smith, 1995) and found that results are consistent within 0.08 GPa. The results of geothermo-barometry are documented in an overview form in Table 2 and as individual P-T values in Fig. 11. Apart from temperature and pressure, we constrain the fO 2 with the formulations given by Ridolfi & Renzulli (2012) for amphibole compositions. The details of all thermo-barometric data are given in electronic appendix 5. Considering these results, the estimated conditions of http://www.petrology.oupjournals.org/ Downloaded from https://academic.oup.com/petrology/advance-article/doi/10.1093/petrology/egab010/6126431 by guest on 09 February 2021 crystallization for the sanidine megacrysts lie between 1 and 0.3 GPa at an oxygen fugacity of logfO 2 between -11 and -13.

Rhyolite-MELTS modeling
Having established that the sanidine megacrysts formed from a rhyodacite magma at around 720 -820 ᵒC and 0.1 -0.3 GPa, we will now try to further constrain temperatures of storage and crystallizations using Rhyolite-MELTS (Gualda et al., 2012). As the starting composition At these pressures, R-MELTS predicts liquidus temperature of about 950 -1,000 ᵒC for different combinations of pressure and water content. Plagioclase (Ab 40 -Ab 45 ) is the first liquidus phase followed by oxides and pyroxene at 940 -990 ᵒC. Sanidine starts to crystalize at 735 -770 ᵒC, (Fig. 12) followed by quartz. A cystallinity of 40 -50 vol%, which is the rheological lock-up condition where the system transforms into a rigid crystal mush (Bachmann & Bergantz, 2008), was reached at 750 -780 ᵒC. Although Gualda et al. (2012) noted that R-MELTS is not very suitable for near-solidus phase equilibria and amphibole is not represented well in the model, the solidus was predicted at a reasonable temperature of ~680 -700 ᵒC.

Links between sanidine-bearing Taápaca hybrid dacites and magma compositions at Taápaca and Parinacota volcanoes
Using again the Ba-Sr relations of Fig. 10, we will further explore the relations between observed magma compositions and their relation to sanidine crystallization. This approach is based on the premise that the endmember compositions of the hybridization process at Taápaca 6, 7 and 8), a full documentation for a total of 24 sanidine crystals mapped is given in supplementary form (electronic appendix 2). At first sight, the variability of zonation patterns and styles is staggering. However, there are certain "motifs" that can be observed in all Bamaps: The most prominent structure are resorption-growth bands (RG-type) that show a typical saw-tooth pattern with an abrupt increase in Ba immediately after re-growth following the resorption event, and a seemingly gradual decline in Ba until the next resorption interface (Fig.   7). These repeated resorption and re-growth bands dominate the zonation pattern across all crystals for 80 -95% of the radial distance between core and rim. The deeper the cut and the more pronounced the resorption interface is, the larger the compositional difference (i.e. the higher the Ba jump; Figs. 6 -8). The Ba-content varies between 0.2 and 2.5 wt% of BaO across the resorption and re-growth bands. In some cases, the resorption interface is highly irregular resulting in a patchy pattern after subsequent growth (Fig. 7c) to these zonations, we observe many irregular glass inclusions in certain growth bands. These are not "normal" inclusions formed during growth but rather from the incipient partial melting of the sanidine just prior to eruption. These inclusions do not relate to the long-term storage conditions of the crystals, which is the focus of this study.
To interpret the growth bands and the origin of the zonations, a few important observations and questions need to be discussed first.

What causes the frequent resorption events and the subsequent Ba-enrichment in post-resorption sanidine growth?
Resorption can be caused by one or a combination of several possible processes (e.g. Ginibre et al., 2007;Shore & Fowler, 1996;Pankhurst et al., 2018). Most likely are (1) direct magma mixing following recharge with high-temperature magmas, (2) a temperature change caused by a recharge without compositional mixing (Ruprecht & Wörner, 2007;Couch et al., 2001), or (3) movement of crystals to an environment of higher temperature or into melt with higher water content (Ginibre et al., 2002(Ginibre et al., , 2004. Mixing of a fresh batch of hot, less evolved and Ba-rich magma with the megacrystic resident rhyodacite magma is a plausible option (Clavero, 2002;Lohnert, 1999) to explain the frequent resorption events and the successive Ba-enrichment recorded in the crystals. After all, the abundant mafic enclaves (up to 6 vol% with 500 -1,500 ppm Ba) in the dacite (with generally ~600 -1,300 ppm Ba) and evidence for mixing between such recharge and a ryhodacite is the fundamental process that formed the Taápaca dacite (see above). However, this mingling/mixing documented by the enclave and the bimodal mineral compositions for plagioclase and amphibole in their host dacite cannot be the actual process that affected the sanidine megacrysts, because (1) the dacite is not the equilibrium melt from which the sanidine formed and (2) this mixing event, that is recorded in the mafic enclaves, is the process that led to the eruption after the sanidine megacryst had already grown. The dacite-forming mixing process is the event that re-activated the megacrysts from prior storage in a rhyodacite magma, but was not responsible for the frequent resorption of these megacrysts which the dacite brought to the surface.
http://www.petrology.oupjournals.org/ Still, repeated mixing events with a hotter and less evolved (therefore, Ba-rich) magma could have affected the rhyodacite magma during its earlier and extended storage history without triggering an eruption. However, such a recharge process would have significantly changed the bulk major element composition of the host magma. This should be recorded in the composition of the growth zones that formed after mixing and resorption. This conflicts with the constant major element composition of the sanidine megacrysts (Na, K) across the resorption zone boundaries (Fig. 6, electronic appendix 1): only Ba (and partly Sr and Ca) show significant variations. Rhyodacite and mafic recharge magmas can also be expected to be different in isotope composition, but δ 18 O values are constant within error at 7.2 ±0.2‰ not only throughout the high-Ba and low-Ba growth zones, but also across multiple crystals ( Fig.   9, Table 1). 87 Sr/ 86 Sr ratios are also found to be uniform at ~0.707 throughout the crystals, regardless of Ba contents (Fig. 9, electronic appendix 4). These observations argue against any compositional mixing with a less evolved mafic magma (or any compositionally different magma for that matter) as a process to explain the resorption and subsequent Ba-rich growth of the megacrystic sanidine crystals.
Alternatively, one could envision mixing between a low-T, megacrystic rhyodacite and a high-T, Ba-rich silicic melt, not involving any mafic recharge. However, to induce heating strong enough to resorb the megacrysts above their liquidus, the high-T silicic melt would have to be hotter than 770 -820 ᵒC as constrained by R-MELTS and the upper limit of thermomertry.
However, such a high-T magma would have to have formed independantly and thus, at such temperatures, would necessarily be compositionally different (less evolved in major elements) compared to the host rhyodacite after it was heated. Specifically, at T >770 -820 ᵒC, minerals like Na-rich plagioclase, K-feldspars and quartz, which were present in the low-T rhyodacite (and now found as inclusions in sanidine), would not be crystallizing (based on R-MELTS modelling) at the relevant pressures. In that case, we would therefore expect significant changes http://www.petrology.oupjournals.org/ Downloaded from https://academic.oup.com/petrology/advance-article/doi/10.1093/petrology/egab010/6126431 by guest on 09 February 2021 in major element composition (e.g. Ab-Or content) in post-resorption growth. But the homogeneity of both sanidine and the abundant plagioclase inclusions in high-and low-Ba zones suggests against this scenario. Similarly, the high crystallinity (at least 40 vol%) of the megacrystic rhyodacite (based on our R-Melts modelling) and the large size of the crystals would also make intra-reservoir transport to any hotter Ba-richer environment difficult, e.g. through convection, plucking from cumulates or dispersion thoughout different melt batches (Ruprecht et al., 2008;Bergantz et al., 2015). Therefore, only heating of the sanidine-bearing rhyodacite magma from an external heat source remains a valid possibility. Therefore, we consider a process such as "self-mixing" (Couch et al., 2001) or "thermal mixing" (Ruprecht & Wörner, 2007) to be the likely explanation for the resorption and inferred heating events in the rhyodacite. In such a scenario, the rhyodacite magma reservoir may have been heated as the consequence of recharge events that were sufficiently weak and/or relatively distant with respect to the location within the rhyodacite reservoir and, therefore, chemical mixing must have been insignificant. Such a scenario is supported by the magnitude of thermal diffusivity (order of 10 -7 m 2 /s; Romine et al., 2012;Hofmeister et al., 2016), which is much faster than chemical diffusivity (by at least 5 -7 logarithmic order faster than Ba or Sr at 800 ᵒC; Zhang et al., 2010). This makes thermal changes much faster and spatially extensive than chemical exchange by mixing or diffusion. A simple calculation shows that heating of the rhyodacite over a distance of a few hundred meters would only take less than 100 years. Gassparging triggered by degassing of an underplating mafic magma and release of bubbles into the reservoir (Bachmann & Bergantz, 2006;Cashman & Giordana, 2014) can also speed up the heat transfer without changing the bulk composition of the rhyodacite.
In any case, we again have to conclude that only the effect of hot recharging magmas without compositional mixing is the likely cause of the observed heating and resorption events.
http://www.petrology.oupjournals.org/ Then, the question arises -how can the thermal effect of a distant hotter recharge explain the post-resorption Ba-enrichment? When such thermal pulses result in partial resorption of crystals, slow diffusing elements that are highly compatible in sanidine (e.g. Ba, Sr; GERM Partition Coefficient Database) will become enriched at the melt boundary layer near the dissolving crystal, while other elements remain relatively constant. Afterwards, when the temperature drops again below sandine-liquidus, growth resumes from this Ba-enriched melt layer. Consequently, new post-resorption growth zones form with a sharp increase in Ba.
A simple modelling of the persistence of such a boundary layer shows that a 0.5 -1 mm-thick Ba-and Sr-rich layer would last for decades to centuries at temperatures as high as 770 -800 °C (see electronic appendix 7) allowing for episodes of enriched regrowth. The conclusion, that changing temperature, not composition, has caused the Ba-zonation patterns in Taápaca sanidine megacrysts, therefore, clearly implies and documents the frequency of the cycles and the extent of heating and cooling during the growth of these megacrysts.
We tested the correlation between the degree of resorption and the "height" of the jump in Ba concentration to the next Ba-enriched growth zone. The degree of resorption was crudely quantified by the curvatue of the resorption interface, measured by the ratio of length of the resorption interface between two points and the actual distance between the same points. When plotted against the relative jump in Ba-content (Fig. 13) across the respective resorption interfaces, this parameter indeed correlates positively with the associated Ba-jumps. This suggests that a higher degree of resorption yields a more Ba-enriched melt near the resorbed crystal, from which re-growth starts after the resorption event. We take this as evidence that the Ba-jump is likely related to an enriched Ba boundary layer after resorption of the sanidine and is a result of the melting itself.

What is the cause for the megacrystic growth?
The reason for the unusually large size (from 1 up to 12 cm) of such megacrysts has been debated for K-feldspar of similar size and zonation patterns observed in granites (Vernon, 1986;Johnson & Glazner, 2010;Glazner & Boudreau, 2011;Holness et al., 2018;Chamber et al., 2020). Higgins (1999Higgins ( , 2011 proposed a process of textural coarsening through temperature cycling abvove and below the K-feldspar liquidus. Small crystals have a higher surface energy per unit volume than larger crystals and thus, to minimize the energy in the system, crystals smaller than a particular size ("critical size", Higgins, 1999) dissolve and contribute to the growth of the lager crystals. This coarsening process can explain the scarcity of sanidine phenocrysts or microcrysts in the dacites. This process is also sensitive to the overall thermal evolution, and thus could explain the differences in crystal sizes between samples. Our new data therefore is fully compatible with the conclusions of Higgins (2011) and Glazner & Boudreau (2011) that temperature cycling is responsible for megacrystic sanidines at Taápaca and we now can even constrain the frequency and degree of the heating/cooling cycles.

What causes the gradual decrease in Ba concentration after the resorption events?
The frequent repetition of high-Ba to low-Ba growth zones with a gradual decrease in Ba following each resorption boundary forms the conspicuous saw-tooth pattern of Ba concentrations (Figs. 8, 14). There are several scenarios to explain this decrease in Ba: Since Ba is compatible in sanidine, the host rhyodacite melt will be continuously depleted in Ba during megacryst growth. The thermometric results and R-MELTS modeling suggest that sanidine growth took place at temperatures <770 ᵒC. At such temperatures, the chemical diffusivity of elements like Ba and Sr (Zhang et al., 2010) becomes too slow to replenish the melt composition near the crystal-melt interface, leading to a Ba-depletion in the crystal-melt http://www.petrology.oupjournals.org/ boundary layer. This process continues until the next resorption event and a sharp increase in Ba concentration occurs due to growth from a Ba-enriched melt after resorption, forming the observed saw-tooth pattern (Figs. 8 and 14). However, without any constraints on the size of the rhyodacite melt reservoir, it is difficult to decide whether this depletion trend reflects only the crystal-melt boundary layer or affects the entire surrounding melt body.
This process of heating-induced resorption followed by cooling, regrowth and compatibe element depletion can also explain the relatively poor correlation of Ba with Sr and Ca: (1) Because of the fast diffusivity of Sr and Ca in melts (by >100 times than Ba; Zhang et al., 2010), any post-resorption boundary layer enrichment will diffuse away much faster, making an enrichment in Sr and Ca less persistent. (2) Additionally, due to relatively smaller partition coefficients of Sr and Ca compared to Ba (GERM Partition Coefficient (Kd) Database and references therein), the degree of enrichment in the newly-grown crystal parts will be less compared to Ba. Both (1) and (2) will result in less prominent variation in Sr and Ca compared to Ba in the crystal. Further, any compositional variation in the crystal will also be more likely to homogenize faster for Sr and Ca compared to Ba because of their ~100 to 1,000 times higher diffusivity in sanidines (Cherniak, 2010). Thus, we observe only poor correlations of Sr and Ca with Ba, which in any case would be more difficult to detect due to the low Ca-and Sr-contents.

What is the cause for the fine-scale oscillatory zonation?
Fine-scale oscillatory or O-type zonation (30 -150 μm width, Fig. 7f) is observed in distinct zones and, at higher frequencies (smaller wavelengths), at the rims of the crystals. Such oscillatory zonation is associated with only small Ba-jumps (from 0.2 -0.3 to a maximum 0.7 wt% BaO). These fine-scale oscillatory growth zones closely trace the crystal shapes and probably reflect incremental diffusion-controlled growth (Haase et al., 1980 Singer et al., 1995). This type of growth implies steady-state cooling, isothermal undercooling, and a stable diffusive boundary layer. Growth may have been associated with minor resorption events as the zone boundaries often have slightly rounded corners and sporadically low amplitude, undulating resorption surfaces. Ginibre et al. (2002Ginibre et al. ( , 2007 and Shore & Fowler (1996) have discussed e.g. 'extrinsic' mechanisms vs 'intrinsic' mechanisms and small-scale boundary layer processes to explain such compositinal oscillations during crystal growth. This process might not be fundamentally different from the process mentioned above explaining the gradual Ba-decrease within each growth zone, but the spatial and temporal scales over which both of these processes operated (and associated temperature fluctuations) must have been very different. Therefore, since their formation does not reflect, and does not allow for estimates of diffusion and residence times of the crystals, we will not consider this type of zonation further.

How can we constrain temperature and time for the growth history of the megacrysts?
If, as argued above, the Ba zonation is primarily formed by only thermal variations, then the zonation patterns, and in particular, the size of Ba-concentration jumps at resorption interfaces, should be a proxy for temperature variations. To test this hypothesis, we correlated the temperatures obtained from amphibole-plagioclase thermometry based on mineral inclusions from individual growth bands (i.e. between major resorption interfaces) to the Bacontent in profiles from core to rim (Fig. 14a). The temperature variation correlates reasonably well with the Ba-content during growth from core to rim, and thus we speculate that even minor Ba-variations may reflect changes in temperature. In Fig. 14b, we show schematically a temperature vs time history during crystal growth (similar to Higgins, 2007). Since there is no growth during periods of melting, we cannot constrain the temperature at the time of resorption http://www.petrology.oupjournals.org/ or how long the resorption period was. Based on previous models (e.g. Ginibre et al., 2004;Cooper & Kent, 2014;Rubin et al., 2017) resorption is thought to be a short event of rapid heating above the liquidus followed by slow cooling and longer periods of (re)growth (Fig.   14b).
Since the zone boundaries are associated with resorption interfaces followed by growth of Ba-rich sanidine the initial compositional profile across the boundary must have been a sharp step profile that was later smoothened by diffusion. Other processes, however, may have modified the initial step profiles such as Ba in-diffusion from surrounding Ba-rich melt. But, based on modeling by Rout & Wörner (2020), this can safely be excluded due to the high diffusivity of Ba in melt compared to sanidine (Zhang et al., 2010). Therefore, the assumption of an initial step-function in Ba concentration across the resorption boundaries satisfies one of the basic criterias and invites the application of diffusion modeling.
However, only compositional gradients in Ba can be modeled due to very high scatter in the x-ray count rates (>35% relative standard deviation) for Ca and Sr and lack of experimental constraints on Ca-diffusivity in sanidine. Zellmer & Clavero (2006) previously applied Sr-diffusion speedometry to a sanidine sample (erupted 14.1 ±1.4 ka) from Taápaca. However, the temperature they used (~875 ᵒC) was more than 100 ᵒC higher than the average storage temperature estimated from our new thermometric calculations as well as the sanidine-liquidus (770 ᵒC, based on R-melts modeling). Such higher temperatures only reflect the higher temperature crystallization of Mg-hastingsites in the hybrid dacite and the mafic enclaves but cannot be applied to diffusion modelling in the sanidine. The http://www.petrology.oupjournals.org/ Downloaded from https://academic.oup.com/petrology/advance-article/doi/10.1093/petrology/egab010/6126431 by guest on 09 February 2021 higher temperature of 875 ᵒC used by Zellmer & Clavero (2006), therefore, underestimates diffusion time-scales by a factor of at least ~100. A problem, however, remains regarding which temperatures to use from the 720°C to 820ᵒC temperature range observed here from the inclusion thermomertry in sanidine.

Temperature constraints for "non-isothermal" diffusion modeling
Temperatures obtained from plagioclase-amphibole-pairs vary significantly (by up to 80 o C) between different growth zones, even within a single crystal, irrespective of the thermometer applied. This clearly indicates non-isothermal conditions and hence calls for a model that takes the temperature variations into account. Therefore, we use a "non-isothermal diffusion modeling" approach to calculate diffusion times across each growth zone boundary using the temperature obtained from the mineral inclusions hosted in the following zone.
We note that the values estimated from different inclusion pairs within a single growth zone are consistent to within ±20 ᵒC for the majority of growth zones. However, in cases where we have a larger range of temperatures within a single zone, we use 0.95 times the highest temperature as the effective temperature as suggested by Chakraborty & Ganguly (1991). Their work has shown that for diffusion processes, there is a simple relationship between the peak temperature (T Pk ) attained in a variety of non-isothermal histories and the characteristic or effective temperature T Ch : T Ch = ~0.95 T Pk in Kelvin. The diffusion times based on modeling isothermal diffusion at T Ch is a good approximation of the actual diffusion time involving nonisothermal variation (Costa et al., 2008).
Our amphibole-plagioclase inclusion temperatures (720 -820 ᵒC) suggest that some of the mineral inclusions might have actually formed at, and thus record magma temperatures above the sanidine liquidus (770 ᵒC). This makes temperatures >770 ᵒC relevant for the http://www.petrology.oupjournals.org/ resorption event, but not for diffusion in the following growth zone. These higher temperatures are, however, very much applicable to diffusion of all interior (i.e. older) boundaries while the rim was melting. Therefore, when the temperature within a growth zone is <770 o C, we have to apply the respective temperature to all prior boundaries. For temperatures >770 o C, we apply this temperature to all but the immediately preceding boundary.

Diffusion coefficients
Diffusion coefficients are calculated from the diffusion temperatures (obtained as described above) using the Arrhenius equation, K)] is the universal gas constant, D 0 (m 2 /sec) is the pre-exponential factor and corresponds to the value of D (m 2 /sec) at infinite temperature and T is temperature in Kelvin.
We use diffusivity data for Ba from Cherniak (2002) who experimentally determined the diffusivity data for Ba in natural sanidine under dry 1 atm conditions. No compositional dependence of Ba diffusion was reported, and no significant diffusion anisotropy was observed between normal to (001) and normal to (010) directions. Activation energy of 455 ±20 kJ/mol and pre-exponential factor of 2.9×10 −1 m 2 /sec (Cherniak, 2002) are used. Although these parameters were determined for the temperature range from 775 to 1124 o C, we extrapolate the Arrhenius relation (Eq 1) to our lower temperatures of 720 to 820 o C. A dependence of Ba diffusion on pressure and fO 2 has not been detected (Cherniak, 2002(Cherniak, , 2010 and thus, these parameters are not essential in our diffusion modeling. http://www.petrology.oupjournals.org/

Non-isothermal diffusion modeling
The Non-Isothermal Diffusion Incremental Step (NIDIS) model that we use here to account for temperature variations during diffusion, was suggested by Petrone et al. (2016Petrone et al. ( , 2018. Rout et al. (2020) modified and simplified the original NIDIS algorithm, which significantly reduces the errors in diffusion time estimates while remaining mathematically unchanged compared to the original algorithm of Petrone et al. (2016). The model calculations start at the rim of a crystal and proceed chronologically backwards towards the core. Each of the calculation steps uses two consecutive diffusion boundaries to estimate the diffusion time.
The calculation steps are described below.
Calculation steps: Each resorption boundary is modeled for different temperatures obtained from amphibole and plagioclase inclusions, as explained above. A least square curve fitting (using MATLAB and ORIGIN) of the diffusion profiles uses the analytical solution for diffusion in a semi-infinite system (Crank, 1975) i.e. Eq 2 ( , ) = where C(x,t) is the Ba concentration at position x after time t (sec). C low and c high are the initial Ba concentrations (observed as the plateaus on both sides of the gradient), D is the temperaturedependent diffusion coefficient (m 2 /sec) and x 0 (m) is the position of the diffusive interface.
The curve fitting of a profile gives the corresponding Dt values.
The Dt value of a diffusion profile that has undergone different temperatures over time t is given as (Crank, 1975;Wilson, 1970): Eq 3 = ∫ 0 ( ) × http://www.petrology.oupjournals.org/ Thus, if the thermal history consists of several smaller isothermal steps, e.g. for temperature of T 1 for time t 1 , T 2 for t 2 and so on, up to T n for t n , then, Eq 4 = 1 1 + 2 2 + … + = ∑ = 1 Where D 1 , D 2 , ...D n are the diffusion coefficients at T 1 , T 2 , ...T n , respectively. Eq 4 means that the overall Dt value of a profile (obtained from curve fitting) will be the summation of individual Dt values of all the subsequent boundaries and the very boundary under consideration. However, note that the individual Dt value and diffusion time of a particular diffusion boundary that was formed after a resorption and re-growth event reflects only the period of diffusion until the formation of the next resorption boundary. Thus, in Eq 4, if Dt belongs to the first boundary that was formed i.e. (Dt) boundary 1 , which is the overall Dt value of boundary 1, then (Dt) boundary 1 = D 1 t 1 + (Dt) boundary 2 . Here, (Dt) boundary 2 = D 2 t 2 + ... + D n t n (from Eq 4). Upon generalization, the overall Dt value of n th boundary, starting from the core, turns out as: (Dt) boundary n = D n t n + (Dt) boundary n+1 . Thus, t n , i.e. the time between the formation of n th and n+1 th boundary can be estimated as: where D n is the diffusion coefficient calculated for the temperature obtained specifically for the n th zone (the zone between n th and n+1 th boundaries). Thus, the difference in the overall Dt values of two adjacent zone boundaries is divided by the diffusion coefficient appropriate specifically for the zone n in between. This way, the diffusion time reflects the time between the formation of the two successive boundaries. Please note that here and onwards, the "diffusion time of a profile or a boundary" only refers to the specific duration of diffusion across that specific boundary until the next boundary forms. Fig. 15 graphically explains the above steps of diffusion modeling using an example of a simple crystal with 3 non-isothermal growth http://www.petrology.oupjournals.org/ zones. Fig. 16 shows an example of non-isothermal diffusion modeling of a natural sanidine crystal based on the application presented in Fig. 15.
Using these steps, the diffusion times for individual profiles are calculated using the different temperature values from thermometry. We use Ba-profiles obtained by x-ray line analysis. For a few crystals, we also use grey scale profiles for the analysis as a proxy for Ba (after an apparent linear relationship was verified; see Fig. 17). The summation of the diffusion times from all boundaries gives the total diffusion time for the entire crystal. This time should be a close estimate of total crystal residence time between the initiation of growth and eruption because heating and resorption events are considered to be short compared to growth and diffusion. The only duration not accounted for in the total diffusion time is the duration of growth of the core before the formation of the first zone boundary.
In order to evaluate the role of different thermometric calibrations, we compared our modeling results based on temperatures using the amphibole-only thermometer by Ridolfi & Renzulli (2012) with diffusion times based on the Holland & Blundy (1994) thermometer. We found that the use of the temperatures from Holland & Blundy (1994) thermometer returns a maximum of 3 to 5 times longer diffusion times compared to the use of Ridolfi & Renzulli (2012) temperatures. Since we are here concerned with estimating the maximum times of diffusion, crystal residence and reactivation prior to eruption, we will use longer diffusion times based on Holland & Blundy (1994) in the following discussions.

Calculated diffusion times
Based on the parameters and constraints developed above, we now present the results for Ba diffusion modeling on major successive resorption interfaces from core to rim for 24 samples of sanidine megacrysts from Taápaca volcano. A compilation of all results of diffusion http://www.petrology.oupjournals.org/ modeling is given in Table 3 (with additional details and extended results in electronic appendix 7). The number of analyzed boundaries per crystal varies between 2 and 6. Some boundaries that do not terminate in compositional plateaus or do not have suitable mineral inclusions in the following growth zone, are not analyzed for diffusion. In such a case, the modeled time of the immediately interior (older) boundary covers that of the following (younger) excluded boundary as well. This is done by using the effective temperature (Chakraborty & Ganguly, 1991, described above) calculated for the preceding (older) as well as all the successive (younger) growth zones with excluded boundaries.
The total residence times obtained from the diffusion modeling for different crystals varies significantly within samples of a particular dome and also between samples from different domes of different stratigraphic stages of Taápaca's eruptive history. For stratigraphic stage II, the total accumulated diffusion times in sanidine megacrysts range between 9 and 499 ky while for stage IV, we estimate between 29 and 396 ky. Using the differential diffusion time and radial distance between the 1st and the last modelled zone boundary, we extract average apparent growth rates of ~16 cm/My or ~5x10 -13 cm/s (with the range being from ~3 to ~32 cm/My or ~10 -13 to ~10 -12 cm/s). The real growth rate, however, must be higher as these crystals suffered frequent events of resorption that indicate interrupting melting periods of unknown duration. The calculated total diffusion times, however, are unaffected by intermittent times of resorption. This is because while the crystal melted on the outside, the internal zone boundaries in the crystal still continue to diffuse and, as described earlier, such scenarios are taken into account in the diffusion modeling. Since we also did not find any correlation between crystal size and residence time, different crystals must have grown at different points of time and within different parts of the crystal-rich rhyodacitic magma reservoir. As a result, different sets of crystals experienced different versions and extents of the thermal history which overall must have followed similar patterns in temperature. Some crystals experienced longer and/or hotter heating events compared to others.
Our estimated residence times from diffusion modeling naturally excludes the time of magma evolution prior to the onset of sanidine growth and also the time between sanidine nucleation and the formation of the first analysed zone boundary (the core-mantle boundary of the crystals). However, the duration of growth of the core until the first analysed diffusion boundary can be approximated from growth rate extracted from the rest of the crystal: Depending on the size of the core this amounts to 5 -30 ky for initial core growth, which does not significantly change the order of magnitude of the total diffusion time for most of the crystals significantly and thus, does not affect our interpretations.
The total residence times of sanidine crystals from eruptive stages II and IV (Table 3) vary significantly. Based on the ages of the different stages, "younger" (29 -77 ky) crystals occur prefereably in younger domes and therefore crystal residence times appear to have decreased over time as the magmatic system matured. By contrast, earlier in the volcano's history, crystals were stored over longer time spans (up to 499 ky). Only sample 11 K1 and 29 B2 from the most recent eruptive stage IV are exceptions with at least 380 ky of storage. This suggests that some relatively older crystals may have been mobilized at a later (younger) stage as well. In any case, crystals were formed and stored for thousands of years in the resident rhyodacitic magma throughout the entire volcanic history of Taápaca volcano, i.e. for more than ~1.5 My.

Mobilization of the sanidine-bearing rhyodacite prior to eruption
http://www.petrology.oupjournals.org/ Diffusion times for individual resorption boundaries, i.e. the duration between consecutive heating events, varies from 1 to ~50 ky (electronic appendix 7). The boundaries in the crystal cores and mantles give longer individual diffusion times (i.e. time until the formation of the next boundary) typically with >30 -50 ky between two successive boundaries. On the other hand, diffusion times for the outermost two to four boundaries closest to the crystal rims are shorter with only 1 to 10 ky (electronic appendix 7). The radial spacing of growth and resorption bands also decreases towards the rims (Figs. 6 -8). This indicates to us that the time intervals between resorption events become shorter towards the rims and implies higher resorption frequencies with time. This increased frequency of heating, and therefore recharge events, suggests that the resident rhyodacite magma reservoir that contained the sanidine megacrysts became more frequently thermally disturbed by recharge about ~3 -11 ky prior to its eventual eruption. This time-scale is obtained by adding the individual diffusion times (1 -10 ky) of the outer, narrowly spaced resorption boundaries. We define this period as the "pre-eruptive mobilization time" (Table 3, electronic appendix 7). The outermost major resorption boundaries give even shorter diffusion times of only 0.4 -3 ky. This implies that the last major heating (and recharge) event that caused the destabilization and the eruption of the sanidine as part of its hybrid host must have happened within this time frame. We note that at the final stage of growth, some sanidine crystals formed a thin overgrowth of only ~1 to 1.5 mm with still narrower growth and resorption bands . Diffusion modeling in this growth zone was, however, impossible due to the lack of temperature constraints and appropriate compositional profiles. This final shift in zonation pattern to even higher frequency must have occurred in less than 0.4 -3 ky and is probably directly linked to the final magma mixing that also formed the hybrid dacite, entrained the sandine megacrysts and led to their eruption to the surface. To further constrain the mixing timescales, we calculated a possible range of temperatures for the mixed product and applied it to model the outermost modellable zone boundary. Using (1) the average temperatures for the rhyodacite (760 ᵒC, Fig. 11) and the basaltic andesite (1015 ᵒC, Fig. 11), and (2)

Cold versus warm storage
Two contrasting models have been proposed to describe magma storage and temperature histories of magmas prior to their eruption. Barboni et al. (2016) and Kaiser et al. (2017) suggest long, tens to hundreds of thousands of years, storage of crystals in a high-temperature meltrich, above-solidus magma reservoir ("warm storage"). Such warm storage implies that magmas remain in a (near-) eruptible state for a long time. Alternatively, Cooper & Kent (2014) and Reid (2008) summarized evidence for near-solidus low-temperature storage ("cold storage", <700 ᵒC) of silicic magmas for thousands to hundreds of thousands of years. Remelting and mobilization of such crystalline mush (>40 -50% crystallinity) is thought to be a common process, driven by mafic recharge, to form, assemble, and erupt silicic magmas (Bachmann & Bergantz, 2008). Regarding the duration of storage of the resident megacrystic rhyodacite, our http://www.petrology.oupjournals.org/ calculations indicate thousands of years of crystal growth, storage and temperature cycling, but we will argue below that our results do not fall clearly into any of these two categories.
Ba diffusion essentially stops for temperatures below the diffusion threshold and thus, 'true' residence times can be significantly longer than those extracted from diffusion modeling.
In such a case, diffusion times will be substantially shorter than the actual crystal "ages" that are recorded in U-Th isochron ages (e.g. Cooper &Kent, 2014 andRubin et al., 2017). We estimate the diffusion threshold for Ba (for details see electronic appendix 8) to be less than 700 ᵒC. Since the amphibole-plagioclase thermometry (Holland & Blundy, 1994) and amphibole thermometry (Ridolfi & Renzulli, 2012) on our samples consistently gave temperatures above 720 ᵒC and the crystals record a frequent change from growth to resorption, we argue that these crystals were never below the diffusion threshold for any significant amount of time. Therefore, our diffusion modeling at these temperatures should record similar residence time scales of tens to hundreds of thousands of years comparable to U-Th isotope isochron ages. This is indeed observed: For a sanidine-bearing dacite (locality TAP 002) with an 40 Ar-39 Ar age of 13 ±8 ka , a pre-eruptive crystal residence time of 158 ±27 ky results from a measured 238 U/ 230 Th mineral isochron age of 171 ±26 ky (Kiebala, 2008). This "crystal age" agrees well with the average diffusion time of 131 ky determined for different sanidine megacrystss from the same dacite deposit (sample 29) at a temperature range of 720 -820 °C.
These consistent results are (i) a further argument to use the longer diffusion times based on the Holland & Blundy (1994) thermometer calibration, and (ii) support our conclusion about the persistence of relatively high storage temperatures. We therefore conclude that the megacrysts never cooled below the diffusion threshold, i.e. ~700 ᵒC for any significant amount of time.
Results from R-MELTS modeling predict the rheological lock-up for the rhyodacite at 750 -780 ᵒC and a solidus at ~680 ᵒC, which is close to the diffusion threshold. Temperatures of 720 -820 ᵒC estimated from amphibole-plagioclase inclusions suggest that the crystals grew http://www.petrology.oupjournals.org/ at at least 40 ᵒC above the magma solidus and were frequently heated by up to 140 ᵒC above the magma solidus (to temperatures above the sanidine-liquidus, i.e. >770 ᵒC, as shown by the abundant resorptions). However, given the relations between crystallinity and temperature, the sanidine crystals would still have been stored in a rigid mush (with ≥40 vol% crystallinity at T <750 -780 ᵒC) for a significant amount of the time. This would make the Taápaca reservoir a distinctive case that frequently transitioned between "cold" (non-eruptible) and "warm" (eruptible) storage (Fig. 18). We define such storage as "Long-term Transitional Temperature Cycling" or LTTC storage.

What triggered the eruptions?
Based on our interpretations and the common occurrence of basaltic andesite magmas (seen as mafic enclaves) throughout the history of Taápaca, we argue that repeated heating events are related to distant and/or weaker recharge of hotter mafic or basaltic andesite magma into a resident colder rhyodacitic crystal-melt mush without resulting in chemical mixing in the magma batch that holds a particular set of sanidine megacrysts. Such earlier recharges obviously did not result in triggering an eruption of that particular magma batch. Similar observations in andesitic systems are consistent with these observations: up to ten resorption events can be recorded during the growth of even small plagioclase phenocrysts in andesites before eruption occurs (Ruprecht & Wörner, 2007). Thus, at Taápaca volcano, after many distant recharge events, only a "critical" (and final series of) recharge event(s) apparently was sufficiently strong and/or sufficiently close to the particular portion of the crystal mush that it eventually did mix and mobilize the mush into an eruptible hybrid dacite magma (Figs. 18 and   19). This scenario was repeated for every set of megacysts over >1.5 My of eruption history.
These "distant" earlier recharge events may very well have mobilized some other portions of http://www.petrology.oupjournals.org/ the rhyodacite reservoir (and their respective crystals/megacrysts) resulting in the eruption of other magma batches at these moments in time. Weaker recharge events may not trigger an eruption at all, but contribute to the overall heat budget of the system to keep it thermally "alive".
The relative uniformity of dacite compositions and their uniform petrography, mineral assemblage, and crystal content throughout the history of Taápaca reflects a typical "breakthrough" threshold where critical mixing occurs. This repeatedly and consequently results in mass proportions of mafic recharge and rhyodacite endmembers that are always similar so that relatively uniform hybrid dacite compositions can reach the surface (Fig. 19).

SUMMARY AND CONCLUSION
We present a model (Fig. 18) for the origin and eruption of sanidine megacrysts at Taápaca volcano. The host dacite formed by mixing between a resident, sanidine-bearing rhyodacite and mafic recharge. We take the Ba-content variation of zoned sanidine and the frequent resorption interfaces as proxy for the thermal variation with time in the rhyodacite magma prior to the formation of the hybrid dacite before every eruption (Fig. 14). The model starts at the onset of crystallization of sanidine and continues for up to ~500 ky of crystal growth and storage. After the formation of the crystal core, the temperature is cycled between near eutectic composition and above sanidine-liquidus, eventually resulting in megacrystic growth (Higgins 2011). Different temperatures and pressures from thermo-barometry, different total diffusion times, and the non-linear correlation between residence time and crystal size suggest that crystals grew in different parts/regimes of the reservoir. Sanidine crystals in distinct magma batches underwent similar histories but on different time scales, at different points in time and locations. This relatively colder period of storage finally ended about 3,000 -11,000 years prior to each eruption when crystals in a particular magma batch slowly started to heat up from more http://www.petrology.oupjournals.org/ frequent recharge events (at intervals of 1,000 -10,000 years). A final pre-eruptive phase of magma mixing occurred from 400 to 3,000 years before each eruption when the recharge frequency significantly increased, leading to a further temperature increase. At this phase, however, the chemical effects of recharge are still not recorded in the sanidine composition.
The fine-scale high frequency zonation at the outermost rims of the crystals marks the arrival of the critical basaltic andesite recharge in a particular part of the reservoir. This finally resulted in overturn, mingling and mixing, formation of the hybrid host dacite, and the onset of an eruption (Fig. 19). This critical stage must have happened after the "heat-up phase" (400 -3,000 years), i.e. within the last few centuries or even decades before the eruption. Such a sequence of events was repeated many times during the history of Taápaca volcano. Each eruption mobilized some parts of the mush reservoir while other parts remained to be erupted at a later stage. Likely, some larger parts of the crystal mush remain(ed) at depth, will not erupt at all and thus, will eventually form a megacrystic granite intrusive rock. The surprising similarities of size, shape, Ba-zonation, growth texture, inclusion trails and host rock composition of sanidines erupted at Taápaca as compared to the megacrysts in granitoids are a strong indication that similar growth processes, duration of crystallization, and temperature cycling also apply to these intrusive rocks. Indeed, models of magma evolution and amalgamation of magma batches over similar time scales result from zircon dating in other volcanic and intrusive rocks (e.g. Coleman et al., 2004;Wotzlaw et al., 2013;Szymanowski et al., 2017;Ratschbacher et al., 2018;Karakas et al., 2019) The model for Taápaca volcano developed here suggests that even relatively small silicic volcanoes can store material well above the solidus at shallow depth for hundreds of thousands of years and can be active in a monotonous style for about 1.5 My. Magmas frequently transition from non-eruptible to eruptiple state above rheological lock-up (c. 40% crystals). This cyclic behaviour between the "cold" and "warm" storage conditions over long http://www.petrology.oupjournals.org/ periods of time, i.e. "Long-term Transitional Temperature Cycling" or LTTC storage, represents a magmatic regime where temperature, crystallinity and eruptability will depend on the frequency and extent of recharge events. This delicate balance, as observed in Taapaca, sets the perfect conditions for the sanidine phenocrysts to develop megracrystic growth that records the complex and long thermal history with frequent resorption and growth cycles. The duration of hundreds to thousands of years between the onset of mobilization and eruption suggests a rather "low and slow" mobilization of magma (Kaiser et al., 2017) in contrast to rapid mobilization in other cases (Druitt et al., 2012;Cooper & Kent, 2014). Long-term temperature cycling in a shallow magma reservoir and K-feldspar megacrystic growth under these conditions are rarely documented in lavas erupted to the surface. However, analogy to and abundance of similar growth patterns observed in megacrystic granites indicate that such conditions could be more common than previously thought in shallow subvolcanic magma systems, but are under-represented in eruptive products. Their reactivation into eruptible magma at Taápaca volcano is a rare case that requires a delicate thermal and compositional balance in a long-lived and stable magmatic system with ephemeral eruptions of uniform magma composition.   The center-cut of this crystal has an area of 58 cm 2 that was mapped for all major elements in its entirety at high resolution (Electronic appendix 2). The errors are within the sizes of the symbols. is a simple 4-stage fractional crystallization model which connects the hybrid mafic magma feeding into the Parinacota magma system to the high-Si rhyolitic compositions from which sanidine has crystallized. It is clear that no natural erupted compositions exist that fall along such a fractionation trend. The orange line (model 2 with two stages 2-1 and 2-2) represents a fractional crystallisation model starting with the least evolved silicic dome lava of Parinacota.
At the end of model 1-1, 45% crystallization occurs with 1% olivine, 8% pyroxene, 32% plagioclase and 3.8% amphibole. At the end of 2-2, additional 22.5% plagioclase, 3% amphibole and 4.5% sanidine have crystallized taking the crystallization to 75%. Similarly, in model 2-1, 12% crystallization occurs with 1.2% pyroxene, 8.4% plagioclase and 2.4% amphibole. By the end of 2-2, total crystallization stands at 85% with 1.6% pyroxene, 53% http://www.petrology.oupjournals.org/ plagioclase, 11% amphibole and 9.3% sanidine. Partition coefficients for Ba and Sr in sanidine D Ba = 25 and D Sr = 3.5 in the model are within the (lower) range given in the literature (e.g. GERM Partition Coefficient (Kd) Database, earthref.org, and references therein). The average uncertainties on the data points are presented by the representative error bars on the Parinacota silicic dome lava matrix (orange hollow circle). The bimodality is evident for amphiboles of different composition in the hybrid dacite. All the temperature data are from amphibole-plagioclase thermometry (Holland & Blundy, 1994). The pressure data for all the Mg-hornblendes are from Al-in amphibole barometry by Mutch et al. (2016). The pressures for the Mg-hastingsites are determined using the barometer by Ridolfi & Renzulli (2012). that does not mix or trigger an eruption but contributes to the overall heat budget of the system to keep it "alive". 2: the "critical" recharge that overcomes the "break-through" threshold, mixes with the host rhyodacite to form the hybrid dacite, which eventually erupts and carries large sanidine crystals to the surface. These "critical" (2) or weaker (1) recharge events of one regime act as the distant events for other regimes which may experience only the thermal effects. 3 and 4: previous or "fossil" recharge events similar to 1 and 2 contributing to the heat budget and forming older domes, respectively.  Table 1. Oxygen isotope data for the sanidine megacrysts and the dacites. All measurements are done using CO 2 laser-assisted irmGCMS (except for 28/1 which is measured using in-situ UV laser ablation).