The Long-Term Life-Cycle of Nevado de Toluca Volcano (Mexico): Insights Into the Origin of Petrologic Modes

The petrologic diversity of volcanic rocks reflects the dynamics of magma reservoirs and the temporal evolution of magma chemistry can provide valuable information for hazard assessment. While some stratovolcanoes monotonously produce intermediate magmas (55–68 wt% SiO2), dominantly erupted magma types (e.g., basaltic andesite, andesite or dacite) frequently differ even between neighboring volcanoes. If such differences arise due to thermal maturation processes over time or are predetermined by other properties of magmatic systems remains poorly understood. This study helps to elucidate the underlying factors modulating the chemistry of the magma preferentially erupted by Nevado de Toluca volcano in Central Mexico. We present a new dataset of bulk-rock and mineral chemistry spanning the entire 1.5 Million years of the volcanos’ eruptive history. The results reveal that Nevado de Toluca dacites and minor andesite originate in a stable configuration of pre-eruptive processes and plumbing system architecture by hybridization between an upper crustal silicic mush and deeper sourced basaltic andesite magmas. Yet, a subtle trend toward increasing silica content with time (2 wt% in 1.5 Ma) and episodicity in magma hybridization conditions are observed. We use thermal simulations of pulsed magma injection to probe the controlling variables on the temporal variation and compositional mode of magma geochemistry. The results show that the subtle temporal trend toward increasing bulk-rock SiO2 content is plausibly explained by slightly dropping recharge rates and continued upper crustal reservoir growth. Our modeling also shows that the dominant composition of eruptible magmas (“petrologic mode”) can shift as a function of magma flux, extrusive:intrusive ratio and temperature of the recharge magma. A comparison of SiO2 whole rock distributions for monotonous Mexican stratovolcanoes and their peripheral cones shows that their petrologic modes vary in concert, indicating that the recharge magma chemistry or temperature is a major control on the preferentially erupted magma composition for these volcanoes.

Differences in the "petrologic mode" occur even between neighboring volcanoes (Hildreth, 2007;Klemetti and Grunder, 2008;Kent et al., 2010;Macias et al., 2017), which indicates that large-scale tectonic parameters (e.g., age or slope of the subducting slab) or crustal thickness are not sufficient alone to explain these variations. For example, sharp differences in the petrologic mode exist between neighboring systems along the Andean chain (e.g., Calbuco-Osorno, El Misti-Ubinas) or in the Trans Mexican Volcanic Belt (e.g., Nevado de Toluca-Popocatépetl). Most agree that local modulations in magma flux, crustal and/or magma properties are some of the likely variables governing the preferential eruption of specific compositions (Carmichael, 2002;Klemetti and Grunder, 2008;Caricchi and Blundy, 2015a;Wörner et al., 2018;Till et al., 2019). The high viscosities of silica-rich magmas and the high densities of basaltic melts may act as physical barriers, preventing magma transport and inhibit volcanic eruptions of such compositions (Stolper and Walker, 1980;Marsh, 1981;Pinel and Jaupart, 2000;Moran et al., 2011). If such barriers are effective in preventing the eruption of very mafic and silicic magma compositions, they may be overcome by mixing/mingling processes induced by mafic magma recharge in crustal mush reservoirs and produce compositionally monotonous intermediate volcanoes (Reubi and Blundy, 2009;Kent et al., 2010;Kent, 2014). Alternatively, magma properties such as the initial water content or nonlinearities in melt fraction-temperature relations can modulate the variety of compositions present in a magmatic system (Melekhova et al., 2013;Nandedkar et al., 2014;Caricchi and Blundy, 2015b;Hartung et al., 2019;Huber et al., 2019).
Thermal modeling of incremental magma reservoir assembly shows that magmatic systems evolve thermally in time with large temperature fluctuations in early stages of temperatures later in the life cycle (i.e., the evolutionary history) of the magmatic system (Petford and Gallagher, 2001;Annen and Sparks, 2002;Dufek and Bergantz, 2005;Annen et al., 2006;Annen, 2009;Gelman et al., 2013;Melekova et al., 2013;Karakas et al., 2017). Thus, progressively more homogeneous magma compositions may be erupted with the progressive thermal maturation of volcanic plumbing systems. Likewise, mechanical modeling predicts an increase of the duration of magma storage during the lifetime of volcanic systems because of the more efficient relaxation of overpressures in thermally primed wall-rocks (Jellinek and DePaolo, 2003;Gregg et al., 2012;de Silva and Gregg, 2014;Karlstrom et al., 2017). This could favor magma homogenization in the plumbing system and lead to more homogeneous erupted compositions. Thermal and chemical heterogeneities in magma reservoirs may decay with time by convective overturn in association with fresh recharge pulses and/or gas exsolution (Ruprecht et al., 2008;Huber et al., 2009, Burgisser andBergantz, 2011;Huber et al., 2012). In highly crystalline reservoirs, melt segregation or reactive flow through vertically extensive crystal piles can occur (Rabionwicz and Vigneresse, 2004;Solano et al., 2012;Cooper et al., 2016;Holness, 2018;Jackson et al., 2018;Bachmann and Huber, 2019;Floess et al., 2019), which may be important in buffering melts to specific compositions governed by the chemistry of cumulate rocks.
Predictions of existing numerical models, relating magmatic processes and the spectrum of chemical compositions erupted by volcanoes can be tested based on long-term eruptive records, but can also help to illuminate physical controls on geochemical distributions and temporal trends. Here we present a dataset of 52 new and 45 previously published , Weber et al., 2020 whole-rock analyses and 6247 new spot analyses on minerals of temporally resolved bulk-rock major and trace element analyses and mineral chemistry spanning the 1.5 Ma life cycle of the dacitic Nevado de Toluca volcano in Central Mexico. We use mineral zonation systematics and thermo-barometric calculations to constrain the origin of magmas at Nevado de Toluca and investigate the temporal evolution of bulk eruptive products and their crystal cargos over the entire lifespan of the volcano. We finally discuss the origin of the observed longterm temporal trends in magma chemistry, as well as the origin of petrologic modes of Mexican stratovolcanoes, based on our dataset, existing data and using thermal modeling to provide a quantitative framework to our interpretations.

