Insights Into Magma Storage Beneath a Frequently Erupting Arc Volcano (Villarrica, Chile) From Unsupervised Machine Learning Analysis of Mineral Compositions

A key method to investigate magma dynamics is the analysis of the crystal cargoes carried by erupted magmas. These cargoes may comprise crystals that crystallize in different parts of the magmatic system (throughout the crust) and/or at different times. While an individual eruption likely provides a partial view of the subvolcanic plumbing system, compiling data from multiple eruptions can build a picture of a whole magmatic system. In this study, we use machine learning techniques to analyze a large (>2,000) compilation of mineral compositions from a highly active arc volcano: Villarrica, Chile. Villarrica's postglacial eruptive activity (14 ka–present) displays large variation in eruptive style (mafic ignimbrites to Hawaiian style effusive eruptions), yet its eruptive products have a near constant basalt‐basaltic andesite bulk‐rock composition. What therefore, is driving explosive eruptions at Villarrica and can differences in storage dynamics be related to eruptive style? Here, we use a hierarchical cluster analysis to detect previously unseen structure in the composition of olivine, plagioclase and clinopyroxene crystals erupted at Villarrica, revealing the presence of compositionally distinct clusters within each crystal population. Using rhyolite‐MELTS thermodynamic modeling, we related these clusters to intensive magmatic variables: temperature, pressure, water content, and oxygen fugacity. Our results provide evidence for the existence of multiple discrete (spatial and temporal) magma reservoirs beneath Villarrica where melts differentiate and mix with incoming more primitive magma. The compositional diversity within an erupted crystal cargo strongly correlates with eruptive intensity, and we postulate that mixing between primitive and differentiated magma drives explosive activity at Villarrica.

Another recent advance in understanding magmatic systems is that crystal cargoes can be rapidly assembled from different parts of a magmatic system after protracted storage (e.g., Bergantz et al., 2015;Cooper & Kent, 2014;Mutch et al., 2019). Together with the TCMS model, they explain how a mineral assemblage erupted during a single event may contain crystals formed from the carrier melt (autocrysts), crystals remobilized from other parts of the magmatic system (antecrysts), those from outside the magmatic system (xenocrysts), and crystals that form due to undercooling upon eruption (microlites; Jerram & Martin, 2008). Therefore, the crystal cargo of an eruptive deposit can be thought of as a snapshot of the underlying magmatic system. While an individual eruption likely provides a partial view of the subvolcanic plumbing, compiling data from multiple eruptions can be used to build up a picture of the whole system. Building this complete picture requires well-characterized magmatic products from as many closely spaced eruptions as possible.
Here, we use established unsupervised machine learning techniques to analyze a large (>2,000) compilation of mineral compositions from a highly active arc volcano: Villarrica, Southern Andes, Chile. We reveal previously unidentified structure in erupted mineral compositions and identify trends throughout Villarrica's eruptive history related to crystal zoning and eruptive style. We utilize thermodynamic modeling to constrain the physical and chemical characteristics of Villarrica's magmatic system to assess the suitability of the TCMS model. Finally we discuss the role of magma mixing in driving explosive behavior at Villarrica volcano and the implications for future eruptions.

