Anatomy of the magmatic plumbing system of Los Humeros Caldera (Mexico): implications for geothermal systems

Understanding the anatomy of magma plumbing systems of active volcanoes is essential not only for unraveling magma dynamics and eruptive behaviors but also to define the geometry, depth, and temperature of the heat sources for geothermal exploration. The Pleistocene–Holocene Los Humeros volcanic complex is part of the eastern TransMexican Volcanic Belt (central Mexico), and it constitutes one of the most important exploited geothermal fields in Mexico with ca. 90 MW of produced electricity. With the aim to decipher the anatomy (geometry and structure) of the magmatic plumbing system feeding the geothermal field at Los Humeros, we carried out a field-based petrological and thermobarometric study of the exposed Holocene lavas. Textural analysis, whole-rock major-element data, and mineral chemistry are integrated with a suite of mineral-liquid thermobarometric models. Our results support a scenario characterized by a heterogeneous multilayered system, comprising a deep (depth of ca. 30 km) basaltic reservoir feeding progressively shallower and smaller discrete magma stagnation layers and batches, up to shallow-crust conditions (depth of ca. 3 km). The evolution of melts in the feeding system is mainly controlled by differentiation processes through fractional crystallization (plagioclase+ clinopyroxene+ olivine+ spinel). We demonstrate the inadequacy of the existing conceptual models, where a single voluminous melt-controlled magma chamber (or “Standard Model”) at shallow depths was proposed for the magmatic plumbing system at Los Humeros. We instead propose a magmatic plumbing system made of multiple, more or less interconnected, magma transport and storage layers within the crust, feeding small (ephemeral) magma chambers at shallow-crustal conditions. This revised scenario provides a new configuration of the heat source feeding the geothermal reservoir at Los Humeros, and it should be taken into account to drive future exploration and exploitation strategies.

be very short, within the span of decades to a few thousands of years for tens to hundreds of cubic kilometers of eruptible magma (e.g., Glazner et al., 2004;Charlier et al., 2007), which are then rapidly evacuated during eruptions of calderaforming ignimbrites (e.g., Begué et al., 2014;Rivera et al., 2014;Wotzlaw et al., 2014;Matthews et al., 2015;Carrasco-Núñez et al., 2018). A key factor in determining the internal architecture of the magmatic systems is the magma intrusion rate. It controls whether successive pulses of magma will coalesce to form progressively larger chambers, as well as the formation of ductile shells surrounding the magma chamber that prevent country rock failure, favoring the inflation of the reservoir (Jellinek and DePaolo, 2003;Annen, 2009). Numerical simulations suggest that caldera systems smaller than 100 km 2 are fed by plumbing systems encapsulated by country rock that remains sufficiently brittle, while larger systems are more ductile, which favors an increase in size (Gregg et al., 2012).
The implications of such innovative conceptual models on the modeling of the heat source in magmatic-bearing geothermal systems are significant. Nonetheless, common numerical modeling of conductive-convective heat transfer in caldera-related geothermal systems has commonly envisaged the classic magma chamber as a single body, chemically stratified, entirely at magmatic temperatures, whose dimensions and depths have been usually constrained by volcanological and petrological data (e.g., Verma, 1985a, b;Wohletz et al., 1999). More complex modeling requires the "unpacking" of the stratigraphy of a volcano by the identification of the various "magma chambers" or magma storage layers that fed the different eruptions in space and time (e.g., Solano et al., 2014;Di Renzo et al., 2016;Cashman et al., 2017;Jackson et al., 2018).
In this paper we present a geothermobarometric study of the post-caldera Pleistocene-Holocene products of the Los Humeros volcanic complex (LHVC), located at the eastern termination of the Neogene-Quaternary Trans-Mexican Volcanic Belt (TMVB) (Fig. 1), with the goal of reconstructing the present-day geometry and structure of the magmatic plumbing system. These data are used to develop a conceptual model for the magmatic heat source of the active and currently exploited geothermal system. Since now, the magmatic heat source for LHVC has been constrained by the geometry of the caldera, the volume and mass balance calculations of the associated ignimbrites Mahood, 1984, 1987;Verma, 1984Verma, , 1985aVerma et al., 1990Verma et al., , 2011Verma and Gomez-Arias, 2013;Verma and Andaverde, 1995), all related to a single magma body. We propose a new and more realistic vision of the magmatic plumbing system, made of multiple interconnected magma stagnation layers within the crust. These new findings must be considered with the developing conceptual geothermal models to improve strategies for exploration and exploitation of the geothermal system within the LHVC. The results and approach presented in this work have also a general value and could represent an efficient strategy to explore and reconstruct the pre-eruptive geometry and the anatomy of active magmatic feeding systems.

Los Humeros volcanic complex
The volcanic evolution of the LHVC consists of three main stages (Carrasco-Núñez et al., 2018): (i) pre-caldera stage, (ii) caldera stage, and (iii) post-caldera stage. The pre-caldera stage is represented by relatively abundant rhyolitic domes, which erupted mainly on the western side of Los Humeros Caldera, with an isolated spot to the south, and in some buried lavas identified in the geothermal well logs (Carrasco-Núñez et al., 2017a). This volcanism has been recently dated by coupled U-Th and 40 Ar/ 39 Ar methods (Carrasco-Núñez et al., 2018), providing ages spanning from 693.0 ± 1.9 ka ( 40 Ar/ 39 Ar, plagioclase) to 270 ± 17 ka (U-Th, zircon), which overlap with the age range obtained from other domes of the western sector outside the caldera, where K-Ar ages (sanidine) of 360 ± 100 and 220 ± 40 ka were obtained (Ferriz and Mahood, 1984). The caldera stage consists of two major caldera-forming events, separated by a large Plinian eruptive episode. The first and largest caldera-forming eruption produced Los Humeros Caldera (18 km in diameter) during the emplacement of the Xaltipan ignimbrite, a rhyolitic, welded to non-welded ash-rich deposit, radially distributed around the caldera. The dense rock equivalent (DRE) volume of this event was estimated at 115 km 3 by Ferriz and Mahood (1984). The age of the Xaltipan ignimbrite was established by whole-rock K-Ar dating at 460 ± 20 ka (plagioclase) and 460 ± 130 ka (biotite) (Ferriz and Mahood, 1984); however, Carrasco-Núñez et al. (2018) based on coupled zircon U-Th dating and the 40 Ar/ 39 Ar method (plagioclase) geochronology provided a younger age of 164.0 ± 4.2 ka.
According to Carrasco-Núñez et al. (2018) during the post-caldera stage, two different eruptive phases occurred (Fig. 1). The first one was a late Pleistocene resurgent phase characterized by the emplacement of felsic domes in the central area at about 44.8 ± 1.7 ka (zircon U-Th ages; Carrasco-Núñez et al., 2018), which is slightly younger than the previously reported whole-rock K-Ar date (60 ± 20 ka, glass; Ferriz and Mahood, 1984). Outside of the caldera, to the north, a rhyolitic dome erupted at 50.7 ± 4.4 ka ( 40 Ar/ 39 Ar, plagioclase; Carrasco-Núñez et al., 2018). This was followed by a sequence of explosive eruptions, producing dacitic pumice fall units (Xoxoctic Tuff, 0.6 km 3 ) and interbedded breccia and pyroclastic flows deposits of the Llano Tuff (Ferriz and Mahood, 1984;Willcox, 2011), with a maximum age of 28.3 ± 1.1 ka (C-14 method; Rojas-Ortega, 2016). The second eruptive phase of the post-caldera stage is a Holocene ring-fracture episode and bimodal activity that occurred towards the south, north, and central parts of Los Humeros Caldera (Carrasco-Núñez et al., 2017a, b). It is character-ized by alternating episodes of effusive and explosive volcanism with a wide range of compositions. The volcanic products span from basalt, basaltic andesitic, trachyandesite, trachyte lava flows; and dacitic, trachydacitic, andesitic, and basaltic pumice and scoria fall deposits erupted by tens of monogenetic eruptive centers located in the LHVC (Ferriz and Mahood, 1984;Dávila-Harris and Carrasco-Núñez, 2014;Norini et al., 2015;Carrasco-Núñez et al., 2017b). Most of the effusive activity initially referred to 40-20 ka (K-Ar method, whole rock; Ferriz and Mahood, 1984), but recent dating reveals that most of this activity is Holocene in age (Carrasco-Núñez et al., 2017b). Trachyandesitic and andesitic basalt lavas erupted to the north of the LHVC at about 8.9 ± 0.03 ka (C-14 method;Carrasco-Núñez et al., 2017b). A rhythmic alternation of contemporaneous bimodal explosive activity produced trachyandesitic and basaltic fall layers grouped as the Cuicuiltic Member erupted at 7.3 ± 0.1 ka (C-14 method;Dávila-Harris and Carrasco-Núñez, 2014). This activity migrated towards the southern caldera rim to form a well-defined lava field. This ring-fracture episode erupted trachyandesite and olivine-bearing basaltic lava flows, at 3.9 ± 0.13 ka (C-14 method;Carrasco-Núñez et al., 2017b). The most recent activity erupted trachytic lava flows near the SW caldera rim, at 2.8 ± 0.03 ka (C-14 method; Carrasco-Núñez et al., 2017b).