GEOLOGICAL BACKGROUND
Nevado de Toluca, also known as Xinantécatl (19°06′30″N; 99°45′30″W; 4,680 m above sea level), is a long-lived stratovolcano located in the central Mexican highland about 80 km southwest of Mexico City (Macías, 2007). It is part of the Trans Mexican Volcanic Belt, a 1,000 km long east-west trending volcanic arc of Miocene to Holocene age that developed in response to the subduction of the oceanic Cocos and Rivera plates beneath the North American continental lithosphere at different angles ( Figure 1; Pardo and Suárez, 1995;Gómez-Tuena et al., 2007a, Gómez-Tuena et al., 2007b, Gómez-Tuena et al., 2018Ferrari et al., 2012).
The volcanic edifice was constructed from various silicic lava flows and domes with intercalated pyroclastic deposits, with the most prominent feature of the volcanos' morphology being a 2.5 × 1.5 km horseshoe shaped amphitheater that opens toward the East in the direction of the City of Toluca (Norini et al., 2004). Geological mapping and dating have shown that volcanism at Nevado de Toluca started in the Early Pleistocene, spanning 1.5 Ma of continuous explosive and effusive activity since then (Garcia-Palomo et al., 2002;Torres-Orozco et al., 2017). The volcanos' history can be divided into several stages based on geomorphology and age relations (Figures 2, 3A;Torres-Orozco et al., 2017). During the early "Old Nevado" stage (1.5-0.13 Ma), lava flows have been erupted to the south and north-east of the current crater area. This early stage is coeval with the extrusion of silicic Domes ("Dome stage") between 1.5 and 0.27 Ma in the periphery of the main volcano. Only a few monogenetic systems have been dated in the area (Bloomfield, 1975;Arce et al., 2013a;Torres-Orozco et al., 2017), but available data point toward a long age range between at least 0.86 Ma to 8.5 ka. The edifice of Nevado de Toluca collapsed several times, which is evident in widely dispersed debris avalanche deposits around the volcano (Capra and Macias, 2000;Caballero and Capra, 2011). The younger eruptive history can be divided into 57-9 ka "Recent Nevado" stage, mostly characterized by silicic lava effusion in the crater area, and a sequence of pyroclastic deposits "Young PD stage" spanning ca. 42-3 ka (Macías et al., 1997;Arce et al., 2006). The young pyroclastic deposits ( Figure 3B) record a complex series of dome destruction events that are preserved in block and ash flow deposits (Garcia-Palomo et al., 2002), as well as by a sequence of at least three Plinian eruptions (Arce et al., 2003;Arce et al., 2005;Arce et al., 2006;Capra et al., 2006).
Geochemical studies, mostly focused on the Plinian deposits, have shown that Nevado de Toluca predominantly erupts subalkaline dacite and minor andesites with calc-alkaline affinity typical for subduction related magmatism Torres-Orozco et al., 2017). Experimental constraints on magma storage pressures range between 150 and 300 MPa at water saturated conditions and eruptive temperatures are typically in excess of 820°C Arce et al., 2013b). Isotopic (Sr, Nd, Pb) and trace element data shows that even though Nevado de Toluca sits on thick continental crust of about 50 km, crustal melting plays only a marginal role in silicic magma genesis (Martínez-Serrano et al., 2004). The shallow part of the magmatic system prior to the Plinian events has evolved as an open system, where fresh pulses of recharge magma enter the preexisting more silicic reservoirs (Smith et al., 2009;Weber et al., 2019).

Bulk Rock Analytical Techniques
Samples were collected in the field to represent the entire eruptive history of the volcano, with particular emphasis on units that have previously been dated by radiocarbon or 40 Ar/ 39 Ar-geochronology (Macías et al., 1997;Torres Orozco et al., 2017). A description of the sampled units and locations, including coordinates, can be found in Supplementary Table 1. All samples were first cleaned from adherent alteration patina, washed with water and dried in an oven at 50°C before further processing. Lava and pumice samples were then reduced in size using a steel jaw-crusher, a hydraulic press and sieved to grain sizes <2 mm. Lava samples were ground to fine powders in an agate mill for 30 min, while pumice samples, due to their higher porosity, were milled for 20 min. In order to monitor the alteration state of the samples, 2 g aliquots were separated from the powders to determine the loss on ignition (LOI). The aliquots were weighted into ceramic crucibles and heated to 1,050°C for 24 h. The LOI was then calculated by weight difference before and after heating. Subsequently, glass beads for X-ray fluorescence analysis (XRF) were prepared. Precisely 6 g of Li tetraborate were mixed for 3 min in a hand mortar with 1.2 g of each of the powdered samples. The so obtained mixtures were then loaded in Pt crucibles and fused at 1,200°C. Major element analysis were carried out at the University of Lausanne using a PANanalytical Axios max X-ray fluorescence spectrometer. BHVO-2 and several in-house standard materials were used for quality control during the analytical session.
Trace element contents in the whole rock samples were obtained on the same fused glass beads by laser ablation inductively coupled plasma mass spectrometry. After XRF analysis, glass beads were crushed and mounted on adhesive tape with the side of the bead not analyzed by XRF facing to the top. We used a quadrupole spectrometer Agilent 7700x coupled to an UP-193FX ArF eximer ablation system at the University of Lausanne. The instrument was optimized by ablating the NIST SRM 612 reference glass in linear scan mode with 20 Hz repetition rate, 75 µm beam size, and 5.0 J cm −2 on-sample energy density. Ablation of the sample glass discs was carried out under a helium atmosphere, using a 12 Hz repetition rate, 5.0 J cm −2 energy density and 150 µm diameter beam size. Relative sensitivity factors were determined by analysis of NIST SRM 612 during the analytical session. On each glass bead, we collected and averaged 3 spot analyses. Data were reduced with the MATLAB based software SILLS (Guillong et al., 2008), using the bulk-rock Ca contents from XRF analysis for quantification. The relative uncertainties (RSE) for the 3-point analysis are for most analyzed elements <2%, with Rare Earth Elements (REE) generally showing RSE <5%. The largest differences for averaged spots were obtained on Mo (RSE: 12.8%), Cu (7.4%) and Ni (6.5%).

In-situ Mineral and Glass Analysis
In this study we present mineral and glass chemistry for 19 eruptions from Nevado de Toluca, which were selected to cover the entire eruptive history of the volcano. Mineral, groundmass and glass compositions were determined on polished 100 µm thick petrographic sections by electron probe microanalyzer. All analyses considered here were carried out, using a JEOL 8200 superprobe at the University of Geneva. Prior to analysis, the petrographic sections were coated with a 20 nm thick layer of carbon to ensure electrical conductivity. For each of the studied eruptions, we analyzed one petrographic section, for which between 8 and 15 BSE images each were obtained in order to reveal textures and internal mineral zonation features. Typically, we analyzed linear rim to core profiles of 10-20 crystals for each mineral phase in the individual petrographic sections. The spacing between points was variable but between 5 and 20 µm for most grains.
Mineral phases (plagioclase, pyroxenes, Fe-Ti oxides and amphiboles) were analyzed using the wavelength-dispersive x-ray spectrometer of the electron probe microanalyzer. Measuring conditions were set to 15 kV acceleration voltage, 20 nA beam current and focused beam size of 2 µm. The peak counting times were set to were 30 s (Si kα TAP), 20 s (K kα PETH), 20 s (Na kα TAP), 30 s (Fe kα LIFH), 30 s (Al kα TAP), 30 s (Ca kα PETJ), 30 s (Ti kα LIFH), 30 s (Mg kα TAP), 30 s (Cr kα PETJ), and 30 s (Mn kα LIFH). Raw analysis for plagioclase and pyroxenes that showed totals outside of the range 100 ± 1.5% were rejected from further analysis. The raw amphibole totals outside the range 95.5-99.5% were filtered out. Compositions that indicated mixed analysis of silicate minerals and glass were also discarded. In total, 2,883 plagioclase, 2,232 pyroxene, and 1,128 amphibole analyses are presented in this study. Matrix and melt inclusion glasses, as well as groundmass analyses, were carried out on the same petrographic sections as the mineral analyses and we used a reduced beam current of 8 nA and 15 kV acceleration voltage. The beam size was set to 10 µm for matrix glass and melt inclusion analysis, while for the mostly microcrystalline groundmass a beam diameter of ∼50 µm was used. In each section, typically 15 groundmass analyses were carried out and averaged. Melt inclusions were hosted in either pyroxene or amphibole. Counting times on the peak position were 30 s (Na kα TAP), 20 s (Si kα TAP), 60 s (S kα PETH), 20 s (K kα PETJ), 30 s (Mg kα TAP), 30 s (Ca kα LIFH), 30 s (Al kα TAP), 30 s (Ti kα PETJ), and 30 s (Fe kα LIFH). Alkali element measurements were performed first in the sequence to mitigate diffusive loss.

