A Recurrent Magmatic Pattern on Observable Timescales Prior to Plinian Eruptions From Nevado de Toluca (Mexico)

Nevado de Toluca is a stratovolcano in the densely populated Mexican Highland and the source of three late Pleistocene Plinian eruptions. While the characterization of the magmatic processes and timescales, leading to these events, is crucial to evaluate potential signals of the reawakening of this volcano, they are not well constrained. Here we present insights gathered from the analysis of mineral zoning. Abrupt changes in plagioclase compositions (anorthite [An] >10 mol. %) record repeated thermochemical disturbance of a shallow magma reservoir prior to the three eruptions investigated. In particular, sudden chemical changes in the outermost crystal rims indicate that the eruptions were triggered by recharge of less evolved magma. Recalculated melt compositions reveal that magma input is chemically heterogeneous. Orthopyroxene compositions span a large range (En53‐En92, Cr2O3 of up to 0.9 wt. %) resulting from magma hybridization, which is consistent with a late‐stage temperature increase recorded by amphibole chemistry. Using diffusion chronometry, we show that this recharge typically occurs decades to centuries before Plinian eruptions and hence on timescales relevant for volcano monitoring. Additionally, modeling of Mg profiles in plagioclase constraints the duration of differentiation cycles to be millennial in order of magnitude. Our study shows that similar processes preceded explosive eruptions from Nevado de Toluca, suggesting that analogous paths of unrest should be considered when evaluating potential future activity. We emphasize that signs of deep crustal activity, such as relatively broad deformation pattern surrounding the volcano, could herald the reactivation of Nevado de Toluca.


