A Pyroxenic View on Magma Hybridization and Crystallization at Popocatépetl Volcano, Mexico

Introduction

Continental arc volcanism is dominated by andesitic and dacitic magmas producing a wide range of eruptive styles (e.g., Reubi and Blundy, 2009; Kent et al., 2010; Ruprecht and Bachmann, 2010). Many studies of arc volcanoes focus on products of large, explosive eruptions, often neglecting effusive and smaller explosive eruptions that dominate the activity of many stratovolcanoes worldwide (Davidson and De Silva, 2000). However, to fully characterize the transcrustal architecture of volcanic plumbing systems (Edmonds et al., 2019) and the dynamics of pre-eruptive processes such as magma mixing, it is crucial to consider entire eruptive histories. Effusive eruptions may offer complementary insight into the controls of magma recharge on eruptive style (Ruprecht and Bachmann, 2010; Koleszar et al., 2012; Bouvet De Maisonneuve et al., 2013), and their intra-eruption chemical variability can reveal heterogeneity and complexity of volcanic plumbing systems at various time- and length-scales (Mangler et al., 2019; Neave, 2020). Moreover, flank vents and cinder cones surrounding stratovolcanoes may erupt magmas that bypassed shallow magma reservoirs, thus sampling more mafic inputs that offer a window into the roots of the plumbing system (e.g., Schaaf et al., 2005; Carmichael et al., 2006; Corsaro et al., 2009; Petrone, 2010; Mangler et al., 2019; Weber et al., 2019).

The Popocatépetl Volcanic Complex (PVC) is an andesitic to dacitic arc volcano with an eruptive history including at least five Plinian and ten effusive eruptions in the last ∼23.5 ka (Espinasa-Pereña and Martín-Del Pozzo, 2006; Mangler et al., 2019). Ongoing eruptive activity commenced in 1995 and is characterized by dome emplacement and destruction cycles that formed and destroyed 85 lava domes by the end of 2019 (Gómez-Vazquez et al., 2016).

Previous work considered Plinian events and present-day activity in isolation from the broader eruptive history, resulting in several competing models of the architecture of the plumbing system and magma petrogenesis. Most studies postulate the presence of one upper crustal evolved magma reservoir periodically recharged by mafic melts, but there is ongoing debate about whether a sustained mafic magma reservoir is present in the shallow plumbing system (Witter et al., 2005), or whether mafic magmas only manifest as short-lived recharge events triggering eruptions (Straub and Martín-Del Pozzo, 2001; Schaaf et al., 2005; Atlas et al., 2006; Cross et al., 2012; Sosa-Ceballos et al., 2012, 2014). Furthermore, both the primary composition of magmas feeding the PVC plumbing system and their fractionation patterns remain poorly constrained. Straub and Martín-Del Pozzo (2001) argue that mantle-derived (basaltic) andesites with ∼55–60 wt.% SiO₂ feed the PVC, whereas Schaaf et al. (2005) suggest feeding by basaltic andesites or basalts similar to the magmas erupted in the neighboring Sierra Chichinautzin Volcanic Field (SCVF). Both studies suggest early crystallization of olivine + Cr-rich spinel at lower crustal levels, however, Straub and Martín-Del Pozzo (2001) propose that olivine ceases to crystallize at shallow depths in favor of clinopyroxene and orthopyroxene, while Schaaf et al. (2005) and Sosa-Ceballos et al. (2012) suggest co-crystallization of orthopyroxene + clinopyroxene + olivine. In addition, the significance of mixing and hybridization between evolved and mafic magmas is recognized in previous studies (e.g., Straub and Martín-Del Pozzo, 2001; Schaaf et al., 2005; Atlas et al., 2006; Sosa-Ceballos et al., 2012; Mangler et al., 2019), however, the mechanisms of hybridization and its effects on crystal cargoes have not yet been studied.

Here, we present a comprehensive petrological examination of effusive and explosive eruptions of the last ∼14.1 ka to test, revise and synthesize existing models of the PVC’s plumbing system. We focus on ortho- and clinopyroxene, because it crystallizes over a wide range of P-T-X conditions, and its compositions and textures record transcrustal magmatic evolution and processes. The petrogenetic archive stored in pyroxene crystals is widely used for thermobarometry (Putirka, 2008; Neave and Putirka, 2017) and, more recently, to reconstruct geochemical evolution of magmas (e.g., Mollo et al., 2018) as well as dynamics and timescales of pre-eruptive processes, such as magma mixing, mush remobilization and cannibalism, and ascent rates (e.g., Petrone et al., 2018; Ubide and Kamber, 2018; Ubide et al., 2019; Di Stefano et al., 2020). In the present study, we use a detailed analysis of ortho- and clinopyroxene populations to reconstruct the architecture and thermal state of the plumbing system, dissect the dynamics of magma mixing and its short- and long-term effects on the crystal cargo, constrain the crystallization sequence, and examine the composition of primary melts feeding the PVC.

Geological Background

The Popocatépetl Volcanic Complex (PVC) is situated ∼70 km SE of Mexico City at the front of the eastern Trans-Mexican Volcanic Belt (TMVB; Figure 1), an E – W trending arc associated with the oblique subduction of the Cocos and Rivera plates beneath the Northern American plate (Figure 1; Ferrari et al., 2012 and references therein). Geophysical evidence suggests that the Cocos plate is located at ca. 80 km depth beneath the PVC, plunging at a 75° angle beneath the Central Mexican arc front (Pardo and Suarez, 1995; Pérez-Campos et al., 2008). The continental crust around the PVC has an inferred thickness of 40–50 km (Ferrari et al., 2012), consisting of Precambrian and Paleozoic terranes (Ortega-Gutierrez et al., 1995; Ferrari et al., 2012) overlain by Mesozoic marine sediments (Fries, 1965) and Tertiary terrigenous and volcanic rocks (Espinasa-Pereña and Martín-Del Pozzo, 2006).

Modern Popocatépetl is constructed on the remnants of at least two older stratovolcanoes dating back to ∼330 ± 10 ka (Cadoux et al., 2011; Sosa-Ceballos et al., 2015). The Nexpayantla volcano (Mooser et al., 1958; gray in Figure 1) collapsed toward the SW between 183 ± 7 ka and 96 ± 8 ka BP (Sosa-Ceballos et al., 2015; Mangler et al., 2019), and the subsequent Ventorrillo volcano (blue colors in Figure 1) collapsed southwards at ∼23.5 ka BP during the Tochimilco Plinian eruption (described as the White Pumice eruption in Siebe et al., 2017). Construction of Popocatépetl volcano began following the Tochimilco eruption with the emplacement of the Chipiquixtle and Ecatzingo lavas on the SW flank, and the San Pedro and La Colonia flows on the SE flank (green colors in Figure 1; Espinasa-Pereña and Martín-Del Pozzo, 2006; Siebe et al., 2017; Mangler et al., 2019). At ∼14.1 ka, the VEI 6 Pumice with Andesite (PwA) Plinian eruption (Mooser, 1967; Sosa-Ceballos et al., 2012) resulted in a NE-directed shift of the central vent location (Siebe, 1997; Ortega-Guerrero and Newton, 1998), and subsequent lavas were emplaced in the north and east of the edifice (yellow – red colors in Figure 1). The majority of effusive eruptions since the PwA Plinian eruption, such as the Las Truchas, Atlimiyaya, Capula and Nealticán lavas, originated from fissures on the NE flank (Espinasa-Pereña and Martín-Del Pozzo, 2006; Mangler et al., 2019). Capula and Nealticán are the youngest volumetrically important lavas at the PVC, which are stratigraphically placed between the VEI 6 Yellow Pumice (YP) eruption at ∼2,150 a BP, and the VEI 4 Pink Pumice (PP) eruption at ∼1,100 a BP (Espinasa-Pereña and Martín-Del Pozzo, 2006; Siebe and Macías, 2006). Since the PP Plinian eruption, dome growth and destruction cycles similar to present-day activity (Gómez-Vazquez et al., 2016) dominated the eruptive style of Popocatépetl volcano (Martín-Del Pozzo et al., 2016).