Los Humeros geothermal system
The LHVC hosts one of the three most important geothermal fields in Mexico, with an installed 93 MW of electric power produced from 20 geothermal wells (Romo-Jones et al., 2017). The existing conceptual models for the Los Humeros geothermal field (LHGF) (see Norini et al., 2015, for a review) stem from the hypothesis of a unique, large, and voluminous cooling magma chamber of 1000-1500 km 3 in volume, residing at a depth of 5 to 10 km from the surface (Verma, 1984(Verma, , 1985a(Verma, , b, 2000Verma et al., 1990Verma et al., , 2011Verma and Gomez-Arias, 2013;Verma and Andaverde, 1995;Carrasco-Núñez et al., 2018) and providing the heat source for the geothermal system (Martínez et al., 1983;Verma, 1983Verma, , 2000Campos-Enríquez and Garduño-Monroy, 1987). However, the LHGF is characterized by a low number of productive geothermal wells (ca. 20 out of 50; Norini et al., 2015;Carrasco-Núñez et al., 2017a). The confined distribution of these productive wells along the NNW-SSE-trending "Maxtaloya-Los Humeros-Loma Blanca" fault system (MHBfs in Fig. 1) cutting across Los Potreros Caldera (e.g., Norini et al., 2015;Carrasco-Núñez et al., 2017a) also corresponds to the almost unique, narrow, and sharp surface thermal anomaly recognized within the caldera (Norini et al., 2015). These observations raise doubts about the existence of a voluminous superficial heat source feeding the LHGF. A revised assessment of the structure of the magmatic plumbing system beneath LHVC is therefore required for a better understanding and exploitation of the geothermal resource.

Materials and methods
In this work we focus on petrographic investigations including a textural and chemical (mineral chemistry and majorelement bulk rock) characterization of the Los Humeros postcaldera-stage (LHPCS) Holocene lavas (Carrasco-Núñez et al., 2017b). These data are used to define the thermobarometric conditions of the magma plumbing system and to test the eventual co-genetic nature of the LHPC melts. Following the recently published geological map (Carrasco-Núñez et al., 2017b) and geochronology (Carrasco-Núñez et al., 2018) of the LHVC, more than 50 samples of the LHPCS lavas were collected in the field with the goal to describe all the compositional variability of erupted products during the LH-PCS activity (Figs. 1,. In the following description of the volcanic units, abbreviations follow Carrasco-Núñez et al. (2017b). The most preserved and representative samples of each LHPCS volcanic unit were then selected (see Fig. 1 for sample location) for bulk and mineral chemistry investigations. With respect to the intracaldera domain (Fig. 2a), we selected lava samples belonging to (i) LH27-1 from the mafic lavas inside the Xalapasco Crater (Qb1), (ii) LH27-2 from the Maxtaloya trachyandesites (Qta4) constituting the rim walls of Xalapasco Crater, (iii) LH4 from San Antonio-Las Chapas lavas (Qta3) outcropping in the Los Humeros town, (iv) LH5-2 from mafic lavas (Qb1) outcropping west to Los Humeros town, (v) LH5-1 from Chicomiapa-Los Parajes felsic lavas (Qt2) outcropping in the northwestern part of Los Potreros Caldera, and (vi) LH6-1 from El Pajaro unit (Qt1) outcropping in the northwestern part of Los Potreros Caldera. In addition to these units, we also selected three more samples (LH13, LH26-1, and LH26-2) from lavas and domes of intermediate compositions, outcropping (Fig. 2bc) in the center of Los Potreros Caldera between Xalapasco Crater and Los Humeros town. These latter lavas and domes are not reported on the published geological map.
The samples were investigated first by optical microscopy and then through back-scattered electron (BSE) imaging for the definition of magmatic fabrics, textures, and constituent mineral assemblages. Mineral chemistry was then defined through electron microprobe analyses (EMPA). Whole-rock (major-element) composition of selected samples was obtained through inductively coupled plasma optical emission (ICP-OE) and X-ray fluorescence (XRF) analyses. Analytical protocols are described in Appendix A. In the following, mineral abbreviations follow Whitney and Evans (2010), whereas types of mineral zoning and textures are after Ginibre et al. (2002), Streck (2008), and Renjith (2014).
Intermediate products fall in the "Basaltic trachyandesite" and "Trachyandesite" fields; these rocks will be referred to hereafter as trachyandesites. The Los Humeros felsic (i.e., SiO 2 > 63 wt %) lava samples fall in the "Trachyte" field. Selected Harker diagrams for major elements are presented in Fig. 3, using SiO 2 wt % as differentiation index. Negative correlations are observed for CaO (Fig. 3c) and Mg# (Fig. 3d), whereas a positive correlation is observed for Na 2 O (Fig. 3c).

Basalts
LHPCS basalts show vesicle-rich (up to 35 vol %) highly porphyritic (phenocrysts up to 50 vol %) textures (Fig. 4ad). Studied basalts do not show presence of fragments from host-rocks or from previous magmatic rocks, therefore can be defined as lithic-free (e.g., Geshi and Oikawa, 2014). The magmatic fabric is fluidal as defined by the alignment of plagioclase laths in the groundmass. Based on the presence of clinopyroxene (Cpx) in the mineral assemblage, basalts can be further subdivided into (i) Cpx-free basalt of the extracaldera Texcal lava flow (LH18) and (ii) Cpx-bearing basalts of the intracaldera lavas at western Los Potreros and at Xalapasco Crater (LH5-2 and LH27-1, respectively).
Cpx-free basalt (LH18) contains euhedral to subhedral olivine (ca. 20 vol %) and euhedral to anhedral plagioclase (ca. 25 vol %-30 vol %) phenocrysts in a holocrystalline groundmass. The latter consists of plagioclase with swallowtail morphology, dendritic to spinifex olivines, and opaque oxides (Fig. 4a, b). Olivine and plagioclase phenocrysts are generally slightly chemically zoned (see below), showing homogeneous cores with normal concordant monotonous zoning texture at outer rim ( Fig. 4a, b). Major phenocrysts of olivine (up to 2.5 mm in size) with Cr-spinel inclusions are observed (Fig. 4a). No pyroxenes are observed in any samples collected from Texcal basalt.
Most of the phenocrysts observed in LHPCS trachyandesites show zoning textures characterized by homogeneous cores surrounded by (i) monotonous zoning at outer rims, (ii) low-amplitude euhedral oscillatory zoning at rims, and (iii) normal concordant step zoning at rims. Homogeneous unzoned clinopyroxene phenocrysts are commonly observed. Major phenocrysts of clinopyroxene (up to 2 mm in size) and plagioclase (up to 2 mm in length) characterized by homogeneous cores and normal concordant monotonous zoning are reported in all studied trachyandesites. We also report the presence of (i) rare clinopyroxene phenocrystals with growth mantle texture, (ii) rare plagioclase phenocrysts with rounded patchy cores, and (iii) very rare clinopyroxenes with homogeneous cores and growth mantle texture at rims. Very rare and large olivine phenocrysts (1.5-2.0 mm in size) presenting resorption patterns at rim and characterized by spinel inclusions are reported in the LH26-1 sample.  (Le Maitre et al., 2002). (c-d) Major-element-selected Harker diagrams for studied LHPCS lavas. The different symbols (circle for basalt, diamond for trachyandesite, and triangle for trachyte) represent the graphic code that will be used coherently throughout the article.

Trachytes
LHPCS trachytes show lithic-free phyric textures, with low porphyritic index (phenocrysts ranging 10 vol %-25 vol %), and fluidal fabrics as shown by iso-orientation of plagioclase, alkali feldspars, and clinopyroxene laths in the groundmass (Fig. 4g, h). They range from vesicle-poor (< 5 vol %) to vesicle-free textures, with the size of vesicles never exceeding 0.05 mm in diameter. The two analyzed trachytic samples (LH5-1 and LH6) are both characterized by the presence of orthopyroxene; however, the two mineral assemblages differ substantially.
The high-SiO 2 (67.58 wt %) LH6 trachyte is made up of sanidine (ca. 10 vol %), plagioclase (ca. 5 vol %-10 vol %), (e-f) BSE images of trachyandesites. It is possible to observe a microcrystalline groundmass where major phenocrysts of Cpx and Pl are dispersed. (g-h) BSE images of trachytes, characterized by a microcrystalline groundmass and Pl + Cpx + Opx phenocryst. Plagioclase phenocrysts show normal monotonous to normal step zoning. Major Cpx phenocrysts present inclusion of Ol + Mt. and orthopyroxene (ca. 10 vol %) phenocrysts in a finegrained trachytic mesostasis. Only major plagioclase and orthopyroxene phenocrysts show core-rim zoning textures with homogeneous cores associated either with normal monotonous zoning or normal low-amplitude oscillatory zoning textures at rims. Dimensions of phenocrysts are comparable to those of LH5-1 trachyte.

Mineral chemistry
Mineral compositions as obtained from electron microprobe analyses and mineral formulae for mineral assemblages of LHPCS lavas are presented in Tables S1, S2, S3, S4, and S5 in the Supplement (for feldspar, clinopyroxene, olivine, orthopyroxene, and spinel and opaque minerals, respectively).

Clinopyroxene
Apart from the LH18 basalt and LH6 trachyte, clinopyroxene is the most abundant mafic phase in analyzed LHPCS samples. It occurs generally as single crystals (Fig. 4e). However, rare crystals showing growth mantle texture are locally reported in trachyandesites. Very rare phenocrysts in trachyandesites show patchy cores. Major clinopyroxene phenocrysts in trachyandesites and trachytes contain inclusions (Fig. 4e, g) of olivine, magnetite, and plagioclase. Polarized light microscopy coupled with BSE images and chemical investigations highlighted the presence of unzoned ( Fig. 4g) and zoned (homogeneous cores surrounded by lowamplitude oscillatory zoning, normal monotonous zoning, or normal step zoning textures at rims; e.g., Fig. 4e) clinopyroxene phenocrysts. Very rare phenocrysts showing growth mantle texture at rim are reported. No evidence of resorption or dissolution textures are observed in analyzed LHPCS samples.
The Cpx population, based on textural observations and mineral chemistry, (Fig. 6a-f) can be classified into five major categories: (i) Cpx1 cluster is represented by homogeneous cores of all zoned phenocrysts in basalts; (ii) Cpx2 population is represented by homogeneous cores of all zoned phenocrysts from trachyandesites and trachytes; (iii) Cpx3 group represents both the unzoned phenocrysts in all studied samples and the rims (low-amplitude oscillatory, normal monotonous, and normal step zoning textures) of Cpx1 and Cpx2 phenocrysts from all studied samples; (iv) Cpx4 population is constituted by microlites and microphenocrystals in groundmass from all analyzed samples; and (v) Cpx5 cluster collects together the emerald-green euhedral to subhedral microlites in groundmass of intracaldera basalts (LH5-2, LH27-1) and the rare green outer rims of major clinopy-roxene (Cpx2) phenocrysts from few trachyandesites (LH15, LH17, LH26-2).
The compositional variation of clinopyroxenes can be summarized in the Na vs. Ti diagram ( Fig. 6e-f). Interestingly, augite-rich (Cpx1, Cpx2, Cpx3, and Cpx4) clinopyroxenes generally show positive correlation and linear distribution characterized by a progressive Ti and Na depletion, from Ti-augite cores (Cpx1) in basalts to DiHd-rich augite (Cpx3, Cpx4) specimens in trachytes. Cpx5, with aegirineaugite and aegirine compositions, diverges from this trend. It shows a negative correlation characterized by a progressive enrichment of Na content, with respect to a general Ti depletion. Aegirine enrichment could be diagnostic of ferric iron (Fe 3+ ) content increasing during the magmatic differentiation, whereas the diopside-hedenbergite enrichment could testify to an increase of ferrous iron (Fe 2+ ) in magma (e.g., Huraiová et al., 2017).

Mineral-liquid thermobarometry
In order to define the thermobaric (T -P ) environmental conditions of the magmatic feeding system of the LHPCS, we integrate thermobarometry models based on olivine (Beattie, 1993;Putirka et al., 2007;Putirka, 2008), orthopyroxene (Putirka, 2008), plagioclase (Putirka, 2005b(Putirka, , 2008, alkali feldspar (Putirka, 2008), and clinopyroxene (Putirka et al., 1996(Putirka et al., , 2003Putirka, 2008;Masotta et al., 2013) chemistry. Due to the paucity or absence of glass, we assume the whole-rock composition as representative of the original liquid (or nominal melt) in equilibrium with phenocrysts (Putirka, 1997(Putirka, , 2008Mordick and Glazner, 2006;Aulinas et al., 2010;Dahren et al., 2012;Barker et al., 2015). We are aware that such a procedure put the focus on early steps of the crystallization history, characterized by high melt to crystal ratios. Relatively late melt compositions, related to the solidification of the groundmass, are not present or can simply not be analyzed. Thermobarometric calculations were developed after the application of mineral-melt equilibrium filters and considering pre-eruptive H 2 O Liq values obtained through the plagioclase-liquid hygrometer model (Eq. 25b in Putirka, 2008). Plagioclase-liquid thermometry and barometry were calculated using Eqs. (24a) and (25a), respectively, of Putirka (2008), mainly based on the Ca/Na distribution between melt and Pl. Alkali-feldspar-liquid thermometry was calculated considering the K-Na exchange, applying Eq. (24b) in Putirka (2008). Olivine-liquid equilibrium thermometry was calculated by integrating the models of Beattie (1993) and Herzberg and O'Hara (2002) with the thermometric Eq. (2) in Putirka et al. (2007). Orthopyroxene-liquid thermometry was calculated by Fe-Mg partitioning following the model of Beattie (1993; in the revised form Eq. 28a in Putirka, 2008). The barometry model of Wood (1974) based on the Na and Al content in Opx, in the revised form Eq. (29a) in Putirka (2008), was applied.

Test for mineral-melt equilibrium
A prerequisite for the application of mineral-liquid thermobarometry models based on mineral-melt equilibrium conditions is to test and verify that the mineral and the chosen liquid compositions represent chemical equilibrium pairs (e.g., Putirka, 2008;Keiding and Sigmarsson, 2012). Petrographic investigations (i.e., polarized light and BSE imaging) and calculation of mineral-liquid partition coefficients were integrated with the aim to select only mineral specimens at equilibrium with the hosting melt (e.g., Putirka, 2008;Keiding and Sigmarsson, 2012).
The predominant euhedral to subhedral habit of crystals is generally considered to be evidence of equilibrium with the surrounding melt (e.g., Keiding and Sigmarsson, 2012). However, mineral-liquid thermobarometric modeling requires paying careful attention to phenocrysts showing strongly zoned textures (patchy, sector, reverse, coarse banding oscillatory zoning), or disequilibrium textures (resorption patterns, dissolution surfaces, reaction rims, and mineral mantles or clots) (e.g., Ginibre et al., 2002;Streck, 2008). These textures imply that core(s) and rim(s), or different portions of the same grain, crystallized and reacted in an evolving liquid with a progressively different composition (e.g., Mordick and Glazner, 2006;Putirka, 2008;Keiding and Sigmarsson, 2012). As defined by Streck (2008), when crystals are complexly zoned, it can be difficult to find criteria to be used for evaluation of crystal populations and their equilibrium with respective hosting melt. However, it is not the case in the studied LHPCS samples, where phenocryst assemblages generally do not show disequilibrium patterns or complexly zoned textures (e.g., Ginibre et al., 2002;Streck, 2008). All microprobe analyses related to those rare crystals presenting morphological evidence of disequilibrium texture, such as patchy zoning, were discarded.
Then, the mineral-liquid equilibria between liquid and previously selected minerals were investigated using (i) the Fe-Mg exchange coefficient, (ii) the An-Ab partitioning coefficient, and (iii) the comparison between observed and predicted normative components of minerals.
The partitioning of An-Ab between mineral and liquid is known as An-Ab exchange coefficient, or K , where Liq is the liquid composition, Pl is the plagioclase composition, and all components are in molar fractions; Carmichael et al., 1977;Holland and Powell, 1992;Putirka et al., 2007;Putirka, 2008;Lange et al., 2009;Keiding and Sigmarsson, 2012;Jeffery et al., 2013;Barker et al., 2015;Waters and Lange, 2015). Figure 10 presents a comparison of measured composition of plagioclase with that calculated from the melt composition, using the thermodynamic model Eq. (31) in Namur et al. (2012). A similar test can be applied for alkali feldspars (Putirka, 2008).

Pre-eruptive H 2 O Liq content estimates
Thermobarometric models for volcanic systems require an initial estimate of the pre-eruptive water concentration (wt %) in melt (H 2 O Liq ). This was determined in this work by using the plagioclase-liquid hygrometer model Eq. (25b) in Putirka (2008). Hygrometry calculations were produced after the application of plagioclase-liquid equilibrium filters. The calculated pre-eruptive H 2 O Liq wt % values (±1σ standard deviation of the weighted mean) are plotted as isolines in Fig. 10. The hygrometer of Putirka (2008) Webster et al., 1999), the anhydrous character is then assumed as a H 2 O Liq < 1 wt % content.
Application of the plagioclase-liquid hygrometer model (Putirka, 2008) defines an anhydrous environment for pressure-temperature calculations in LHPCS basalts. Hydrous conditions are required for evolved LHPCS melts and in particular for trachytic lavas, where the effect of 1 wt % H 2 O is expected to generate a temperature decrease of ca. −40 • C and a pressure increase of ca. +1.0 kbar in geothermometers and geobarometers, respectively (Putirka, 2008;Keiding and Sigmarsson, 2012).

Trachyandesites
Based on the Opx presence or absence criterion, two populations of trachyandesites have been discriminated in this study.

Major-element mass balance modeling
Based on the documenting of textural evidence -(i) Cpxbearing basalts are mainly characterized by euhedral olivine and plagioclase and subhedral-anhedral clinopyroxene, indicating crystallization of olivine and plagioclase prior to clinopyroxene (e.g., Bindeman and Bailey, 1999) and (ii) all LHPCS volcanic rocks do not show disequilibrium textures (such as fine-sieve textures, resorption surface, crystal clots, disequilibrium growth-mantel, reverse zoning, reaction rims, breakdown mantle, and dissolution; e.g., Streck, 2008) typical of assimilation and fractional crystallization (AFC) mixing processes -we suggest that the studied LHPCS volcanic rocks represent cogenetic melts, belonging to the same line of descent, excluding major mass change due to assimilation and mixing (AFC-mixing) processes. In order to test this hypothesis, we applied fractional crystallization (FC) modeling (e.g., White et al., 2009;Moghadam et al., 2016;. The FC modeling is focused on the following hypotheses: (i) a direct cogenetic relationship between all LH-PCS basalts and (ii) a common genesis for all LHPCS trachyandesites and trachytes through differentiation via fractional crystallization starting from the same basaltic parental melt.
If Parent melt (for RFC) or Cumulate (for crystal accumulation) compositions are assumed as matrix b, and the FC model is solved for b, then b = Liquid (Daughter or Melt) + Minerals (fractionating or accumulated assemblage). If compositions of Liquid and Minerals are known (matrix A), it is possible to estimate, by least squares approximation, their proportion (in matrix c). The similarity of b (matrix c multiplied with matrix A) to b (real value) is quantified with the sum of the square of the residuals ( r 2 ) as (1) RFC and Cumulate model results are considered acceptable when r 2 < 1.0. Proportion of Liquid (Daughter or Melt) is expressed with F in matrix c.
Major-element mass balance models are calculated in the system SiO 2 -TiO 2 -Al 2 O 3 -FeO * -MnO-MgO-CaO-Na 2 O-K 2 O. The LH5-2 Cpx-bearing basalt, with the lowest SiO 2 and the highest MgO contents, was selected as a possible source for all pyroxene-bearing trachyandesites and trachytes. The fractional crystallization hypothesis is then tested for all the studied LHPCS rocks and considering the magmatic mineralogy made of An-rich plagioclase, Ti-rich clinopyroxene, Mg-rich olivine, and spinel. The same mineral assemblage was used then to verify the cogenetic relationship between studied LHPCS basalts through progressive crystal accumulation. All calculations were managed with Microsoft Office Excel 2019. Results of FC models are presented in Table S6 in the Supplement.
The results obtained from FC models thus indicate that the LHPCS volcanic rocks are genetically linked melts, due to crystal accumulation (basalts) and fractional crystallization (intermediate and felsic rocks) of a Pl + Cpx + Ol + Spl mineral assemblage. Trachyandesites and trachytes represent different degrees of fractionation (RFC values in the range 45 %-74 %) starting from a Cpx-bearing basaltic source. Cpx-bearing basalts are interpreted as the result of crystallization and accumulation of Cpx, together with Pl + Ol + Spl, in a pristine Cpx-free basaltic melt. Results from FC models also confirm the possibility to produce hydrous felsic melts starting from a nominal anhydrous (H 2 O < 1 wt %; e.g., Webster et al., 1999) mafic parental melt. Integrating FC model and hygrometer (Putirka, 2008) results, LHPCS trachytes show H 2 O ca. 1.4 wt %-2.0 wt % and represent the ca. 25 wt % fractionated residual melt from a parental basaltic source characterized by H 2 O in the range 0.3 wt %-0.5 wt %.

Magma evolution beneath Los Humeros
The conceptual model of the present-day LHPCS magmatic plumbing system beneath the Los Humeros Caldera is presented in Fig. 12. Based on textural observations, mineral chemistry, and thermobaric estimates the early hightemperature (1230-1270 • C) stage of LHPCS magma evolution is represented by high-anorthite plagioclase phenocrysts and Mg-rich olivine (X Fo = 80 %-85 %) crystallizing in the deep (ca. 8 kbar) basaltic reservoir. Where these magmas erupted directly, they formed Cpx-free Ol-basalt lava flows such as the Texcal lava flow (LH18). This scenario, for LH18 basalt sample, is confirmed by (i) olivine and plagioclase with homogeneous cores and normal monotonous zoning textures at rims, indicating a fast growth during ascent of magma (e.g., Streck, 2008); (ii) olivine with spinifex, dendritic, and skeletal textures, interpreted as supercooling mineral texture largely resulting from rapid olivinesupersaturated magma rise from deeper level during the eruption (e.g., Donaldson, 1974;Nakagawa et al., 1998;Fowler et al., 2002;Dahren et al., 2012;Welsch et al., 2013); and (iii) plagioclase specimens with swallow-tailed crystal morphology, interpreted as rapid plagioclase growth due to undercooling related to the eruption process (e.g., Renjith, 2014).
On the other hand, a permanence of these basaltic melts in the deep reservoir, together with a temperature decrease of ca. 100 • C, can lead to clinopyroxene appearance and/or crystallization in the system (e.g., Grove, 2000) and its progressive accumulation in the phenocryst assemblage. This hypothesis is supported by Cpx-Liq thermometry models for Cpx1 (Ti-rich augites in basalts), indicating Cpx appearance at ca. 7-8 kbar and 1150 • C (mean values) and by FC models indicating a Pl + Cpx + Ol + Spl crystal accumulation up to 15 wt %-17 wt % in the pristine basaltic melt to produce the Cpx-bearing basalts.
The shallowest magma stagnation level (< 3 kbar; mean 1.5 kbar) has been here interpreted as a complex magma Figure 12. Schematic representation (not to scale) of the magmatic plumbing system feeding LHPCS activity, beneath Los Humeros Caldera, as derived by pressure-temperature estimates obtained from mineral-liquid thermobarometry models. The conceptual model is integrated with the crustal structure of the study area as derived by the measure No. 10 of the Crust 1.0 global Model (Davies, 2013). The gray shaded field indicates the depth and thickness of the existing conceptual model of a single, huge, classical magma chamber proposed by Verma (1985a, b) and mainly related to the Los Humeros caldera-stage activity. plexus constituted by a system of small magma volumes, distributed in locally interconnected pockets and batches. In this plexus, mafic and intermediate magmas shortly stall prior to erupt. More evolved melts reside for a relatively longer time, enough to crystallize orthopyroxene and to enable the escape of part of the exsolved volatiles (e.g., Sparks et al., 1998;Feng and Zhu, 2018;Clarke et al., 2007), as suggested by phenocryst textures and compositions and by poor-vesicle textures observed in Opx-trachyte samples (LH5-1, LH6-1).
Compositional reverse zoning, associated with disequilibrium textures and dissolution or resorption patterns in phenocrysts, are widely considered indicators of both magma replenishment or assimilation processes (e.g., Wright and Fiske, 1971;Duda and Schminkcke, 1985;Clague et al., 1995;Yang et al., 1999;Klügel et al., 2000;Zhu and Ogasawara, 2004;Stroncik et al., 2009;Ubide et al., 2014;Viccaro et al., 2015;Gernon et al., 2016;Feng and Zhu, 2018). In the case of LHVC, almost all investigated LH-PCS samples, from basalts to trachytes, contain mainly phenocrysts with homogeneous cores and low-amplitude oscillatory, normal monotonous zoned rims (Pl + Ol + Cpx), or unzoned homogeneous phenocrysts (as in the case of Cpx3 and Opx). Rare specimens not suitable for mineral-liquid thermobarometry, such as plagioclase and clinopyroxene with patchy cores or olivine xenocrysts, are reported. The general absence of disequilibrium textures and dissolution patterns in studied LHPCS samples, is therefore interpreted as a lack of evidence of major mixing/recharge and/or assimi-146 F. Lucci et al.: Anatomy of the magmatic plumbing system of Los Humeros Caldera lation processes acting in the plumbing system (e.g., Cashman et al., 2017, and references therein). This hypothesis is in line with the results obtained from tests for mineralmelt equilibria. The Rhodes diagram (Rhodes et al., 1979;Putirka, 2008) for olivine compositions (Fig. 8a) highlights a progressive decrease in Mg# Liq from basalts to trachytes coupled with a general absence of xenocrystal and/or antecrystal cargo. This behavior is compatible with a complete removal from the melt of previously crystallized Mg-olivine (Roeder and Emslie, 1970;Dungan et al., 1978;Rhodes et al., 1979;Putirka, 2008;Melluso et al., 2014). All LHPCS melts (from basalts to trachytes) invariably show suites of olivines with maximum forsterite (Fo) contents in equilibrium with the respective whole rocks, and vertical trends consistent with closed-system melt differentiation (Roeder and Emslie, 1970;Rhodes et al., 1979;Putirka, 2008;Melluso et al., 2014). Similar behavior is obtained for orthopyroxene (Fig. 8b), where again the Rhodes test highlights (i) absence of antecrystals and (ii) Opx suites progressively and normally Fe enriched from trachyandesites to trachytes. The absence of clinopyroxene clots and overgrowth mantle textures on orthopyroxene crystals again excludes the occurrence of magma mixing and/or recharge processes (Laumonier et al., 2014;Neave et al., 2014;Zhang et al., 2015;Feng and Zhu, 2018). Such interpretation is supported also by field observations, where the interbedded basaltic andesite and trachydacite fall deposits of the ca. 7 ka Cuicuiltic Member show no evidence of magma mixing (Dávila-Harris and Carrasco-Núñez, 2014).
An-Ab partition coefficients (e.g., Putirka, 2008;Jeffery et al., 2013) show a comparable scenario (Fig. 10) in which (i) the LHPCS basalts are characterized by suites of plagioclases with maximum anorthite (An) contents in equilibrium with the respective whole rocks, and a progressive An Pl decrease consistent with closed-system differentiation, and (ii) the progressive decrease in predicted An Liq from basalt to trachyte is compatible with evolved melt differentiation via fractional crystallization. The LHPCS intermediate and evolved products show plagioclase phenocrysts characterized by An-rich homogeneous cores (An 70 %-85 %), with compositions comparable to those of basalts. These Anrich cores can be crystallized in two possible scenarios. The first one is related to the H 2 O content in magma. Increasing the water content in melt strongly favors crystallization of An-richer plagioclase. A water content rise from 0.5 wt % to 2.0 wt % could lead to an increase in the An component up to 6-8 mol % (Bindeman and Bailey, 1999;Sano and Yamashita, 2004;Ushioda et al., 2014). In this view, the Anrich plagioclase in intermediate and felsic rocks can be interpreted as the response to the increasing water content in the fractionated melt. The second scenario implies that Anrich plagioclase taps a more primitive stage of basalt segregation. Since plagioclase phenocrysts with An in the range 65 %-81 % are commonly found in LHPCS basalts, the Anrich plagioclase cores in trachyandesites and trachytes could represent either antecrystals derived from crystallization of early sills in the magmatic reservoir system (sensu Jackson et al., 2018) or crystallization products in an earlier stage of the trachyandesite and trachyte segregation from the basaltic reservoir (e.g., Bindeman and Bailey, 1999;Kinman and Neal, 2006). We suggest that both scenarios concurred for the genesis of An-rich phenocrysts in trachyandesites and trachytes. Note that when An-rich plagioclase crystals are found (in mafic and intermediate rocks with Pl + Ol + Cpx assemblages), it implies that no significant clinopyroxene crystallization has occurred prior to the anorthitic plagioclase (Bindeman and Bailey, 1999).
With respect to plagioclase, a similar behavior is observed also for clinopyroxene and in particular for Cpx1 and Cpx2 (clinopyroxene cores in basalts and in trachyandesites + trachytes, respectively) populations. Since these mineral cores (Pl, Cpx1, and Cpx2) generally present normal growth rims (i.e., Ab-rich Pl and Cpx3), we suggest that stagnation levels at both intermediate and shallower depths underwent crystallization in a closed system. Otherwise, features such as (i) diffused reverse zoning textures; (ii) hightemperature crystal clots, mantling, and overgrowth textures; (iii) disequilibrium and dissolution textures (e.g., Stroncik et al., 2009;Cashman et al., 2017;Feng and Zhu, 2018, and references therein) should be widely observed, but this is not the case in the studied LHPCS lavas.

The magma plumbing system
The petrological archive constituted by the LHPCS lavas, spanning from transitional and alkali basalts to trachytes, describes the Holocene activity of the LHVC. Harker diagrams for major-element bulk compositions of the LHPCS lavas are characterized by linear trends (Fig. 3b-d) comparable to those expected for cogenetic melts (e.g., Giordano et al., 2012). Major-element FC modeling confirms the hypothesis of a common genesis for the LHPCS volcanic rocks through crystal fractionation or accumulation processes of the same mineral assemblage (Pl + Cpx + Ol + Spl). Furthermore, textural observations and results from FC models permit us to exclude mass-change or mass-addition processes driven by AFC-mixing processes.
Results obtained from the application of different and independent thermobarometry models (Fig. 11) confirm the working hypothesis of a complex magmatic plumbing system rather than a single (i.e., Standard Model) magma chamber (e.g., Keiding and Sigmarsson, 2012;Cashman and Giordano, 2014;Cashman et al., 2017;Feng and Zhu, 2018) developed beneath the active Los Humeros Caldera and feeding the LHPCS volcanism.
With the aim to propose an updated and realistic conceptual model of the present-day main storage zones and magma plumbing system within the crust below Los Humeros Caldera, we integrate pressure-temperature estimates acquired in this study with the existing data related to the crustal structure and corresponding physical parameters of the study area. The resulting model is shown in Fig. 12.
The density of the TMVB crust shows a large range between 1800 kg m −3 for unconsolidated sediments via about 3000 kg m −3 for the lower crust to 3300 kg m −3 for the upper mantle (Dziewonski and Anderson, 1981;Campos-Enríquez and Sánchez-Zamora, 2000;Davies, 2013). The available compilation of crustal data for LHVC is recovered by the measure No. 10 of the Crust 1.0 global model (Dziewonski and Anderson, 1981;Davies, 2013). The measure No. 10 (yellow star in Fig. 1) is located within the study area at the southern termination of the Tepeyahualco lava flow and describes a crust made of five main seismic layers (Fig. 12): (i) upper sediments (thickness: 1 km, density 2110 kg m −3 ), (ii) middle sediments (thickness: 0.5 km, density 2370 kg m −3 ), (iii) upper crust (thickness: 13.6 km, density 2740 kg m −3 ), (iv) middle crust (thickness: 15.3 km, density 2830 kg m −3 ), and (v) lower crust (thickness: 13.6 km, density 2920 kg m −3 ). Inferred (seismic) Moho depth is reported at −41.7 km with an upper mantle density of 3310 kg m −3 (Dziewonski and Anderson, 1981;Davies, 2013). Here we use a five-tiered density model as derived from the Crust 1.0 global model to convert obtained pressure estimates to crustal depths below LHVC.
The thermobarometry models applied to the LHPCS lavas define a broad region of crystallization between ca. 0 and 30 km in depth that can be described with a quadrimodal distribution of pressure values (Fig. 12). This allows us to propose a complex polybaric continuous heterogenous multilayered transport and storage magmatic system.
A second magma transport and storage system can be recognized at depths of 15-30 km (ca. 4.5-7.8 kbar), in continuity with the deeper basaltic reservoir and distributed along the whole middle crust thickness, as recorded by the wide range of pressure estimates obtained from plagioclase (X An = 40 %-70 %) and Cpx2 clinopyroxene cores (Mg# 59-84; TiO 2 mean value 0.99 wt %). Thermometry models based on plagioclase, Cpx2 clinopyroxene, and olivine show convergence for hydrous temperature values in the range of 979-1263 • C. Thermobarometry models, together with textures and petrographic relations in all analyzed trachyandesite and trachyte samples, suggest that all plagioclase, all Cpx2 clinopyroxene phenocrysts, and part of microlites grew in this second storage system. In particular, it is possible to observe two main crystallization temperature conditions: (i) at ca. 1190 • C (weighted mean value, MSWD = 2.2, n = 205) plagioclase phenocryst crystallization in trachyandesite melts is observed, whereas (ii) at the lower temperature of ca. 1070 • C (weighted mean value, MSWD = 1.7, n = 155) the crystallization of all olivine, all Cpx2 phenocrysts, and plagioclase phenocrysts in trachytes is reported. We interpret the common Pl + Cpx2 phenocryst-forming barometric conditions as the evidence of a growth-dominated regime within this second magma storage zone (e.g., Barclay et al., 1998;Humphreys et al., 2006;Scott et al., 2012). The smaller crystals (microcrystals and microlites) represent the nucleationdominated regime (Scott et al., 2012) that can be associated with ascent-related decompression of melts at shallower levels (e.g., Cashman, 1992;Cashman and Blundy, 2000;Humphreys et al., 2009).
The third melt storage zone occurs at shallower depths of ca. 10-15 km, possibly corresponding to the transition between middle and upper crusts, as indicated by convergence of barometric estimates (weighted mean value of 3.9 ± 0.2 kbar (±1σ ), MSWD = 0.80, n = 203; P ranging ca. 1-7 kbar) obtained from Cpx3 clinopyroxene (i.e., unzoned phenocrysts and overgrowth or rims around earlierformed Cpx1 and Cpx2 cores) populations. For this third storage zone, the Cpx3-Liq thermometry model indicates a mean temperature of 1040 • C (weighted mean value, MSWD = 2.6, n = 203; T ranging ca. 940-1210 • C), comparable to those calculated for Ol + Cpx2 assemblages in the previously described second and deeper stagnation system. The obtained pressure estimates for the second and the third storage systems are compatible with multiple magma storage pockets in which melts of comparable compositions ascend slowly enough for phenocrysts to form (e.g., Scott et al., 2012) and start cooling before the final ascent to shallower conditions (e.g., Dahren et al., 2012;Chadwick et al., 2013;Gardner et al., 2013;Jeffery et al., 2013;Preece et al., 2013;Troll et al., 2013). Taking into account the textures and the chemistry of Cpx3 clinopyroxene phenocrysts, the obtained thermobarometric estimates could be interpreted as the pressure-temperature environment of last major levels of magma stagnation and fractionation (Putirka, 1997;Klügel et al., 2005;Galipp et al., 2006;Stroncik et al., 2009).
On the other hand, the Holocene eruptive phase of the LHPCS is characterized by a bimodal volcanism (Carrasco-Núñez et al., 2017a, b, 2018, typified by alternating episodes of effusive and explosive activity with a wide compositional range of volcanic products, spanning from basaltic to trachytic lava flows and mafic to felsic pumice and scoria fall deposits, erupted by tens of monogenetic centers located in the LHVC (e.g., Norini et al., 2015;Carrasco-Núñez et al., 2017a, b, 2018. The LHPCS volcanic activity is characterized by spatially distributed small volumes of erupted material (ca. 6 km 3 of mafic lavas, 10 km 3 of intermediate and felsic lava, and 1 km 3 of mafic and felsic tephra; Carrasco-Núñez and Branney, 2005). Furthermore, key elements, such as the lithic-free character of the LHPCS volcanic products, their overall textures, and chemistry of the constituent mineral assemblages, coupled with the results from RFC models, suggest that LHPCS magmatism is characterized by batches of magma evolving in a nearly closed system, unaffected by magmatic assimilation and mixing and/or recharge processes. In particular, the almost complete lack of magma mixing and/or recharge events (e.g., Lee et al., 2014) is confirmed by the absence of the typical expected mineral textures (e.g., Streck, 2008;Renjith, 2014) such as (i) fine-sieve textures and resorption surfaces due to reaction with a more primitive magma; (ii) glomerocryst forming due to the recrystallization or suturing at rim of resorbed crystals; (iii) reverse zoning textures due to compositional inversion in an open or recharged system; and (iv) reaction rims, breakdown mantles, and crystal clots due to the disequilibrium-triggered recrystallization into a new set of minerals.
The existing literature, focused on magma recharge processes (e.g., De Paolo, 1981;Hofmann, 2012;O'Neill and Jenner, 2012;Lee et al., 2014), highlights that a high evacuation and/or eruption efficiency would shorten the residence time of magma in the storage chamber, and it would reduce the effect of crystallization on modifying the magma composition (Lee et al., 2014). Moreover, in the case of eruption and/or evacuation rates higher than the recharge rates (e.g., Lee et al., 2014), it is possible to hypothesize a magmatic system dominated by ephemeral closed-system magma batches not affected by major mixing processes prior to their evacuation and/or eruption (e.g., De Paolo, 1981;Hofmann, 2012;O'Neill and Jenner, 2012;Lee et al., 2014). This scenario best approximates the characteristics observed for all the Holocene LHPCS magmatic products. In addition, the lack of liquid-dominated zone(s) (e.g., Bachmann and Bergantz, 2008), where mixing could occur (e.g., Cashman and Giordano, 2014), suggests that the remnants of the huge magma chamber of the LH caldera stage are now completely solidified and crosscut by the uprising LHPCS mafic and felsic magmas. This scenario is also coherent with the postcaldera eruption behavior observed in other volcanic complexes, such Ischia (e.g., Casalini et al., 2017), and it is consistent with the recent literature proposing complex magma chamber reservoirs made up of multiple discrete melt pockets with no mass exchange and reactivated shortly before eruption (e.g., Cashman and Giordano, 2014;Cashman et al., 2017;Casalini et al., 2017).
Thermobarometric estimates obtained in this study, combined with the existing literature and integrated with information from the crustal structure beneath Los Humeros Caldera, therefore permit us to discard the Standard Model of the huge voluminous chamber in favor of a more feasible conceptual model characterized by a polybaric magmatic plumbing system made up of multiple, more or less interconnected, magma transport and storage layers, i.e., transient batches and ponds of different magmas, localized beneath Los Humeros nested caldera and feeding the Holocene activity of the LHVC. In particular, our results indicate that magma transport and storage levels beneath Los Humeros Caldera are vertically distributed across the whole crust from ca. 30 to 3 km (from the lower to the very upper crust) with density contrasts between the different crustal layers acting as a controlling parameter for ascending or stalling magmas (e.g., Dahren et al., 2012), reflecting the buoyant magma compositions and the melt fractions (e.g., Cashman et al., 2017;Jackson et al., 2018). Moreover, it is possible to propose that each of these crust or density boundaries have determined lateral transport and grow magma stagnation pockets (e.g., Dahren et al., 2012;Jackson et al., 2018). At depths < 5 km, buoyant magmas and fractionated melts (from mafic to felsic) ascending from all the lower storage zones are stalled once more. The shallowest complex multistorage system is interpreted as a plexus of scattered, more or less interconnected, ephemeral small-volume batches and pockets of melts, without any defined spatial distribution, as confirmed by field locations of the studied LHPCS lavas eruptive centers.
A shallow storage zone presenting magmas with heterogenous compositions (from mafic to felsic) has been already proposed by Dávila-Harris and Carrasco-Núñez (2014) to explain the eruptive history of the intracaldera Cuicuiltic Mem-ber that was produced by the coeval eruption of mafic and felsic unmixed magmas. However, a shallow ponding system characterized by heterogeneous composition of magmas involved beneath Los Humeros Caldera is not an exceptional case. Examples of shallow heterogeneous reservoirs beneath active volcanic complexes are widely reported (e.g., Kratzmann et al., 2009;Sigmarsson et al., 2011;Keiding and Sigmarsson, 2012).
Our results also agree with the work of Créon et al. (2018), where calculated fluid saturation depths derived for melt inclusions in post-caldera lavas indicate different magmaponding levels within a range of depths between 5 and 13 km, together with a possible deeper reservoir (26-32 km) and a final shallow stagnation level at ca. 1.5-3.0 km.

Implications for the active geothermal system
The geothermal activity of a volcanic complex is expected to be the result of stagnation and cooling of magmas in the shallower storage zone (e.g., Gunnarsson and Aradóttir, 2015), where classic conductive models are adopted to model the heat source, mainly controlled by age and volume of the magmatic system (Smith and Shaw, 1975;Cathles et al., 1997;Duffield and Sass, 2003;Gunnarsson and Aradóttir, 2015;Carrasco-Núñez et al., 2018). As widely demonstrated (e.g., Smith and Shaw, 1975;Cathles et al., 1997), a very large intrusion would produce a long-lived hydrothermal and geothermal system. Many numerical models (e.g., Cathles et al., 1997) suggest that, in the most favorable conditions, a voluminous (> 2000 km 3 ) intrusion or chamber of mafic melt would be able to sustain a convective geothermal system up to 800 kyr. On the other hand, very small mafic sill and dike intrusions (< 10 km 3 ) would produce very localized thermal anomalies and could cool down to the solidus temperature in less than 0.1 kyr (Nabelek et al., 2012) and definitively cool in ca. 1 kyr (e.g., Cathles et al., 1997). Convection, due to hydrothermal fluid circulation, increases the cooling rate of a magmatic intrusion (Cathles et al., 1997).
The present geothermal activity of LHVC is characterized by a limited NNW-SSE nonhomogeneous areal distribution within the Los Potreros nested caldera (e.g., Norini et al., 2015;Urbani et al., 2019). Based on (i) the young age (Upper Pleistocene-Holocene) of most of the LHPCS volcanic activity; (ii) the relatively small erupted volumes of the LHPCS lavas, in particular of those erupted within the Los Potreros Caldera; and (iii) the existence of a shallower magmatic plexus characterized by heterogeneous unmixed magmas (this study), we therefore discard the hypothesis of a single, large, and voluminous shallow magmatic chamber homogenously distributed beneath the caldera in favor of a more feasible scenario characterized by an upper crustal plexus made of small single-charge ephemeral pockets of different magmas localized beneath Los Humeros nested caldera, very close or within the Los Humeros exploited geothermal field. In this scenario, every LHPCS magma pocket and cryptodome within the Los Humeros Caldera (see Urbani et al., 2019) could be interpreted as a scattered and localized short-lived (ca. 0.1-1 kyr; Cathles et al., 1997) heat source, whereas the cooling and solidified remnants of the huge magma chamber of the caldera stage could still represent a background positive thermal anomaly affecting the volcanic field.
Our reconstruction of the Los Humeros heat source therefore suggests the possible existence of a wide backgroundpositive thermal anomaly associated with the cooling, solidified remnants of the voluminous magma chamber of the caldera stage, in juxtaposition to scattered high-frequency heat sources related to the very shallow intrusive complex that makes up the surficial (upper crustal) plexus of the LH-PCS magmatic plumbing system.
In the light of our results, a revision or update to the heat source feeding the Los Humeros geothermal system is needed to produce correct and up-to-date geothermal potential estimates of the geothermal field and to develop efficient geothermal exploration and exploitation strategies.

Conclusions
In this study we propose an integrated field-based petrographic and mineralogical approach to unravel the evolution and configuration of the present-day magmatic plumbing system feeding the post-caldera-stage activity of LHVC. The main results of this study can be summarized as follows: i. The Rayleigh fractional crystallization (RFC) models demonstrate that all LHPCS magmas, from basalts to trachytes, belong to the same line of descent and evolve through a progressive fractionation of the Pl + Cpx + Ol + Spl mineral assemblage.
ii. A complex polybaric magmatic transport and storage system, characterized by multiple magma levels more or less interconnected in space and time, has been recognized based on application of mineral-melt thermobarometry models.
iii. A deep mafic reservoir (at ca. 30 km depth) is identified by the Pl + Ol assemblage in basalts. Intermediate magma storage systems (in the whole middle crust) are described by the composition of the Cpx phenocrysts, whereas a shallow magmatic stagnation system (ca. 1.5 kbar; 3-5 km depth) is defined by crystallization of Cpx microlites (aegirine clinopyroxenes in basalt) and, in particular, by Opx growth in most evolved melts. All the Cpx-bearing lavas are produced by progressive differentiation via polybaric fractional crystallization during magma ascent through the plumbing system. iv. The chemical composition of the main phases (Ol, Pl, Cpx, Opx), together with results from FC modeling, does not support a magmatic feeding system dominated by magma mixing and magma replenishment. They are instead compatible with a plumbing system dominated by discrete levels, pockets, and batches of melts.
v. The thermobarometric results indicate that the configuration of the magmatic plumbing system is vertically extensive across the entire crust, with a deeper residence zone for basalts at ca. 8 kbar (ca. 30-33 km depth). A complex zone, from middle-(6-4 kbar) to upper-crust (0.5 kbar) depths, where rapidly ascending basalts stall before their eruption, is proposed. This complex zone also corresponds to depths where smaller batches of mafic magmas, at times interconnected with the lower feeding zone, differentiate to trachyandesites and trachytes.
vi. The main outcome for the modeling of the magmatic heat source of the LHVC geothermal system is the inadequacy of conservative conceptual models based on the classical melt-dominated, single, long-lived, and voluminous magma chamber (i.e., Standard Model) in favor of an innovative and more realistic vision of the magmatic plumbing system made of multiple, more or less interconnected, magma transport and storage layers within the crust, feeding small (ephemeral) magma pockets at shallow-crust conditions.
vii. The proposed model for the magmatic plumbing system at LHVC provides a new configuration of the heat source feeding the present geothermal reservoir that must be taken into account for geothermal exploration and exploitation purposes.  Whitney and Evans (2010).

A2 Bulk major-element geochemistry
After washing in distilled water, samples were grounded in an agate mill and pre-contaminated with an aliquot of sample. Whole-rock major-element concentrations (four samples) were measured at the Activation Laboratories (Ontario, Canada), through inductively coupled plasma (ICP) optical emission (OE). For major elements, the uncertainty (1σ ) is estimated better than 2 % for values higher than 5 wt % and better than 5 % in the range 0.1 wt %-5 wt %. Additional samples (nine samples) were analyzed by X-ray fluorescence (XRF) using a ZSX Primus II (Rigaku Co., Japan) at Nagoya University, Japan. Loss on ignition (LOI) was measured from the sample powder weight in a quartz glass beaker in the oven at 950 • C for 5 h. XRF-analyses were carried out following the procedure presented in Azizi et al. (2015Azizi et al. ( , 2018a. For major elements, the uncertainty (1σ ) is estimated better than 1 % for values higher than 10 wt % and better than 5 % in the range 0.1 wt %-10 wt %.

A3 Mineral chemistry
Polished thin sections (13 samples) selected for petrography investigations, were then studied for mineral chemistry, and ca. 2400 analyses of mineral phases were obtained with a Cameca SX100 electron microprobe (EMP) at the Institut für Anorganische Chemie, Universität Stuttgart, Germany. Operating conditions were 15 kV and 10 to 15 nA, counting times of 20 s both for peak and background. Spot sizes were 1-10 µm depending on the phases analyzed. Compositions were determined relative to natural and synthetic standards.
A set of reference materials (i.e., natural and synthetic oxides and minerals) was used for routine calibration and instrument stability monitoring. In particular, we used (i) Si, Ca: natural wollastonite (P&H Developments); (ii) Si, Fe: natural fayalite USNM 85276 (Jarosewich et al., 1980); (iii) K: natural orthoclase (P&H Developments); (iv) Na: natural pure albite from Crete (Greece); (v) Al: synthetic corundum (P&H Developments); (vi) Mg: synthetic periclase (P&H Developments); (vii) Mn: natural rhodonite (P&H Developments); (viii) Ti: synthetic rutile (P&H Developments); (ix) Cr: synthetic chromium oxide (P&H De-velopments). Repeated analyses of the standards (Table S7) resulted in one-sigma (1σ ) standard deviations close to the ones calculated from counting statistics. For the major minerals, calculated 1σ (%) precisions are (i) better than 1.5 % for Si; (ii) better than 2 % for Al; and (iii) 1 % to 5 % for Ca, Mg, Fe, Mn, Ti, and Cr, applying the abovementioned applied conditions. For Na and K, calculated 1σ (%) precisions are below 5 % for analyses of feldspars and Aeg-rich clinopyroxene. The 1σ accuracy is estimated to be up to 3 times larger than the precision because additional effects such as uncertainty of the mass absorption coefficients that are used for the matrix correction of the microprobe raw data or instability of the beam may play a role.
Data availability. Data presented in this paper are all available in the Supplement.
Author contributions. FL, FR, GCN, and GG conceived the initial idea of this study. FL, FR, GG, GCN, and SU participated in the fieldwork. FL conducted petrographic investigations and performed the thermobarometry modelling. FL and TT performed and processed BSE and EMP (electron microprobe) analyses. HA and YA performed and processed whole-rock geochemistry analyses. FL and JCW performed fractional crystallization (FC) and massbalance modelling. FL prepared the paper with contributions from all co-authors.
Competing interests. The authors declare that they have no conflict of interest.
Acknowledgements. The authors are grateful to the editor (Johan Lissenberg), to Chiara M. Petrone, and to an anonymous reviewer for their helpful and constructive comments that deeply contributed to improving the article. The authors wish to thank the Comisión Federal de Electricidad (CFE, Mexico) for their assistance and support. This paper presents results of the GEMex project, funded by the European Union's Horizon 2020 program for Research and Innovation under grant agreement no. 727550 (scientific responsibility Guido Giordano), and by the Mexican Energy Sustainability Fund CONACYT-SENER, project 2015-04-268074 (WP 4.5, scientific responsibility Gerardo Carrasco-Núñez). More information can be found on the GEMex website: http://www. gemex-h2020.eu (last access: 14 January 2020). Authors would like to thank Gianluca Norini for useful discussions in the field. Special thanks to Javier Hernández, Jaime Cavazos, Francisco Fernández, and Alessandra Pensa for their support in the fieldwork and logistics. The grant to the Department of Science, Roma Tre University (MIUR-Italy Dipartimenti di Eccellenza, ARTICOLO 1, COMMI 314-337 LEGGE 232/2016) is gratefully acknowledged. Gerardo Carrasco-Núñez is grateful to the PASPA-DGAPA program (UNAM) for support during his sabbatical stay at the University of Roma Tre (Rome, Italy).
Financial support. This research has been supported by the GEMex project, funded by the European Union's Horizon 2020 program for Research and Innovation (grant no. 727550) and the Mexican Energy Sustainability Fund CONACYT-SENER (project 2015-04-268074).
Review statement. This paper was edited by Johan Lissenberg and reviewed by Chiara Maria Petrone and one anonymous referee.