Thermobarometric Calculations
The mineral, glass and groundmass compositions analyzed in this study were used to estimate the pressures, temperatures, and oxygen fugacity of Nevado de Toluca magmas using a range of geothermobarometric techniques. Amphibole crystallization temperatures and pressures were calculated for euhedral outermost amphibole rim compositions using Equation 7b from Putirka (2016) by iterative solution, which requires input of amphibole and liquid composition, as well as an estimate of water content. We used average groundmass for each of the studied eruptions to approximate the input liquid composition. In order to assess equilibrium between the assumed liquid and outermost rim amphibole, we discarded amphibole-liquid pairs with mineral-melt exchange coefficient K D (Fe t -Mg) outside the range 0.17 and 0.39 (Putirka, 2016). 39% of the analyzed amphibole rims were found to be in equilibrium. The impact of the assumed liquid H 2 O content was assessed by performing calculations with 4, 6, and 8 wt% H 2 O. However, the resulting average differences in pressures and temperature calculations of 0.4 kbar and 0.3°C are much smaller than the uncertainties of amphibole thermobarometry (±4 kbar and ±30°C for single spot analyses; Putirka, 2016). Amphibole crystallization temperatures can also be estimated independent of liquid composition and pressure, not necessarily with worse precision and accuracy when compared to amphibole-liquid thermometers (Ridolfi and Renzulli, 2012;Putirka, 2016). We calculated liquid-independent amphibole temperatures using Eqn. 5 presented in Putirka (2016) on outermost rim compositions in order to compare them to the amphibole-liquid pair estimates. This thermometer has been calibrated over a temperature range applicable to intermediate magmas and produces calibration data with standard error of estimate (SEE) of ±53°C (Putirka, 2016). Plagioclase-amphibole temperatures were calculated using the calibration based on the edenite + albite richterite + anorthite exchange reaction on analyses corresponding to average outermost rim composition within 2SD in samples that showed euhedral crystal textures (Holland and Blundy, 1994). Several samples contain both clino-and orthopyroxenes, but equilibrium between outermost rim compositions based on a instead of an Kd (Fe-Mg) of 1.09 ± 0.14 was only achieved in the monogenetic cone samples (Tenango, Sabanillas). We thus calculated 2-pyroxene temperatures only on these samples, using Eq. 36 of Putirka (2008). The composition of coexisting titanomagnetite and ilmenite was used to constrain temperatures and the oxygen fugacity for various eruptions of Nevado de Toluca. We used the model of Ghiorso and Evans (2008) to calculate the temperatures of Fe-Ti exchange between touching equilibrium pairs (Bacon and Hirschmann, 1988) of titanomagnetite and ilmenite. Oxygen fugacity is reported as values relative to Nickel-Nickel-Oxide buffer.

Thermal Modeling
We used heat conduction modeling, in order to simulate the temporal evolution of temperatures in crustal rocks subjected to repeated magma injection. An explicit finite difference scheme was used to solve the axisymmetric formulation of the heat equation, which can be written as: Frontiers in Earth Science | www.frontiersin.org September 2020 | Volume 8 | Article 563303 6 where t is the time, T is the temperature, z is the vertical coordinate, r is the radial distance from the axis of rotational symmetry, k is the thermal conductivity, L is the latent heat of crystallization, ρ is the density, c is the specific heat and X c is the fraction of crystals. As we do not aim to model the particular phase equilibria of Nevado de Toluca in this study, but to investigate the general controls on compositional modes, we used a hypothetical linear melt fraction temperature relation in the presented model calculations. While the specifics of melt-fraction temperature relations impact on compositional distributions (Caricchi and Blundy, 2015a), it does not impact the direction of compositional change as a function of magma recharge or extraction rate (Weber et al., 2020) and is therefore considered of secondary importance here. A temperature dependent thermal conductivity (k) was implemented in the model based on the experiments of Whittington et al. (2009) for average crust. Zero heat flux was imposed in all simulations in the direction perpendicular to all lateral boundaries, apart from the surface, where temperature was fixed at 8°C. The numerical code used in this study has been benchmarked against existing models (Caricchi et al., 2014;Caricchi et al., 2016), which used a different numerical method to solve heat conduction problems.
Magma injection was modeled at an initial crustal intrusion depth of 5 km by accretion of basaltic sills, intruded at their liquidus temperature of 1,170°C, which causes downward advection of crustal rocks. All presented simulations were run with a linear initial geothermal gradient of 25°C/km. In order to investigate the impact of changing extrusive:intrusive (E:I) ratio of a magmatic systems compositional diversity, we implemented heat and mass removal in the model. In such simulations, sills were removed at a rate that is equivalent to E:I ratio of 0.4 and 0.8 of the final intrusion volume built by a specific injection rate. This was accomplished by removing sills matching the previous injection site and dimension, which ensures removal of magma at the highest temperature in the reservoir, by advecting the temperature field upwards. We performed petrological modeling, in order to invert the computed temperature distributions into potentially eruptible magma chemistry. All petrologic calculations were carried out using the hypothetical liquid line of descent presented in Supplementary Table 7. For a  (Sun and McDonough, 1989). (D) Chondrite normalized Rare Earth Elements (McDonough and Sun, 1995).
Frontiers in Earth Science | www.frontiersin.org September 2020 | Volume 8 | Article 563303 7 detailed description of the calculation procedure, the reader is referred to Weber et al. (2020).