Unsupervised Machine Learning Applied to Geochemical Data
The term "unsupervised machine learning" describes a class of algorithms that are designed to find patterns in multidimensional data. They are termed "unsupervised" as they do not require any prior knowledge of the relationships between data. Unsupervised machine learning techniques can be used to gain insight into the structure of multivariate compositional data (data that describe quantity relative to a whole, i.e., close to a constant value such as 100%), which is common in geochemistry (e.g., Chiasera & Cortés, 2011;X. Liu et al., 2020;Templ et al., 2008). The main advantage of these methods is that they allow the data to be interrogated in a multivariate sense and can highlight trends that otherwise would be difficult to identify using traditional Harker-style bivariate plots (Cortés, 2009). Cluster analysis is an unsupervised machine learning technique that attempts to group data by some measure of similarity. The most similar data points are iteratively combined, until all data points belong to a single cluster. Hierarchical cluster analysis describes the relationship of all the data points to each other during this process, producing a hierarchy of data clusters organized by their similarity.
The majority of past studies that utilize cluster analysis have used it to characterize the composition of volcanic products via the analysis of whole-rock data sets, both for individual volcanoes (e.g., Mt Etna, Italy [Corsaro et al., 2013]; Izu-Oshima, Japan [Kuritani et al., 2018]) and regional volcanism (e.g., the Virunga Volcanic Province [Barette et al., 2017]). Others have applied the cluster analysis to individual phases from a single volcano for example, silicate melt inclusions (Hamada et al., 2020); tephra glass (E. J. ; and minerals (Caricchi et al., 2020;Cortés et al., 2007;Gleeson et al., 2021). Some, but not all, of these studies have recognized that compositional data require special mathematical treatment prior to analysis. In this study we use a log-ratio approach (i.e., Aitchison, 1986) to transform the compositional data into a compatible Euclidean geometry.

Geological Background
Villarrica (39.5°S, and 71.9°W) is a Quaternary stratovolcano in the Andean volcanic arc, and part of the Chilean Central Southern Volcanic Zone (Figure 1a). The volcano is Chile's most active, with more than 50 recorded eruptions since 1558 (Petit-Breuilh, 2004), and is considered by the Chilean geological survey (SERNAGEOMIN), as the most hazardous of the 92 geologically active Chilean volcanoes (SERNAGEOMIN, 2019). Villarrica's postglacial (14 ka-present) eruptive activity displays a wide range in eruptive intensity and magnitude. This includes two major eruptive events, that generated the Licán (ca. 13.9 ka BP) and Pucón (ca. 3.7 ka BP) mafic ignimbrites, with estimated volumes of 10 and 5 km 3 (non-DRE), respectively (Lohmar et al., 2012;Silva Parejas et al., 2010). In contrast, historic eruptions (1900-present) have ranged from Hawaiian to violent Strombolian and are dominated by effusive lava flows (Pizarro et al., 2019). However, the paroxysmal March 2015 eruption, which lasted just 30 min, was characterized by a 1.5 km high fire fountain (Romero et al., 2018) and is the most recent demonstration of Villarrica's explosive potential. Despite this variety in postglacial eruptive style, Villarrica's volcanic 10.1029/2022GC010333 3 of 24 products have a limited compositional range: 98% of the juvenile whole-rock compositions collated in this study have 52-57 wt% SiO 2 ( Figure 1). However, whole-rock data may not fully reflect the compositional variety of a magmatic system as heterogeneity may be present at a scale smaller than that of the whole rock sample (Pichavant et al., 2007). Therefore in this study, we focus on the main mineral phases erupted at Villarrica, whose compositions give insight into magma dynamics and the physical conditions of the magmatic system.
To ensure that only phases relevant to the magmatic system were considered, any analyses labeled as xenolith or microlite by the original authors were removed from the data sets. The analyses were then manually screened for errors for example, misclassified, mixed-phase and poor-quality analyses. Any analysis without an analytical total between 98 wt% and 102 wt% was removed. The cations per formula unit (cfu) were calculated for each mineral analysis based on stoichiometry, using 4, 32, and 6 oxygens for olivine, plagioclase, and clinopyroxene, respectively. Total iron was assumed to be entirely FeO for plagioclase and olivine analyses. The method of Droop (1987) was used to calculate the proportion of Fe 2+ to Fe 3+ in clinopyroxene analyses. Clinopyroxene analyses with a total cfu outside 4.00 ± 0.02 were removed. After screening, 2,267 analyses (out of an initial 2,611) were deemed suitable, and these are broken down by eruption in Table 1.

Data Quality and Limitations
The database used in this study contains a large (>2,200) total number of analyses, but the number of analyses per eruption is much smaller: only the 2015 eruption has nearly 100 analyses for either of the two most modally abundant minerals, olivine and plagioclase (Table 1). This has implications for relating compositions to textural features for example, crystal zoning. Cheng et al. (2017) suggests that 100 or more analyses are required to characterize the compositional and textural features of a complexly zoned mineral population. However, we find no relationship between zoning and compositional clustering.
All eruptive deposits have not been sampled and analyzed equally. Table 1 shows that there are a higher number of analyses of high-volume eruptions (the Licán and Pucón ignimbrites), historic eruptions, and those with easily accessible deposits. The main route up Villarrica, the Pucón Ski Center Road, cuts both the Pucón ignimbrite and the Chaimilla Fall Deposit and provides easy access to historic eruptions. Conversely, the Los Nevados and Chaillupén parasitic cones, and the Dacitic Domes are further from established roads and tracks. While this bias prevents us from commenting on individual eruptions that are poorly studied, a first order understanding of Villarrica's magmatic system is still possible.
Finally, there are likely systematic errors related to the different analytical equipment in different labs, which are in turn calibrated with different standards. This results in different analytical uncertainties for each of the past studies complicating direct comparison. There is also little overlap between the eruptions sampled by different studies. Without this, there is no way to quantify these potential inter-lab biases. Despite this, over half the data was obtained from a single lab, which reduces the chance of this having an effect (Lohmar, 2008). Furthermore, as discussed later in the paper, we find no systematic bias in the clustering that results from the source of the data.

Compositional Data Constraints and the Log-Ratio Transformation
Compositional data carry only relative information, subject to nonnegative and constant-sum constraints (100 wt% etc.). These constraints mean that often used statistical methods that assume that data are normally distributed (i.e., unconstrained to a constant value and varying from −∞ to +∞), are not directly applicable. In his seminal work, Chayes (1971) established that it is the ratios between compositional data that measure variability rather than absolute differences. To circumvent this problem, Aitchison (1986) introduced several log-ratio transformations: linear transformations that map compositional data into an unconstrained Euclidean space.
There are several proposed log-ratio transformations: the additive log-ratio transformation (Aitchison, 1986), the centered log-ratio transformation (Aitchison, 1986), and the isometric log-ratio transformation (Egozcue et al., 2003). In this study, the mineral data sets were transformed using the isometric log-ratio (ilr) transformation (Equation 1), implemented using the Pyrolite python library (M. Williams et al., 2020 Moreno and Clavero (2006). b Moreno (1993). c Clavero-Ribes (1996). d Lara and Clavero (2004  was chosen as it ensures the transformed data have a nonsingular covariance matrix and preserves the geometric properties of the raw compositional data (Egozcue et al., 2003): where, x is a compositional analysis, i is a specific part, D is the number of parts (elements analyzed), and g(x i ) is the geometric mean of the parts of x: A requirement of using any log-ratio transformation is that the data cannot contain zeros (Cortés et al., 2007;Fry et al., 2000). Zeros in compositional data can be structural (e.g., K in clinopyroxene), below detection of the analytical technique, or simply not analyzed (Fry et al., 2000). We removed structural zeros by only considering elements that reasonably exist in a mineral's structure. To deal with zeros resulting from detection limits, a detection limit of 0.05 wt% was assumed for all elements. The detection limits were generalized because detection limits were not reported in all cases. Zeros resulting from detection limits were replaced using the Multiplicative Replacement method of Martín-Fernández et al. (2000) and Fry et al. (2000), which is equivalent to distributing a detection limit threshold evenly among the below detection zeros. Not all elements (especially minor elements) were analyzed in every study. Elements that were measured for less than half of each of the mineral data sets were not used in our analysis. A table containing the elements used in the cluster analysis for each of the mineral data sets can be found in Supporting Information S1.
Variables with low abundances but high relative variances (often minor and trace elements) have high log-ratio variances and therefore dominate any analysis of the complete log-ratio transformed data set (Baxter et al., 2005;Greenacre, 2019;Greenacre & Lewi, 2009). To prevent this, we normalized the log-transformed data set using the column (part) medians and standard deviations: Where X′ is the normalized data set, x ij the ith part of the jth analysis, is the column median, and σ i is the column standard deviation. The median was chosen over the arithmetic mean as the transformed data sets were nonnormal: all four ilr-transformed data sets failed the Henze-Zirkler multivariate normality test with a specificity of 0.05 (Henze & Zirkler, 1990). Relationships between the transformed compositional data were then explored using hierarchical clustering methods over a Euclidean space.

Hierarchical Cluster Analysis
Hierarchical clustering algorithms attempt to group data into clusters by some measure of similarity. Agglomerative clustering methods start by grouping the two most similar data points into a cluster, then treating them as a single data point. The most similar data points or clusters are then iteratively combined until only a single cluster containing all the data remains. To cluster the data, a measure of dissimilarity (often called distance) must be chosen. A popular measure is the Euclidean distance which is the equivalent of Pythagoras's Theorem but over more than two dimensions: Where d E is the Euclidean distance matrix, i is the number of parts, and x and y are the points considered. Next, a linkage method must be chosen, that is, a method describing how the distances between points and clusters are used to relate them. We used ward clustering (Ward, 1963), which minimizes the in-cluster variance, over alternatives such as single-linkage or average-linkage as it does not suffer from chaining (W. T. Williams & Lambert, 1966;Wishart, 1969). The hierarchical clustering algorithm was implemented using the scikit-learn python library (Pedregosa et al., 2011). Dendrograms depicting the resultant hierarchy were used to determine a suitable number of clusters independently for each mineral data set. Choosing the number of clusters that best represent the data set is a balance between showing global versus local variability: a small number of clusters will highlight the largest differences in the data at the cost of potentially obscuring more subtle differences, and vice versa for a large number of clusters. Figure 2 shows that the olivine and plagioclase data sets are well described by two clusters until almost half of their respective maximum distances. However, to maximize our ability to detect more subtle compositional changes between clusters, we chose the next lowest number of clusters that well describe them. This was chosen by finding the maximum separation (second highest for olivine and plagioclase) between branches on the dendrogram. This results in three, five and four clusters for olivine, plagioclase and clinopyroxene, respectively.
We assessed the robustness of the identified clusters by repeatedly (1000x) resampling half of each mineral's data set and performing hierarchical clustering. The clusters from subsampled data were consistent with those identified in the complete data set, demonstrating that they are not strongly dependent on the size of our complete database. Further details can be found in Supporting Information S1.

Multivariate Outlier Detection
Each cluster contained points normally distributed about its center; therefore outliers are, according to the empirical rule for normal distributions, those values located at distances larger than three standard deviation from its center. To identify potential outliers, the Mahalanobis distance, that is, the distance of a given data point x and a distribution (Mahalanobis, 1936), was calculated for each cluster: where D M is the Mahalanobis distance, is a matrix of data points in each cluster, and and C are the location and covariance estimators. The location and covariance estimators were robustly calculated using the Minimum Covariance Determinant estimator (MCD) (as in Filzmoser and Hron [2008]), implemented using the FastMCD algorithm (Rousseeuw & Driessen, 1999) available in the scikit-learn python package (Pedregosa et al., 2011). The Mahalanobis distances for each cluster can be approximated by a χ 2 distribution (Rousseeuw & Zomeren, 1990). A critical Mahalanobis distance was determined for each cluster, above which a point is considered an outlier: which is the square root of the upper-α quartile of the χ 2 distribution with p degrees of freedom, which were 5, 7, and 10 for olivine, plagioclase and clinopyroxene, respectively. The typical choice for (1-α) 0.975 (Rousseeuw & Zomeren, 1990) was used that is, the outliers will be contained in the upper 2.5% of the χ 2 distribution. Any identified outliers were removed from the data sets.

Estimating Cluster Centers
To characterize each of the identified clusters, a measure of central tendency is required. For compositional data, the best unbiased estimator of the expected value is defined by the closed geometric mean (Cortés et al., 2007, and references within). To ensure that the center of the cluster is a valid mineral composition, the closed geometric mean of each cluster was calculated, and the nearest (according to a Euclidean distance matrix) data point to each mean was used to represent each cluster.

Thermodynamic Modeling
The current consensus, especially in arc settings, is that the majority of erupted minerals are remobilized antecrysts from a mushy magmatic system . This means it is unlikely an erupted mineral is in equilibrium with its carrier liquid, making traditional mineral-liquid thermobarometry impossible. For the same reason, mineral-mineral thermobarometry (e.g., olivine-augite) may not be accurate, even if both minerals appear to be in chemical equilibrium, unless textural evidence also suggests they are in equilibrium. To circumvent these issues, we compare measured mineral compositions to many thermodynamic simulations of a simple Villarrica-like system, performed at a range of intensive variables. This avoids the need for mineral-liquid and mineral-mineral equilibrium.
Thermodynamic modeling was carried out to investigate the likely provenance of the identified mineral clusters with respect to the physical and chemical conditions of magma evolution at Villarrica. We used the rhyolite-MELTS v1.2.0. algorithm (Ghiorso & Gualda, 2015;Gualda et al., 2012), via the alphaMELTS 2.0 front end (originally Adiabat_1ph [Smith & Asimow, 2005]), to model crystal fractionation from a primitive input melt. The initial bulk composition of each model was estimated using the most primitive glass composition reported at Villarrica: a postentrapment-crystallization corrected, olivine-hosted melt inclusion (VL15B-ol5-inc1, Mg# 65 at QFM+1) from an upper unit of the Chaimilla Fall Deposit (Pioli et al., 2015). Isobaric fractional crystallization models were performed with temperature decreasing in 1°C increments from above the liquidus, at a range of pressures (25, 50, 75 and 100 MPa and then increments of 100 up to 700 MPa), variable initial water contents (0.0-6.0 wt% in increments of 0.5), and oxygen fugacity buffers (fO 2 ) QFM-0.5 (0.5 units below the Quartz-Fayalite-Magnetite buffer) to QFM+2.25, in increments of 0.25. A total of 1,560 rhyolite-MELTS simulations were performed, one for each pressure, H 2 O and fO 2 permutation.
We identified the best-fit intensive variables by comparing the composition of minerals crystallized in each rhyolite-MELTS simulation to those in each identified cluster. We calculated the cfu, and then ilr-transformed the composition of olivine, plagioclase and clinopyroxene crystallized at each temperature step for all 1,560 rhyolite-MELTS simulations by the same procedure as for the measured compositions. As detection limits are not defined sensu stricto in MELTS, any zeros in the simulated compositions were replaced with 0.0001 cfu. With both the simulated and measured compositions in the same Euclidean space, we used the Euclidean distance between each as a measure of similarity (Equation 4). The best fitting intensive variables were those that produce the "closest" compositions to those measured: those that have the smallest distance.
To back out the best-fit intensive variables for each cluster, we calculated the weighted arithmetic mean and weighted standard deviation of the best-fitting T, P, fO 2 , and H 2 O from all measured compositions that comprise the cluster. The weights used were the inverse of the Euclidean distance: the measured compositions that were best reproduced by a simulation are given a higher weight. The weighted average and standard deviation for each mineral cluster are given in Table 3 and plotted in Figure 7.
We used a Monte-Carlo approach to estimate the effects of the analytical uncertainty in the initial bulk-composition inputted into the Rhyolite-MELTS simulations and the analytical uncertainty in the clustered mineral compositions. The maximum resultant uncertainties are comparable to those of mineral-liquid or mineral-mineral thermobarometric methods: ca. ±50°C, ±200 MPa, ±0.5 log units, and ±1.5 wt% H 2 O. This is to be expected as both MELTS and thermobarometers are calibrated on experimental data sets from the Library of Experimental Phase Relations. Further details can be found in Supporting Information S1.

Results
The compositions of minerals erupted in Villarrica's crystal cargoes are diverse, especially in comparison to whole-rock compositions (Figure 1c). Olivine compositions range from Fo 50-87 , plagioclase An 50-96 , and clinopyroxene Mg 57-92 . Despite this range in mineral compositions, hierarchical cluster analysis and outlier detection were performed successfully on the three mineral data sets. We identified 4, 94 and 15 outliers for olivine, plagioclase and clinopyroxene, which corresponds to 0.6%, 8.4%, and 3.5% of each respective data set. The representative compositions (centers) of each identified mineral cluster are shown in Table 2. The clusters are ordered from most primitive to most evolved (e.g., Ol 1 is the most primitive olivine cluster). Violin and box-plot diagrams showing the compositional differences for each mineral's clusters are included in Supporting Information S1. We compared the cluster designations with each composition's source study and found no obvious dependence. Cluster designations appear to be purely compositional.
To visualize the distribution of the identified clusters and verify that they were compositionally distinct, two techniques were chosen to reduce the dimensionality of the ilr-transformed data. First, a principal component analysis (PCA) was performed and the resulting first and second principal components plotted (Figures 3a-3c). Then a t-distributed Stochastic Neighbor Embedding (t-SNE) algorithm (Van der Maaten & Hinton, 2008) was implemented with a perplexity value of 40. 10,000 iterations were used to ensure that the projection was stable (Figures 3d-3f). Both the PCA and t-SNE projections show that the clusters identified by the hierarchical cluster analysis are well defined and have good separation.
A major advantage of the cluster analysis and the multivariate methods used in this study is that they allow patterns to be identified in multivariate space and reveal trends that are hard to detect using traditional bivariate plots. A good example of this is shown by the olivine data set. Figure 4 shows probability density functions (PDFs) for each of the mineral data sets and the identified clusters generated using Kernel Density Estimation (Silverman, 1986). The molar forsterite content (X Fo = 100 × Mg/(Mg + Fe 2+ + Mn)), molar anorthite content (X An = 100 × Ca/(Ca + Na + K) in moles), and molar magnesium number (Mg # = 100 × Mg/(Mg + Fe 2+ )) were calculated, for olivine, plagioclase and clinopyroxene, respectively. The resulting PDFs were multiplied by the number of data points in each cluster, making them analogous to nonnormalized histograms. The distribution of molar Fo for the complete olivine data set might suggest the existence of two clusters with Fo 76 and Fo 84 . However our clustering approach, which considers the concentration of multiple elements (including minor elements), detects three clusters which could not be identified looking at the PDF of the complete data set alone. This demonstrates how the cluster analysis can identify and reveal otherwise hidden structure in mineral composition data sets by detecting similarities in all elements analyzed.
Olivine compositions were grouped into three clusters, whose cfu contents show the expected proportionality with Fo (Table 2). Fe 2+ , Si and Mn are all inversely proportional to Fo. However, Ca varies independently of Fo and the other elements, and therefore has potential to give insight into differing parental magma minor-element compositions and/or differing H 2 O-contents (Gavrilenko et al., 2016;Kamenetsky et al., 2006).
Plagioclase compositions were grouped into five clusters. These show expected trends with increasing An-content: increasing Ca and Al; and decreasing Si, Na, and K (Table 2). However, Mg and Fe covary independently of the other elements and An-content. That they covary implies that the variations reflect one or more physical processes and are not due to the potentially high analytical uncertainties associated with Fe and Mg in plagioclase (Ginibre & Wörner, 2007). Variations in Fe and Mg in plagioclase independent of other elements, have been attributed to changes in Fe and Mg in the parental melt (Ginibre & Wörner, 2007;Singer et al., 1995), but might reflect melt-plagioclase disequilibrium caused by rapid growth (Ginibre & Wörner, 2007;Mollo et al., 2011;Singer et al., 1995).
Clinopyroxene compositions were grouped into four clusters. However, the representative compositions show more complex trends than for plagioclase and olivine ( Note. The concentration of FeO and Fe 2 O 3 in clinopyroxene was calculated with the method of Droop (1987). Clusters are ordered from the most to the least primitive. Blank cells denote elements not used in the cluster analysis. b.d., below detection limit; -, not calculated; Fo, molar forsterite content; An, molar anorthite content; Mg # , magnesium number. be the result of the higher-complexity of the clinopyroxene mineral structure, compared to feldspar and olivine. However, the identified clusters still appear to show some expected correlations with Mg # , with the most primitive Cpx 1 having the lowest incompatible elements, for example, Mn and Ti, and Cpx 4 having higher incompatible elements and Ca. Cpx 3 has much lower calculated Fe 3+ content compared to the other three clusters, which could reflect significantly different oxidation conditions during crystallization.

Crystal Zoning
Original studies classified each mineral analysis as unzoned and zoned-cores, intermediates, or rims. However, this information was not reported for a significant number of analyses (olivine: 20%, plagioclase: 32%, clinopyroxene: 35%), and there is no way to judge the accuracy of the designations. Therefore it was assumed that each of the source studies categorized their analyses in the same manner. Figure 5 shows the proportion of crystal zones in each cluster for each category. There is no clear relationship between the cluster distribution and zoning for any of the three minerals.  Figure 6 shows the proportion of phenocrysts from each cluster within individual eruptive deposits sorted by age. There is a stark contrast between the composition of crystal cargoes of historic eruptions (erupted in the last 100 years) and prehistoric eruptions (ca. 95-1.6 ka). Generally, the crystal cargoes of historic eruptions contain a smaller variety of crystal compositions. This is especially pronounced for feldspar and clinopyroxene compositions (Figures 6b and 6c). Both feldspar and clinopyroxene compositions become markedly more homogeneous through to the present day: erupted plagioclase compositions are dominated by Pl 4, and clinopyroxene compositions are dominated by Cpx 1 . Generally, historic eruptions contain a higher proportion of Ol 2 than prehistoric eruptions, but olivine compositions are more variable than both plagioclase and clinopyroxene ( Figure 6).

Eruption Age
Of the three minerals considered in this study, olivine cluster proportions vary the most from eruption to eruption. For example, the Licán ignimbrite has a large proportion of primitive olivine (Ol 1 ) whereas the Pucón ignimbrite Clinopyroxene cluster membership broadly corresponds to that of olivine and plagioclase: the greatest amount of variation is present in the oldest eruptions. Clinopyroxene becomes both more scarce and compositionally homogeneous in later eruptions: cluster Cpx 2 , Cpx 3 , and Cpx 4 are most abundant in the prehistoric eruptions while Cpx 1 is noticeably more abundant in the historic eruptions. Very little was erupted after the 1948 eruption.

Eruptive Center Location
Villarrica's flanks host over 30 parasitic or adventitious cones which form two groups: the Los Nevados group to the northeast and the Chaillupén group to the south (Moreno & Clavero, 2006) (Figure 1b). Other than the scoria and lavas of these two groups of cones, all other eruptive deposits are thought to originate from Villarrica's main vent, although the location of this main vent has changed through time (Moreno & Clavero, 2006). The cluster membership for the two groups of parasitic cones are shown in Figure 6 (labeled LNPC and CPC, respectively). The proportion of phenocrysts in each cluster for the two groups of cones are distinct. For example, the Chaillupén cones contain the most primitive olivine (Ol 1 ), and none of the most evolved clinopyroxene (Cpx 4 ), unlike the Los Nevados cones. This could be because the two sets of cones tap different parts of Villarrica's magmatic system, or they sample the system at different points in time. Compared to the eruptions from the central vent, the two sets of parasitic cones are more diverse than the historic eruptions, and most similar to the prehistoric eruptions. This is especially apparent when . Plots of probability density functions (PDF) generated using Kernel Density Estimators of each mineral subset and its clusters. Each PDF is scaled by the number of data points it has generated from, analogous to a histogram. Comparisons of the distribution of identified clusters with each respective complete data set highlights the ability of the cluster analysis to identify otherwise hidden structure in mineral compositions.
comparing the olivine and clinopyroxene cluster memberships. More detailed studies of the two groups of cones are needed to constrain their ages and their petrological differences to each other and products from the main vent (e.g., Robidoux et al., 2021).

Eruption Volume and Intensity
Villarrica's postglacial eruptions show large variation in volume and intensity; most notable are the differences between the explosive high-volume mafic ignimbrites (Licán and Pucón) with the smaller-volume, comparatively effusive, historic eruptions. The large volume of the ignimbrites and their association with caldera collapse suggests that a larger portion of the magmatic system was evacuated during those eruptions. This might be expected to produce products with a higher compositional variety, for example, both evolved and primitive compositions, especially if the magmatic system were to contain multiple somewhat isolated bodies that differentiated independently. This is indeed reflected in the cluster memberships of the historic versus older eruptions. Generally, the more explosive, older eruptions have a higher diversity in mineral compositions, especially the Licán and Pucón ignimbrites. Of the historic eruptions the intense March 2015 eruption, that produced a 1.5 km high fire fountain, erupted more diverse olivine compositions than any other historic eruption.

Thermodynamic Modeling
The representative compositions identified by hierarchical clustering were well reproduced by the relatively simple model setup, shown by the low minimum distances (Table 3). However, the most primitive clusters (Ol 1 and Pl 1 ) have higher minimum distances than the other olivine and plagioclase clusters. This suggests that the initial bulk composition used, a primitive melt inclusion found at Villarrica, was not sufficiently primitive to reproduce exactly these compositions. However the calculated distances are not so different as to suggest that the models do not provide a reasonable indication of the conditions of crystallization.
The best fitting pressure and temperature conditions from the rhyolite-MELTS simulations agree with the broad estimates available from previous thermobarometric studies (Figure 7). Pressures from both our simulations and past thermobarometry imply that polybaric crystallization, extending to at least the mid to lower crust (ca. 600 MPa), is required to produce the variety of erupted mineral compositions at Villarrica. The predicted temperatures from rhyolite-MELTS (800-1100°C) are less than those calculated by thermobarometry (1050-1250°C). This discrepancy is likely due to the high (>2 wt%) water contents of majority of the best fitting rhyolite-MELTS simulations as the olivine-augite geothermometer used by past studies was not calibrated on experiments with high water contents, and therefore likely overestimated temperatures (Loucks, 1996).
The most primitive olivine (Ol 1 ) formed at relatively high temperatures and pressures (1093°C, 382 MPa). Furthermore, Ol 1 has both the highest fO 2 and water content (QFM+0.73, 4.65 wt%). This is in agreement with experimental and natural studies of olivine crystallization from basaltic melts (Feig et al., 2010;Gavrilenko et al., 2016). This suggests that it formed from relatively primitive undegassed melts with a relatively low residence time in the crust prior to crystallization. In contrast, the most evolved olivine cluster (Ol 4 ) formed at lower temperature, pressure, fO 2 , and water content. The two most primitive plagioclase clusters (Pl 1 and Pl 2 ) have high water contents (>4.0 wt%) and relatively high temperatures (°C), as expected from experimental studies (Panjasawatwong et al., 1995;Takagi et al., 2005;Waters & Lange, 2015). In contrast, the most evolved plagioclase (Pl 5 ) has a significantly lower crystallization temperature (978°C), and a notably higher pressure and fO 2 (336 MPa, QFM0.98). All representative clinopyroxene compositions are best fit at lower pressures (−100 MPa), and at a variety of temperatures (970-1041°C) and water contents (0-2.0 wt%).

Vertically and Laterally Extensive Magma Processing
Several models of volcanic systems that integrate geophysical, geochemical, and petrological information have concluded that magma processing beneath arc volcanoes likely takes place throughout the crust in transcrustal magma systems (Annen et al., 2006(Annen et al., , 2015Cashman et al., 2017;Hildreth & Moorbath, 1988). There are multiple lines of evidence that suggest that this is the case at Villarrica volcano. We have compared representative mineral compositions identified via clustering, with thermodynamic fractional crystallization models of a simple Villarrica-like system. The comparisons suggest that polybaric crystallization, at pressures up to ca. 600 MPa, are required to produce the mineral compositions erupted at Villarrica (Figure 7). This is supported by thermobarometry from past studies that focused on individual eruptions (Figure 7). The majority of pressures calculated from thermobarometry are in the shallow crust (0-150 MPa), but pressures calculated from olivine-augite pairs erupted during historic eruptions extend to 700 MPa. Additionally, Kapinos et al. (2016) detected a low conductivity zone beneath Villarrica volcano that extends from 19 to 50 km, beyond the base of the crust (deeper than any of the crystallization pressures identified in our study).

Mineral
Cluster Note. One weighted standard deviation of the best fitting conditions per cluster is shown in brackets. Oxygen fugacity is shown in log units relative to the QMF buffer. T, temperature; P, pressure; Min Dist, the minimum distance between ilr-transformed simulated and measured compositions.

Table 3
The

Weighted Average of Best Fitting Conditions for All Compositions in Each Identified Mineral Cluster Compared to Rhyolite-MELTS Thermodynamic Models
In addition to vertical connectivity, there is also evidence of lateral connectivity within Villarrica's magmatic system (e.g., Ebmeier et al., 2018;Lerner et al., 2020) ( Figure S1 in Supporting Information S1). Nested calderas, kilometers in diameter, were produced during large eruptions at Villarrica and surround the current central vent: Caldera 1 (ca. 100ka), Caldera 2 (ca. 14ka), and Caldera 3 (ca. 3.7ka; Moreno & Clavero, 2006). Calderas 1 and 2 cover an area of 6.5 by 4.2 km (long axes), and Caldera 3 is roughly circular and 2 km in diameter. Domes and sills erupted on caldera walls demonstrate lateral transport of magma during large eruptions (Moreno & Clavero, 2006). After the March 2015 eruption, Delgado et al. (2017) detected reinflation at depths of 4.2 km in part of Villarrica's magmatic system, ca. 5 km SE of the central vent near the edge of Calderas 1 and 2. Further afield still, are the two groups of parasitic cones on Villarrica's flanks. The Los Nevados group contain fissures and cones that extend to ca. 10 km to the NE and the Chaillupén group extend ca. 12 km to the south of the present central vent. Recently, Pavez et al. (2020) detected a low conductivity zone (ca. 4 km deep) associated with the Los Nevados group. Combined, all these lines of evidence demonstrate that Villarrica's magmatic system is spatially extensive and imply that different parts of the system are tapped to accumulate the crystal cargo of each eruption. Figure 7. Weighted average of best fitting temperature and pressures from rhyolite-MELTS simulations for each of the identified mineral clusters. Gray error bars show ± one weighted standard deviation of best fit pressures and temperatures per cluster. pTb, compilation of past thermobarometry from past studies that correspond to a shallow, mid and deep-crustal storage (Lohmar, 2008;Lohmar et al., 2012;Morgado et al., 2015;Pioli et al., 2015;Pizarro et al., 2019;Witter et al., 2004). The error bars correspond to the range of uncertainty for each result. The deep region is based off a single olivine-augite pair and therefore is less certain than the shallow and mid-crustal results. The corresponding depth is calculated assuming a constant density of 2.7 gcm −3 . The black dashed line is the calculated depth (ca. 4.2 km) of reinflation after the March 2015 eruption (Delgado et al., 2017). The gray shading corresponds to an area of high conductivity (ca. 19-50 km) (Kapinos et al., 2016).

The Role of Magma Mixing
The mixing of compositionally distinct magmas has been suggested as the trigger for eruptions at many arc volcanoes globally (e.g., Bouvet de Maisonneuve et al., 2012;Cassidy et al., 2015;Kahl et al., 2011;Ruprecht et al., 2012), to the extent that it is considered ubiquitous . Here, we define magma mixing as the physical interaction of magmas that have different chemical composition and/or intensive variables. It does not imply total homogenization of the two magmas.
The compositional variety observed within erupted crystal cargoes, combined with evidence for spatially extensive magma processing, suggest that Villarrica's magmatic system comprises multiple, variably evolved reservoirs distributed throughout the crust. This is supported by widely varying best fitting intensive variables from rhyolite-MELTS thermodynamic simulations (Table 3). In addition to the range of best-fit pressures, there is a large range of best fitting crystallization temperatures that reproduce crystal compositions (850-1100°C). There is also a large variation in best-fitting water contents from anhydrous to 5.5 wt%. This implies that degassing plays a large role in driving crystallization in Villarrica's magmatic system, as has been suggested for other arc volcanoes (e.g., Blundy et al., 2006;Bouvet de Maisonneuve et al., 2012). The large range of best fitting fO 2 values (−0.5 to 2.00 ΔQFM) suggests that crystallization is occurring from melts that have undergone different amounts of degassing and fractional crystallization (Carmichael, 1991;Lindsley & Frost, 1992;Sato, 1978). The presence of minerals produced by different intensive variables in a single crystal cargo means there must be physical interaction between reservoirs in Villarrica's magmatic system Reservoirs that are infrequently disturbed by ascending primitive magma, are able to cool and differentiate via fractional crystallization to produce evolved crystal compositions for example, Pl 5 (An 64 ). Prior to eruption, ascending primitive melt interacts with multiple reservoirs, accumulating differing minerals. This mixing of magmas with different compositions produces new compositions, which may be recorded in zoned crystals ( Figure 5). These antecrysts are accumulated as melts ascend resulting in the variety of erupted mineral compositions in a single crystal cargo ( Figure 6).
The important role of magma mixing at Villarrica is supported by our observation that the proportion of identified clusters do not show clear trends when they are compared to crystal zoning ( Figure 5). If fractional crystallization is the dominant process during magma processing in Villarrica's magmatic system, we would expect mineral cores to be dominated by primitive compositions and rims by more evolved compositions. However, both primitive and evolved clusters are present as cores, intermediate, and rim zones of all three of the commonly erupted minerals. This implies that magma mixing plays a role in assembling crystal cargoes at Villarrica (e.g., Ruprecht et al., 2012;Streck, 2008).
Furthermore, textural evidence for magma mixing is present in almost all Villarrica's eruptive deposits. The Dacitic Dome contains reverse and oscillatory-zoned plagioclase, olivine with fayalitic rims (Fo 58 ) and reversezoned clinopyroxene (Lohmar, 2008). The Licán ignimbrite has both reverse and oscillatory-zoned plagioclase, reverse-zoned clinopyroxene, and orthopyroxene rims surrounding olivine crystals (Lohmar et al., 2012). A pre-Pucón surge deposit contains banded pumice, the pale bands have a bulk SiO 2 of 63 wt% (Lohmar, 2008;Moreno et al., 1994). The Pucón ignimbrite contains plagioclase crystals with resorbed cores and rims, sieve textures, and low-An microlites. Clinopyroxene is often reverse-zoned (Lohmar, 2008). Additionally, what is thought to be a dacitic enclave in a basaltic scoria bomb of the Pucón ignimbrite has been observed (McCurry & Schmidt, 2001;McCurry et al., 2004). The Chaimilla Fall Deposit contains complexly zoned plagioclase crystals with both reverse rims and widespread evidence of resorption (Pioli et al., 2015). Multiple historic eruptions also contain similar indicators: An-poor plagioclase is present as both rims of oscillatory-zoned crystals and crystals with resorption textures. Both olivine and clinopyroxene crystals display resorption textures (Morgado et al., 2015;Pizarro et al., 2019).
The ubiquity of evidence for magma mixing in Villarrica's eruptive products implies complex magma dynamics involving multiple variably evolved reservoirs. These may be intermittently connected and tapped prior to eruption resulting in the variable crystal cargoes of eruptive deposits.

Linking Crystal Cargo and Whole-Rock Composition to Eruptive Behavior
The complex temporal trends in the compositions of erupted crystal cargoes show strong correlations with both whole-rock composition and eruptive style. There is a stark contrast between the composition of crystal cargoes of historic eruptions and prehistoric eruptions ( Figure 6). Generally, the crystal cargoes of historic eruptions contain a smaller variety of crystal compositions; this is especially pronounced for feldspar and clinopyroxene compositions. This disparity strongly correlates with whole-rock compositions (Figure 6d): historic eruptions have average bulk SiO 2 and MgO contents of 54 and 5 wt%, respectively, whereas prehistoric eruptions range from 54 to 65 wt% SiO 2 and 1-5 wt% MgO. Both the composition of erupted crystal cargoes and whole-rock compositions strongly correlate with eruptive style: the most differentiated eruptive deposits contain the more variable crystal cargoes and are the most explosive and high-volume.

Dacitic Dome
The Dacitic Dome (c.a. 95 ka) is the most differentiated of Villarrica's products analyzed to date, with substantially higher bulk SiO 2 content than any subsequent eruption (65 wt% vs. 55 wt%, Figure 1). However the magma that formed the Dacitic Dome is almost certainly a mix of evolved and more primitive magmas (e.g., Eichelberger et al., 2006;Reubi & Blundy, 2009), supported by textural observations described above (Lohmar, 2008). The dacitic composition may have formed from mafic magmas ascending through evolved mushes (e.g., Reubi & Blundy, 2009) or from differentiated melts ascending through mafic mushes (e.g., Kent et al., 2010). Of the three eruptions that contain it, the crystal cargo of the Dacitic Dome contains the highest proportion of An-poor plagioclase (Pl 5 , Figure 6b), along with its high bulk-SiO 2 content.

The Licán and Pucón Mafic Ignimbrites
The Licán (ca. 14 ka) and Pucón (ca. 3.7 ka) mafic ignimbrites are deposits from the two most explosive, high-volume eruptions at Villarrica (Lohmar et al., 2012;Silva Parejas, 2008). Both have high, but variable, bulk-SiO 2 contents ( Figure 6d) compared to historic eruptions. They are also the last two eruptions whose crystal cargoes contain the evolved Pl 5 . Their whole-rock chemistry correlates with the proportion of Pl 5 in their crystal cargoes: the Licán has lower MgO and slightly higher SiO 2 than the Pucón, and has a higher proportion of Pl 5 (Figures 6b  and 6d).
Rhyolite-MELTS models show that Pl 5 likely crystallized from an evolved melt (ca. 60 wt% SiO 2 ), after substantial differentiation via crystal fractionation, resulting in a bulk crystal fraction of at least 50%. The eruption of a reservoir with such a high crystal fraction requires liberation via mush disaggregation to reduce its viscosity (Sparks & Marshall, 1986). This may result from a combination of chemical or thermal mixing, and/or volatile fluxing from an incoming more primitive basaltic melt (e.g., Bachmann & Bergantz, 2006;Bergantz et al., 2015;Bouvet de Maisonneuve et al., 2012;Pistone et al., 2017;Zellmer et al., 2016). Upon mixing with the primitive magma, plagioclase crystals from the evolved reservoir would become reversely zoned, with low-An cores and higher-An rims. In turn, primitive plagioclase carried by this primitive melt would grow low-An rims upon mixing (as suggested in Lohmar, 2008;Lohmar et al., 2012). Therefore mixing between an evolved and primitive magma explains how Pl 5 can exist as cores, intermediates, and rims ( Figure 5).
The correlation of whole-rock compositions with the proportion of low-An Pl 5 in erupted crystal cargoes, combined with textural data and thermodynamic modeling, supports the triggering of these mafic ignimbrite generating eruptions via destabilization of a differentiated (dacitic) mush by an influx of primitive magma.

The March 2015 Fire Fountain
In contrast to the ignimbrite-forming eruptions, historic eruptions have been dominated by effusive activity, punctuated by fire fountaining which appears to be related to lahar generation (e.g., 1908, 1948-49, 1963, 1964, 1971(Lara & Clavero, 2004). The paroxysmal eruption of March 2015 is significant in that it was the most intense historic eruption at Villarrica, producing a fire fountain that was 1.5 km in height but lasted just 30 min (Romero et al., 2018). The composition of olivine erupted in 2015 are more varied than both all preceding historic eruptions, and in spatter erupted later that year ( Figure 6). This is in stark contrast to the March 2015's otherwise homogeneous cargo, which consists of only intermediate plagioclase clusters Pl 3 and Pl 4 and few clinopyroxene crystals. This implies that the mafic portion of the March 2015's crystal cargo was accumulated by primitive melt as it ascended through Villarrica's magmatic system. Plagioclase may have formed in response to degassing upon ascent (as in Blundy et al. [2006]) as with the other historic eruptions. The low-SiO 2 and high-MgO bulk composition of erupted spatter suggests that the intensity of the eruption is unlikely to be caused by the same mechanism as the mafic ignimbrites. Instead high primary volatile contents likely drove fast magma ascent, resulting in the vigorous fountaining behavior (Allison et al., 2021;Barth et al., 2019;La Spina et al., 2021) and allowing it to punch through existing reservoirs and assemble its more varied crystal cargo.

A Model for Generating Explosive Eruptions at Villarrica Volcano
The eruption of the Dacitic Dome (ca. 95 ka) demonstrates that evolved portions of Villarrica's magmatic system existed prior to the most explosive postglacial eruptions known at Villarrica, the Licán ignimbrite. Therefore portions of Villarrica's magmatic system were differentiated beyond the typical whole-rock compositions of erupted historic products (Figure 1). The existence of this evolved (dacitic) reservoir has implications for magma dynamics. It has been shown that felsic mushes are sufficiently viscous that they are noneruptible (Marsh, 2002;Sparks & Marshall, 1986). In combination with their low relative density, they can act as a density filter, preventing more dense mafic melts from ascending to the surface (Kent et al., 2010). Therefore the establishment of a significant volume of evolved mush within the magmatic system could dampen the eruption of mafic magma and increase the overall volume of evolved melt until a sufficiently large volume of ascending primitive melt is able to destabilize it (Sparks & Marshall, 1986).
The eruption of large-volume eruptive deposits, such as Licán and Pucón ignimbrites, requires significant volumes of mobile magma to exist within the magmatic system (Druitt & Sparks, 1984). If the majority of Villarrica's magmatic system is composed of near-solidus mush (as implied by the TCMS model), this would require a large volume of magma to destabilize it (Marsh, 2002;Sparks & Marshall, 1986). The composition of minerals (Ol 1 and Pl 1 , Table 2) and melt inclusions at Villarrica suggest primitive magmas were present within its system (Pioli et al., 2015). The best fitting conditions for the most evolved plagioclase (Pl 5 ) are in the middle crust, at low temperatures and after significant crystal fractionation (Figure 7). Therefore we propose that the trigger of the Licán and Pucón ignimbrites was a large influx of primitive magma, which mixed with an evolved reservoir (Figure 8a). The remnants of the mush that this primitive magma interacted with is shown by the presence of evolved plagioclase, Pl 5 in both ignimbrites' crystal cargo. In the case of the Licán ignimbrite, the trigger may have been facilitated by widespread deglaciation in the Southern Andes around 14 ka ago, altering the stress state of the crust to facilitate magma ascent (e.g., Jellinek et al., 2004;Rawson et al., 2016;Watt et al., 2013;Wilson & Russell, 2020). After the eruption of the Licán ignimbrite, another evolved reservoir of similar composition likely accumulated over a period of ca. 10 ka. The existence of an evolved reservoir prior to the Pucón Ignimbrite (ca. 3.7 ka) is suggested by the mineralogy of a Pre-Pucón pyroclastic surge deposit Lohmar (2008), which contains a large proportion of evolved (63 wt% SiO 2 ) pumice. The Pucón ignimbrite forming eruption may then have been triggered in a similar fashion to the Licán eruption, that is, destabilization of an evolved reservoir by an influx of primitive magma.
These evolved reservoirs were likely depleted by mixing with primitive magma prior to the ignimbrite-forming eruptions, as subsequently erupted crystal cargoes (3.7 ka to present) contain no trace of Pl 5 (Figure 6). Without the density filter provided by the evolved reservoirs, the most recently erupted magmas were able to ascend through Villarrica's magmatic system faster, perhaps only mixing in the shallowest parts, resulting in less differentiated erupted compositions (Figure 8b). The absence of mixing between highly differentiated and primitive compositions results in substantially lower explosivity. Instead recent lava fountaining is likely driven by rapid ascent, enabled by high volatile-contents (e.g., the March 2015 eruption). The near continuous activity at Villarrica's lava lake over the last ca. 40 years demonstrates the ability of magma to ascend relatively unhindered through the magmatic system. Absence of activity at the summit and/or the eruption of crystals in equilibrium with differentiated melt compositions (e.g., Pl 5 ) may signal the onset of evolved reservoir development, and signal an increased likelihood of high-volume explosive eruption occurring.

Conclusions
Use of the multivariate cluster analysis has allowed us to identify previously unidentified structure in the composition of minerals erupted throughout Villarrica's eruptive history. Comparisons of identified representative compositions with >1,500 Rhyolite-MELTS thermodynamic simulations show that magma processing at Villarrica takes place at a range of pressures, temperatures, water contents, and oxygen fugacities. These thermodynamic models strongly suggest that magma storage beneath the volcano is characterized by a series of ephemerally connected, variably evolved mush dominated sills. Prior to eruption, magma ascends through this transcrustal magmatic system, accumulating different antecrysts, which are erupted as variable crystal cargoes.
We have identified temporal trends in the composition of erupted crystal cargoes that correlate with trends in whole-rock composition and eruption style. We propose that prior to high-volume, ignimbrite-forming explosive mafic eruptions (the Licán and Pucón ignimbrites) much of Villarrica's magmatic system had differentiated by fractional crystallization. These midcrustal-evolved reservoirs acted as density filters, suppressing the eruption of small volumes of primitive magma. Only when a sufficiently large influx of primitive magma entered the system was this evolved portion destabilized. Magma mixing ensued, producing the complex textures observed by past studies, and resulting in the relatively primitive, basaltic-andesite compositions of the two ignimbrites. The only remnants of the evolved reservoir are low-anorthite plagioclase, which is present as both cores and rims of more primitive crystals.
The lack of evidence for an evolved reservoir after these high-volume explosive eruptions explains the homogeneous, more primitive, and whole-rock compositions of subsequent eruptions. Without the density filter of the evolved reservoirs, more primitive melts ascend relatively uninhibited, mixing and degassing in the shallow subsurface before eruption. Further petrological work utilizing melt inclusion compositions, trace-element data, and diffusion chronometry will allow further investigation of magma dynamics at Villarrica.