Interplinian effusive activity at Popocatépetl volcano, Mexico: New insights into evolution and dynamics of the plumbing system

Effusive eruptions dominate the eruptive record of many arc volcanoes and may hold crucial information about their plumbing systems, yet they are underrepresented in geochemical and petrological studies. Here, we present whole rock major and trace element data as well as Sr–Nd–Hf isotopic compositions for 14 lava flows and four Plinian eruptions of the Popocatepetl Volcanic Complex (PVC) in the last ~23.5 ka, allowing the first comprehensive geochemical characterisation of the dynamics and evolution of its plumbing system. Lavas and pumices of the PVC are andesites–dacites with a narrow compositional range showing no first-order geochemical trends in the last ~23.5 ka. Trace element and isotope ratios show that PVC magmas are derived from a depleted mantle source with a component of subducted sediments. Assimilation-fractional crystallisation models show that magma compositions are modified to varying degrees by assimilation of lower and upper crust en route to the surface. In the shallow plumbing system, geochemically distinct magmas coexist and undergo extensive mixing and hybridisation, thus buffering erupted whole rock compositions. Only few flank eruptions sample more primitive magmas from deeper reservoirs that circumvented the shallow plumbing system. Some Plinian eruptions caused compositional shifts reflecting reconfigurations of the plumbing system, which also affected subsequent effusive eruptions. Our study thus shows that the geochemical variability of PVC magmas in the last ~23.5 ka is dominated by crustal processes, and magma hybridisation is the primary mechanism to produce the buffered whole rock compositions of the PVC.