Bulk-Rock and Glass Geochemistry
We present and analyze a combined dataset of 97 bulk rock major and trace element analysis for the long-term eruptive history of Nevado de Toluca volcano. In agreement with previous studies of the long-term history of the volcano, the major element geochemistry is focused on dacite throughout the history of the volcano with an average SiO 2 content of 63.69 (±2.68 1SD) wt% SiO 2 ( Figure 4A). Andesites are less frequently erupted from Nevado de Toluca but occur in all eruptive stages. The most mafic erupted bulk compositions with SiO 2 contents around 55 wt% are observed in monogenetic cones, extending into the basaltic andesite field, also in agreement Frontiers in Earth Science | www.frontiersin.org September 2020 | Volume 8 | Article 563303 8 with previous work on peripheral cones surrounding Nevado de Toluca (Martínez-Serrano et al., 2004;Torres-Orozco et al., 2017). Samples erupted from the main volcano show typical trends of increasing incompatible major element components (e.g., K 2 O, Figure 4A) and decreasing trends of compatible elements, when plotted against an index of differentiation (e.g., MgO or SiO 2 ). Different eruptive stages show slight differences in slope and variance in bivariate plots of MgO vs. FeO total ( Figure 4B). Melt inclusions, hosted in amphibole and orthopyroxene, are consistently rhyolitic in composition with average silica contents of 75.20 (±1.63 1SD).
Primitive mantle normalized trace element pattern of all analyzed samples show typical features of subduction related magmatism such as negative Nb-Ta anomalies and enrichment of large ion lithophile elements (LILE) like Cs, Ba or Sr ( Figure 4C). Rare Earth Element contents in the bulk-rock samples normalized to chondritic compositions are enriched in the light REE over heavy REE. Based on trace element variation pattern, three compositional groups can be identified (Figures 4C,D): 1) A main calc-alkaline suite, which comprises >90% of the analyzed samples. This group contains samples from all eruptive stages, including some monogenetic cones. 2) A subset of monogenetic cones, however, shows enrichment of light REE, less pronounced Nb-Ta negative anomalies, and greater enrichment in Sr and Ba relative to the dominant trace element group. Chemical differences between these samples and the main volcano are also evident in major elements, in particular higher K 2 O contents ( Figure 4A). 3) The third trace element group is characterized by depletion relative to group 1 of mid to heavy REE elements (Sm-Lu) with the exception of Eu. These distinct compositions are only evident in three samples from silicic Domes surrounding the main volcano and are consistent with an important role of amphibole fractionation during magma evolution (Davidson et al., 2007) before plagioclase crystallization.

Petrography
We collected petrographic images and determined the modes of macro-crystal and groundmass phases in 17 eruptions from Nevado de Toluca volcano and two peripheral cones, spanning its entire eruptive history of 1.5 Ma ( Figure 3C; Supplementary Figure 1). A representative selection of BSE images ( Figures 5A-F) shows that the predominant phase assemblages are comprised of plagioclase, orthopyroxene, amphibole and Fe-Ti oxides. Occasionally, quartz, clinopyroxene and biotite are observed in the eruptive products. The mineral textures are complex but similar in all studied eruptions from the main volcano. Plagioclase textures are characterized by multiple resorption and overgrowth zones, while pyroxenes and amphiboles show strong, but generally binary core to rim zonation. Reaction coronas surrounding amphibole are observed in some petrographic sections ( Figure 5E) and are particularly evident in two of the early history samples ("U5 Raíces," "U14 Tepehuisco"), which lack strongly zoned orthopyroxenes. The analyzed monogenetic eruptive products contain two-pyroxene macro-crystal assemblages and a groundmass dominated by plagioclase microlites. Crystallinities of the eruptive products from the main volcano range from 20 to 45 vol% with average of 36 vol%, while crystal contents of the monogenetic eruptions are below 12 vol% (Supplementary Figure 1). Average fractions of phenocrystals are 0.7 for plagioclase, 0.14 for amphiboles and 0.08 for orthopyroxenes. Samples from the early eruptive stages (i.e., "Old Nevado" and "Domes") typically contain clinopyroxene more frequently compared to the younger eruptive units (Supplementary Figure 1).

Mineral Chemistry
Major and minor elements in silicate and oxide mineral phases from 19 eruptions have been analyzed. Anorthite contents in plagioclase crystals range from 13 to 69 mol% with average value of 40 mol% and inter quartile range (IQR) between 34 and 46 mol%. FeO total contents in plagioclase show a large range between 0.05 and 0.85 wt% and are positively-skewed with average of 0.20 wt% and IQR between 0.14 and 0.31 wt%. Anorthite contents for the different eruptive stages of the main volcano vary only slightly with maximum difference of 2 mol% for average values and 3 mol% for IQRs. Greater variation between the eruptive stages of Nevado de Toluca is evident in FeO total contents, which show average values of 0.31 and 0.22 for the "Old Nevado" and "Domes" stages, while the younger history ("Young PD" and "Recent Nevado") units have lower average in FeO total contents of 0.17 and 0.15 wt%, respectively.
Amphibole crystals in the studied eruptions are Ca-rich and can be classified as magnesio-hornblende and pargasite. Individual eruptions show large differences in amphibole compositions both in central tendency and variance for several elements (Supplementary Figure 2). For example, Al T per formula unit (p.f.u.) varies between 1.01 and 1.95 with average of 1.41. Much of this range is for instance captured by amphibole from the 31 ka (cal. BP) BAF28 eruption (average: 1.39 Al T p.f.u., IQR: 1.31-1.41 Al T p.f.u.), while the ca. 15 ka earlier Pink Pumice eruption contains amphiboles of rather restricted composition with average of 1.15 and IQR of 1.12-1.20 Al T p.f.u. In accordance with zonation pattern in BSE images ( Figures 5A-E), chemical compositions of Nevado de Toluca amphiboles are often bimodal with different compositions of crystal cores and rims (Supplementary Figure 2). Systematic differences between the eruptive stages and amphibole compositions are, however, not evident.
Orthopyroxenes are present in most of the studied eruptions and range between 39 and 88 mol% in enstatite (En) and from below detection limit to 0.90 wt% in Cr 2 O 3 content. Average En and Cr 2 O 3 contents in the monogenetic cone samples of 84 mol% and 0.25 wt%, respectively, are higher compared to eruptions from the main volcano with average En of 66 mol% and Cr 2 O 3 of 0.1 wt%. However, high En and Cr 2 O 3 in heavily zoned crystal cores (c.f., Figure 5D) are a characteristic feature of almost all sampled silicic eruptions from the main volcano as well. Clinopyroxenes occur in 9 of the 19 studied eruptions and in all eruptive stages except the "Young PD" sequence (Supplementary Figure 1). Cpx compositions cluster around average values of 45 mol% and IQR of 41-49 mol% in En. All cpx can be classified as augitic with the exception of the "Tepehuisco Frontiers in Earth Science | www.frontiersin.org September 2020 | Volume 8 | Article 563303 9 Dome" eruption that contains diopsitic cpx. Similar to opx, cpx crystal show strong compositional zoning with Cr 2 O 3 contents of up to 1.17 wt%.