Introduction
Many of the largest eruptions that occurred over the last 40 years were produced by volcanoes, which previously did not receive much scientific attention due to their extended periods of dormancy (Carn et al., 2009;Scaillet & Evans, 1999;Luhr et al., 1984). Thus, the study of the long-term eruptive history of volcanoes capable of explosive eruptions and the comparison between such systems is essential to anticipate their potential reawakening (e.g., Bacon & Lanphere, 2006 ;Óladóttir et al., 2008 ;Siebe et al., 1996 ;Singer et al., 2008). While erupted volume and age estimates are essential in this respect, petrological and geochemical studies can provide critical information on the processes leading to an eruption (e.g., Befus & Gardner, 2016 ;Caricchi & Blundy, 2015 ;Costa et al., 2004 ;Muir et al., 2014), as well as their duration (e.g. Costa et al., 2008, Druitt et al., 2012, Hartley et al., 2016, Morgan et al., 2006, Ruprecht & Plank, 2013, Saunders et al., 2012, Zellmer et al., 1999. The chemical zoning of minerals is an archive of the variations of the chemical and physical conditions in subvolcanic magma reservoirs. Plagioclase crystals have proven to be useful recorders of magmatic environments for various reasons, including its ubiquitous occurrence in almost all magma types, compositional sensitivity to key magmatic system variables such as P H2O -T-X (Cashman & Blundy, 2013), and the preservation of original chemical zoning due to sluggish coupled diffusion of CaAl-NaSi (Grove et al., 1984). However, plagioclase are also complex because of their dependence of trace element partitioning on their major element compositions (Bindeman et al., 1998;Dohmen & Blundy, 2014). A number of studies have proposed that combined textural and compositional variation in plagioclase can be used to track diverse magmatic processes such as recharge of fresh magma batches into crustal reservoirs (e.g., Shcherbakov et al., 2011;Ruprecht & Wörner, 2007;Humphreys et al., 2006;Ginibre et al., 2002;Singer et al., 1995), thermally driven convection and overturn (Couch et al., 2001;Triebold et al., 2006), melt segregation in silicic crystal mushes (Andersen et al., 2018;Singer et al., 2016), or H 2 O undersaturated magma ascent (Blundy & Cashman, 2001;Nelson & Montana, 1992;Viccaro et al., 2010).
Orthopyroxene zoning is widely employed as diffusion chronometer to estimate the timing of magma hybridization (Allan et al., 2013;Allan et al., 2017;Barker et al., 2016;Chamberlain et al., 2014;Cooper et al., 2017;Dohmen et al., 2017;Singer et al., 2016) and to directly compare petrological information to monitoring data sets at active volcanoes (Kilgour et al., 2014;Saunders et al., 2012). While mineral studies have documented the processes leading to specific eruptions, petrological investigations at the mineral scale for sequences of volcanic eruptions are so far scarce. However, such information is required to identify temporal variations of preeruptive magmatic conditions, which can be a powerful tool to interpret the evolution of monitoring signals (Andersen et al., 2018;Singer et al., 2014;Stock et al., 2016;Weber & Castro, 2017).
In this contribution, we focus on a sequence of three Plinian eruptions from Nevado de Toluca volcano in the central Trans Mexican Volcanic Belt (Arce et al., , 2003Bloomfield et al., 1977;Capra et al., 2006;Macías et al., 1997). Using textural and chemical information stored in plagioclase, amphibole, and orthopyroxene crystals, we constrain the preeruptive architecture of magma storage and dynamics preceding these event and assess the role of magma heterogeneity in petrogenesis at Nevado de Toluca. Further, we determine the timescale of preeruptive magmatic processes using diffusion modeling of Mg in plagioclase and Fe-Mg exchange in orthopyroxene. We finally discuss the implications of our study for the interpretation of monitoring signals preceding large explosive volcanic eruptions of Nevado de Toluca.

Petrology and Eruptive History of Nevado de Toluca
Nevado de Toluca is a long-lived stratovolcano in the central Mexican highland about 80 km southwest of Mexico City (Figure 1). It is part of the Trans Mexican Volcanic Belt, a 1,000-km-long, east-west trending volcanic arc that developed in response to subduction of the oceanic Cocos and Rivera plates beneath the North American continental lithosphere (Ferrari et al., 2012;Pardo & Suárez, 1995). Detailed stratigraphy and radiometric dating shows that volcanic activity at Nevado de Toluca started in the early Pleistocene at 1.5 Ma and was frequently interrupted by long periods of dormancy between eruptions and sector collapses (Arce et al., , 2003Bloomfield & Valastro, 1974;Capra et al., 2006;Capra & Macıas, 2000;García-Palomo et al., 2002;Torres-Orozco et al., 2017).
The most recent eruptive history (since 42 ka) is dominated by several dome collapse events that are preserved in block-and-ash flow deposits (García-Palomo et al., 2002) and by a sequence of three Plinian eruptions. The oldest and least voluminous of these Plinian events is the 0.85-km 3 dense rock equivalent (DRE) Lower Toluca Pumice (LTP) that was dated at 26 ka cal. 14 C years BP (Capra et al., 2006). In the following 13 ka, Nevado de Toluca was again dominated by dome formation and destruction, notably at about 16,500 cal. 14 C years BP when a massive block and ash flow of 0.5 km 3 was produced . Activity resumed approximately 14,000 years ago with the 1.8-km 3 DRE Middle Toluca Pumice (MTP) Plinian eruption, the deposits of which are distributed mainly to the northeast of the volcano (Arce et al., 2005). The largest and last Plinian-type eruption is the 8km 3 DRE Upper Toluca Pumice (UTP), which occurred around 12.5 ka. The deposits of this eruption are distributed in the area today occupied by Mexico City within the 10-cm isopach (Arce et al., 2003). The last eruptive activity at Nevado de Toluca was described as a phreatomagmatic surge and dated by Macías et al. (1997) at 3,768 years BP (calibrated 14 C age).
Magmas erupted by Nevado de Toluca throughout its eruptive history are almost exclusively dacitic and only in minor proportions andesitic (Torres-Orozco et al., 2017). Several petrological studies have been conducted on the Plinian deposits Arce et al., 2013;Smith et al., 2009) and are summarized here. Petrographic analysis of the Plinian pumice shows that UTP and MTP deposits are identical in terms of their mineral phase assemblage pl > opx > hbl ≫ ilm + mt. The subordinate content of opx and the presence of schists lithic fragments distinguish the LTP from MTP and UTP (Arce et al., 2013). Estimates of preeruptive conditions (P-T-X-fO 2 ) for all three eruptions have been obtained by experimental phase equilibria and twooxide thermometry. Preeruptive temperatures of~868°C, 802 ± 7°C, and 824 ± 12°C were estimated by two-oxide thermometry for the LTP, MTP, and UTP, respectively Arce et al., 2013). Experimental results suggest similar storage pressures for the three eruptions of 150 to 200 MPa and H 2 O saturated conditions (Arce et al., 2013. The shallow part of the magmatic system prior to the UTP eruption evolved as an open system, where fresh pulses of magma entered the preexisting reservoirs (Smith et al., 2009). Martínez-Serrano et al. (2004 used isotopic and trace element evidence to show that crustal melting plays only a marginal role in the petrogenesis of the silicic magmas.

Samples and Analytical Methods
We present 11 new whole rock major and trace element analysis for the three Plinian eruptions (UTP, MTP, and LTP) and for several peripheral cones surrounding the stratovolcano (Figure 1). The latter are scarce in the area and were chosen based on recent geomorphological mapping (Torres-Orozco et al., 2017) to probe the potential mafic input into the plumbing system of Nevado de Toluca, which erupts monotonously dacitic magma ( Figure 1). The samples were processed according to standard techniques, and major element analysis was conducted using the PANanalytical Axios max X-ray fluorescence spectrometer on fused glass discs at the University of Lausanne. Whole rock trace element contents were determined on the same fused glass beads used for major element analysis. The measurements were carried out by laser ablation-inductively coupled plasma-mass spectrometry (LA-ICP-MS) using a quadrupole spectrometer Agilent 7700x coupled to an UP-193FX ArF excimer ablation system at the University of Lausanne. Major and minor element contents in plagioclase, orthopyroxene, amphibole, matrix, and plagioclase-hosted melt inclusions were determined using a JEOL 8200 electron probe microanalyzer at the University of Geneva. We analyzed trace element contents along profiles from core to rim by (LA-ICP-MS) using a sector-field spectrometer Element XR coupled to a RESOlution ArF excimer ablation system equipped with an S155 two-volume ablation cell at the University of Lausanne. Detailed analytical conditions and modeling procedures are presented in the supporting information extended methods file.  Arce et al., 2005Arce et al., , 2003. Inset on the top left shows the location of Neogene to recent Trans-Mexican Volcanic Belt in Mexico (tan shaded area) modified after Ferrari et al. (2012). Dashed lines in the inset represent the inferred depth of the subducting slab (Pardo & Suárez, 1995), and triangles mark the location of active volcanic systems: Ce: Ceboruco, Tq: Tequila, Co: Colima, MG: Michoacán-Guanajuato, NT: Nevado de Toluca, P: Popocatépetl, Or: Pico de Orizaba. The <35-ka eruptive history and dense rock equivalent volumes of the studied Plinian eruptions are shown in the inset on the lower right side (Arce et al., , 2003Capra et al., 2006). Colored triangles and circles correspond to the studied eruptions, and white triangles mark dome collapse events and other pyroclastic activity (García-Palomo et al., 2002;Macías et al., 1997).

Whole Rock and Glass Geochemistry
Major and trace element contents in whole rock and matrix glass for the three studied Plinian eruptions (UTP, MTP, and LTP), as well as for six peripheral cones surrounding the volcano, are displayed in Figure 2. In agreement with previous analyses Capra et al., 2006), the UTP and MTP eruptions are of dacitic whole rock composition with a median SiO 2 contents (wt. %) of 65.73 ± 0.35 (1 sigma) and 65.55 ± 0.33. The products of LTP are slightly more mafic and more variable ranging from the andesite to dacite field with an average composition of 62.55 wt. % and a standard deviation of 1.94 wt.%. The peripheral cones erupted magmas of diverse chemistry with one sample showing elevated total alkali (K 2 O + Na 2 O) contents of 7.32 wt. % (Figure 2a  lithophile element (LILE) contents. While the general shape of the trace element pattern for the Plinian eruptions is nearly identical, differences exist in the degree of enrichment relative to primitive mantle. Samples from the MTP eruption are typically less enriched in incompatible elements such as Th (median [Th] PM = 31.45) or REE compared to the LTP (median [Th] PM = 51.04), which shows the highest degree of enrichment for most trace elements, even though the SiO 2 content of this eruption is the less evolved. The trace element contents of the UTP eruption are intermediate between MTP and LTP, but exceptions exist for some of the LILE [Cs and Rb], which are more enriched in the UTP and MTP by a factor of 1.5 to 1.7 on average for Cs and 1.3 to 1.4 for Rb.
The peripheral cones exhibit similar trace element patterns that are, however, more diverse compared to the Plinian eruptions. In particular, the LILEs are highly variable with normalized Ba contents of 64.82-122.35 and [Sr] PM = 29.60-67.02. Nb-Ta contents are comparable to the Plinian samples for most monogenetic cones except the more alkali-rich samples that show less pronounced negative anomalies. A further characteristic is the greater diversity in LREE/HREE ratios, which result from the higher degrees of enrichment in LREE in the cone samples.

Mineral Textures and Compositions
Plagioclase crystals with sodic cores are the most common textural type in the studied samples and comprise between 36% and 48% of the crystals within the sample population for the three eruptions (Figures S1 and S2). Typically, the transition between the sodic core and mantle is characterized by rounded or irregular wavy resorption surfaces filled with more calcic plagioclase. Normal zoned or calcic core crystals comprise 24% to 37% of the sample population. In stark contrast to that, sieve textured cores appear as complex intercalation of distinct sodic and calcic patches and melt inclusions of variable sizes between few and 100 μm. Especially in larger grains, sieve textured cores can account for about 90% of the crystal volume. The mantle and rims of crystals show oscillatory zoning and resorption zones. Major resorption horizons show wavy embayment and sudden increase in An content of 10-20 mol. %, often followed by progressive decrease. These are ubiquitous and define repetitive sawtooth-like patterns in An content.
The relationship between plagioclase anorthite content and textural features is displayed in Figure 3. The range of An contents is very similar for the three eruptions and varies mostly between 25 and 65 mol. %. Median compositions are 44 mol. % for the UTP and 41 mol. % for the MTP. The LTP plagioclase is overall more calcic with a median value of 48 mol. %. Generally, the different textural zoning features display a wide range of An contents that almost spans the entire variety of recorded compositions for the three eruptions. However, the outermost rim composition is defined by a sharp peak at An 47 and An 53 for the UTP and LTP, respectively, indicating that a large subset of the crystal population in these eruptions recorded similar conditions shortly before eruption. In contrast, MTP rims show a broad bimodal distribution with maxima at An 43 and An 48 . The compositions of prominent resorption zone and sieve-texture-plagioclase for all three eruptions are similar to the rims. Core and mantle compositions are more sodic compared to resorption zone overgrowth and rims in the case of the UTP and LTP eruptions. Plagioclase from MTP samples show a broad compositional overlap between cores and rims, even if mantles and cores tend to be slightly more anorthitic.
Orthopyroxene (opx) crystals are present in the three studied eruptions and show considerable textural and chemical diversity. On the base of textures (Figure 4), we identify two main types of opx crystals based on contrasting brightness in BSE images. Relatively homogeneous ones (e.g., Figure 4e) were only observed in the UTP samples. These crystals, which comprise about 30% of the sample population for UTP, typically show a slight enrichment in En in the outermost rim, as indicated by darker colors in BSE images, and faint sector zoning. The second textural type is composed of crystals that show high contrast between low En rim and high En core compositions, which frequently show wavy to rounded boundaries indicative of resorption and overgrowth (Figures 4b-4d and 4f). In some cases, pervasive finger-like resorption and overgrowth pattern are present. Opx crystals in the upper part of the MTP deposit are similar to the high En cores but rimmed by hornblende instead of low En opx. Some high En cores show internal reverse zonation pattern (Figures 4a and 4d).
Major element analysis of opx are shown in Figure 5. En content spans a large range of 53.2 to 91.5 mol. %. Cr 2 O 3 contents range from undetectable corresponding to En 55 -En 65 to up to about 1-wt. % Cr 2 O 3 for En > 80 mol.% (Figures 5a and 5b). The compositional variability of opx for the studied eruptions spans a large fraction of the compositional diversity observed in opx in different volcanic arcs, extending into the compositional field of peridotites. Based on their major element compositions, the opx can be divided into three distinct groups, which can be present in a single crystal ( Figure 4f). The first includes pyroxenes with En contents <66 mol. % and low Cr 2 O 3 contents, which correspond to homogeneous textural Type 1 crystals of the UTP eruption, and bright rims of textural Type 2. This group can be clearly identified in the UTP and MTP eruptions, between which it is distinguishable by slightly higher En and MnO contents in UTP with respect to MTP (Figures 4c-4f). Group one cannot be clearly identified in LTP as <En 70 compositions are more scattered compared to the other two eruptions. The second group is characterized by intermediate En contents roughly between 70 and 80 mol. %. These crystals, which are present in the UTP and LTP but not in the MTP samples, show elevated Cr 2 O 3 contents of up to 0.58 wt. % and higher TiO 2 and Al 2 O 3 (wt. %) contents compared to zones of group one. The third group is only well defined for MTP and LTP samples and encompasses crystal zones with >En 90 and typically high Cr 2 O 3 contents, while Al 2 O 3 and TiO 2 (wt. %) are lower compared to group two compositions.
Amphibole phenocrysts (n = 65) are omnipresent in the studied deposits. Crystals are typically euhedral to subhedral in shape and often show sieve textures with abundant melt and mineral inclusions in the crystal cores ( Figure S4). Texturally simpler amphiboles are also present, but the transition between these features is gradual. LTP amphiboles show a greater proportion of sieve textures compared to the MTP and UTP. Notably, several crystals (4/22 UTP, 5/21 MTP, and 10/22 LTP) show compositionally distinct overgrowth bands on the outermost rims. Chemically, the amphiboles are Ca rich with Ca ranging between 1.34 and  (Ulmer et al., 2018;Blatter et al., 2013;Müntener et al., 2001;Blatter & Carmichael, 2001 1.81 a.p.f.u. (based on 23 oxygens) and can be classified mainly as pargasite with minor Mg-hornblende and edenite. Si, which is a strong function of temperature in amphibole (Putirka, 2016), varies between 6.10 and 7.34 a.p.f.u and shows linear inverse correlations with Al (R 2 = 0.86), Na (R 2 = 0.64), and Ti (R 2 = 0.66). Exchange component calculations show that compositional variability in these crystals is mainly driven by Fe-Mg exchange, with contributions from edenite and tremolite exchange. Amphibole compositions from the three studied eruptions show similar compositional diversity in most of the analyzed elements. An exception are variations in K content (Figure 6a), for which different trends can be distinguished between the three studied eruptions, but also within the amphibole population of a single eruption (MTP).

Architecture of the Plumbing System
Based on experimental evidence, it is well established that dacitic to andesitic eruptions from Nevado de Toluca, similar to many other stratovolcanoes, are fed from a shallow crustal storage region (Castro et al., 2013;Costa et al., 2004;Cottrell et al., 1999;Muir et al., 2014;Scaillet et al., 2008;Weber & Castro, 2017), in this case located between 4.5-and 6-km depth (Arce et al., 2013. Amphibole barometry provides a further test of where magmas accumulate beneath Nevado de Toluca prior to the Plinian eruptions ( Figure 6). The Al-in hornblende barometer has been widely used to constrain magma storage pressures (e.g. , Mutch et al., 2016, Ridolfi & Renzulli, 2012. However, a recent evaluation of the available data used to calibrate the barometric models shows that this approach faces severe challenges (Putirka, 2016), as it yields uncertainties of ±4 kbar for single spot analysis. Averaging of individual analysis from many amphibole crystals can improve the resolution, reducing the uncertainty to ±1 kbar (Putirka, 2016). Figure 6b shows results of thermobarometric calculations for the eruptions under consideration. Averaged pressure estimates are 2 kbar for the UTP and MTP eruption, while LTP amphiboles record higher pressures of 2.8 kbar. These findings are in good agreement with experimental reconstruction (Arce et al., 2013 and independently validate that silicic magmas are stored in and erupted from a shallow feeding system beneath the volcano. At Nevado de Toluca, the presence of opx consistently spanning a large range of compositions in the Plinian eruptive products provides a unique opportunity to constrain the chemistry of magma input and together with amphibole barometry to gather useful information on the architecture of the magmatic plumbing system. Constraining the conditions at which magma evolves in deeper portions of the systems, where more primitive melts might stall and differentiate, is difficult. Such difficulties arise especially because of overprinting effects such as shallow crustal magma mixing and homogenization Burgisser & Bergantz, 2011;Charlier et al., 2007;Huber et al., 2009;Reubi & Blundy, 2009), resorption and peritectic reactions consuming mafic mineral phases (Coombs & Gardner, 2004;Erdmann et al., 2014;Kushiro, 1969;Zellmer et al., 2016), or silicic magma petrogenesis at deep crustal levels (Annen et al., 2005;Annen et al., 2008;Jackson et al., 2018;Solano et al., 2012).
Fe 2+ -Mg exchange between opx and basaltic liquids is mainly a function of melt composition, which in turn might respond to changes in temperature, pressure, or variability in water content through destabilization of other phases (Frey & Lange, 2011;Gaetani & Grove, 1998;Putirka, 2008;Waters & Lange, 2017). Using phase equilibria experiments, Waters and Lange (2017) showed that in rhyolitic magmas, Fe-rich opx (low En content) is stabilized by high water contents, while high-Mg opx crystallize from low water content melts. The most Mg-rich opx crystals were measured in LTP, which is consistent with the lower water content of the magma from this eruption as calculated from the plagioclase-liquid hygrometer of Waters and Lange (2015;4.5-5 wt. % for LTP with respect to 5.5-6 wt.% for the UTP and MTP). This effect may arise due to preferential complexing of dissolved hydroxyl groups with Mg 2+ relative to Fe 2+ , lowering the activity of the enstatite component (Waters & Lange, 2017). Like for basaltic magmas, no significant effects of temperature and pressure on the Kd[Fe 2+ -Mg] opx were observed, at least in the range of experimental conditions (ΔT = 65°C; ΔP = 105 MPa). In summary, the major element contents in opx are mainly a function of melt chemistry and water content in the case of evolved magmas (Putirka, 2008). Thus, we suggest that the three groups of opx compositions crystallize from melts of different compositions (i.e., degree of chemical evolution). Zones from group one with median compositions En = 59.3 mol. % are consistent with derivation from the host dacites, which is corroborated by their presence as euhedral phenocrysts in UTP and equilibrium shaped subhedral to euhedral crystal rims in all three eruptions. While a clear cluster of opx compositions with En < 70 mol. % is missing for the LTP eruptions, the abundant but more scattered En 55-70 compositions are not surprising taking into account that the whole rock data span a compositional range between andesite and dacite. Comparison of the En-Cr 2 O 3 systematics with opx from different volcanic arcs and peridotites from the GEOROC database (http://georoc.mpch-mainz.gwdg.de/georoc/) shows that the crystals from group three of LTP and MTP resemble pyroxenes from mantle equilibrated primitive magmas (Figure 4a). Curvilinear trends in crystals of group three, descending from high to low Cr 2 O 3 and slightly lower En contents (Figures 4 and 5), indicate that the compositional evolution of these grains is controlled by a melt crystallizing decreasing fractions of phases, such as olivine and spinel, with high Cr partition coefficient. Such trends resemble the shape of the volcanic arc opx rather than mantle peridotite data (Figure 5a). These evidences exclude that the high En-Cr 2 O 3 opxs are actually mantle xenocrystals.
Comparison of the natural opx with experimental phase equilibria of primitive arc magmas at various pressures (Blatter et al., 2013;Blatter & Carmichael, 1998;Blatter & Carmichael, 2001;Müntener et al., 2001 ;Ulmer et al., 2018) can help to further constrain the origin of these crystals, as they may have crystallized over a considerable pressure range. Enstatite contents of >80 mol. % are reproduced both by crystallization between 400 and 50 MPa (Blatter et al., 2013;Blatter & Carmichael, 2001) and at lower crustal pressures of 1-1.2 GPa (Müntener et al., 2001;Ulmer et al., 2018). Cr 2 O 3 contents >0.6 wt. % in experimental opx also crystallize from melts with compositions similar to some of the most primitive natural compositions. For both shallow and deep crystallization, Cr 2 O 3 and En contents are decreasing with temperature, which is related to the progressive depletion of the melt mainly due to olivine or Cr-spinel precipitation. The high Cr 2 O 3 content of opx in the natural samples thus reflect its near-liquidus crystallization from a liquid with initial Fe/Mg ratios low enough to stabilize high-En opx (Blatter et al., 2013;Nandedkar et al., 2014). Basaltic andesite, which typically has opx as the liquidus phase (Müntener et al., 2001;Ulmer et al., 2018) and is commonly observed in the central part of the Trans Mexican Volcanic Belt (Blatter et al., 2007;Ferrari et al., 2012), is a likely candidate for the origin of group three opx. Good correspondence between experimental and natural data can be observed for contents of various major elements (e.g., TiO 2 , MnO, and CaO) at En contents typical of group three (>80 mol. %), as well as for the evolutionary trends of these elements with decreasing En content and presumably decreasing temperature ( Figure 5).  (Blatter et al., 2013;Blatter & Carmichael, 2001. This mismatch in Al 2 O 3 contents suggests that opx from group two originated at relatively shallow crustal levels where plagioclase is saturated along with orthopyroxene. For example, opx En-Al 2 O 3 contents that compare well to group two compositions were experimentally grown by Blatter et al. (2013) by crystallization of a basaltic liquid with 2 wt. % initial water content at 400 MPa and 975°C in presence of plagioclase. Alternatively, group two opx can also be generated at lower pressure and water-saturated conditions (P H2O between 100 and 141 MPa at 909-950°C).
Chemical profiles were collected for selected crystals from the three eruptions ( Figure 4). Group three opx in the LTP samples show intracrystal core variation that is inconsistent with closed system crystallization, such as increasing En and Cr contents toward the outer portions of the crystal core ( Figure 4a). Interestingly, the resorption and overgrowth zones, which formed by interaction with the erupted host magma, are only a few micrometers wide. This indicates that these cores, which are in disequilibrium with the host melt, became opx-saturated and grew a rim shortly before eruption. In contrast, crystals from group two show broad overgrowth bands of several tens of micrometers, which shows that they were either interacting with the dacitic magma for longer time or that the magma was opx saturated early on after entrainment. Similar degrees of interaction are preserved in the zoning pattern of crystals from the MTP eruption, which sometimes show evidences of high-pressure crystallization in the innermost cores (i.e., elevated Al contents in 104 out of 3,332 spot analyses). Opx from groups two and three are not simply related by equilibrium or fractional crystallization but seem to be the product of the interaction between magmas of different chemistry, as exemplified by sharp increases of Cr and Ti contents along the crystal stratigraphy. Late stage increases in En are present in many outermost rims within group one crystals. While a decrease of water content driven by magma ascent could increase En content (Waters & Lange, 2017), the coupled increase of En, Ca, and Al rather suggest that the interaction with a different magma or chemical changes due to mineral breakdown (Allan et al., 2013) is responsible for this chemical signal.
The wide compositional array preserved in opx shows that the plumbing system of Nevado de Toluca consists of storage regions that are spatially separated and, given the spread of pressures obtained by amphibole barometry (Figure 6b), may be located at different depths. Magma hybridization is often inferred from the presence of disequilibrium crystal populations and zonation patterns (e.g., Ganne et al., 2018;Kent et al., 2010;Humphreys et al., 2009;Klemetti & Grunder, 2008;Ruprecht & Wörner, 2007;Browne et al., 2006) but could also be explained by entrainment of xenocrysts or antecrysts from previously fully or partly solidified intrusions.

Magma Recharge
Plagioclase textures and chemistry have proven useful to trace recharge in magma reservoirs (e.g., Cashman & Blundy, 2013;Humphreys et al., 2006;Davidson & Tepley, 1998;Singer et al., 1995). The An content in plagioclase varies as function of temperature, PH 2 O, and melt composition (Cashman & Blundy, 2013), and the combination of major, minor, and trace elements measurements can be used to disentangle these 10.1029/2019JB017640 Journal of Geophysical Research: Solid Earth effects (Ruprecht & Wörner, 2007). Minor and trace element partitioning are typically mainly dependent on melt composition and to a lesser extent on variables such as temperature or PH 2 O (Bindeman et al., 1998;Dohmen & Blundy, 2014). Previous work on plagioclase from the UTP eruption (Smith et al., 2009) has revealed late stage sudden increases in An content (>10 mol. %) for many crystals coupled with increases in Mg content, indicative of magma recharge shortly prior to eruption. Our results confirm these findings and moreover show repetitive cycles of magma injection, which affect plagioclase crystals especially in their outer portions. Chemical profiles in plagioclase from Plinian eruptions show good correspondence between abrupt changes in An content and elements such as Fe, Mg, Ti, and Sr (dashed lines in chemical profiles of Figure 7). These correlations, together with the large abundance of reverse zoning and stepped sawtooth textures, point toward an open magmatic system that undergoes differentiation and repeated recharge cycles that eventually trigger an eruption.
An alternative hypothesis is that the magma was initially H 2 O undersaturated and that late stage volatile saturation triggered the eruptions, as gas volume expansion can be expected to create large overpressures in magmatic reservoirs (Tait et al., 1989). Plagioclase-liquid hygrometry shows that preeruptive H 2 O contents for the three eruptions are likely between 5 and 6 wt. %, which when compared to the 2 ± 1 kbar (UTP and MTP) 2.8 ± 1 kbar (LTP) constrained by amphibole barometry, and experimental phase equilibria (Arce et al., 2013, indicate that the Toluca magmas could have been close to or beyond volatile saturation. In either case, a key observation is that the ubiquitous resorption and overgrowth textures found in plagioclase and orthopyroxene are corresponding to chemical changes that are not easily explained by changes in H 2 O content (Figures 3b and 7). Thus, we suggest that compositional mixing rather than changes in volatile content are recorded by mineral chemistry in the period immediately preceding eruption.
Similarities and differences in the preeruptive record of magmatic processes for the three Plinian eruptions are evident from the systematics of FeO t and An (Figure 3b). Fe partitioning in plagioclase is a function of melt composition and to a much lesser extent temperature (Bindeman et al., 1998). As Fe occurs in two oxidation states in most crustal magmas (Carmichael, 1991), a strong dependence of Fe partitioning on oxygen fugacity has been shown experimentally (Wilke & Behrens, 1999). The partition coefficient of Fe in plagioclase increases from 0.08 to 0.55 from reducing to oxidizing conditions in the range of log fO 2 = −15.78 to −7.27 (Wilke & Behrens, 1999), showing that Fe is mostly incorporated into plagioclase by the substitution of Al 3+ by Fe 3+ (Hofmeister & Rossman, 1984). Likewise, the partition coefficient of Eu in plagioclase increases from 0.095 to 1.81 in the same fO 2 range (Wilke & Behrens, 1999). This implies that if changes in fO 2 drive chemical variation in plagioclase, correlation between redox sensitive elements should be observed. Fe contents in the plagioclase from Plinian eruptions correlate well with Ti and Mg and not obviously with fO 2 -dependent Eu ( Figure S5), which indicates that Fe content in plagioclase from the three Plinian eruptions is mainly controlled by melt composition rather than oxygen fugacity.
Plagioclase core and mantle compositions from the three Plinian eruptions typically show a weak positive An-FeO correlation mostly between 28 and 55 mol. % An and 0.1 to 0.2 wt. % FeO (Figure 3d-3f). Individual crystals record cycles of both increasing and decreasing An content, frequently linked to prominent resorption textures, which is inconsistent with simple cooling but points toward an open dynamic system in which crystals experience repeated changes in intensive parameters. The relatively low slope in FeO content indicates that these changes are not controlled by interaction with a significantly more mafic melt but can rather be explained by changes in temperature or volatile content (Ginibre et al., 2002;Ruprecht & Wörner, 2007). Petrological experiments at 200 MPa on similar compositions show that changes in An content of the same order of magnitude as those observed in our study can be translated into temperature changes of 60°C or changes in xH 2 O fluid between 1 and 0.25 (Riker et al., 2015;Rutherford & Devine, 2008). The presence of similar groups of plagioclase showing An-FeO correlation in cores from the studied eruptions indicates that the magmatic systems of Nevado de Toluca was in a similar background state of activity prior to these Plinian events. Plagioclase crystallization occurred within a reservoir of rather evolved melt, likely highly crystallized, that was repeatedly disturbed by thermal perturbations and/or invasion of CO 2 -rich fluids (Caricchi et al., 2018), resulting in the oscillating An and Fe contents in crystal cores.
The An-FeO systematics reveals a further magmatic pattern prior to each Plinian eruption. Plagioclase mantles and outermost crystal rim compositions show elevated FeO contents after prominent resorption zones with steep increase of FeO at relatively high An content (Figures 3d-3f). Such pattern are consistent with 10.1029/2019JB017640 Journal of Geophysical Research: Solid Earth injection of more mafic magma into a dacitic magma reservoir (Ruprecht & Wörner, 2007), as also testified by opx chemistry. Interestingly, the distributions of An and FeO of the outermost plagioclase rims differ for the three eruptions. UTP and LTP rims are characterized by relatively homogeneous An content and a bimodal distribution of FeO content (Figures 6b and 6f). The differences in An and FeO for these two groups of rim compositions are not consistent with contamination by microinclusions of glass or Fe-Ti oxides ( Figures S5 and S11). Such distributions can be interpreted as resulting from magma homogenization prior to eruption with differing fractions of crystals that were in close proximity to the recharge magma, or reflect plagioclase growth from the different magmatic end-members that are evident at least in case of the LTP deposit, which is zoned from andesite to dacite (Arce et al., 2013). In contrast to these two eruptions, plagioclases from MTP show more heterogeneous An and FeO, which could reflect incomplete hybridization of potentially smaller mafic magma input into a heterogeneous reservoir at lower temperature.

Melt Evolution and Diversity
Trace elements in plagioclase show two different behaviors (Figures 7 and S6). Both Sr and Ti are positively correlated with An ( Figure S6). Except for the crystal cores, also, Mg is positively correlated with An, which can be accounted for by the partial diffusive reequilibration of Mg in the plagioclase cores ( Figure S6; Cherniak, 2010), before further mafic recharge resulted in growth of a new equilibrium rim. On the other hand, incompatible elements such as Ba and La correlate well among each other along the profiles but show decoupled behavior from the other trace elements (Figure 7). While the decrease of Sr, Ti, and Mg with An can be generally accounted for by the combined effects of crystallization of the phase assemblage and two end-member mixing of melts at different stages of chemical evolution, such an explanation is not consistent with the observed pattern for Ba and La. An alternative hypothesis to explain the different behavior of trace elements in plagioclase is that the recharge melts are chemically heterogeneous, which seem supported by the chemical heterogeneity of magma erupted by peripheral cones distributed around Nevado de Toluca (Figure 2b). Using An-dependent partition coefficients (Bindeman et al., 1998;Dohmen & Blundy, 2014), we recalculated the melt content of La, Ba, and Sr from their contents in plagioclase (Figure 8). Ba contents calculated using the Bindeman et al. (1998) partition coefficients are slightly lower than observed in UTP glass (Smith et al., 2009), while the values obtained using the model of Dohmen and Blundy (2014) are similar and higher than measured matrix glass contents ( Figure S7). Recalculated La and Sr contents are comparable to the matrix glass measurements (Smith et al., 2009) using both models. Low An contents are associated with low Ba melt contents for the UTP and MTP, while at high An the content of Ba in LTP is significantly lower than in plagioclases from UTP and MTP (Figure 8). Due to the incompatible behavior of Ba in the Nevado de Toluca magmas, such trends cannot be explained by melt evolution and are more consistent with plagioclase crystallization from magma produced by hybridization between a resident felsic magma and recharge magmas of distinct compositions. As illustrated by the wide spread in recalculated La-Ba and Sr-Ba melt contents, at least three general mixing components can be identified (Figures 8c and 8d), with the silicic end-member corresponding to the low  Figure S7.
Ba-La-Sr magma, toward which all compositions converge. High-Ba and low-La and low-Ba and high-La contents characterize the two mafic components. Hybridization of heterogeneous melts and convergence of these toward a common silicic end-member appears to be a viable process to explain trace element diversity in plagioclase and is consistent with the findings that monogenetic cone magmas near the central volcano are heterogeneous, while Nevado de Toluca produces relatively chemically homogeneous magmas. Interestingly, similar contrasts in magma geochemistry between monogenetic centers and the main edifice, as well as hybridization between these heterogeneous melts prior to Plinian eruptions, have been observed at Volcán Colima (Crummy et al., 2014), suggesting that the mixing of mafic magmas produced by different degrees of partial melting of the mantle with more felsic resident magmas occurs at regional scale.

Timescales of Preeruptive Processes
The diffusive relaxation of Mg in plagioclase and the Fe-Mg exchange in orthopyroxene can be used to calculate timescales of preeruptive crystal residence and interaction between magmas (Anderson et al., 2018;Chamberlain et al., 2014;Ruprecht & Cooper, 2012;Druitt et al., 2012;Costa et al., 2010). The textural complexity of plagioclase crystals from the Plinian eruptions of Nevado de Toluca indicates that the thermal and therefore diffusive history of these grains is complicated. However, diffusion timescales for plagioclase residence at magmatic conditions can be obtained from crystals that show a mismatch between observed and recalculated initial Mg profiles ( Figure 9).
Timescales of Mg relaxation in plagioclase were calculated for three normally zoned crystal cores and at the rim for one grain (Figure 9). Comparison between modeled curves and Mg content shows that the plagioclase cores resided at elevated temperatures for millennia (between 2 and 9 ka). As the thermal history of these grains is not known a priori, we used both the temperatures of preeruptive two-oxide equilibration (813-824°C; Arce et al., 2006) and the An-T relationship (777-784°C; Figure S8), constrained from outermost rim plagioclase composition and thermometry as input (Arce et al., 2013. The plagioclase core profiles are normally zoned with only minor oscillatory zoning or resorption horizons and decreasing anorthite content toward their outer portions. This part of the profile represents one major differentiation cycle prior to rejuvenation, which is marked by sharp boundaries and increase in An content in the outer portions of the crystals (Figure 7). We suggest that our calculated core diffusion timescales provide an order of magnitude estimate for the duration of a magma differentiation cycle prior to rejuvenation of the magmatic system by injection of more mafic magma. Mg contents in the crystal rims show good correspondence between measured and calculated initial values indicating that these were not much impacted by diffusion as already shown by Smith et al. (2009). One single profile could be modeled by diffusion over decades to centuries between the formation of the chemical boundary and eruption ( Figure 9c). Diffusion modeling on Toluca plagioclase is challenging, as most crystal cores are either fully equilibrated or show reverse zoning in the core with no observed mismatch between equilibrium and calculated initial concentrations. Nevertheless, the calculated diffusion timescales suggest that Nevado de Toluca is underlain by a long-lived mushy magmatic system.
Opx diffusion timescales are summarized in Figure 10c using two different formulations for the diffusion coefficient and preeruptive temperatures as constrained by two-oxide thermometry (Dohmen et al., 2016;Allan et al., 2018, Allan et al., 2013. Yearly to decadal time estimates were obtained by using a corrected version (Allan et al., 2018;Fabbro et al., 2018) of the formulation presented in Allan et al. (2013). Diffusion timescales modeled using the Dohmen et al. (2016) equation for diffusion coefficient are systematically longer. It is not clear which formulation can be considered to be more realistic in our case. Independently of the formulation considered, the duration of preeruptive processes differs for the different eruptions. MTP diffusion timescales are typically longer compared to the UTP, indicating that more time passed between the magma input event that caused the drastic compositional zoning in orthopyroxene and eruption. The fastest diffusion timescales were obtained from two LTP crystals with 1.1-3.3 and 1.5-2 years. As discussed in the previous section, the large compositional boundaries in orthopyroxene likely reflect magma mixing events. Magma mixing is more consistent with the range of diffusion timescales we obtain, compared to entrainment opx xenocrysts from the wall rock, as the latter can be expected to operate continuously over much longer timescales.
Uncertainties from temperature variations in the opx diffusion modeling are accounted for by running models within the temperature range given by Fe-Ti oxide thermometry ( Fresh incoming magma will typically heat up the resident magma depending on the mass and temperature of the recharge magma. Amphibole crystals from all studied eruptions show an increase in rim ward temperature, yielding final outermost rim temperatures, which are 10-40°C higher but similar to the Fe-Ti oxide constraints within error. Thus, our opx timescales represent maximum values. The differences in diffusion times we find can result either by injection of magma in a thermally zoned magma reservoir or by different distances between the locus of injection and the crystals on which the profiles were collected. Plagioclase and opx diffusion modeling suggest that injection of mafic to intermediate magma into the upper crustal silicic storage region beneath Nevado de Toluca triggered the Plinian eruptions and occurred between a few years and up to several decades before.

Implications for Volcano Monitoring
Extrapolating petrological data of ancient or historic eruptions can be useful to anticipate the potential future behavior of volcanoes (Weber & Castro, 2017;Stock et al., 2017;Muir et al., 2014;Singer et al., 2014). Geochemical patterns in plagioclase and orthopyroxene crystals record recurrent cycles of fresh magma injection and differentiation prior to each of the last three Plinian eruptions of Nevado de Toluca. XMg versus distance (μm) for crystals shown in (a). White dots are measured profiles. Dashed lines indicate initial profiles and blue lines correspond to best fit diffusion model. (c) Ranked time in years of the modeled crystal boundaries is shown for the three studied eruptions. Each vertical pair with darker and brighter colors corresponds to the same modeled profile using different formulation of the diffusion coefficient. The lower time estimates (brighter colors) were obtained using Allan et al. (2013Allan et al. ( , 2018, and greater time for the same crystal (darker colors) were modeled using the equation of Dohmen et al. (2016). Modeled crystal boundaries were interfaces between group one and group two orthopyroxene (LTP and UTP) and between group one and group three for the MTP samples. Two crystals with variation in En content within group one were modeled for the UTP eruption, which are marked by arrows.
Inter-Plinian Block-and-Ash flows are compositionally very similar to the Plinian eruptions and are likely sourced from the same plumbing system as the Plinian events . The consistency in preeruptive processes and timescales between these events, separated in time by millennia (Figure 1), would suggest that similar magmatic events will precede a future Plinian eruption of Nevado de Toluca.
The surprisingly mafic orthopyroxene compositions found in all studied eruptions stress the important role of deep magmatic processes on the reactivation of the subvolcanic magma reservoir. Magma injection at deep crustal level could be indicated by deep seismicity and broad deformation pattern surrounding a volcano and registered by satellite or ground-based geodetic techniques such as InSAR, leveling, or GPS (e.g., Magee et al., 2018). Numerical modeling shows that such broad deformation patterns are mainly controlled by chamber depth, crustal heterogeneity, and rheology (Gottsmann & Odbert, 2014;Hautmann et al., 2010). While such pattern have been associated with the eruption of mafic magmas (Fedotov et al., 2010;Sturkell et al., 2013), we emphasize that at Nevado de Toluca signs of deep crustal unrest might herald potential Plinian eruptions. Given the volcano's current long dormancy of over 3 ka, which is comparable to the order of magnitude of timescales recorded in plagioclase core diffusion, it is likely that renewed magma input into the shallow magma reservoir is required for Plinian or other eruptive activity to resume. Extended periods of dormancy lasting many centuries to millennia are frequently observed at Arc volcanoes (Hildreth, 2007), which can hamper recognition of the volcanic hazard potential by authorities and the local population. Our data show that Nevado de Toluca is a long-lived mushy system that can transition into a critical state on sufficiently short timescales relevant for volcanic hazard assessment.
Considering the potential effect of magma input on magma reservoir pressurization and failure, we calculated the approximate volume of magma to be injected to generate sufficient overpressure to trigger a Plinian eruption of Nevado de Toluca (Figure 11; Cañon-Tapia, 2014). The critical pressure at which failure of a reservoir occurs can be calculated as a function of material properties such as tensile strength and Poisson's ratio of the rock, which was varied between 0.5-6 MPa and 0.25-0.5, respectively (Gudmundsson, 2011;Perras & Diederichs, 2014). Overpressure generated by magma injection was modeled for several volume ratios between input (ΔV) and reservoir volume (V), using a bulk modulus of magma between 1.15 · 10 10 and 2 · 10 10 Pa and the pressure derivative between 4 and 7 (Canon-Tapia, 2014). ΔV/ Figure 11. Comparison of overpressure generated due to magma injection and the depth range of the silicic Nevado de Toluca reservoir. Critical pressures (black lines) for the rupture of a magma chamber and pressure inside a magma chamber due to injection of more magma were calculated as function of depth, using the analytical equations of Canon-Tapia (2014). Dashed curves show the pressure due to magma injection for two different tensile strength of the crust (red: 0.5 MPa; tan: 6 MPa) and are labeled for the ratio of volume change to magma chamber volume (ΔV/V). The gray shaded area marks the depth range of the silicic magma reservoir that fed Plinian eruptions from Nevado de Toluca (Arce et al., 2013. V ratios between 0.002 and 0.01 generate sufficient overpressure to induce failure of a magma reservoir at depth corresponding to the feeding reservoir of the Plinian eruptions of Nevado de Toluca. Combining these volume fractions with the erupted magma volumes for the three eruptions (Figure 1), recharge volumes of 0.0017-0.0085 km 3 for the LTP, 0.0036-0.018 km 3 for the MTP, and 0.016-0.08 km 3 are required to trigger the UTP eruption by magma injection. Small volcanic eruptions, occurring with recurrence intervals of days to weeks on Earth, produce bulk volumes very similar to the calculated recharge trigger threshold between 0.001 and 0.01 km 3 (Pyle, 2015), indicating that even a relatively small injection of magma in the shallow reservoir may be sufficient to trigger Plinian events at Nevado de Toluca. Using median opx diffusion timescales for the UTP and LTP of 18 years, and 44 years for the MTP yield magma injection rates between 9.44 × 10 −5 km 3 /year and 8.89 × 10 −4 km 3 /year required which are likely episodic rather than continuous injections, as such low rates might be not sustainable due to conduit freezing, to reach the critical volumes and trigger the eruptions. The estimates volumes and injection rates are minima, as our calculations do not consider inelastic effects or magma compressibility (e.g., Kilbride et al., 2016). Nevertheless, volume changes of this magnitude measured during unrest are resolvable from modeling of volcano deformation data derived by InSAR (Jay et al., 2014) or GPS (Sigmundsson et al., 2010) and should be considered as potential indicator of a future Plinian eruption, especially if following the signs of deep crustal magmatic activity.

Conclusions
Detailed chemical and textural analyses of plagioclase, amphibole, and orthopyroxene crystals for three major Plinian eruptions of Nevado de Toluca volcano in Central Mexico reveal a recurring pattern of magmatic processes prior to these events. The major and trace element zoning in plagioclase is consistent with cycles of magma differentiation separated by injection of magma from depth. Recharge magmas are chemically heterogeneous basalts or basaltic andesite judging from the compositions of opx cores that partly equilibrated with the host dacite prior to each of the Plinian eruptions. Distinct groups of opx compositions together with amphibole barometry further reveal that the plumbing system of Nevado de Toluca consists of a complex system of separated magma storage possibly located at different crustal levels. Plinian eruptions are preceded by homogenization processes. Diffusion modeling of Mg in plagioclase indicates that differentiation cycles between input of mafic magma at shallow depth occur over millennia, while Fe-Mg interdiffusion in orthopyroxene constraints the timing of preeruptive magma injection in the range of years to less than two centuries. Based on our petrological results and modeling, we conclude that future voluminous explosive eruptions of Nevado de Toluca could potentially be foreseen by geophysical monitoring. the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (Grant agreement 677493-FEVER). We are grateful for the logistic support offered by UNAM during the field campaigns. All data supporting the analysis and conclusions are available in electronic form in the online supporting information file. We thank Stephen Parman for editorial handling. Comments by Phillip Ruprecht and an anonymous reviewer are appreciated as they helped to significantly improve the quality of this contribution. The data used in this study are available in the EarthChem data repository DOI: 10.1594/IEDA/111395, URL: https:// doi.org/10.1594/IEDA/111395.