Arc volcanoes are characterised by a variety of eruptive styles, ranging from effusive to violent explosive episodes. Stratovolcanoes in arc settings commonly exhibit several types of activity during an eruptive cycle or even a single eruption (e.g. Novarupta volcano, Interplinian effusive activity at Popocatépetl volcano Mangler et al., 20191912: Adams et al. [2006). Understanding the transition between eruptive styles is a major challenge in modern volcanology, both from a scientific and hazard mitigation perspective. Extensive research has deepened our knowledge of how the interplay of volatile and crystal contents, shallow magma degassing, rheological properties, magma flow rate and conduit geometry impacts the style of eruptive activity [e.g. Wilson et al. 1980;Woods and Koyaguchi 1994;Parfitt and Wilson 1995;Dingwell 1996;Gonnermann and Manga 2003;Ripepe et al. 2005;Edmonds and Herd 2007;Platz et al. 2007;Castro and Gardner 2008;Degruyter et al. 2012]. Such investigations are generally focussed on explosive activity and deposits, often neglecting the voluminous effusive eruptions that dominate the activity of many stratovolcanoes worldwide [Davidson, De Silva, et al. 2000]. However, effusive eruptions are key to understanding the evolution of volcanic complexes and may hold important information about the build-up to and aftermath of explosive eruptions [e.g. Ruprecht and Bachmann 2010;Koleszar et al. 2012;De Maisonneuve et al. 2012;Tepley et al. 2013].
Popocatépetl is an andesitic-dacitic arc stratovolcano with highly variable eruptive activity, situated in a densely populated area in central Mexico (Figure 1). The current eruptive episode started in December 1994, after a 67-year repose period, and has since emplaced and destroyed at least 38 domes [Gómez-Vazquez et al. 2016]. Present-day activity is associated with mild to moderately explosive Vulcanian activity, and ashfall frequently reaches Mexico City and the city of Puebla ( Figure 1). However, past activity includes at least six Plinian eruptions in the last ∼23.5 ka and the emplacement of at least ten voluminous lava flows. Previous work has focussed on the present-day activity as well as on the Plinian eruptions, whereas only scarce and unsystematic data are available on past interplinian activity. In this study, we present a comprehensive geochemical characterisation of both effusive and explosive eruptions in the last ∼23.5 ka, which allows us to fully reconstruct the evolution of Popocatépetl's plumbing system, and to identify magmatic processes related to the transition from effusive to explosive activity.

Eruptive history of Popocatépetl volcano
Popocatépetl is a Quaternary volcanic complex in Central Mexico, at the front of the Central Trans-Mexican Volcanic Belt (TMVB), an E-W trending volcanic arc generated by the oblique subduction of the Cocos and Rivera plates under the North American plate [Ferrari et al. 2012]. It is the southernmost volcanic centre in a N-S oriented, southward-younging volcanic chain known as the Sierra Nevada, comprising Tláloc, Iztac-cíhuatl and Popocatépetl volcanoes (Figure 1). Modern Popocatépetl is constructed on the remnants of at least two earlier volcanic edifices that were partially destroyed by sector collapses [Mooser 1958;Robin and Boudal 1987;Espinasa-Pereña and Martín-Del Pozzo 2006;Sosa-Ceballos et al. 2015;Siebe et al. 2017]. In this work, Popocatépetl will be used to refer to the modern cone, whereas the term Popocatépetl Volcanic Complex (PVC) refers to the entire volcanic centre including the remnants of previous edifices.
Previous studies relate the earliest preserved volcanism at the PVC to the Nexpayantla edifice (∼330 ka-183±7 ka BP; Sosa-Ceballos et al. [2015]), the remnants of which constitute the NW flank of the modern edifice [Espinasa-Pereña and Martín-Del Pozzo 2006;Cadoux et al. 2011;Sosa-Ceballos et al. 2015]. After the sector collapse of the SSW flank of Nexpayantla volcano, construction of Ventorrillo volcano began inside the collapse caldera, with the oldest dated rocks yielding an 40 Ar/ 39 Ar-age of 96±8 ka BP [Sosa-Ceballos et al. 2015]. Ventorrillo lavas and pyroclastic deposits make up large parts of the proximal cone W of the current vent [Espinasa-Pereña and Martín-Del Pozzo 2006;Sosa-Ceballos et al. 2015], while other Ventorrillo lavas have travelled up to 16 km from the vent in the northeast (Xalitzintla) and east (Ixtepec and Metepec) (Figure 2). The southern sector of Ventorrillo partly collapsed at ∼23.5 ka BP, triggering the Tochimilco Pumice Plinian eruption (TP), which was assigned a VEI of 5 (Volcanic Explosivity Index; Newhall and Self [1982]) and described in detail as the White Pumice Plinian eruption by Siebe et al. [2017].
The collapse of the Ventorrillo edifice marked the onset of modern Popocatépetl, construction of which was characterised by voluminous effusive activity punctuated by at least five Plinian eruptions as well as mild explosive periods such as the present-day activity [Siebe et al. 1996;Panfil et al. 1999;Siebe and Macías 2006;Plunket and Uruñuela 2008;Arana-Salinas et al. 2010;Sosa-Ceballos et al. 2012;Siebe et al. 2017]. Early lava flows were emplaced within the collapse caldera of the Ventorrillo volcano towards the south and southeast [ Figure 2; Siebe et al. 2017]. The subsequent Pumice with Andesite Plinian eruption (PwA) at ∼14.1 ka BP is considered the largest explosive eruption of modern Popocatépetl [VEI 6;Siebe et al. 1999;Espinasa-Pereña and Martín-Del Pozzo 2006;Sosa-Ceballos et al. 2012;Siebe et al. 2017]. The PwA produced a combined 5.0 km 3 of tephra [Siebe et al. 2017] during an initial stage of Vulcanian to subplinian eruptions that escalated into the cataclysmic main stage [Sosa-Ceballos et al. 2012]. After the PwA Plinian eruption, the primary direction of lava emplacement rotated counterclockwise to the NE, where the youngest large-scale lavas are found ( Figure 2).
In addition to summit eruptions producing largely proximal lavas, flank eruptions associated with a SW-NE trending fault system account for most of the recent lavas, especially in the distal sectors. For instance, the youngest prominent lavas, Nealticán and Capula, originated from a fissure on the NE flank, and the Nealticán flow reached more than 17 km in length ( Figure 2). These two flows are bracketed by the two most recent Plinian eruptions, 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 [Siebe et al. 1995;Siebe et al. 1996;Panfil et al. 1999;Siebe and Macías 2006;Plunket and Uruñuela 2008]. The YP Plinian eruption consisted of an initial 'blast', followed by the main Plinian pumice fall and subsequent ignimbrite deposition, producing 3.2 km 3 of tephra [Siebe et al. 1995;Panfil et al. 1999;Siebe and Macías 2006]. The PP eruption produced <1.69 km 3 tephra during three pulses of pumice fall, accompanied by surges, ashfall and lahars [Siebe et al. 1995;Siebe et al. 1996;Panfil et al. 1999;Siebe and Macías 2006]. Historical accounts show that Popocatépetl's activity since the PP eruption has been dominated by prolonged dome-building and -destruction cycles similar to the eruption ongoing since 1994 [Martín-Del Pozzo et al. 2016]. Vulcanian and Strombolian eruptions related to this activity have produced ash, lapilli and ballistics covering the upper flanks of Popocatépetl [Espinasa-Pereña and Martín-Del Pozzo 2006], and occasional larger eruptions (VEI 3) in the 1500s as well as in 1664, 1997CE have produced pumice [Martín-Del Pozzo et al. 2016].
for major, trace and radiogenic isotope compositions [Schaaf et al. 2005]. Major and some trace element data are available for the Metepec [Robin 1984;Conte et al. 2004;Espinasa-Pereña and Martín-Del Pozzo 2006] and Buenavista [Nixon 1989] lavas. Boudal and Robin [1988] present some major and trace element data for lavas from the east flank, but they lack GPS coordinates that would allow discrimination of lavas. To the best of our knowledge, no geochemical data have been published on the SW Ventorrillo, Xalitzintla, Ixtepec, La Colonia, Ecatzingo, Las Truchas, Olga and Capula lavas.

Sampling strategy
In view of the limited exposure of lavas at the densely vegetated PVC, and the preservation issues associated with minor explosive eruptions, every effort has been made to produce the most detailed and comprehensive record of eruptive activity at the PVC in the past ∼23.5 ka. Ten lavas of the most recent Popocatépetl edifice and four lavas of the preceding Ventorrillo edifice were identified, mapped and sampled in three field campaigns during 2014-2016 ( Figure 2). Stratigraphic relationships between the units were established and a geological map produced based on field observations in conjunction with a detailed survey of aerial photographs, digital elevation models, and published stratigraphic columns [ Figure 2; Espinasa-Pereña and Siebe et al. 2017]. Sample locations are shown in Figure 2 and GPS coordinates are given in Supplementary Data 1 (Table S1.1). Whenever possible, each effusive unit was sampled at several locations covering different lobes, however, inaccessible terrain and thick sediment cover prevented extensive sampling in many cases. Excellent exposure of the most recent Nealticán lavas allowed for detailed sampling of discrete proximal and distal lobes, which are interpreted to represent early (POP-101), main  and late (POP-99) stages of the prolonged fissure eruption ( Figure 2).
In addition to the suite of 37 lava samples from 19 locations around the volcano, pumice samples of the four most prominent Plinian eruptions of the last ∼23.5 ka were collected (Figure 2). The TP Plinian eruption and the entire eruptive sequence of the PwA Plinian eruption were sampled at single locations, and samples of the PP and YP Plinian eruptions were collected at four different localities, respectively. In addition, two post-PP pumice samples (<1100 ys BP) and a ballistic block related to present-day activity were collected ( Figure 2). The Ochre Pumice Plinian eruption, dated at ∼5 ka BP and described by Arana-Salinas et al. [2010], was not Presses universitaires de �rasbourg
Finally, five crustal xenoliths were sampled, including two dm-sized granodiorites from the PwA deposits, as well as a skarn, a fine-grained Tertiary rhyolite with contact reactions, and a ∼1 m 3 block composed of monomineralic anhydrite from the Nealticán lava.

Analytical methods
Care was taken in the field to sample fresh lava, avoiding stained surfaces and rocks showing signs of alteration such as incipient pyroxene and plagioclase weathering. Whenever possible, large pumice clasts (>5 cm) were sampled and weathered crust removed. Samples were cleaned in MQ water (Milli-Q, resistivity >18.2 MΩ.cm at 25 • C) and powdered to particle sizes <2 µm using a Fritsch Pulverisette jaw crusher and a Retsch PM100 planetary ball mill with a 250 mL agate drum at the Natural History Museum, London. A detailed discussion of the analytical methods used in this study can be found in Appendix A, and data on the accuracy and repeatability of respective methods are compiled in Supplementary Data 1.
Powdered bulk rock samples were analysed for major element compositions at Activation Laboratories Ltd., Ancaster, Canada, using X-ray Fluorescence analysis. United States Geological Survey (USGS) reference materials BCR-2, AGV-2 and BHVO-2 were processed as unknowns, showing a relative deviation from recommended values of <5 % and a repeatability of <5.5 % 2 SD for major elements (see Supplementary Data 1: Table S1.2). Trace element concentrations were determined on an Agilent 8800 Triple Quadrupole inductively coupled plasma mass spectrometer (ICP-MS) at the Open University, Milton Keynes, and on an Agilent 7700x ICP-MS at the Natural History Museum, London. Replicate measurements of USGS reference material BCR-2 (n = 10) indicates a repeatability of <7 % 2 SD and a relative deviation from recommended values of <5 % for most elements (see see Supplementary Data 1: Table S1.3).
Sample digestion, chromatographic separation and analysis of Sr and Nd isotopic compositions of all samples, as well as Hf isotopic compositions on a representative subset of 21 samples were determined at the MAGIC laboratories at Imperial College London. Strontium isotopic compositions were analysed on a Thermo Finnigan Triton thermal ionization mass spectrometer. The instrumental precision within a single sample run is ≤6 ppm 2 SE. Aliquots of USGS reference material BCR-2 digested with the samples and measured during the analysis period yield an average of 87 Sr/ 86 Sr = 0.705008 ± 7 (2 SD; n = 14; see see Supplementary Data 1:  [Weis et al. 2007]. Reported 143 Nd/ 144 Nd ratios are averages of multiple sample measurements, and the 2 SD uncertainties given are either the standard deviations of these averages, or of the BCR-2 averages analysed with the samples, whichever is larger.

Major and trace elements
Whole rock major and trace element concentrations of PVC lavas, pumices and xenoliths are compiled in Supplementary Data 1 (Table S1.1). Lavas and pumices from the Popocatépetl and Ventorrillo edifices are andesites and dacites (Figure 3), defining a narrow subalkaline range with dominantly intermediate silica contents (SiO 2 = 61-65 wt.%). Lower silica contents observed for some pumices coincide with low alkali metal concentrations and high LOI values, and are interpreted to reflect secondary alteration processes. The majority of effusive and explosive samples fall within the published range of Plinian eruptions of the past ∼23.5 ka [dashed field in Figure 3; Siebe et al. 1999;Schaaf et al. 2005;Larocque et al. 2008;Arana-Salinas et al. 2010;Sosa-Ceballos et al. 2012;Sosa-Ceballos et al. 2014]. This compositional range is significantly smaller than published data for the present-day activity [white field in Figure    pared to primitive mantle [Sun and McDonough 1989], a steep LREE slope (La N /Sm N = 2.1-2.8), pronounced negative anomalies for the high field strength elements Nb, Ta and Ti, and a strong positive Pb-anomaly. However, Zr and Hf do not display negative anomalies, and no significant Eu anomaly is observed (Eu/Eu* = 0.83-1.01). The HREE show flat patterns (Tb N /Yb N = 1.3-1.7), with the steepest slopes for the PwA and PP explosive eruptions. Compositional differences between different units of the Popocatépetl and Ventorrillo edifices are subtle and manifest themselves in absolute enrichment levels rather than in changes in overall trace element patterns. The most enriched samples are the PwA pumices and the units since the YP Plinian eruption ∼2150 a BP, which are on average 15 % more enriched than the lavas erupted between the PwA and the YP Plinian eruptions. Overall, trace element enrichments and patterns of the PVC resemble those of the bulk continental crust (Figure 4). Table 1 presents Sr-Nd-Hf isotopic compositions of the Popocatépetl and Ventorrillo edifices. Both lavas and pumices fall within the mantle array for Sr and Nd isotope compositions (inset in Figure 5), and most overlap with literature data for Plinian eruptions of the PVC  Figure 4: Trace element concentrations normalised to primitive mantle [Sun and McDonough 1989] for major Plinian eruptions and the respective interplinian averages, based on all effusive samples. Present-day activity data from Schaaf et al. [2005]; Torres-Alvarado et al.

Sr-Nd-Hf isotopes
[2011]. Inset shows PVC range compared to bulk continental crust [Rudnick and Gao 2003] and representative data for high-Nb and calc-alkaline (CA) lavas from the neighbouring Sierra Chichinautzin Volcanic Field (SCVF) [Straub et al. 2015].
[dashed field in Figure 5; Siebe et al. 1999;Schaaf et al. 2005;Arana-Salinas et al. 2010]. Chipiquixtle, the most primitive unit in terms of major and trace elements, shows the most radiogenic 143 Nd/ 144 Nd ratios and least radiogenic 87 Sr/ 86 Sr ratios of this study, while literature data of present-day activity extends to even more radiogenic Nd isotopic compositions [white field in Figure 5; Schaaf et al. 2005 [Siebe et al. 1999;Schaaf et al. 2005;Arana-Salinas et al. 2010] are shown. Inset shows observed range from PVC compared to high-Nb and calc-alkaline series of the SCVF [Straub et al. 2015].
Nd isotope compositions, the YP samples display a bimodal distribution.

Discussion
In the following sections, we use our new chemical constraints to investigate the petrogenesis of the PVC. First, source characteristics of the sub-arc mantle are constrained. We then examine magmatic evolution during the last ∼23.5 ka to identify potential geochemical trends over time, changes related to transitions from effusive to explosive activity and vice versa, and differences between summit and flank eruptions. New data for xenoliths sampled from PVC lavas and pumices are used to revisit the potential for lower and upper crustal contamination to modify PVC magmas. Finally, we discuss our findings in the context of magma mixing and hybridisation to constrain the architecture and dynamics of the plumbing system of the PVC.

Mantle source characteristics
The TMVB produces volcanic rocks with starkly contrasting geochemical signatures in close spatial association, which is widely attributed to a heterogeneous mantle source [e.g. Luhr 1997; Wallace and Carmichael 1999;Ferrari et al. 2001;Petrone et al. 2003]. In the Sierra Chichinautzin Volcanic Field (SCVF) west of the PVC (Figure 1), calc-alkaline monogenetic centres coexist with OIB-type, high-Nb lavas on a km-scale [Figures 4, 5 and 6;Verma 1999;Siebe 2004;Schaaf et al. 2005;Straub et al. 2013]. This spatial association has been argued to reflect small-scale mantle   Figure 6: εHf vs. εNd of selected Ventorrillo and Popocatépetl lavas and pumices. Present-day terrestrial array from Vervoort et al. [2011], SCVF data from Straub et al. [2015], PVC data from Straub et al. [2015] and Cai et al. [2014]. heterogeneities, with depleted and undepleted mantle domains of <km size in the mantle wedge [e.g. Siebe 2004]. To explain this heterogeneity, Wallace and Carmichael [1999] proposed corner flow as a mechanism to transport unmodified back-arc mantle into the depleted melting region beneath central Mexico, while Ferrari [2004] invoked propagation of a slab detachment to cause upwelling of deep-seated mantle material.
While OIB-type lavas are thought to derive from primitive mantle domains [Gomez-Tuena et al. 2006;Cai et al. 2014;Straub et al. 2015], several contrasting petrogenetic models have been proposed to explain geochemical features of calc-alkaline lavas of the central TMVB.
• Partial melting of sub-arc mantle and veins of pyroxenites formed by mantle reaction with hydrous slab components [e.g. Luhr 1997;Straub et al. 2008;Straub et al. 2011;Straub et al. 2013]; • Melts of the subducting slab and overlying sedimentary strata contribute to partial melts of variably depleted mantle [Gomez-Tuena et al. 2006;Cai et al. 2014]; • Subduction-eroded fore-arc granodioritic basement and altered oceanic crust are added to the depleted mantle in the form of 'slab diapirs' before melting [Straub et al. 2015].
Our new geochemical dataset allows an evaluation of magma petrogenesis at the PVC with respect to existing models of the sub-arc mantle beneath the Central TMVB.
Compared to the bimodal compositional trends of the SCVF (calk-alkaline vs. OIB-type), PVC rocks show little variability in terms of major and trace element compositions as well as Sr-Nd-Hf isotopes (Figures 4,

Presses universitaires de �rasbourg
Page 51 5 and 6). Trace element and Sr-Nd-Hf isotopic ratios of the PVC overlap with the calc-alkaline series of the SCVF (Figures 4 and 5). The PVC rock suite shows no affinity to OIB-type, high-Nb characteristics, which are defined by less enriched LILE and more enriched HREE compared to the PVC (Figure 4) and systematically lower radiogenic 87 Sr/ 86 Sr ratios (inset in Figure 5). Hence, PVC rocks are akin to the calcalkaline series of the SCVF, with no indication for significant mantle source heterogeneity or primitive mantle influence. Indeed, subchondritic Nb/Ta ratios (average Nb/Ta = 13) at Nb = 4-9 ppm indicate a depleted mantle source ( Figure 7A), which has also been proposed for calc-alkaline rocks of the SCVF [Straub et al. 2015]. The chemical characteristics of our sample suite are thus consistent with previous studies of the PVC suggesting a sub-arc background mantle depleted by multiple melting events based on relatively low HFSE abundances and depleted Sr-Nd isotopic signatures [Straub and Martín-Del Pozzo 2001;Schaaf et al. 2005;Straub et al. 2015]. In terms of additions to the depleted mantle beneath the PVC, Straub and Martín-Del Pozzo [2001] invoked melts of subducted slab components, and Gomez-Tuena et al. [2006] extended this hypothesis to propose that direct slab melts dominate the trace element and Sr-Nd-Pb characteristics of Central TMVB rocks, including the PVC. However, a significant contribution of slab melts to PVC magmas is unlikely considering the low Sr/Y ratios <40 observed (Y = ∼20 ppm). Moreover, ƒa correlation between Sr/Y and 87 Sr/ 86 Sr ratios towards the Sr isotopic ratio of Pacific MORB would be expected if subducted oceanic crust played a key role in the petrogenesis of PVC magmas [Gomez-Tuena et al. 2006], however no such correlation is apparent. Furthermore, the slab melt model of Gomez-Tuena et al. [2006] predicts a correlation between SiO 2 concentrations and Nb/Ta ratios, which is not observed at the PVC. Similarly, low Ba/Th ratios argue against significant contributions of altered oceanic crust to PVC rocks [Elliott 2003], in contrast to other Quaternary TMVB magmas, which show elevated Ba/Th ratios ( Figure 7B).
Significant contributions of slab melt or altered oceanic crust to the depleted mantle beneath the PVC are thus considered unlikely. Instead, enriched LREE (La/Sm N = 2-3, Figure 7B) and a strong positive Pb anomaly (Figure 4) suggest sediment input to the depleted mantle wedge beneath the Central TMVB. This is consistent with the findings of Schaaf et al. [2005], who interpret high Pb concentration anomalies in PVC rocks to indicate 'fluxing' of background mantle with subducted sediment, similar to the pyroxenite vein model proposed for the TMVB by Straub et al. [2008], Straub et al. [2011], Straub et al. [2013], and Straub et al. [2013]. The potential addition of subduction-eroded granodiorite to the sub-arc mantle beneath the PVC as slab diapirs [Straub et al. 2015] is more difficult to evaluate, because of significant compositional overlap be- The overall compositional range of PVC rocks is narrow compared to the geochemical variability of the neighbouring SCVF, and there is no evidence for a heterogeneous mantle source such as proposed for the latter. To further constrain the stability of the mantle source, the geochemical evolution of the PVC will be reconstructed using our comprehensive dataset of effusive and explosive eruptions in the last ∼23.5 ka.

Temporal geochemical evolution of PVC magmas
Long-lived arc volcanoes commonly show significant changes in magma compositions over time. This can express itself either as a secular evolution related to a gradual change in the magmatic source (e.g. Cayambe: Presses universitaires de �rasbourg Page 54  Tormey et al. [1995]), or as abrupt changes associated with temporal or permanent reconfigurations of the plumbing system due to large explosive eruptions and/or sector collapse events (e.g. Stromboli: Francalanci et al. [1989]; Colima: Crummy et al. [2014]).
The PVC shows no systematic first-order secular geochemical evolution since >23.5 ka BP with respect to major and trace element and Sr-Nd-Hf isotope compositions ( Figure 8). Effusive and explosive eruptions define long-term narrow compositional ranges, highlighted as the 'main PVC range' in Figure 8, which includes both the oldest lavas and the present-day activity. This substantiates the hypothesis of a homogeneous and stable mantle source, and it further suggests that the plumbing system has buffered the compositions of effusive and explosive eruptions to some degree. However, significant compositional variability is evident both within the main PVC range and beyond ( Figure 8). Some of these excursions manifest as isolated flank eruptions sampling more primitive compositions (i.e. Chipiquixtle and Atlimyaya flank eruptions; Figure 8), which might suggest that these eruptions largely bypassed the shallow plumbing system and thus preserve distinct magmatic signatures. On the other hand, Plinian eruptions can cause longerlasting compositional shifts (i.e. PwA and YP; Figure 8), which dominate the geochemistry of subsequent effusive eruptions, thus shaping the geochemical evolution of the PVC. Subtle changes in isotope compositions of effusive eruptions before and after the 14.1 ka BP PwA Plinian eruption (yellow in Figure 8) indicate that the PwA caused a reconfiguration of the plumbing system, resulting in a compositional shift of erupted magmas. Indeed, the initial stage of the PwA Plinian eruption (samples POP-01-04, Table 1; Figure 9) is characterised by more radiogenic 87 Sr/ 86 Sr ratios of ∼0.70444 and the highest MgO recorded for explosive eruptions at the PVC, while the cataclysmic main stage (samples POP-05-06) is associated with less radiogenic 87 Sr/ 86 Sr ratios of ∼0.70424 and more evolved compositions (Figure 9). Neodymium isotopic ratios show a similar albeit less pronounced shift at the transition from initial to main stage of the PwA eruption (Figure 9). This is consistent with isotopic data for minerals and glass reported by Sosa-Ceballos et al. [2014], who suggest that during the early stages of the PwA Plinian eruption a deeper, less evolved magma reservoir with more radiogenic Sr compositions was tapped, and later stages erupted shallower magmas with less radiogenic Sr isotopic signatures. Our geochemical data for the PwA eruption supports this hypothesis, and Sr isotope data for pre-and post-PwA effusive eruptions further suggest that this shift in the feeding reservoir was longlasting: Sr isotopic signatures of pre-PwA lavas (e.g. La Colonia and Ecatzingo) are similar to those of the initial  stages of the Plinian eruption, whereas post-PwA lavas scatter around the less radiogenic 87 Sr/ 86 Sr ratios of the main PwA stage (e.g. Las Truchas and Olga; Figure 8). The only outliers to this pattern are the Chipiquixtle and Atlimiyaya flank eruptions, which are inferred to have circumvented the PVC's shallow plumbing system, and the Buenavista flank eruption (marked with a 'B' in Figure 8), which has been argued to be connected to the magmatic system of the neighbouring Iztaccíhu- atl volcano rather than directly to the PVC [Nixon 1989]. We therefore suggest that the magma reservoir feeding pre-PwA effusive eruptions was largely emptied in the initial stage of the PwA eruption, and that post-PwA lavas were instead predominantly fed by the magma reservoir tapped during the main stage of the PwA Plinian eruption.
The absence of systematic geochemical trends in lavas erupted between the PwA and the YP Plinian eruptions is considered to reflect the relative stability of the plumbing system between 14.1 ka and 2150 a BP. A tendency towards lower MgO concentrations during this time period (Figure 8) may indicate some degree of melt fractionation, and further secondary variation may be due to recharge with more primitive magma batches and small-scale heterogeneities within the feeding magma reservoir. However, during the YP eruption, Sr isotopic compositions pinpoint a major shift in the main magma input, from unradiogenic pre-YP 87 Sr/ 86 Sr ratios to some of the most radiogenic values observed at the PVC (∆ 87 Sr/ 86 Sr = 388 ppm; Figure 8). A reciprocal pattern is observed for 143 Nd/ 144 Nd ratios, where the system is reset to the least radiogenic values of the PVC (∆εNd = −1.4). Notably, the distinct shift in isotopic signatures is not reproduced in major and trace element compositions ( Figure 8). Like the PwA, the YP Plinian eruption caused a compositional shift that reverberated in the subsequent effusive activity: the Capula flank eruption as well as a proximal lobe of the following Nealticán eruption (interpreted to represent its early stages; POP-101; Figure 2) show major, trace and isotopic compositions similar to the YP eruption ( Figure 8). By contrast, distal lava lobes of the Nealticán eruption (inferred to represent later stages; Figure 2) trace the gradual evolution from extreme YP/Capula isotopic signatures back to 'main PVC range' characteristics (Figure 8). This may indicate that during the YP eruption, an evolved magma reservoir with distinct geochemical characteristics was tapped, which also fed the subsequent Capula eruption, before being exhausted and hybridised with the main feeding reservoir during the Nealticán effusive eruption. As a result of this hybridisation process, 'main PVC range' compositions were restored by the time of the PP Plinian eruption. Pumices of the PP eruption show compositions similar to pre-YP lavas, with the exception of Sr isotope ratios, which are more radiogenic in the PP pumices ( Figure 8). Presentday activity is characterised by significant geochemical variability [Figures 3, 4 and 8;Straub and Martín-Del Pozzo 2001;Schaaf et al. 2005;Witter et al. 2005;Torres-Alvarado et al. 2011], which is inferred to reflect Presses universitaires de �rasbourg
Lavas and pumices of the PVC in the past ∼23.5 ka show no first-order geochemical evolution, pointing to a relatively homogeneous mantle source and a feeder system buffering erupted compositions to a certain extent. However, secondary compositional variations over time suggest a complex plumbing system with several, compositionally distinct magma reservoirs (Figures 8 and 9). To explain this compositional variability, crustal processes need to be invoked. Fractional crystallisation and magma mixing can account for variations in major and trace element compositions, but they do not cause isotopic fractionation. Therefore, while both processes are relevant for the PVC [e.g. Schaaf et al. 2005], neither can explain the observed isotopic range by itself. However, crustal assimilation can introduce isotopic variability. To constrain the role of crustal assimilation at the PVC, the diverse suite of xenoliths collected from the PwA and Nealticán units is used to model assimilation-fractional crystallisation (AFC) processes after DePaolo [1981].

Crustal assimilation modelling
Geophysical evidence suggests a 40-50 km thick continental crust beneath the PVC [Ferrari et al. 2012], comprised of Precambrian and Paleozoic terranes [Ortega-Gutierrez et al. 1995;Ferrari et al. 2012], Mesozoic limestones, evaporites, shales and sandstones, and Tertiary terrigenous and volcanic rocks [Fries 1965;Goff et al. 1998;Siebe 2004;Schaaf et al. 2005;Espinasa-Pereña and Martín-Del Pozzo 2006;Siebe and Macías 2006;Sosa-Ceballos et al. 2014]. Despite the thick continental crust underlying the Central TMVB, crustal assimilation is widely considered negligible for the mantle-derived magmas of the region [e.g. Martínez-Serrano et al. 2004;Straub et al. 2013;Cai et al. 2014;Straub et al. 2015]. Correlations of isotopic ratios with fractionation indices, which are typically attributed to AFC processes, have been described for high-Nb lavas of the SCVF [e.g. Straub et al. 2013]. However, these correlations have been interpreted to reflect addition of a silicic slab component to the mantle wedge rather than crustal assimilation [e.g. Straub et al. 2013;Straub et al. 2014;Straub et al. 2015;Cai et al. 2014]. On the other hand, the abundance of granodiorite, rhyolite, skarn and marble xenoliths in both pumices and lavas of the PVC (up to 55 vol.% in the main PwA stage; Sosa-Ceballos et al. [2012]), some of which are partially digested, prompted the examination of potential crustal overprints in PVC magmas. Goff et al. [1998] and Goff et al. [2001] suggest carbonate and evaporite assimilation to account for elevated CO 2 and SO 2 degassing during present-day activity. Furthermore, the PVC's elevated 87 Sr/ 86 Sr ratios at similar 143 Nd/ 144 Nd ratios compared to the neighbouring SCVF high-Nb   lavas ( Figure 5) have been interpreted to be caused by assimilation of local carbonates [Siebe 2004;Schaaf et al. 2005]. However, Sosa-Ceballos et al. [2014] rule out significant upper crustal assimilation at the PVC as their AFC models of skarn and sandstone assimilation fail to reproduce the compositional range of the PVC. They acknowledge that minor carbonate assimilation is likely, but argue that isotopic compositions are similarly enriched at Nevado de Toluca, which does not overlie shallow carbonates [Sosa-Ceballos et al. 2014]. Instead, they propose partial melting of variably sediment-modified mantle and subsequent assimilation of lower crustal rocks to explain the isotopic range of the PVC [Schaaf et al. 1994].
Variations of Sr and Nd isotopic compositions with MgO ( Figure 10 Capula flow, some Nealticán samples and one Ventorrillo lava at the evolved end (blue field in Figure 10). In addition, the Atlimiyaya flank eruption consistently forms part of the trend defined by Chipiquixtle-YP. Samples of the Chipiquixtle-YP trend overlap with the main PVC range for MgO concentrations of 4-5 %, and some units within this range cannot be unequivocally attributed to the one or the other group, including the PwA pumices ( Figure 10). However, all the eruptions that exceed the main PVC range in Figure 10 form part of the Chipiquixtle-YP trend.
The correlation of isotopic ratios with fractionation indices displayed by the Chipiquixtle-YP trend has not previously been described for rocks of the PVC and suggests that geochemical variations are dominated by crustal contamination. To test whether assimilation of local crust can produce the Chipiquixtle-YP trend, and to evaluate the role of crustal contamination in the main PVC range, a suite of xenoliths collected from the PVC (Table 1, Supplementary Data 1: Table S1.1) was used to model AFC processes after DePaolo [1981]. Assimilation-fractional crystallisation models were produced for Sr, Nd and Hf isotopic ratios, La/Nb ratios and Mg concentrations. The least evolved sample of the dataset, POP-95 (Chipiquixtle), was used as the magmatic endmember. A skarn (POP-97sk) and a Tertiary rhyolite (POP-97sst) xenolith sampled from the Nealticán lava were selected as shallow crustal assimilants (Table 1, Supplementary Data 1:  Table S1.1), since alongside comagmatic granodiorites, these lithologies are the most common xenoliths in PVC lavas and pumices [Schaaf et al. 2005;Sosa-Ceballos et al. 2014]. To assess the prospect of lower crustal assimilation, central Mexican metasedimentary granulite sample JC 101 from Schaaf et al. [1994] was used. Since no Hf isotope data is available for JC 101, Hf isotopic ratios were modelled using metapelitic granulite sample OAX-5-86 Vervoort et al. [2000]. Modelling parameters are given in Table 2, and the results shown in Figure 11. Full models are given in Supplementary Data 2.

Lower crustal assimilation
Assimilation of lower crustal rocks can reproduce the Chipiquixtle-YP trend in major and trace element, and isotope-isotope space (blue line in Figure 11). Moderate crystallisation (∼10-15 %) at an assimilation rate three times lower than that of fractional crystallisation (r = 0.30) can reproduce the Chipiquixtle-YP compositional range. Minor discrepancies between the models and the data likely reflect the uncertainty related to the composition of the lower crustal endmember. The lower crust beneath Mexico is composed of a variety of igneous and metasedimentary rocks [Ruiz et al. 1988;Roberts and Ruiz 1989;Schaaf et al. 1994;Straub et al. 2015], and considerable compositional variation is expected. Since no constraints for the composition of   Schaaf et al. [1994]; c from Vervoort et al. [2000]. r = ratio of assimilation to crystal fractionation. D = Bulk distribution coefficients calculated for given shallow and lower crustal phenocryst assemblages from Luhr and Carmichael [1980], Mahood   the lower crust beneath the PVC are available, the central Mexican granulite JC 101 from Schaaf et al. [1994] used here for modelling of Sr-Nd isotopic ratios, Mg and La/Nb ratios (Table 2) is considered the best estimate for the lower crust beneath the central TMVB.
The model produces a self-consistent explanation for major, trace and isotopic data at reasonable degrees of crystallisation. Magmas of the Chipiquixtle-YP trend are thus proposed to show lower crustal signatures, acquired in reservoirs of stalling magmas in the lower crust and/or during their ascent to the surface. Variable amounts of assimilation could be related to differences in the extent and timescale of exposure to the lower crust.

Upper crustal assimilation
The shallow crustal basement beneath the PVC is represented in the models by rhyolite and skarn xenoliths POP-97sst and POP-97sk. Modelling of 20-45 % crystallisation at minimal assimilation rates (r = 0.05) produces trends that largely coincide with the compositional variation displayed by the main PVC range (yellow and orange lines in Figure 11). This implies that assimilation of shallow basement may indeed play a role in the late-stage petrogenesis of PVC rocks, potentially modifying whole rock geochemical characteristics of PVC lavas and pumices. Granodioritic xenoliths show significant compositional overlap with the PVC andesites and dacites (stars in Figure 11), allowing for the assimilation of large quantities of granodioritic material without produc-

Presses universitaires de �rasbourg
Page 59 Interplinian effusive activity at Popocatépetl volcano Mangler et al., 2019 ing any distinct geochemical variation. In accordance with Schaaf et al. [2005], such granodiorites are interpreted as comagmatic rocks representing solidified domains of the igneous body beneath the PVC, e.g. cool margins of magma reservoirs or wall rock in the conduit. The volumetric predominance of such granodiorites over other types of xenoliths in PVC rocks indicates that melt reservoirs in the shallow plumbing system are predominantly hosted in solidified comagmatic rocks rather than in rocks of the crustal basement, thus limiting the scope for PVC magmas to extensively assimilate country rock. Substantial assimilation of local basement might only occur in association with mafic injections into the system, when heat, melt and gas addition pressurises and remobilises the magma reservoir sufficiently to fracture and incorporate wall rock [e.g. Beard et al. 2005;Glazner 2007;Andrews et al. 2008;Huber et al. 2011]. Uptake of shallow crustal xenoliths by magmas en route to the surface is considered an additional important factor. Both processes imply pre-eruptive disaggregation and bulk assimilation as the dominant form of interaction between magmas and local basement, suggesting only a minor role for assimilation by melting of xenoliths. Field observations of m-sized xenoliths in PVC rocks showing minimal evidence for digestion support this model (Appendix B).
The low r-value of 0.05 for the rhyolite and skarn assimilation models (Table 2; Figure 11) is consistent with the field evidence for limited interaction of PVC magmas with the surrounding country rocks. Assuming interaction with local basement is dominated by bulk assimilation rather than by melting, the modelled fractionation values shown in Figure 11 need to be considered semi-quantitative, since the AFC model of De-Paolo [1981] does not strictly account for bulk assimilation. Instead, the models presented here illustrate that even minor shallow assimilation of local country rocks can explain some of the geochemical characteristics of the main PVC range. In fact, higher assimilation rates of carbonate rocks would result in silicaundersaturated magmas with exceedingly high CaO concentrations [Mollo et al. 2010;Di Rocco et al. 2012], which are not observed at the PVC.
Finally, the local upper crust is composed of a variety of lithologies, which have undergone variable degrees of contact metamorphism with the magma reservoirs feeding the PVC [Schaaf et al. 2005]. This is evident from the wide range of PVC xenolith compositions sampled by Schaaf et al. [2005], Sosa-Ceballos et al. [2014], and during this study. Constraints on whether crustal assimilation is a significant process rely heavily on the specific xenolith used and may lead to conflicting results. The dataset and models presented here demonstrate that limited assimilation of skarn and rhyolite xenoliths can reproduce key geochemical trends observed in PVC rocks, and shallow crustal assimilation is thus considered a significant process since at least ∼23.5 ka.

Magma mixing
Assimiliation-fractional crystallisation of PVC magmas in the lower and upper crust can explain the compositional range of PVC lavas and pumices ( Figure 11). Geochemical variations over time (Figure 8) suggest that the PVC plumbing system hosts several magma reservoirs with distinct crustal signatures, which intermittently interact with each other. There is abundant petrological evidence for magma mixing and mingling at the PVC since at least ∼23.5 ka and up to the presentday activity [e.g. Straub and Martín-Del Pozzo 2001;Martín-Del Pozzo et al. 2003;Schaaf et al. 2005;Witter et al. 2005;Sosa-Ceballos et al. 2014], and field observations confirm that even the compositionally most primitive Chipiquixtle flow features mm-to m-thick bands of intermingled melts (Appendix B). This suggests that all PVC magmas were subject to pre-eruptive mixing and mingling processes, and the whole rock chemistry of PVC lavas and pumices can be explained by the mixing of magmas with distinct crustal signatures, producing melts of hybrid compositions (green field in Figure 11). Indeed, samples of the YP Plinian eruption and the subsequent Capula and Nealticán effusive eruptions trace such a hybridisation process from lower crustal signatures (YP) to main PVC range characteristics (Nealticán and PP; Figure 11). With these considerations in mind, the architecture of the PVC plumbing system and its dynamics over the past ∼23.5 ka are discussed below and represented schematically in Figure 12.

The shallow plumbing system
Melt inclusion and vesicle size distributions point toward a network of conduits and interfingering dykes at shallow depths beneath the PVC [Atlas et al. 2006;Roberge et al. 2009;Cross et al. 2012]. Based on seismic velocity models, Kuznetsov and Koulakov [2014] suggest that the modern Popocatépetl edifice is made up of solidified magmatic rocks perforated by fluid-and melt-filled fractures, cracks and pores down to a depth of at least 4 km below present sea level (i.e. ∼9.5 km depth beneath the summit crater). This is supported by abundant granodiorite xenoliths in PVC rocks [Sosa-Ceballos et al. 2012], which indicates the presence of a solidified comagmatic body hosting the active magma reservoir(s) [Schaaf et al. 2005]. It is proposed that this granodioritic body accumulated and increased in volume throughout the eruptive history of the PVC (dark grey in Figure 12). It could thus accommodate significant volumes of both crystal mush and eruptible magma in a dynamic, interconnected network, facilitating hybridisation between geochemically distinct magma batches (yellow/orange in Figure 12) and/or magma injections from deeper crustal levels (purple in Figure 12). The plumbing system beneath the PVC is thus envisioned as a network of transient melt reservoirs of varying compositions and crystal mush, hosted in a voluminous comagmatic body. This is consistent with emerging views of magma reservoirs as extensive, complex and highly dynamic systems of mush, melt and fluids rather than stagnant, stable magma 'chambers' [Christopher et al. 2015;Bachmann and Huber 2016;Druitt et al. 2016;Cashman et al. 2017;Magee et al. 2018;Pankhurst et al. 2018].
The upper crustal plumbing system described above is the locus of magma hybridisation at the PVC. Magma injections from lower crustal levels are suggested to feed the shallow plumbing system and mix with resident magma reservoirs, producing the hybrid melts of the 'main PVC range' that have dominated volcanism at the PVC since >23.5 ka (Figures 11 and 12A). Extensive magma mixing precludes detailed constraints on geochemical characteristics of the parental melts. However, the more primitive whole rock compositions of the Chipiquixtle and Atlimiyaya flank lavas suggest that magmas feeding these eruptions largely bypassed the shallow plumbing system and underwent limited hybridisation, as schematically shown for the Chipiquixtle eruption in Figure 12A.

The primitive endmember
The Chipiquixtle lava is the most primitive unit of the Popocatépetl edifice reported to date, and it may thus provide insights into the composition of parental melts of the PVC. Its whole rock composition is similar to some lavas of the neighbouring SCVF (Figure 13), and the general geochemical affinity of SCVF calc-alkaline rocks to the PVC at generally more primitive compositions (e.g. Figure 5) has prompted suggestions of a genetic relationship [Straub and Martín-Del Pozzo 2001;Schaaf et al. 2005]. More primitive compositions of the monogenetic SCVF magmas (i.e. basalts and basaltic andesites) may be related to shorter crustal residence times compared to PVC magmas. Indeed, prolonged crustal residence of SCVF magmas would facilitate crustal assimilation, which could produce melt compositions similar to PVC rocks ( Figure 13). AFC modelling using a calc-alkaline SCVF andesite as the magmatic endmember (CH-09-2: Straub et al. [2015]), and crustal assimilants and parameters identical to the models shown in Figure 11, can reproduce the vectors of both the Chipiquixtle-YP array and the main PVC range (Figure 13). This supports a genetic relationship between calc-alkaline SCVF and PVC rocks, and we suggest that primitive calc-alkaline rocks erupted in the SCVF approximate the parental melts feeding the PVC plumbing system. While such magmas pass through the crust rapidly and with little modification in the SCVF [e.g. Straub et al. 2013;Straub et al. 2015], they are significantly overprinted by AFC and hybridisation processes in the complex plumbing system be-neath the PVC (Figures 10-13).

An eruptive record dominated by magma hybridisation
The YP Plinian eruption and subsequent effusive eruptions illustrate the interplay of crustal assimilation at variable depths and magma hybridisation that dominates PVC whole rock compositions. Magmas stalling and fractionating at lower crustal levels are predicted to show evolved compositions coupled with lower crustal signatures (Figure 10). A magma with such characteristics was tapped during the YP Plinian eruption (Figure 8), and its lower crustal signatures continue to manifest during the subsequent Capula and Nealticán flank eruptions. However, the Nealticán lavas trace the waning influence of the lower crustal signature magma in a gradual hybridisation process with a 'main PVC range' magma, which eventually results in the eruption of main PVC range compositions during the PP Plinian eruption at ∼1100 a BP (Figures 8,10,11 and 12C). Hybridisation between deep-seated magmas and melts of the upper crustal plumbing system appears to largely buffer compositions of the main PVC range. However, subtle compositional differences are preserved within the shallow plumbing system, as shown by the two-stage PwA Plinian eruption (Figures 9  and 12B). Both phases erupted pumices with hybrid crustal signatures (green field in Figure 11) representative of the main PVC range (Figure 8), however they show systematic compositional differences (Figure 9). This indicates that both magmas involved in the eruption were shallow hybrids of slightly different compositions ( Figure 12B), which were mingled during the transition from the initial to the main stage of the PwA [Sosa-Ceballos et al. 2012]. Whole rock geochemical data further suggest that the PwA Plinian eruption caused a long-lasting change in the configuration of the shallow plumbing system: The shallow evolved reservoir primarily feeding pre-PwA lavas was largely exhausted during the PwA, and post-PwA effusive eruptions were dominated by the second reservoir tapped during the PwA Plinian eruption (Figures 8, 9 and 12). A significant impact on the plumbing system is consistent with the inferred magnitude of the PwA. With a volume of 1.8 km 3 DRE [Siebe et al. 2017], it is the largest eruption recorded at the PVC, which requires large portions of the plumbing system to be remobilised and tapped during the event. Furthermore, the high energy release resulted in significant fracturing and uptake of wall rock ( Figure 12B In summary, as shown graphically in Figure 11, the interplay of crustal assimilation at various depths on the one hand, and hybridisation between shallow   [Straub et al. 2015] for comparison. AFC model calculated using calc-alkaline SCVF sample CH-09-2 with same assimilants and parameters as in Figure 11A-C. See text for explanation. stalling reservoirs and fresh injections on the other hand, can reproduce the entire range of compositions displayed by PVC lavas and pumices. The long-lived upper crustal plumbing system beneath the PVC produced a voluminous environment of solidified comagmatic rocks hosting distinct, interconnected melt reservoirs. This facilitates extensive magma hybridisation within the shallow plumbing system, which buffers most compositions erupted by the PVC. Hybridisation textures found in present-day samples [e.g. Straub and Martín-Del Pozzo 2001;Schaaf et al. 2005;Witter et al. 2005;Roberge et al. 2009] and compositional variability exceeding the range displayed by pumices and lavas erupted in more than 23.5 ka (Figures 3-5), suggest similar dynamics for the present-day plumbing system of the PVC.

Conclusions
Research on active stratovolcanoes tends to be focussed on present-day activity and past cataclysmic eruptions, Presses universitaires de �rasbourg while past effusive and minor explosive eruptions are often neglected. This results in an incomplete understanding of volcanic systems and their dynamics, as Plinian eruptions merely provide isolated snapshots in time that lack potentially crucial information about their cause and effects on the plumbing system. Here, we present a comprehensive geochemical study of effusive and explosive eruptions of the last ∼23.5 ka at the PVC, to highlight mantle source characteristics, reconstruct the evolution and dynamics of the plumbing system, and show that the compositional range of the PVC rock suite is dominated by crustal processes such as assimilation-fractional crystallisation and magma mixing.
PVC magmas originate from a mantle source depleted by multiple melting events and modified by addition of subducted sediment and/or eroded forearc granodiorites. The PVC rock suite shows a restricted range for major and trace element, and Sr-Nd-Hf isotope compositions, with no long-term geochemical trends or systematic changes leading up to Plinian eruptions. Secondary geochemical variability of lavas and pumices reflects the evolution and dynamics of the plumbing system in the past ∼23.5 ka: The sector collapse of the Ventorrillo edifice and the associated VEI 5 TP Plinian eruption ∼23.5 ka BP had no significant effect on the plumbing system of the PVC, as indicated by unchanged magma compositions. During the VEI 6 PwA Plinian eruption ∼14.1 ka BP, the hitherto dominant magma reservoir was exhausted, and a new, geochemically distinct melt domain was tapped. After the PwA eruption, the plumbing system remained relatively steady until ∼2150 a BP, when the VEI 6 YP Plinian eruption tapped a magma reservoir with lower crustal signatures. The subsequent Capula and Nealticán effusive eruptions record a gradual hybridisation, which resulted in the reset to average PVC compositions for the VEI 4 PP Plinian eruption ∼1100 a BP. Throughout the past ∼23.5 ka BP, effusive flank eruptions produced more primitive magmas from deeper reservoirs, which experienced limited hybridisation by circumventing the shallow plumbing system. The compositional range of lavas and pumices erupted at the PVC is the result of crustal processes rather than mantle source heterogeneity. Shallow fractional crystallisation and limited assimilation of local upper crustal host rocks (e.g. skarn and older volcanics) can account for the main compositional range defining most PVC rocks. In addition, lower crustal assimilation can explain the distinct geochemical characteristics displayed by some pumices of the YP Plinian eruption ∼2150 a BP and by several effusive flank eruptions. Geochemical data and textural observations suggest that such geochemically distinct magmas undergo extensive mixing and mingling in the shallow plumbing system. All PVC magmas are subject to pre-eruptive hybridisation processes, albeit to different degrees. The relatively narrow compositional range of the majority of PVC rocks indicates a high degree of hybridisation. This illustrates the efficiency with which a volumetrically dominant, long-lived plumbing system such as beneath the PVC can buffer magma compositions through AFC processes combined with hybridisation between shallow stalling reservoirs and fresh injections. and HNO 3 added. The bombs were then sealed and placed in an oven at 160 • C for five days. This initial treatment replaced the first digestion step described for Sr-Nd. After oven digestion, subsequent treatment of the samples was as described for Sr-Nd digestions. Hafnium isolation was achieved in a two-stage column chemistry following the procedure of Bast et al. [2015]. Hafnium isotopic compositions were analysed on a Nu Instruments Nu Plasma HR multiple collector ICP-MS. Samples were introduced as 30 ppb Hf solutions in distilled 0.5M HNO 3 -0.3M HF and analysed in one block of 60 cycles. Wash time between each sample was 2 × 160 seconds to ensure complete Hf washout back to background levels between samples and standards. Mass fractionation was corrected online using the exponential law and 179 Hf/ 177 Hf = 0.7325, and potential mass interferences for Yb and W were monitored and corrected using natural abundances ratios 176 Yb/ 172 Yb = 0.58341, 174 Yb/ 172 Yb = 1.45990 and 180 W/ 182 W = 0.00450. Measured 176 Hf/ 177 Hf ratios were corrected for the accepted value for Hf standard JMC-475 of 176 Hf/ 177 Hf = 0.282160 [Nowell et al. 1998]. The instrumental precision within a single Hf sample run is ≤6 ppm 2 SE (Table 1). Two separate BCR-2 digestions were analysed, yielding an average of 176 Hf/ 177 Hf = 0.282866 ± 7 (2 SD; n = 5; see Supplementary Data 1: Table S1.4), in agreement with the recommended value of 176 Hf/ 177 Hf = 0.282870 ± 8 (2 SD; Weis et al. [2007]). B Appendix 2: Field evidence for magma mingling and xenoliths  Presses universitaires de �rasbourg