Petrogenesis and Deep Structure of the Plumbing System
Intermediate and silicic arc magmas are widely believed to originate by a combination of processes including fractional crystallization (Grove et al., 2005;Lee and Bachmann, 2014;Müntener and Ulmer, 2018), reactive melt flow through porous cumulates (Solano et al., 2012;Cooper et al., 2016;Jackson et al., 2018), assimilation of crustal rocks (Hildreth and Moorbath, 1988), and hybridization of evolved and mafic magmas (Nixon and Pearce, 1987;Eichelberger et al., 2000;Tepley III et al., 2000;Witter et al., 2005;Reubi and Blundy, 2009;Straub et al., 2011;Kent, 2014), with the relative importance of these processes being still widely discussed and possibly changing from one system to another. Plagioclase and orthopyroxene crystals from three Plinian eruptions in the younger (<26 ka) history of Nevado de Toluca have previously been used to show that interaction of a resident shallow crustal silicic magma with two distinct more mafic melts is the dominant mechanism that produced the Plinian eruptive products (Smith et al., 2009;Weber et al., 2019). Our much extended dataset of mineral and bulk-rock geochemistry for 19 additional eruptions, spanning the entire 1.5 Ma history of the volcano, allows us to evaluate the configuration of the plumbing system and pre-eruptive conditions through time.
The long-term mineral record of Nevado de Toluca reveals that strongly zoned pyroxene crystals, similar to those found in the Plinian deposits , are a characteristic feature in eruptive products throughout the volcano's history ( Figure 6). Most remarkable about these crystals are elevated Cr 2 O 3 contents of up to 1.2 wt% in clinopyroxenes and up to 0.91 wt% in orthopyroxenes at En fractions typically between 0.4 and 0.5 for cpx and 0.7 and 0.9 for opx. BSE images and chemical profiles reveal that these compositions are only present in crystal cores that have undergone extensive resorption and overgrowth with a significantly more evolved composition and Cr 2 O 3 contents often below detection limit (Figure 7). The data also show that both cpx and opx evolve along two different major trajectories, which can be distinguished by the different En, and MnO at comparable Cr 2 O 3 content ( Figure 6A). Based on the En-Cr systematics in opx, 3 groups of data can be distinguished: 1) High Cr (>0.1 to 0.91 wt%) compositions at high En (80-90) contents (cluster 1). 2) Fairly high Cr (up to 0.54 wt%) at  intermediate En70-80 and 3) Low Cr contents (frequently below detection limit) at low En contents between ∼50 and 70 mol%. Differences in En contents between the two major Cr 2 O 3 trajectories (clusters 1 and 2) are likely representing 2 distinct and spatially separated melts, as the Fe-Mg exchange between mafic melt and pyroxene is mostly dependent on melt composition (Gaetani and Grove, 1998;Hirschmann et al., 2008;Putirka, 2008;Waters and Lange, 2017). As shown in Weber et al. (2019), the highest Cr pyroxene cores from Nevado de Toluca are overlapping in composition with pyroxenes from mantle peridotite, indicating that these have grown from primitive melts in their early stages of evolution and therefore provide information about the deep structure of the plumbing system. Within a radius of 80 km to Nevado de Toluca, a large number of primitive mantle melts have been identified in the neighboring Valle del Bravo to the west (Blatter and Carmichael, 1998;Wallace and Carmichael, 1999;Blatter et al., 2007) and in the Sierra Chichinautzin Volcanic Field to the east Gomez-Tuena et al., 2007b). As recently compiled by Schmidt and Jagoutz (2017), the great majority of these melts are high-Mg basaltic andesites with SiO 2 contents mostly between 52-56 wt% and Cr contents of 200-500 μg/g, while calc-alkaline arc basalts are rare in the area. Pyroxene compositions in these melts are fully consistent with the interpretation that the high Cr pyroxene cores from Nevado de Toluca are derived from primitive basaltic andesite melts that may either originate by direct mantle melting (Straub et al., 2011;Straub et al., 2014) or derive from crystallization of primary calc-alkaline basaltic mantle melts at upper mantle pressures (Schmidt and Jagoutz, 2017). Interestingly, pyroxenes in the analyzed Frontiers in Earth Science | www.frontiersin.org September 2020 | Volume 8 | Article 563303 monogenetic basaltic andesite to andesite samples ("Tenango", "Sabanillas") are very similar in composition to the mafic pyroxene cores found in the dacite samples from the main volcano ( Figure 6). It is therefore plausible that the shallow crustal plumbing system of Nevado de Toluca is recharged by deep basaltic andesite to andesite magmas with at least two prominent compositions. Two-pyroxene thermometry (Putirka, 2008; Eqn. 36) on average outermost pyroxene compositions indicate similar temperatures of 1,030 ± 58°C ("Sabanillas") and 1,050 ± 58°C ("Tenango") for the 2 monogenetic eruptions, indicating that the temperatures of magmas feeding the shallow storage region are >1,000°C. The presence of primitive pyroxene cores, likely originating from mantle or lower crust, in shallowly sourced silicic eruptions over the long-term history of Nevado de Toluca is consistent with a conceptual model in which magmatic systems are transiently interconnected through the entire Earth crust (Ruprecht and Plank, 2013;Cashman et al., 2017). The occurrence of mafic pyroxene cores in almost all of the analyzed samples testifies that the plumbing system of Nevado de Toluca has been operating in a relatively stable configuration and by similar pre-eruptive processes throughout its entire history. While no systematic differences are evident in mafic pyroxene cores between explosive and effusive eruptions, the record reveals that the early stages of activity ("Old Nevado" and "Domes") differ from the more recent history ("Young PD" and "Recent Nevado") by the more frequent occurrence of high Cr cpx and less well defined mafic clusters in opx ( Figures 6B,C). These variations in phase abundances and clustering of the data may reflect differences in the composition of recharge magmas with time due to changes in the melt source. However, experimental phase equilibria of basaltic andesite magmas and thermal modeling offer a more nuanced framework to explain these data. High pressure (i.e., 1 GPa) fractional and equilibrium crystallization experiments on hydrous basaltic andesite starting compositions by Ulmer et al. (2018) show that opx crystallizes as the liquidus phase in such melts and is joined by cpx at slightly lower T, while amphibole is stabilized at even lower T at the expense of opx in a peritectic reaction. Increasing temperatures of the mafic recharge magma with time could therefore explain the occurrence of the early stage mafic amphibole and lack of opx ("U5 Raíces," "U14 Tepehuisco"), which are indicative of lower recharge magma T, as well as the decreasing abundances of mafic cpx through time. This scenario is consistent with thermal modeling of deep crustal magma reservoir assembly, which typically show trends of increasing temperature with time due to the thermal maturation of the crust (Annen and Sparks, 2002;Annen et al., 2006;Weber et al., 2020).