Previous petrological studies of pumices erupted over the past 23.5 ka describe a compositionally and texturally diverse crystal cargo comprising plagioclase (An_21–84) + orthopyroxene (enstatite; Mg# = 100*${\displaystyle \frac{Mg}{Mg+F{e}^{2+}+Mn}}(a.p.f.u.)=$ 54 – 84) + clinopyroxene (augite; Mg# = 52–91) ± amphibole ± olivine (Fo_72–90) ± Fe-Ti-oxides (Straub and Martín-Del Pozzo, 2001; Schaaf et al., 2005; Sosa-Ceballos et al., 2012, 2014). Mineralogical data for lavas of the PVC is relatively scarce (Schaaf et al., 2005; Sosa-Ceballos et al., 2015), and petrological constraints on the dynamics and evolution of the PVC’s plumbing system are therefore hampered by an incomplete dataset that does not account for temporal evolution and different eruptive styles.

Mangler et al. (2019) presented whole-rock major and trace element data as well as Sr–Nd–Hf isotopic ratios of explosive and effusive activity at the PVC in the past ∼23.5 ka. The study found restricted compositional ranges, which the authors ascribed to a high degree of magma hybridization (i.e., homogenization of compositionally distinct melts and crystal cargoes) in the shallow plumbing system beneath the PVC (Mangler et al., 2019). Some flank eruptions probe more mafic magmas with distinct isotopic signatures consistent with lower crustal assimilation (e.g., Chipiquixtle and Atlimiyaya; Mangler et al., 2019), offering insights into the deeper portions of the PVC plumbing system. Building on the work of Mangler et al. (2019), and to test existing models of the PVC plumbing system and petrogenesis, the present study examines petrographic and mineralogical data for a subset of their sample suite.

Materials, Methods and Nomenclature

Five effusive and two explosive units of the last ∼14.1 ka BP were selected for analysis. Sample locations are shown in Figure 1 and coordinates are given in Mangler et al. (2019). The Las Truchas, Atlimiyaya and Olga lavas represent effusive activity between the PwA and YP Plinian eruptions (∼14.1–2.15 ka BP), and the Capula and Nealticán lavas are bracketed by the YP and PP Plinian eruptions (∼2.15–1.1 ka BP; Figure 1). In addition, the pre-PwA Chipiquixtle flank eruption (56.9 wt.% SiO₂) was chosen to investigate less evolved magma influx to the PVC (Mangler et al., 2019).

High resolution backscattered electron (BSE) images were acquired at the Imaging and Analysis Centre at the Natural History Museum, London, using a FEI Quanta field emission gun scanning electron microscope at 15 kV acceleration voltage, 15 mm working distance, and 10–30 μs dwell time. Chemical compositions of minerals, groundmass and melt inclusions were determined at the Natural History Museum on a CAMECA SX-100 electron microprobe using a 20 kV, 20 nA, focussed beam (∼2 μm) for olivine, pyroxene and Fe-Ti-oxides, a defocussed 10 μm diameter beam for plagioclase, and a 20 μm diameter beam for glass and bulk analyses of micro- to cryptocrystalline groundmass. Peak counting times were 20 s for major elements and 30 s for minor elements, except for Na and K, which were measured for 10 s in the beginning of each analysis to minimize element migration. Background counting times were 50% of peak counting times. The accuracy of the method is given by relative deviations of generally <5% from recommended values for USGS reference basaltic glass BCR-2G, and the repeatability of the method is generally better <10% (2σ) (Supplementary Table 1).

Recent advances in our understanding of volcanic plumbing systems as transcrustal mush systems (e.g., Cashman et al., 2017) have highlighted the need to clarify the terms used to describe these systems. In this work, we distinguish between magma reservoirs as physical entities of “partially or wholly molten rock with varying proportions of melt, crystals and exsolved volatiles” (Edmonds et al., 2019), and magmatic environments (ME) as virtual entities defined by a set of P, T, X, fO₂, and fH₂O conditions (Kahl et al., 2011, 2017). This definition allows to resolve compositional heterogeneities as well as transient processes in magma reservoirs.

To avoid ambiguity when describing crystal phases and their complex relationship with their host melts, some authors (Thomson and Maclennan, 2012; Neave et al., 2013) have recently re-introduced the term macrocryst (originally used for kimberlites by Skinner and Clement, 1979) instead of the genetically connotated terms phenocryst, antecryst, and xenocryst. We agree that these terms are misleading and recognize the need for a unifying, unambiguous definition. However, the term macrocryst is problematic in itself as it implies large crystals, a subjective term which lacks a unified understanding of what size this corresponds to. Lower size thresholds used to define macrocrysts in recent studies range from 0.2 to 5 mm (Le Maitre et al., 1989; Thomson and Maclennan, 2012; Neave et al., 2013; Óskarsson et al., 2017; Bennett et al., 2019). This range reflects variations in crystal size distributions for different magmatic systems and demonstrates the ambiguity introduced by terminology based on relative crystal size. To avoid any ambiguity, we use an alternative nomenclature based on the simple distinction between crystals that crystallized in the magmatic system, and those that did not: xenocrysts, as defined by Sollas (1892), are crystals of “foreign derivation” and can usually be easily identified as such under the microscope or using microanalytical techniques; philocrysts are crystals that formed in the volcanic plumbing system, regardless of their individual histories. The broad petrogenetic classification into philocrysts and xenocrysts facilitates unambiguous, accurate and consistent description of crystal cargoes in all volcanic settings.

Results

Petrography

Mineral assemblages, textures and compositional ranges are generally similar for both effusive and explosive eruptions, and differences between units predominantly manifest in variable modal proportions of philocrysts and groundmass (Table 1).

Lavas and pumices are characterized by a philocryst assemblage of plagioclase + orthopyroxene + clinopyroxene with modal proportions of [pl]:[opx]:[cpx] = [55–70%]:[7–30%]:[4–17%]. Minor philocryst phases are olivine (<1–15% modal), which is commonly replaced by orthopyroxene – oxide intergrowths, as well as titanomagnetite and ilmenite (<<5% modal). Amphibole is absent in all studied samples except for single relic crystals in the Atlimiyaya and Las Truchas lavas.

Philocrysts account for 10–45 vol.% of lava samples and exhibit a seriate texture in a micro- to cryptocrystalline groundmass of varying vesicularity (<1–25 vol.%: Figure 2A). The groundmass is basaltic andesitic to rhyolitic in composition (Figure 3A) and comprises micro- to nanocrystals of plagioclase, orthopyroxene and Fe-Ti-oxides, as well as varying amounts of interstitial glass. Pumices are characterized by vesicularities of 65–75 vol.%, with ≤10 vol.% philocrysts in a groundmass of light to dark brown glass with dispersed Fe-Ti-oxides and domains of cryptocrystalline material (Figures 2B–D).

A proportion of philocrysts (1–25% modal) are present as glomerocrysts <1 to >10 mm in size. Glomerocrysts consist of euhedral to subhedral pyroxene philocrysts with minor plagioclase in an evolved interstitial groundmass comprised of anhedral plagioclase, apatite, Fe-Ti-oxides and glass (Figure 2C).

Phase Assemblage

Plagioclase

Plagioclase crystals range in size from <1 μm to ≤5 mm and are predominantly euhedral with elongate to tabular habits (Figures 2A,B,D). Plagioclase is dominantly andesine and labradorite with An₃₀ – An₆₅, extending to ≤An₈₀ only in the Las Truchas lava (Figure 3B; representative mineral compositions given in Supplementary Table 2). The overall range is consistent with literature data for plagioclase compositions (Figure 3B; Straub and Martín-Del Pozzo, 2001; Witter et al., 2005; Sosa-Ceballos et al., 2014). Unresorbed crystal cores generally show compositions <An₅₀ (Figure 3C) and are dominated by high frequency (<10 μm), low amplitude (ΔAn of <<5) oscillatory zoning (Figure 2D). Most crystals record several reversely zoned bands with sharp increases in An content (ΔAn = 10–20; Figure 3C) and subsequent normal zoning (Figure 2D). Reverse zoning is associated with abundant disequilibrium textures such as surficial and pervasive resorption producing rounded or embayed dissolution surfaces, sieved cores and bands (Figures 2A,B,D). Outermost rims of plagioclase crystals are commonly reversely zoned with > An₅₀. Melt and Fe-Ti-oxide inclusions, as well as rare apatite and corroded clinopyroxene and orthopyroxene inclusions are present in plagioclase crystals of all sizes (Figure 2D).

Orthopyroxene

Orthopyroxene forms euhedral to subhedral crystals with stubby to elongate habit, and shows a seriate distribution with crystal size rarely exceeding 1 mm (Figures 2A,C,E). Philocrysts are clinoenstatites (1–4 wt.% CaO) with a range of En₅₄–En₉₀ (Figure 3D), corresponding to Mg# of 56–90. Unusually magnesian orthopyroxene up to En₉₀ (Mg# 90) in the Olga, YP, Nealtican and PP units exceed the range reported in literature (Mg# = 54–84; Figure 3D) and are amongst the most enstatite compositions reported for arc magmas in general (see arc orthopyroxene compilation in Weber et al., 2019). The pre-PwA Chipiquixtle lava, distinguished by its least evolved whole-rock composition and highest abundance of olivine (Table 1), notably shows lower En contents and a more restricted range (En₆₁–En₇₆; Mg# 65–78; Figure 3D), a pattern reproduced by the olivine-rich Atlimiyaya and Las Truchas lavas (En₆₂–En₇₇; Mg# 64–80). Normal and reverse zoning is common in orthopyroxene and associated with varying degrees of diffusion or surficial and pervasive resorption (Figure 2E). Inclusion systematics in orthopyroxene are complex and include sub- to anhedral titanomagnetite and ilmenite (<1–100 μm), apatite laths of 1–30 μm length (Figure 2E), pyrrhotite globules of up to 10 μm, and rare cubic Cr-spinels.

Clinopyroxene

Clinopyroxene is present as euhedral to subhedral philocrysts 250 μm −1 mm in size, with a lower aspect ratio than orthopyroxene. The majority of clinopyroxenes are augites with Wo₃₉–Wo₄₅ and Mg# similar to orthopyroxenes (Mg# 69–91), consistent with literature data for the PwA eruption and present-day activity (Figure 3E). Clinopyroxene crystals tend to exhibit more strongly embayed crystal surfaces than orthopyroxene, and rounded or patchy rather than sieved crystal cores; some crystals show oscillatory zoning (Figure 2F). Mineral inclusions in crystal cores mirror those of orthopyroxene, although inclusions are generally rarer in clinopyroxene.

A compositionally distinct group of low-Ca clinopyroxene crystals with systematically less calcic compositions (Wo₃₁–Wo₄₄) at similar Mg# (71–88) is found together with abundant olivine philocrysts in the Chipiquixtle, Atlimiyaya and Las Truchas flank eruptions, and sparsely in the Nealticán lavas, PP and YP pumices. Low-Ca pyroxene crystals exhibit prismatic habits and complex zoning patterns dominated by resorption and oscillations. Mineral inclusions are sparse and restricted to Fe-Ti-Cr oxides.

Olivine

While not a major philocryst phase in most units, olivine constitutes up to ∼15% of the philocryst assemblage in the low-Ca clinopyroxene bearing Chipiquixtle, Las Truchas and Atlimiyaya lavas and is the second most abundant phase after plagioclase in the former two. Olivine forms euhedral to anhedral philocrysts <50 μm to 1 mm in size with a compositional range of Fo₅₆–Fo₈₉ (literature: 70–91; Figure 3F). Pristine olivine shows homogeneous cores with Fo₈₂–Fo₈₉, up to 5400 ppm Ni, and abundant Cr-spinel inclusions of typical cubic shape (Figure 2G). Such pristine olivine is abundant in the Chipiquixtle, Las Truchas and Atlimiyaya lavas, whereas in other units, olivine is rarer and exhibits pervasive disequilibrium textures such as extensive diffusion, coarse sieving, and rounded and embayed crystal faces. Many philocrysts are partly or entirely dissolved and replaced by orthopyroxenes or a vermicular intergrowth of orthopyroxene and Fe-Ti-(Cr-) oxides (Figure 2H). Clusters of such orthopyroxene–oxide intergrowths (± peripheral clinopyroxene and plagioclase) can be found on μm- to mm-scale in all units, even in those where olivine is scarce (i.e., Olga lava and all explosive units).

Pyroxene Populations

Ortho- and clinopyroxene philocrysts show exceptional compositional diversity extending to highly magnesian compositions (Mg# 54–90 for orthopyroxene and 69–91 for clinopyroxene). The modal abundance and diversity of pyroxene philocrysts provide the most comprehensive petrogenetic record at the PVC, which can be used to trace pre-eruptive magmatic evolution, and to reconstruct chemical and physical processes of mixing, mingling and hybridization. We couple compositional data with textural characteristics and the mineral inclusion record of ortho- and clinopyroxene to define distinct crystal populations that reflect changing MEs in the plumbing system.

We first distinguish two Types of ortho- and clinopyroxene, based on their broadly bimodal core compositions: an evolved Type 1 and a mafic Type 2 (Figure 4B). Evolved Type 1 crystals are characterized by core compositions of Mg# <71 for orthopyroxene and <78 for clinopyroxene, and are further subdivided into Type 1a and Type 1b based on distinct mineral inclusions found in crystal cores (Figure 4A). Type 1a (T1a) crystals are characterized by abundant apatite, Fe-Ti-oxide and sulfide inclusions; Type 1b (T1b) cores lack apatite inclusions and rarely show Fe-Ti-oxide or sulfide inclusions. T1b crystal cores further tend to show compositional gradients preserving slightly elevated Mg, Ca, and Cr-concentrations (green symbols in Figure 4B), whereas T1a cores are compositionally more homogeneous (yellow symbols in Figure 4B).

Type 2 pyroxene (T2) is characterized by Mg-rich core compositions of Mg# >78 for orthopyroxene and >82 for clinopyroxene, and Cr₂O₃-concentrations of >0.25 wt.% (blue symbols in Figure 4B). Mineral inclusions present in T2 crystals are Cr-spinel (Figure 4A) and rare sulfide globules. Aspect ratios of T2 crystals are generally higher than those of T1a pyroxene (Figure 4A), and they sometimes show sector zoning or skeletal growth, indicating rapid growth rates at varying degrees of undercooling (Lofgren, 1980; Watson and Liang, 1995; Brophy et al., 1999; Neave et al., 2019; Ubide et al., 2019; Masotta et al., 2020).

Type 1b ortho- and clinopyroxenes are most abundant in PVC lavas and pumices, making up 44–49% of erupted pyroxenes (Figure 4C). Type T1a and T2 crystals account for 33–40% and 16–20% of erupted pyroxene crystals, respectively. By contrast, olivine-rich flank eruptions Chipiquixtle, Atlimiyaya and Las Truchas show a notable absence of T2 crystals and instead carry chemically and texturally distinct low-Ca clinopyroxene philocrysts. This suggests a distinct petrogenetic evolution for these flank lavas, which will be discussed separate from other PVC units.

Crystal cores of T1a, T1b, and T2 ortho- and clinopyroxene are mantled by at least one band or rim (Figure 4A); entirely homogeneous crystals without any signs of compositional zoning are extremely rare. Bands and rims span the entire range of compositions observed for ortho- and clinopyroxene, but commonly show Mg# intermediate between T1a and T2 compositions (Figure 5). Intermediate band and rim compositions produce reverse zoning in T1a and T1b crystals (Figures 2E,F, 4A) and normal zoning in T2 crystals (Figures 2C, 4A). Reverse zoning in evolved crystals is commonly associated with diffusion as well as surficial and pervasive resorption, i.e., rounded or embayed surfaces for both pyroxenes, with additional sieving for orthopyroxene (Figures 2E,F, 4A). By contrast, normal zoning in mafic crystals is rarely associated with resorption textures, but merely with varying extents of diffusive replacement (Figures 2C, 4A). Single crystals can show several generations of overgrowths (Figure 2E).

Discussion

Evolved and Mafic Magmatic Environments in the PVC Plumbing System

The three pyroxene populations identified in PVC magmas (T1a, T1b, and T2) are abundant throughout the last ∼14.1 ka and across eruptive styles (Figure 4), as exemplified by ortho- and clinopyroxene compositions in the three most recent large eruptions – the YP and PP Plinian eruptions and the intercalated Nealticán flank eruption (Figure 5). Pyroxene populations suggest the sustained presence of at least two distinct magmatic environments in the plumbing system over >14.1 ka: a mafic magmatic environment crystallizing Cr-spinel-bearing T2 pyroxene cores with Mg# >78 for ortho-, and Mg# >82 for clinopyroxene (blue in Figures 4, 5); and an evolved magmatic environment crystallizing T1a pyroxene cores with Mg# 60–71 for ortho-, and Mg# 69–79 for clinopyroxene (yellow in Figures 4, 5) under apatite and Fe-Ti-oxide saturated conditions, as evidenced by abundant apatite and Fe-Ti-oxide inclusions (Figure 4A).

Disequilibrium textures in T1 and T2 crystal cores imply that mixing between mafic and evolved magmas is ubiquitous, consistent with previous studies stressing the importance of magma mixing at the PVC (e.g., Boudal, 1985; Straub and Martín-Del Pozzo, 2001; Schaaf et al., 2005; Mangler et al., 2019). Some T1a crystal cores exhibit extensive sieving associated with elevated Mg# (Figure 2E; white symbols in Figure 5), suggesting partial dissolution in a hotter, more mafic melt (Edwards and Russell, 1996). Conversely, some T2 cores show textural and compositional evidence for extensive diffusion (Figures 2C, 5), suggesting diffusive re-equilibration in a more evolved ME. Not all T1 crystals show strong disequilibrium textures, but all T2 crystals show varying extents of diffusion, which indicates that mafic recharge into the evolved ME (rather than vice versa) controls magma petrogenesis at the PVC. This implicitly places the evolved ME at shallower depths than the mafic ME.

Type 1b pyroxene cores are of particular interest in the context of magma mixing. Compositionally, T1b crystal cores largely overlap with evolved T1a compositions, but they lack apatite inclusions, which indicates that they did not crystallize from the evolved ME. A tendency toward elevated core Mg# (Figures 4, 5) and lower aspect ratios (Figure 4A) exhibited by some T1b crystals point toward a possible hybrid origin for T1b crystals related to mixing between the evolved and mafic magmas (green in Figures 4, 5; see section “Long-Term Diffusive Hybridization of the Crystal Cargo: The Case of T1b Pyroxene”).

Clear evidence for the sustained presence of evolved and mafic magma environments crystallizing ortho- and clinopyroxene in the PVC plumbing system is in contrast to most existing models of the PVC plumbing system, which postulate mafic magmas to be confined to discrete recharge events rather than representing sustained magma bodies (Straub and Martín-Del Pozzo, 2001; Schaaf et al., 2005; Atlas et al., 2006; Cross et al., 2012; Sosa-Ceballos et al., 2012, 2014). We suggest that this is partly due to incomplete datasets resulting from a research focus on explosive eruptions. For example, Straub and Martín-Del Pozzo (2001) report less magnesium-rich orthopyroxene (Mg# <84), and they do not take into account evidence for the compositional gap between mafic and evolved pyroxene, and therefore assume crystallization from a single magma reservoir. On the other hand, Sosa-Ceballos et al. (2012) describe distinct low- and high-Mg# pyroxene populations in PwA pumice samples and relate them to silicic and mafic magmas in the plumbing system, but they do not investigate the nature of the mafic magma any further. Only Witter et al. (2005) postulate the presence of a mafic reservoir in the shallow plumbing system, consistent with our findings.

Dynamics of Magma Mixing and Hybridization

Compositional zoning in ortho- and clinopyroxene offers constraints on the dynamics of magma mixing and hybridization. Pyroxene bands and rims cover the entire range from evolved to mafic compositions for all pyroxene types (Figure 5). Most bands are of intermediate compositions with elevated Mg, Cr, Ca and Al concentrations compared to evolved crystals, producing reverse zoning in T1a and T1b crystals, and normal zoning in T2 crystals (Figures 4, 5). Discrete compositional changes for both slow- (e.g., Cr) and fast-diffusing elements (e.g., Mg) between cores and bands indicate that zoning is controlled by abrupt changes in melt composition (Figure 6; ΔMg# ≤27 for orthopyroxene and ΔMg# ≤13 for clinopyroxene) rather than to degassing-related changes in water concentrations (Frey and Lange, 2011; Waters and Lange, 2017). Such features are consistent with injections of mafic magma into the evolved ME. Initial chaotic mixing results in the juxtaposition of mafic and evolved philocrysts within small-scale transient melt domains of variable, intermediate compositions and temperatures (Petrone et al., 2018). Pyroxene crystals immersed in such transient hybrid magmas crystallize overgrowths of intermediate compositions reflecting the rapidly changing magmatic environment (Wallace and Bergantz, 2005; Bergantz et al., 2015; Petrone et al., 2018). The sharp compositional boundaries between zones (Figure 6) suggest that this process is near-instantaneous.

Intermediate zones sometimes constitute the outermost rims of pyroxene crystals, but more commonly, they are overgrown by more evolved rims with compositions akin to T1a cores. Transitions from intermediate into evolved zones are texturally discrete and recorded by both slow- and fast-diffusing element concentrations (Figure 6), indicating another near-instantaneous process. We suggest that these evolved zones reflect the rapid re-establishment of physico-chemical conditions and architecture of the evolved magma reservoir(s) through efficient melt hybridization (Petrone et al., 2018) and mush reorganization following the recharge event. Resettling and incorporation of philocrysts into colder, degassed parts of the magma body (i.e., the crystal mush) could facilitate near-instantaneous growth of more evolved zones (Humphreys, 2009; Tornare and Bussy, 2014). It should be noted that a rapid return to pre-injection magmatic conditions requires a high ratio of evolved over mafic magma volumes and/or efficient transit of the recharge melt through the reservoir.

Glomerocrysts in PVC lavas and pumices (Figure 2C) provide evidence for long-term crystal storage as a magma mush. Notably, both evolved and mafic pyroxene with variable zoning and diffusion patterns may coexist in a single glomerocryst. Juxtaposition of crystals with different growth histories is commonly inferred to be the result of repeated magma recharge and crystal mush recycling (e.g., Davidson et al., 2007; Cashman and Blundy, 2013; Ubide et al., 2014). In accord with this view, we posit that mafic injections at the PVC disrupt and remobilize the crystal mush in the evolved ME, bringing different populations and generations of philocrysts into direct contact. Numerical models of granular interactions and fluid flow in response to magma recharge into a crystal mush (Bergantz et al., 2015) predict rotating motions in the “mixing bowl” to effectively mix different regions of the mush, which could explain the random dispersal of T1a, T2, and T1b crystals in a texturally and compositionally heterogeneous crystal mush. Some T2 and T1 crystals exhibit more than one reverse zoning pattern (Figure 2E), suggesting that magma mixing and mush reconfiguration is a common process at the PVC.

Magma Storage Temperatures

Geothermometers assume compositional equilibrium between a pair of minerals, or between a mineral and its host melt (e.g., Putirka, 2008). However, compositional and textural disequilibria dominate the philocryst assemblage of the PVC. For example, 73–100% of pyroxene rims (n = 154) and 46–67% of plagioclase rims (n = 79) in PP and YP pumice samples are not in Fe-Mg/An-Ab exchange equilibrium with their host glass (Putirka, 2008; Neave and Putirka, 2017; Supplementary Figures 1, 2). Abundant textural evidence for mixing-induced crystal-melt disequilibrium in PVC lavas and pumices further demonstrates the need to establish a robust thermometer that does not depend on melt compositions. Touching pyroxene pairs in chemical equilibrium are rare in PVC lavas and pumices as a result of frequent disruption and reassembly of the crystal mush. However, almost identical Mg# ranges and trace element differentiation patterns shown by ortho- and clinopyroxene of all units suggest that they are co-crystallizing phases in both the mafic and evolved magmas (Figure 4B and Supplementary Figure 3 for Cr; Supplementary Figures 4, 5 for Ti and Al). We use this chemical evidence for simultaneous pyroxene crystallization in both mafic and evolved MEs to develop a new approach to implementing 2-pyroxene thermometry. Clino- and orthopyroxene record the same magmatic processes at an offset of ΔMg# = 3–6 (Figure 4B), translating into a K_D(Fe–Mg)^cpx–opx = 1.0–1.10, which coincides with the Fe-Mg equilibrium exchange coefficient K_D(Fe–Mg)^cpx–opx = 1.09 ± 0.14 given by Putirka (2008). Equilibrium Fe-Mg exchange conditions between ortho- and clinopyroxene are thus fulfilled if the entire pyroxene dataset is considered, and temperature estimates for the PVC plumbing system can be obtained by testing the entire suite of orthopyroxene measurements (n = 329) for equilibrium against the entire suite of clinopyroxene data (m = 330). Recent work has shown that simple equilibrium criteria such as Fe-Mg exchange coefficients may produce inaccurate thermobarometric results as they are relatively insensitive to kinetic effects during crystallization (e.g., Mollo et al., 2010, 2013; Neave et al., 2019). Crystal-melt partition coefficients vary as a function of melt water content and cooling rate (Mollo et al., 2010; Neave et al., 2019), which may, for example, lead to temperature overestimations for rapidly cooled clinopyroxene-melt pairs (Mollo et al., 2013). Similar dynamic crystallization effects could also affect 2-pyroxene thermobarometry. However, pyroxene textures do not generally indicate large degrees of undercooling in either the evolved or the mafic ME of the PVC. Kinetic effects would most likely affect rim crystallization during magma mixing, but temperatures derived from such pyroxene rims with our method are consistent with those of pyroxene cores. We therefore consider our approach to produce robust constraints on temperatures and pressures of PVC magmatic environments.

Average temperatures for PVC magmas were calculated as a function of orthopyroxene Mg# (Figure 7). Each orthopyroxene was matched with the entire clinopyroxene dataset (e.g., opx₁ ↔ cpx_1–m, opx₂ ↔ cpx_1–m, …, opx_n ↔ cpx_1–m) and m temperatures were calculated for each of the n orthopyroxenes solving equation 36 of Putirka (2008; standard error of estimate = 45°C). Temperatures for opx–cpx pairs that did not fulfill equilibrium criteria were discarded. An average temperature ± 1 SD was then calculated from equilibrium temperatures for each orthopyroxene data point. Individual results were then binned by Mg# and average temperatures for each Mg# calculated. For example, the average temperature for orthopyroxenes with Mg# 80 was calculated from 319 individual temperatures of orthopyroxenes with Mg# 79.5–80.4. Temperatures thus derived for each Mg# (Figure 7) are averages of between 1 and 1990 individual equilibrium temperatures, with 1SD uncertainties of typically 23–31°C (total range including averages based on <100 equilibrations = 6–53°C). Combined with the compositional ranges defined for T2 and T1 pyroxenes, this allows temperatures of the mafic and evolved MEs to be constrained (Figure 7).

Mafic T2 pyroxene (orthopyroxene Mg# >∼78) shows a temperature range of 1000–1115°C (1027 ± 27°C at Mg# 80 and 1074 ± 42°C at Mg# 88) and an average of ∼1040°C. By contrast, evolved T1a pyroxene (orthopyroxene Mg# ≤71) yields temperatures of 940–995°C at an average of ∼970°C. Hybrid compositions between T1a and T2 (orthopyroxene Mg# 71–78) record a sharp decrease in temperature relative to more mafic crystals and show a temperature range of 960–1030°C and an average of ∼990°C (Figure 7).

Equilibrium pressures were calculated analogously to temperatures using equation 39 from Putirka (2008), yielding pressures between 2 and 6 kbar for both T1 and T2 pyroxenes, with uncertainties similar to the standard error of estimate (±2.8 kbar) associated to the method (Supplementary Figure 6). Calculated pressures for the evolved and mafic MEs are therefore indistinguishable from each other and constrain mid- to upper crustal storage depths.

Results of our 2-pyroxene thermobarometry constrain temperatures of 1000–1115°C for the mafic ME, and 940–995°C for the evolved ME, and mid- to upper crustal pressures for both. Temperature and pressure ranges for the evolved ME are consistent with previous estimates for evolved magmas of 0.3–3 kbar and 870–1000°C, with an average of 950°C, based on Fe-Ti-oxide pairs (Schaaf et al., 2005; Witter et al., 2005; Sosa-Ceballos et al., 2012), single pyroxene pairs (Straub and Martín-Del Pozzo, 2001; Witter et al., 2005; Sosa-Ceballos et al., 2012), H₂O and CO₂ concentrations in melt inclusions in pyroxene and plagioclase (Atlas et al., 2006) and apatite saturation temperatures (Witter et al., 2005). Minimum temperature estimates for the evolved ME are >150°C above solidus, implying hot rather than cold storage (cf. Cooper and Kent, 2014; Barboni et al., 2016; Petrone and Mangler, in press) in the PVC’s evolved ME for >14.1 ka. Our results further provide the first evidence for the presence of a sustained mafic magmatic environment at depths < 6 kbar and temperatures of 1000–1115°C. Prolonged storage of mafic magma in the shallow PVC plumbing system has previously only been suggested by Witter et al. (2005), who propose pressures of 1–4 kbar and temperatures of 1050–1300°C based on MELTS modeling. Finally, most pyroxene zones crystallized at hybrid temperatures of 960–1000°C (Supplementary Figure 7). This suggests that mixing of mafic and evolved magmas promotes pyroxene crystallization in the form of bands and rims around crystal cores (Figure 5).

Long-Term Diffusive Hybridization of the Crystal Cargo: The Case of T1b Pyroxene

Mafic T2 crystals injected into the evolved ME are exposed to new, cooler and more evolved conditions, resulting in the crystallization of a more evolved rim. Upon prolonged storage in the evolved ME, initially discrete boundaries between mafic crystal cores and evolved growth rims develop into wide diffusive fronts, gradually replacing mafic cores with more evolved compositions (Figures 6C,D). The width of such diffusive zones in T2 crystals of any given lava or pumice of the PVC varies considerably, from sharp transitions (Figure 6C) up to extensive zones >30 μm (Figure 6D), which corresponds to highly variable crystal residence times of days to centuries in the evolved ME for T2 crystals (Petrone and Mangler, in press). Even more prolonged storage of a T2 pyroxene crystal in the evolved ME may therefore result in near-complete re-equilibration to evolved Mg#, and its mafic origin may no longer be discernible based on non-quantitative imaging techniques or major element compositions. Indeed, lack of apatite inclusions and a tendency toward elevated Mg#, Cr (Figures 4A,B), Ca and Al concentrations (Supplementary Figure 8) in T1b crystal cores is consistent with crystallization from a more mafic magmatic environment. We suggest that T1b pyroxenes initially crystallized as T2 crystals in the mafic ME and were subsequently injected into and stored in the evolved ME long enough for Fe-Mg (but not Cr and Al) diffusive equilibration to obscure their origin.

To estimate timescales of diffusive re-equilibration of a T2 crystal in the evolved ME, we modeled Fe-Mg diffusion in the T2 crystal shown in Figure 4A at average evolved conditions (970°C and fO₂ = NNO + 0.7) until it reached a core-rim Mg# gradient comparable to the T1b crystal shown in Figure 4A (ΔMg# <9; see Petrone et al., 2016 for method). The resulting model yields a minimum timescale of ∼500–1000 years at 970°C to produce such a T1b-like gradient; lower long-term storage temperatures would increase the time needed. This shows that residence times of at least several hundreds to thousands of years in the evolved ME are required to transform a T2 into a T1b crystal, depending on the size of the crystal and its specific thermal history. Formation of T1b crystals by prolonged diffusional re-equilibration of T2 crystals therefore represents a magma mixing-related hybridization process, which acts over much longer periods than the near-instantaneous growth of bands and rims from hybrid melts. Considering the mafic origin of T1b crystals, it follows that 60–67% of pyroxene philocrysts in lavas and pumices of the PVC were of mafic origin (Figure 4C, T1b + T2), and 44–49% (T1b) had been stored for hundreds to thousands of years in the evolved ME prior to their eruption. The crystal mush feeding eruptions of Popocatépetl is therefore dominated by pyroxene crystals of mafic origin.

It cannot be excluded that some T1b crystals have nucleated during mixing in an intermediate melt and represent true hybrid crystals rather than overprinted mafic philocrysts. However, models of nucleation and crystal growth rates during magma mixing by Hort (1998) suggest that relatively small thermal perturbations are compensated predominantly by overgrowth on pre-existing crystals. Hence, nucleation of primary hybrid T1b pyroxene crystals is considered a minor process.

Crystallization Sequence of the PVC

Pyroxene trace element patterns provide a detailed record of crystallization in both the evolved and the mafic magma environments. For example, the steep decline of Cr₂O₃ concentrations in T2 pyroxene (blue in Figure 4B) reflect Cr-spinel crystallization in the mafic magma, and the kinks at the T2–T1 transition (Mg# 78 for orthopyroxene, Mg# 82 for clinopyroxene) show that Cr-spinel ceases to crystallize upon injection into the evolved ME (Figure 4B). Similarly, TiO₂ concentration patterns in pyroxene trace Fe-Ti oxide crystallization in the PVC plumbing system (Supplementary Figure 4). Abundant inclusions of Cr-spinel, Fe-Ti-oxides and other mineral phases provide an alternative crystallization record of mafic and evolved PVC magmas. A survey of the Mg# of pyroxene cores, bands and rims, combined with the types of inclusions present in respective zones, shows consistent patterns in ortho- and clinopyroxene (Figure 8), which can be used to constrain the crystallization sequence of PVC magmas in the mid- to upper crust.

Cr-spinel inclusions are typically present in orthopyroxene with Mg# >77 and clinopyroxene with Mg# >83, coinciding with the respective lower compositional boundaries for T2 pyroxene (Figure 8). This shows that Cr-spinel only crystallizes in the mafic magma, as suggested by Cr₂O₃-concentrations in pyroxene. Cr-spinel inclusions in pyroxenes with lower Mg# are exclusively associated with T1b pyroxene crystals, which represent former T2 crystals that were diffusively overprinted during storage in the evolved ME.

Sulfide globules are found in all types of pyroxene (Mg# 87–65 for ortho-, Mg# 86–71 for clinopyroxene; Figure 8), indicating that PVC magmas become sulfur saturated at an early stage (Larocque et al., 1998). This is consistent with observations of large SO₂ emission rates (≤13 kt/d; Delgado-Granados et al., 2001) at present-day Popocatépetl, which has been associated with deep degassing of mafic magmas (Roberge et al., 2009).

Ilmenite and magnetite inclusions are found in evolved orthopyroxene (Mg# ≤76) and clinopyroxene (Mg# ≤80: Figure 8). Fe-Ti-oxides are not stable in the mafic ME and only crystallize upon mafic recharge into the evolved ME, likely controlled by abrupt changes in melt compositions, oxidation state and water content during magma hybridization (e.g., Jenner et al., 2010). The onset of apatite crystallization coincides with ilmenite and magnetite crystallization (Figure 8), which may indicate that apatite saturation is triggered by magnetite fractionation (Jenner et al., 2010). Most apatite inclusions are found in T1a orthopyroxene crystals with Mg# ≤71 and T1a clinopyroxene with Mg# ≤78, corresponding to temperatures of 975 ± 23°C (Figure 7). This is consistent with apatite saturation temperatures of 945–959°C (Witter et al., 2005) and suggests that the peak of apatite crystallization is only reached after post-injection re-establishment of evolved magmatic conditions. The presence of apatite and Fe-Ti oxide inclusions are therefore robust indicators for the evolved ME.

Plagioclase inclusions are too rare in pyroxene to robustly reconstruct its stability field. However, combined textural and compositional data for YP, Nealticán and PP plagioclase philocrysts shows that most unresorbed plagioclase cores are characterized by <An₅₀, whereas bands and rims associated with surficial or pervasive resorption invariably show >An₄₅ and commonly extend to An₆₀ and higher (Figure 3C). Crystal cores with >An₅₀ commonly show extensive sieving and patchy textures (white symbols in Figure 3C). This suggests that most plagioclase crystallized in the evolved ME, with An-rich bands and rims as well as sieved cores forming as a result of mafic injections. Evidence for plagioclase crystallization in the mafic magmatic environment is rare (Figures 3B,C), indicating that plagioclase is generally suppressed in the mafic magma due to higher pressures and/or water activities (Blatter et al., 2013). Depressurization and degassing during mafic recharge events would increase the stability field of plagioclase and facilitate the crystallization of An-rich bands and rims.

Olivine dissolution and overgrowth by high-Mg# peritectic orthopyroxene (Zellmer et al., 2016) is prevalent in most PVC rocks (Figure 2H), indicating that olivine is generally not stable in shallow PVC magmas, including the mafic ME. Moreover, the conspicuous absence of T2 pyroxenes in the most olivine-rich units (Chipiquixtle, Atlimiyaya and Las Truchas; Figure 3D) suggests that olivine does not crystallize from the same melt as mafic T2 pyroxene. This observation contradicts the petrogenetic model of Witter et al. (2005), who propose crystallization of olivine + high-Mg# orthopyroxene + high-Mg# clinopyroxene from a mafic magma body at 1–4 kbar. Textural evidence rules out olivine crystallization in the evolved or mafic MEs of the mid- to upper crust, and a distinct, deeper magmatic origin of olivine is required.

Collectively, these observations constrain the crystallization sequence of both the mafic and evolved magmatic environments in the mid- to upper crust. The mafic ME crystallizes Cr-spinel + T2 orthopyroxene (Mg# 78–90) + T2 clinopyroxene (Mg# 82–91) ± sulfide ± plagioclase at temperatures of 1000–1115°C, and the evolved ME crystallizes plagioclase (An_28–50) + T1a orthopyroxene (Mg# 71–56) + T1a clinopyroxene (Mg# 64–79) + apatite + Fe-Ti-oxides ± sulfide at 940–995°C. Frequent injections of mid- to upper crustal mafic magma into the evolved ME trigger extensive crystallization of pyroxene and plagioclase in the form of hybrid bands around pre-existing mafic and evolved crystal cores at 960–1030°C. Onset of major apatite, magnetite and ilmenite crystallization at orthopyroxene Mg# 71 (corresponding to 975 ± 23°C) coincides with the subsequent return to pre-injection evolved magmatic conditions. Upon further magma cooling, clinopyroxene becomes unstable, as shown by the lack of groundmass clinopyroxene, or clinopyroxene in equilibrium with orthopyroxene with Mg# <60.

The phase assemblage and crystal compositional ranges of the evolved ME are reproduced by MELTS models (rhyolite-MELTS_v1.1x, Gualda et al., 2012; Ghiorso and Gualda, 2015) of a typical PVC andesite (POP-97; Table 2) under storage conditions inferred for the evolved reservoir (Atlas et al., 2006; Sosa-Ceballos et al., 2014; this work). A representative model of simultaneous cooling (1030–940°C) and degassing (2.6–0.8 wt.% H₂O) at 3 kbar and fO₂ = NNO + 0.7 demonstrates the good fit of modeled orthopyroxene, clinopyroxene and plagioclase with the observed mineralogy (yellow, orange and red lines in Figure 9).

The phase assemblage and compositional ranges of the mid- to upper crustal mafic magma have also been modeled in MELTS using the least evolved PVC whole-rock composition, the magnesian andesitic Chipiquixtle lava (POP-95; Table 2), at conditions inferred for the mafic ME (T = 1000–1115°C; P = 2 – 6 kbar; fO₂ = NNO + 0.7; H₂O = 2.2–3 wt.%). The equilibrium phase assemblage for the modeled temperature range is orthopyroxene (Mg# 76–86) + clinopyroxene (Mg# 79–86), in excellent agreement with the observed phase assemblage and compositions; a representative model (P = 3.5 kbar; H₂O = 3 wt.%) is shown in Figure 9 (blue lines). In addition to orthopyroxene + clinopyroxene, plagioclase (An_49–67) appears at lower H₂O concentrations and temperatures (<1045°C at 2.2 wt.% H₂O), indicating that plagioclase crystallization may be favoured in degassed, cooler magma domains within the mafic ME, consistent with experimental findings (e.g., Blatter et al., 2013). Olivine is not stable alongside clino- and orthopyroxene throughout the modeled range, in agreement with textural observations. However, MELTS models predict olivine crystallization at the expense of clino- and orthopyroxene at high temperatures >1075°C, low pressures <2.5 kbar and water contents >2.3 wt.% H₂O (Supplementary Datasheet 1), and experimental results for water-saturated andesites indicate that olivine stability extends to lower pressures and higher temperatures (Blatter and Carmichael, 1998). Such hydrous, hot conditions favoring shallow olivine crystallization could be given for small batches of rapidly ascending and erupting mafic magma, as suggested for the present-day activity (Atlas et al., 2006; Roberge et al., 2009), but it cannot account for large, Cr-spinel-bearing forsteritic olivine crystals frequently found in PVC rocks (Figures 2G,H).

Primitive Magmas Feeding the PVC: Constraints From Olivine + Low-Ca Clinopyroxene in Flank Eruptions

The good fit of the Chipiquixtle whole-rock composition as an endmember for T2 pyroxene crystallization models at mid- to upper crustal conditions is in stark contrast to its observed phase assemblage (Table 1), which is distinguished by the presence of pristine olivine as well as low-Ca clinopyroxene, while T2 pyroxene is essentially absent. This phase assemblage is also observed in other flank lavas (i.e., Atlimiyaya and Las Truchas; Figures 1, 3D, 10) and offers an opportunity to constrain the nature of mafic magmas feeding the PVC. Mangler et al. (2019) identified distinct geochemical signatures of these lavas consistent with lower crustal contamination, proposing that the erupted magmas largely circumvented the shallow plumbing system of the PVC. The lack of T2 pyroxene supports the hypothesis that the magmas did not traverse the shallow mafic ME, while the presence of T1 pyroxene (Figure 3D) suggest some degree of mixing with the evolved ME. We therefore posit that olivine and low-Ca clinopyroxene crystallized from a deep-seated, primitive magma that bypassed the mafic ME in the mid- to upper crust, but interacted with an evolved magma at shallow depths. Crystallization of olivine + Cr-spinel + clinopyroxene in the lower crust beneath the PVC has been suggested by Straub and Martín-Del Pozzo (2001), and the distinct phase assemblage of Chipiquixtle, Atlimiyaya and Las Truchas lavas (Figure 10) provides the first robust evidence for this model. The olivine + spinel + low-Ca clinopyroxene assemblage can therefore be used to test the viability of a range of primitive endmembers that have been suggested for the PVC (Straub and Martín-Del Pozzo, 2001; Schaaf et al., 2005; Mangler et al., 2019).

Primitive magmas feeding the PVC have been widely suggested to be mantle-derived, high-Mg melts, but their exact compositions remain controversial. Straub and Martín-Del Pozzo (2001) argue that primary (basaltic) andesites with 55–62 wt.% SiO₂ feed the PVC, whereas Schaaf et al. (2005) and Mangler et al. (2019) suggest basalts or basaltic andesites similar to the magmas erupted in the neighboring SCVF. To examine the viability of these potential primitive magmas, we used the Chipiquixtle magnesian andesite (POP-95; 57 wt.% SiO₂), a SCVF calc-alkaline basaltic andesite (53 wt.% SiO2; Straub et al., 2015; Table 2), and a SCVF basalt (49 wt.% SiO2; Straub et al., 2015; Table 2) to test whether they could crystallize olivine + Cr-spinel + low-Ca clinopyroxene under plausible crustal conditions. Each endmember was equilibrated in MELTS over a wide range of P-T-H₂O-fO₂ conditions informed by literature constraints on the primitive melt: temperature (1080–1220°C) and pressure (1–12 kbar) ranges were chosen to cover estimates by Straub and Martín-Del Pozzo (2001) for olivine – clinopyroxene crystallization (1085–1170°C; 3–12 kbar), and their fO₂ value of FMQ −0.2 was adopted. Water contents used (0–6 wt.%) are consistent with reported H₂O-concentrations in olivine melt inclusions of 0.3–5.1 wt.% (Roberge et al., 2009). For every P-T-H₂O combination, the modeled equilibrium assemblage was assessed against the observed phase assemblage, and only low-Ca clinopyroxene with compositions matching the characteristics defined in Figure 10 (Mg# 71–88; CaO/(FeO + MgO) = 0.7–0.97; Al₂O₃ ≤6.5 wt.%; TiO₂ <1.4 wt.%) were allowed. Model results are summarized in Supplementary File 1, and the compositions of modeled low-Ca clinopyroxene are shown in Supplementary Figure 9. The SCVF basalt crystallizes olivine + Cr-spinel + low-Ca clinopyroxene at lower to upper crustal conditions (2–8 kbar, 1080–1220°C, 0–5 wt.% H2O; Figure 11), in agreement with experimental data for oxidized, hydrous basalt (Figure 3B in Blatter et al., 2013). Equilibrium temperatures increase and water contents decrease as pressure decreases (Figure 11). In contrast, models of the SCVF basaltic andesite only match the phase assemblage at a restricted upper crustal range, consistent with experimental results by Blatter and Carmichael (1998), and Chipiquixtle andesite models do not return olivine + Cr-spinel + low-Ca clinopyroxene under any P-T-H₂O-fO₂ conditions (Supplementary File 1). These contrasts in phase assemblages for basalt, basaltic andesite and andesite models are significant even considering the uncertainties related to (Cr-)spinel stability models in MELTS due to a lack of data on the cation-ordering state of (Cr-)spinels under relevant P-T-X conditions and resulting uncertainties in spinel molar volumes (Hamecher et al., 2013).

MELTS modeling results clearly support a hydrous, high-Mg basaltic magma similar to the most primitive melts of the neighboring SCVF to crystallize the olivine + Cr-spinel + low-Ca clinopyroxene assemblage sampled by the Chipiquixtle, Atlimiyaya and Las Truchas. The wide range of equilibrium pressures (Figure 11) is consistent with polybaric fractionation during ascent from the mantle proposed by Straub and Martín-Del Pozzo (2001). However, these authors propose olivine crystallization to cease at lower crustal levels followed by polybaric crystallization of clinopyroxene, whereas our findings suggest simultaneous crystallization of all three phases over a range of crustal depths. Low-Ca clinopyroxene is not only present in flank eruptions clustered around the PwA Plinian eruption ∼14.1 ka BP, but also in the most recent YP and PP pumices (2150 and 1100 a BP), and Straub and Martín-Del Pozzo (2001) report clinopyroxene compositions with low-Ca characteristics in present-day ejecta. This indicates that basaltic injections into mid- to upper crustal magma reservoirs are a common feature at the PVC since at least 14.1 ka. This prescribes a deeper origin for the olivine + Cr-spinel + low-Ca clinopyroxene phase assemblage, which is at odds with the low-pressure crystallization window (<3 kbar) for the SCVF basaltic andesite. Neither the basaltic andesites suggested by Schaaf et al. (2005), nor the magnesian andesites with ∼55–60 wt.% suggested by Straub and Martín-Del Pozzo (2001) can therefore be considered parental melts to olivine + Cr-spinel + low-Ca clinopyroxene assemblage observed in PVC magmas. We instead argue that the most primitive philocrysts of the PVC plumbing system are inherited from a basaltic feeder melt. Basaltic melt inclusions in PVC olivine (Roberge et al., 2009) provide supporting evidence for primary basaltic input into the PVC plumbing system. However, our MELTS models suggest that none of the tested mafic endmembers are related to evolved PVC magmas by straightforward fractionation alone, since resulting magmas would be more K-rich than the PVC rock suite, as suggested by Straub and Martín-Del Pozzo (2001). Another magma source is therefore required, and it is conceivable that the PVC is fed by a wide variety of high-Mg mantle-derived magmas, ranging from basaltic to andesitic compositions similar to those erupted in the SCVF (Schaaf et al., 2005; Straub et al., 2015).

A Modified Petrogenetic Framework for the PVC

Similar textural characteristics and compositional ranges of PVC lavas and pumices erupted during the last ∼14.1 ka suggest that magma petrogenesis at Popocatépetl has been controlled by similar processes and rates throughout this period. Magma mixing is the dominant process in the late-stage petrogenesis of both explosive and effusive PVC rocks. Philocryst populations and assemblages in lavas and pumices point toward at least three different magmatic environments. In the following, we present our petrogenetic model for the PVC based on a crystallization sequence calibrated against orthopyroxene–clinopyroxene compositions and temperatures derived from 2-pyroxene-thermometry (Figure 12), as well as MELTS modeling. A schematic of the magmatic environments of the PVC is shown in Figure 13.

(1) Cr-spinel-bearing olivine + Cr-spinel + low-Ca clinopyroxene philocrysts are inherited from a basaltic mantle melt that feeds the mid- to upper crustal plumbing system of the PVC. Primary basalts are Mg-rich and hydrous (≤5 wt. H₂O) and crystallize Cr- spinel + forsteritic olivine + low-Ca clinopyroxene at lower- to mid-crustal conditions (1080–1220°C, <8 kbar). Previous studies concluded primary melts of the PVC to be high-Mg basaltic andesites (Schaaf et al., 2005) or andesites (Straub and Martín-Del Pozzo, 2001), which we consider to complement rather than contradict our findings. We propose that the PVC is fed by a wide range of primary magmas of basaltic to andesitic composition, similar to the range erupted at the neighboring SCVF (Schaaf et al., 2005; Straub et al., 2013, 2015). Geochemical evidence for lower crustal contamination in the PVC lavas with the most primitive signatures (Mangler et al., 2019) suggests stalling of primitive melts in the lower crust rather than rapid ascent, as proposed by Straub and Martín-Del Pozzo (2001).

(2) A mafic magmatic environment at mid- to upper crustal pressures (2–6 kbar) crystallizes T2 orthopyroxene (Mg# 78–90) + T2 clinopyroxene (Mg# 82–91) + Cr-spinel ± sulfide over a temperature range of 1000–1115°C. Crystallization of An-rich plagioclase may commence at the lower end of the temperature range but is rarely represented in erupted rocks. The mafic magma is inferred to be a magnesian andesite (∼57 wt.% SiO₂) similar to the Chipiquixtle lava, with higher silica contents than suggested by Witter et al. (2005; 53 wt.%). Forsteritic olivine is inherited from primary basalts periodically recharging the magnesian andesite, where it is subject to breakdown and replacement by high-Mg orthopyroxene and oxides (Coombs and Gardner, 2004; Zellmer et al., 2016). Recharge with primary melts may trigger ascent and injection of the shallow mafic magma into the evolved ME. Shallow olivine crystallization (Roberge et al., 2009) may be associated with such rapid ascent of andesitic magma batches and subsequent eruption.

(3) A shallow evolved magmatic environment (0.3–3 kbar; Atlas et al., 2006) crystallizes plagioclase (An_28–50) + T1a orthopyroxene (Mg# 71–56) + T1a clinopyroxene (Mg# 64–79), as well as apatite + Fe-Ti-oxides + sulfide, which are found both as abundant inclusions in pyroxene and dispersed in the groundmass. Inferred temperatures for the evolved ME range from 940 to 995°C, with an average thermal state of ∼970°C, consistent with other geothermometers (Schaaf et al., 2005; Witter et al., 2005; Sosa-Ceballos et al., 2012), and indicating long-term hot storage of PVC magmas. Erupted composition generally show 61–63 wt.% SiO₂, however, groundmass and glass compositions range from ∼55 to 72 wt.% SiO₂, representing the endmembers of mixing between shallow mafic and evolved melts.

Injections of the shallow mafic magma into the evolved ME induce compositional and thermal hybridization, causing enhanced pyroxene crystallization at hybrid temperatures of 960 – 1030°C, predominantly as hybrid bands around pre-existing T1 and T2 pyroxenes. In addition to sustained crystallization of orthopyroxene + clinopyroxene ± sulfides, plagioclase stabilizes during magma hybridization and crystallizes calcic rims (≤An₆₅) around more evolved cores (≤An₅₀) from the evolved ME (Figure 3C). Apatite, magnetite and ilmenite start crystallizing slightly after initial chaotic magma mixing as the system gradually cools down. The onset of apatite crystallization occurs at ≤975°C (opx Mg# 71), marking the return of the magma to pre-injection evolved magma conditions (Figure 12). As the magma continues to fractionate and cool, clinopyroxene becomes unstable, as shown by the lack of clinopyroxene in equilibrium with orthopyroxene of Mg# <60 (Figure 12).

Abundant glomerocrysts in lavas and pumices provide evidence for long-term storage conditions in the evolved ME (Figure 2C). After a mafic recharge event, pyroxene philocrysts resettle to form a texturally and compositionally heterogeneous crystal mush, while the majority of plagioclase remains dispersed in the melt due to their lower density. Crystals may be stored for up to millennia in the crystal mush, causing diffusive re-equilibration of mafic T2 pyroxene to form T1b pyroxene crystals. Between 60 and 67% of philocrysts in the studied lavas and pumices originate from the mid-to upper crustal mafic ME (T2 + T1b crystals), and 69–74% of these crystals were stored in the evolved crystal mush for hundreds to thousands of years before being erupted (T1b crystals; Figure 4C). Mafic recharge into the evolved ME is a recurring process since at least the PwA Plinian eruption ∼14.1 ka, effectuating long-term hot storage of magmas in the PVC plumbing system (cf. Barboni et al., 2016). Each magma injection remobilizes the crystal mush, and magma mixing serves as an agent of magma hybridization on various length- and timescales, and potentially as an eruption trigger (Anderson, 1976; Sparks et al., 1977; Mangler et al., 2019). Literature data indicates that similar processes dominate the present-day activity of Popocatépetl volcano (Straub and Martín-Del Pozzo, 2001).

Conclusion

This study examines the petrogenesis of the Popocatépetl Volcanic Complex through the lens of ortho- and clinopyroxene and for the first time integrates data from both explosive and effusive eruptions of the last ∼14.1 ka. Compositional ranges of philocrysts and groundmass are similar in both lavas and pumices, reflecting a steady-state plumbing system governed by frequent mixing between at least three magma environments defined by distinct pyroxene populations. We use the co-crystallization of ortho- and clinopyroxene throughout the shallow petrogenesis to develop novel approaches to constraining the magmatic temperatures and crystallization history of the PVC, resulting in a temperature- and composition-calibrated crystallization sequence. Primitive basaltic magma crystallizes olivine + Cr-spinel + low-Ca clinopyroxene at 1080–1220°C and feeds the mid- to upper crustal plumbing system, which hosts a mafic (magnesian andesitic) magmatic environment crystallizing high-Mg (T2) orthopyroxene + T2 clinopyroxene + Cr-spinel ± sulfide at 1000–1115°C and 2–6 kbar, and an evolved (up to rhyolitic) magmatic environment crystallizing plagioclase + evolved (T1a) orthopyroxene + T1a clinopyroxene + magnetite + ilmenite + apatite ± sulfide at long-term storage temperatures of ∼970°C and pressures of 0.3–3 kbar.

Complex textural and compositional disequilibrium characteristics of pyroxene and plagioclase philocrysts in PVC rocks reflect repeated mixing events between the three magmatic environments in the shallow crust. Each injection of mafic magma into the evolved magmatic environment results in transient magmas of hybrid compositions and temperatures (960–1030°C), which crystallize intermediate pyroxene bands and rims before the system returns to pre-injection conditions and crystals settle into a heterogeneous crystal mush. Up to 67% of erupted pyroxene philocrysts are of mafic origin, of which up to 74% were stored in the crystal mush of the evolved magmatic environment for centuries to millennia, resulting in re-equilibration to more evolved compositions (T1b pyroxene). Mafic injections into the evolved magmatic environment frequently remobilize the crystal mush and are likely to trigger eruptions of the hybrid andesites of the PVC.

The presence of mantle-derived basaltic magmas in the lower crust feeding a mafic magmatic environment in the mid-to upper crust and more evolved magmas with up to rhyolitic melt compositions at shallow depth beneath the PVC, are consistent with the current and emerging understanding of volcanic plumbing systems as a transcrustal mush system (e.g., Marsh, 2006; Cashman et al., 2017; Edmonds et al., 2019; Sparks et al., 2019). Our detailed study of pyroxene philocrysts provides new constraints on short- and long-term dynamics of magma mixing and hybridization, as well as on the formation and composition of crystal mush in a typical arc magmatic system. Thermal constraints imply long-term hot (i.e., above-solidus) storage of PVC magmas at all crustal levels (cf. Barboni et al., 2016; Petrone and Mangler, in press) fueled by frequent mafic injections from deeper levels, and argue against the emerging notion of cold (near- to sub-solidus) conditions dominating magma storage in the crust (e.g., Cooper and Kent, 2014). While there is clear evidence for mush formation in the evolved magmatic environment, we suggest that eruptible magma was present beneath the PVC at all times throughout >14.1 ka. We suggest that magma petrogenesis at arc volcanoes of comparable size and composition as the PVC may be dominated by similar hot storage conditions and short- to long-term magma hybridization dynamics.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

MM, CP, JP, and HD-G collected the samples. MM, CP, and SH acquired and processed the imaging and microanalytical data. MM led the data analysis and interpretation with contributions from CP and JP. MM wrote the manuscript in close collaboration with CP and JP. All authors contributed to discussions about the conceptual framework.

Funding

This work was supported by a Janet Watson Ph.D. Scholarship from the Department of Earth Science and Engineering, Imperial College London to MM, NERC Grant NE/M014584/1 to CP, RS-Newton International Exchanges Grant IE140605 to CP, and a NHM Collection Enhancement Fund to CP.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank Guillem Gisbert-Pinto for assistance in the field, Enrica Bonato, Tony Wighton, and Callum Hatch for thin section preparations, and John Spratt and Tomasz Goral for assistance with EPMA and SEM analyses. We thank the reviewers for their insightful comments, which helped improve the original version of this manuscript, and the editor for her careful handling.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2020.00362/full#supplementary-material

References