The Role of Shallow Crustal Hybridization
Pyroxene and plagioclase crystals provide an extended record of magma hybridization processes in the upper crustal plumbing system of Nevado de Toluca over the entire history of the volcano. Interaction of mafic and silicic magma upon eruption is recorded texturally and compositionally between low En-Cr pyroxene rims and the mafic cores (Figures 6, 7). The An-FeO relation of plagioclase crystals provides a complementary record of magma hybridization. As shown in Figure 8, plagioclase rims are systematically higher in FeO and An contents compared to crystal cores in most of the studied samples, indicating that the plumbing system was subjected to the episodic input of more mafic magma prior to almost each eruption over 1.5 Ma. Plagioclase cores in the analyzed samples typically record a range of An contents at low FeO (Figure 8), which can be interpreted to result from variable thermal conditions during long-term magma storage rather than pre-eruptive compositional mixing (Couch et al., 2001;Ruprecht and Wörner, 2007).
Petrologic experiments on the Plinian deposits of Nevado de Toluca revealed that these eruptions were fed from an upper crustal reservoir at about 2 kbar , Arce et al., 2013a, Arce et al., 2013b, which is in agreement with most magma storage estimates of arc volcanoes worldwide . Crystallization pressures calculated between amphibole rims and average groundmass are presented in Figure 9 and are generally in agreement with experimental estimates but suffer from large uncertainties (4 kbar for individual analyses; Putirka, 2016) and only approximately known liquid compositions, which may be particularly problematic due to evidence for magma hybridization. Nevertheless, the agreement in mineral records, erupted bulk-rock compositions and barometry indicates that dacites and andesites are erupted from an upper crustal reservoir over the eruptive history of Nevado de Toluca. We used a range of mineral  (Holland and Blundy, 1994) of average outermost rim compositions are shown as red rectangles with 2 standard deviation error bars. Amphibole-liquid temperatures (Putirka, 2016) were calculated by using outmost amphibole compositions from euhedral crystals and average groundmass analyses of the corresponding thin section (green triangles). Amphibole temperatures were also calculated using a calibration independent of melt chemistry (Putirka, 2016;blue diamonds), which for most samples shows agreement within errors to amphibole-liquid temperatures. Fe-Ti oxide temperatures and oxygen fugacity (relative to Nickel-Nickel-Oxide) are shown as orange circles and purple squares respectively and were calculated using Ghiorso and Evans (2008). Boxplots of pressure (kbar) were calculated using outermost amphibole rims and average groundmass compositions for different samples following Putirka (2016). Calculation details and associated uncertainties are discussed in the main text. Crystal fractions shown on the right side were estimated by point counting.
Frontiers in Earth Science | www.frontiersin.org September 2020 | Volume 8 | Article 563303 thermometry techniques to constrain pre-eruptive temperatures ( Figure 9). The T estimates derived from the different methods are variable but for most samples agree within their uncertainty, which likely in part reflects the mentioned issue of melt composition approximation by average groundmass. Still, the thermometers show rather consistent results in terms of different average temperatures for individual eruptions, which typically range between 800 and 950°C. As Nevado de Toluca is a system markedly impacted by hybridization processes, these temperatures cannot be interpreted to represent long-term magma storage conditions, but reflect the equilibration of crystals following a recharge event. Thus, temperature differences between eruptions can be associated with variations in the recharge temperature and the relative volumes of mafic and silicic magma interacting prior to these events. While this relation can be exemplified by the lower temperatures and less evidence for mixing in the "Pink Pumice" eruption compared to the following "BAF37" and "BAF28" eruptions ( Figure 9), no systematic relation between eruption type and intensive parameters (P-T-fO 2 -Crystallinity) of the magmatic system are evident from the data. Altogether, textural observations and geochemical data indicate that andesites and dacites at Nevado de Toluca are consistently produced by hybridization and mingling of mafic recharge magmas and a silicic reservoir resident in the upper crust over the entire history of the volcano. In the following, using thermal modeling results, we will discuss how the evolution of the thermal and chemical architecture of the magmatic system can control the subtle but measurable evolution of the chemical composition of the erupted magmas at Nevado de Toluca. This is of interest, because establishing the long-term pattern of chemical evolution of this volcano may serve to evaluate its potential future activity.

Bulk-Rock, Magma Recharge, and Hybridization History
Our analyses confirm that the Nevado de Toluca has erupted dacites and andesites throughout its life-cycle. Within the andesitic Frontiers in Earth Science | www.frontiersin.org September 2020 | Volume 8 | Article 563303 to dacitic compositional range, only few major and trace elements show minor but consistent variations in time ( Figure 10). The average SiO 2 content of erupted bulk-rock compositions increases by about 2 wt% over the 1.5 Ma history of the volcano, indicating slightly more evolved magma geochemistry in time ( Figure 10A). This is consistent with the drop of the average V (compatible) contents from 83.2 μg/g for eruptions older than 1,200 ka to 70.0 μg/g in the later eruptive history ( Figure 10B) and the increase of average Nb (incompatible) content from 3.8 to 4.6 μg/g over time ( Figure 10C). The long-term temporal evolution of mineral chemistry at Nevado de Toluca provides further insights into its life-cycle. The An contents in plagioclase, which are mostly sensitive to temperature, melt composition, and P H2O (Cashman and Frontiers in Earth Science | www.frontiersin.org September 2020 | Volume 8 | Article 563303 16 Blundy, 2013) record similar conditions for most of the studied eruptions ( Figure 11A). Variation mostly between An fractions of 0.25 and 0.55 are thus a reflection of the range of thermo-chemical conditions that the system experiences through time. Distributions of An contents in crystal cores are bimodal in many of the studied eruptions and rims show higher An content in average compared to crystal cores. Such distributions are in agreement with complex resorption and overgrowth textures observed in plagioclase and are best explained by repeated heating and cooling cycles (Smith et al., 2009;Weber et al., 2019), which are ubiquitously observed in intermediate and silicic arc magmas (Tepley III et al., 2000;Davidson et al., 2005;Humphreys et al., 2006Humphreys et al., , 2010Zellmer and Turner, 2007;Edmonds et al., 2010;Ruprecht et al., 2012). While An fractions record similar ranges through time, FeO contents are clearly higher in plagioclase crystal rims and interiors from older eruptions >1,100 ka ( Figure 11B). As discussed before, elevated FeO contents in the crystals likely reflect more mafic conditions. As Fe is mostly incorporated into plagioclase as ferric iron (Wilke and Behrens, 1999), changing fO 2 conditions may also be invoked to explain the temporal trend and higher contents in crystal rims. However, the fO 2 estimates do not show any correlation with differences in FeO contents of plagioclase (Figures 9, 11). We therefore conclude that the high FeO contents in early erupted plagioclase have grown from more mafic melts or a higher contribution of the mafic melt to the hybrid rocks. As An contents record similar ranges over time, while FeO decreases (Figure 11), either lower P H2O or temperature would be required in the early phase to compensate the effect of more mafic melt chemistry and thus keep An contents constant through time. As discussed before, a lower recharge magma temperature is consistent with the occurrence of mafic amphibole and cpx  in the early stages and thermal modeling of incremental hot-zone assembly (Annen et al., 2006;Weber et al., 2020). Interestingly, the plagioclase chemistry shows episodes with a stronger and weaker mafic magma component, indicated by higher or lower FeO contents in crystal rims. These chemical variations are broadly consistent with the thermometry data (Figure 9), corroborating the interpretation that these variations result from differences in recharge volumes or temperatures. This interpretation is generally consistent with Martínez-Serrano et al. (2004), who proposed that temporal variations in bulk-rock and Sr-isotopes are due to input of fresh mafic magma into the system. In any case, a detailed analysis of higher resolution stratigraphic records would be required to fully resolve the nature of these episodic changes. Amphibole compositions in the sequence of studied eruptions show large differences in variance and systematic changes with time ( Figures 12A,B). This is exemplified in Ti (C-site) p.f.u. and Al IV p.f.u., both proxies for temperature variation (Putirka, 2016), which is higher in early erupted amphiboles (>1,000 ka) and decreases progressively toward lower values. Similar to the plagioclase data, amphibole records chemical variations, such as increasing and decreasing Ti contents ( Figure 12A), which are likely reflecting the extent of replenishment and differentiation in the subvolcanic reservoir. Using the empirical major element inversion method from Zhang et al. (2017), we calculate melt compositions in equilibrium with the amphiboles ( Figures  12C,D). The recalculated melt compositions are in good agreement with the bulk-rock and glass geochemical trends of Nevado de Toluca (Supplementary Figure 3). Amphibole equilibrium melt compositions span a large range of SiO 2 contents between 55 and 80 wt%, revealing a large diversity of chemical compositions present in the plumbing systems, which are not observed in melt inclusion or glass analyses. The amphibole equilibrium melts also show more mafic conditions in the early history of the volcano, which is particularly evident in recalculated SiO 2 and FeO melt contents of the "Raices" and "Tepehuisco" eruptions, indicating that amphiboles in these dacites are derived from basaltic andesite to andesitic liquids. This is consistent with the presence of reaction coronas surrounding such grains ( Figure 5E). High pressure phase equilibria experiments on several primitive basaltic andesite compositions typical of arc magmas at 1 GPa by Ulmer et al. (2018) shows that amphibole crystallizes from such liquids <1,050°C, while at higher temperatures opx and cpx assemblages dominate. The absence of mafic amphiboles in the later history of Nevado de Toluca may thus indicate either 1) complete digestion of mafic amphibole in the silicic magma of later eruptive stages, 2) a more evolved recharge magma composition in later stages of activity, in which mafic amphibole is not stable (e.g., Andesite >1,000°C; Blatter and Carmichael 2001) or 3) increasing temperature of the recharge magma with time beyond the stability field of mafic amphibole. While the second hypothesis is inconsistent with the presence of ultramafic pyroxene cores, the latter is consistent with thermal modeling of pulsed magma injection, which shows increasing temperatures with time during the build-up of reservoirs in mid to lower crustal hot zones due to progressive heat-up of the surrounding wall-rocks (Annen and Sparks, 2002;Annen et al., 2006;Weber et al., 2020). The temporal variation of En contents in orthopyroxenes shows that most eruptions record a large range with slight tendency toward lower variance in the early history of the volcano ( Figure 13). As these variations are associated with textural evidence for the interaction of contrasting magma compositions (Figures 5, 7) and chemical evidence for hybridization of a mafic high Cr and silicic low Cr melt, we suggest that this variation is resulting from differences in melt chemistry of the resident silicic magma and mafic recharge magma. Orthopyroxene crystals from the early history of Nevado de Toluca show a tendency toward lower Cr contents and variability when compared to later eruptive stages, but clear evidence for mafic magma input in the early stages is present in the form of abundant high Cr cpx and mafic amphibole cores. As discussed before this observation is consistent with an increasing temperature of the recharge magma with time, but may also reflect a change in recharge magma composition. In the following we discuss and identify the most likely controls on the long-term trends in bulkrock and mineral chemistry based on thermal considerations.

Controls on Long-Term Geochemical Trends
The long-term geochemical trend toward more silicic compositions observed in bulk-rock and mineral chemistry may be explained by a range of changing conditions either related to the rates of magma recharge, the relation of eruptive to intrusive fluxes, and the modalities and rates of magma reservoir growth. We use thermal modeling of incremental upper crustal magma assembly and withdrawal in order to illustrate and test these processes ( Figure 14). As shown in Figures 14A,B, for an equivalent recharge magma composition and temperature, higher magma injection rates will result in higher average temperature and therefore a more mafic eruptible magma composition of the reservoir. From this follows that the slight trend toward more silicic magma chemistry observed at Nevado de Toluca may be explained by a drop in average recharge rates through time. Such a process may be related to progressive decrease of fertility of the mantle source, if the timescales of melt production and extraction are faster than supply of fertile mantle by convection. Mechanical considerations provide another mechanism that could potentially result in a drop of recharge rates with time. Continuous magma injection into the crust leads to progressive heat-up of wall-rocks surrounding magma reservoirs, which favors the viscous relaxation of overpressures generated by magma injection in lower and upper crustal reservoirs (Jellinek and DePaolo, 2003;Gregg et al., 2012;Karlstrom et al., 2017). This implies a greater potential for magma storage rather than transfer with increasing lifetime of a magmatic system and therefore result in more homogeneous and more silicic erupted compositions with time. The temporal evolution of erupted volumes from Nevado de Toluca shows higher volumetric eruptive fluxes in the early history of the volcano compared to later evolutionary stages (Torres-Orozco et al., 2017), which is consistent with a decrease in magma fluxes with time.
Thermal modeling of pulsed magma injection and withdrawal also shows that the E:I ratio impacts the thermal and compositional evolution of crustal magma reservoirs ( Figures  14A,B). An increasing eruption efficiency (higher E:I ratio) removes more heat from the system and therefore shifts the composition toward more silica rich. A systematic increase in eruption efficiency may thus be invoked to explain the long-term FIGURE 15 | SiO 2 (wt%) density distributions of whole-rock data from Nevado de Toluca (gray), Popocatépetl volcano (orange), Citlaltépetl ( Pico de Orizaba) (gold) and the Colima Volcanic Complex (cyan) with their peripheral cones. Data for Colima (N 322 analysis), Popocatépetl (N 220) and Citlaltépetl (N 99) were minded from the GEOROC database (http://georoc. mpch-mainz.gwdg.de/georoc). Care was taken to exclude oversampling of prominent eruptions. Note that the position of the volcanos' mode and the composition of more mafic cones vary in concert.
Frontiers in Earth Science | www.frontiersin.org September 2020 | Volume 8 | Article 563303 geochemical trend of Nevado de Toluca. In fact, many volcanoes show drastic and non-linear changes in volumetric eruption rates throughout their histories (e.g., Hildreth, 2007). If such differences reflect changing E:I ratios or are related to other internal or external magmatic system variables is, however, not well understood. Yet, an increasing E:I ratio is an unlikely explanation for the geochemical trend, as it is inconsistent with 1) the evolving mechanical properties of magmatic systems through time and 2) the temporal record of erupted volumes from Nevado de Toluca, which indicate higher eruptive flux in the early history of the volcano (Torres-Orozco et al., 2017). The conditions of magma reservoir growth offer a further potential explanation for the long-term geochemical trend observed at Nevado de Toluca. Thermal modeling shows that early injected magma batches cool so rapidly below the solidus temperature that no coherent reservoir exists, until after some incubation time a magma reservoir starts to grow (Annen et al., 2006). From this follows that early injected recharge magmas have less potential to interact with resident extractable silicic mush compared to later intrusions. This process can be illustrated by plotting the volume fractions of resident and recharge magma with time ( Figure 14C). Individual recharge pulses likely remain within a certain volumetric range, while the upper crustal silicic reservoir grows in size. Later injections can therefore be expected to be buffered by the composition of the silicic mush, which will shift the erupted compositions toward more evolved chemistry in time, in particular for systems dominated by magma hybridization such as Nevado de Toluca. Altogether, we propose that the faint long-term compositional shift toward higher silica contents at Nevado de Toluca is controlled by slightly dropping recharge rates, accompanied by hotter recharge magma with time, which is consistent with the effects of thermal and mechanical maturation of the plumbing system (Karlstrom et al., 2017;Weber et al., 2020), and continued reservoir growth, increasing the volume fractions of silicic magma relative to the incoming more mafic recharge melts ( Figure 14C).

Petrologic Modes of Mexican Stratovolcanoes
Our 1.5 Ma dataset from Nevado de Toluca shows that tuning of erupted magma geochemistry to monotonous compositions can be a primordial and persistent property of magmatic systems that do not undergo substantial temporal evolution. Besides Nevado de Toluca, several volcanic systems in the TMVB erupt magmas of geochemically restricted variety (Figures 1, 15). While compositionally resolved long-term eruptive histories are sparse, the available data indicates that preferential eruption of specific intermediate compositions is a long-term property of various Mexican stratovolcanoes. The Colima Volcanic Complex has consistently produced basaltic andesite and mafic andesite magma over the last 30 ka (Luhr and Carmichael, 1980;Crummy et al., 2014;Crummy et al., 2019) and while dacites have been described in a transitional phase of the early history of Nevado de Colima, the bulk of the over several hundred thousand years built Colima Volcanic Complex is composed of mafic andesitic rocks with mode around 60 wt% silica (Robin et al., 1987). Popocatépetl volcano and its precursory structure El Ventorrillo were built over about 310 ka and dominantly erupted andesitic magmas (Schaaf et al., 2005;Sosa-Ceballos et al., 2015;Mangler et al., 2019), which are on average more evolved compared to the Colima Volcanic Complex (Figure 15). Citlaltépetl volcano (i.e., Pico de Orizaba) in the eastern TMVB shows a less well defined compositional mode in the andesite to dacite spectrum which defines the bulk of the volcano, but has generally produced a range of geochemical diversity from basalt to rhyolite compositions in its about 500 ka history (Schaaf and Carrasco-Núñez, 2010). Given lifetimes of several hundred thousand years or longer, the compositional distributions reflect the long-term behavior of these systems.
Thermal modeling results show that differences in recharge rate, temperature and E:I ratio can shift the compositional spectrum of potentially eruptible magma compositions in crustal reservoirs ( Figure 14B). The preferentially erupted magma of volcanoes may thus indicate either differences in magma flux, eruption efficiency or temperature and composition of the recharge magma. The Colima Volcanic Complex, the most mafic of the systems considered here, could thus be built by higher recharge rates compared to Popocatépetl volcano, which accordingly may evolve by higher rates compared to Nevado de Toluca (Figures 14, 15). Large differences in E:I ratio would be required to explain the shift in erupted compositions between the TMVB volcanoes which is, however, an unlikely scenario, as it implies that the eruptibility of viscous dacite magmas is greater than that of more mafic compositions (Takeuchi, 2011). Interestingly, the compositional mode of the summit volcanoes and peripheral monogenetic cones varies in concert for the systems shown in Figure 15, suggesting that differences in the temperature or composition of the recharge magmas are likely controlling the petrochemistry of the most frequently erupted products ("petrologic mode"). This would require that the composition of monogenetic cones in each case approximates the chemistry of recharge magmas feeding the main volcanoes.
All systems considered in Figure 15 show evidence for hybridization of mafic and more evolved magma prior to eruption (Cantagrel et al., 1984;Robin and Potrel, 1993;Schaaf et al., 2005;Witter et al., 2005;Smith et al., 2009;Sosa-Ceballos et al., 2014;Macias et al., 2017;Crummy et al., 2019;Mangler et al., 2019;Weber et al., 2019), but the composition of the mafic mixing endmember is challenging to constrain. In the case of Nevado de Toluca, equivalence between monogenetic centers and feeding recharge magmas is indicated by the similarity of characteristic mafic pyroxene cores ( Figure 6). Monogenetic cones surrounding the Colima Volcanic Complex, are composed of shoshonites (trachybasalt) and calc-alkaline basalts (Carmichael et al., 2006). Based on petrogenetic considerations, the later have been proposed to represent the parental magmas of the Colima andesite suite (Luhr and Carmichael, 1980) and are very similar in composition to the proposed gabbroic mixing endmember of Reubi et al. (2019). A contribution of alkaline mafic magmas to the andesite erupted at Colima has also been documented prior to some eruptions (Crummy et al., 2014;Crummy et al., 2019). Notably, basaltic magmas, equivalent to the calc-alkaline cones, have also been explosively erupted from the main volcano (Crummy et al., 2014;Crummy et al., 2019). Constraints on the mafic recharge magma of Popocatépetl volcano were derived from the presence of olivine crystals (Fo87-90) and their melt inclusions (50-52 wt% SiO 2 ) in recent effusive and explosive products (Roberge et al., 2009), which are identical to olivine compositions from mafic cones with bulk-rock SiO 2 contents of 50-54 wt% in the vicinity of the volcano, belonging to the Chichinautzin Volcanic Field (Schaaf et al., 2005). Further evidence for the equivalence of basaltic to basaltic andesite compositions of the Chichinautzin Volcanic Field and the recharge magmas feeding the magmatic system of Popocatépetl volcano has been presented based on mineral chemistry (Mangler et al., 2020). Overall, the constraints at different major Mexican volcanoes are consistent with the hypothesis that peripheral cone compositions are equivalent to the mafic hybridization endmember that contributes to the production of andesitic to dacitic volcanic rocks in the TMVB. The systematic shift in the recharge and summit cone geochemistry may therefore reflect that compositionally monotonous stratovolcanoes are built by similar recharge rates but with different compositions and temperatures of the feeding magmas, indicating that the petrologic mode for different volcanoes is already set in the deep crust or Earth mantle.

CONCLUSIONS
Our analysis of the long-term volcanic record of Nevado de Toluca volcano in Central Mexico reveals that: (1) Dacite and andesite petrogenesis is dominated throughout the entire life cycle of the volcano by hybridization processes of deep basaltic-andesite recharge magma and an upper crustal silicic mush. This process is captured in compositional zoning of all major mineral phases. (2) The compositional tuning to preferentially erupted dacite compositions ("petrologic mode") is a primordial and resilient feature throughout the volcanos' history.
(3) Minor long-term temporal trends toward more silicic compositions are observed from bulk-rock compositions and mineral chemistry and show that these data are best explained by increasing interaction of mafic recharge magma with growing volumes of silicic mush through time and slightly dropping magma supply rates with increasing temperature of the recharge magma. (4) Very mafic clino-and orthopyroxene compositions indicate that Nevado de Toluca is fed from 2 spatially separated magma sources either in the upper mantle and/or lower crust, at least in the younger history (<59 ka) of the volcano. (5) A comparison of petrologically monotonous stratovolcanoes in the TMVB shows that the mode of the preferentially erupted composition in the main volcano and in peripheral cones varies in concert. This points toward an important role of recharge magma composition and/or temperature in controlling the compositional mode of Mexican stratovolcanoes.

DATA AVAILABILITY STATEMENT
All datasets presented in this study are included in the article/ Supplementary Material.

AUTHOR CONTRIBUTIONS
GW and LC conceived this study. GW carried out the geochemical analysis and numerical modeling. Field work was carried out by GW, LC, and JA. The initial manuscript draft was prepared by GW. All authors contributed to data interpretation and writing of the final manuscript.