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 Trans-Mexican 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.

1 Introduction

Recent views on the structure of volcanic plumbing systems have moved from the “Standard Model” (sensu Gualda and Ghiorso, 2013) of a single bowl-shaped magma chamber where all petrological processes of differentiation and assimilation occur (e.g., Hildreth, 1979, 1981; Hildreth and Wilson, 2007) to more complex arrays of stratified and variably interconnected transient magma accumulation zones, set in largely crystallized and vertically extensive mush zones (e.g., Bachman and Bergantz, 2004, 2008; Cashman and Giordano, 2014; Cashman et al., 2017). Furthermore, the time required for the assembly of large magma chambers is now believed to 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 caldera-forming 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² 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).

A key to deciphering where magmas are stored, and therefore the anatomy of a magmatic plumbing system, is the understanding of pre-eruptive processes such as mineral crystallization and the migration and stagnation of melts prior to their eruption (Feng and Zhu, 2018; Putirka, 2008; Keiding and Sigmarsson, 2012; Scott et al., 2012; Barker et al., 2015; Jeffery et al., 2013; Cashman and Giordano, 2014; Pamukcu et al., 2015; Lucci et al., 2018). Early, segregated minerals reflect the magmatic environment (i.e., pressure–temperature, magma and fluid compositions, oxidation state), and thus their growth, texture and chemistry provide an important archive of information (Ginibre et al., 2002; Feng and Zhu, 2018; Ginibre et al., 2007; Streck, 2008; Giuffrida and Viccaro, 2017; Viccaro et al., 2016; Putirka, 2008; Lucci et al., 2018). Accordingly, petrographic observations and mineral chemistry of primary minerals, integrated with opportunely selected thermobarometry models (e.g., Putirka, 2008; Masotta et al., 2013), could lead to the comprehension and reconstruction of the magmatic storage and/or feeding systems of the erupted products (Feng and Zhu, 2018; Giuffrida and Viccaro, 2017; Elardo and Shearer, 2014; Petrone et al., 2016; Zheng et al., 2016; Eskandari et al., 2018; Shane and Coote, 2018; Scott et al., 2012; Stroncik et al., 2009; Barker et al., 2015; Jeffery et al., 2013; Keiding and Sigmarsson, 2012).

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 (Ferriz and Mahood, 1984, 1987; Verma, 1984, 1985a, b; Verma et al., 1990, 2011; Verma 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.

Figure 1 Geological context. (a) The Los Humeros volcanic complex (LHVC, blue dot) with respect to the Trans-Mexican Volcanic Belt (TMVB). (b) Shaded relief image obtained from 15 m resolution digital elevation model (DEM) of the LHVC. Volcanic products of the Los Humeros post-caldera stage are redrawn from Carrasco-Núñez et al. (2017b). The description of the volcanic units, their names, and abbreviations follow Carrasco-Núñez et al. (2017b). The map shows locations (white dots) and volcanological significance of the samples used in this study. The yellow star indicates the locality of measure No. 10 of the Crust 1.0 global model (Dziewonski and Anderson, 1981; Davies, 2013).

2 Geological setting

2.1 Regional geology

LHVC is the largest and easternmost Quaternary caldera (Fig. 1) of the 1200 km long and active Trans-Mexican Volcanic Belt (TMVB), generated since ca. 20 Ma by the subduction of the Cocos plate beneath central Mexico (e.g., Demant, 1978; Ferrari et al., 1999, 2012; Gómez-Tuena et al., 2003, 2007a, b, 2018; Norini et al., 2015). The eastern sector of the TMVB, where LHVC is located, is characterized by monogenetic volcanism, scattered basaltic cinder, and scoria cones, maar volcanoes of basaltic and rhyolitic composition, large rhyolitic domes and major stratovolcanoes such as Pico de Orizaba (or Citlaltépetl) and Cofre de Perote (e.g., Yáñez and García, 1982; Negendak et al., 1985; Carrasco-Núñez et al., 2010, 2012a).

The Paleozoic to Mesozoic crystalline basement of eastern TMVB is exposed along the Teziutlán Massif (Viniegra-Osorio, 1965; Yáñez and García, 1982; Ferriz and Mahood, 1984). The crystalline basement is partially covered by a thick, highly deformed Mesozoic sedimentary succession, part of the Sierra Madre Oriental fold and thrust, structured during the Late Cretaceous–Paleocene Laramide orogeny (e.g., Campos-Enríquez and Garduño-Monroy, 1987; Suter, 1987; Fitz-Díaz et al., 2018). Oligocene to Miocene granodiorite and syenite intrusions (whole-rock K–Ar ages spanning 15–31 Ma; Yáñez and García, 1982) are randomly exposed in the region. Miocene volcanism consists of andesites of the Cerro Grande volcanic complex (Gómez-Tuena and Carrasco-Núñez, 2000), with an age of ca. 9–11 Ma (K–Ar method, whole rock; Carrasco-Núñez et al., 1997) and the Cuyoaco andesite dated at ca. 11 Ma (K–Ar method, whole rock; Yáñez and García, 1982) to the west of LHVC, which can be correlated with the Alseseca andesite (Yáñez and García, 1982) exposed to the north. Pliocene andesitic volcanism (Ferriz and Mahood, 1984; Yáñez and García, 1982) is represented by the Teziutlán andesite, which, recently dated by the ⁴⁰Ar∕³⁹Ar method at 2.61±0.43 to 1.46±0.31 Ma (Carrasco-Núñez et al., 2017a), correlates with most of the thick andesitic successions of the subsurface geology of LHVC.

2.2 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 ⁴⁰Ar∕³⁹Ar methods (Carrasco-Núñez et al., 2018), providing ages spanning from 693.0±1.9 ka (⁴⁰Ar∕³⁹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³ 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 ⁴⁰Ar∕³⁹Ar method (plagioclase) geochronology provided a younger age of 164.0±4.2 ka.

Following this catastrophic event an eruptive pause occurred, resuming with a sequence of intermittent Plinian episodes at 70±23 ka (⁴⁰Ar∕³⁹Ar method on plagioclase; Carrasco-Núñez et al., 2018), separated by short gaps marked by thin paleosoils. The deposits consist of thick (1–6 m), coarse, pumice-rich well-sorted, massive, and diffuse-stratified layers, rhyodacitic in composition, which are grouped as the Faby Tuff (Ferriz and Mahood 1984; Willcox, 2011). The second caldera-forming episode produced the 9–10 km large Los Potreros Caldera, which is associated with the emplacement of the compositionally zoned andesitic-rhyodacitic-rhyolitic Zaragoza ignimbrite (Carrasco-Núñez et al., 2012b). This is an intraplinian pyroclastic flow deposit, with an estimated volume of ca. 15 km³ DRE (Carrasco-Núñez and Branney, 2005). Previous ages of this unit were reported at 100 ka (K–Ar dating, plagioclase; Ferriz and Mahood, 1984) and at 140±24 ka (⁴⁰Ar∕³⁹Ar method, plagioclase; Willcox, 2011). However, a new younger ⁴⁰Ar∕³⁹Ar plagioclase age of 69±16 ka for the Zaragoza ignimbrite was recently obtained, supported by the field evidence that the Zaragoza ignimbrite overlies a rhyodacitic lava flow dated at 74.2±4.5 ka (zircon U–Th dating) (Carrasco-Núñez et al., 2018).

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 (⁴⁰Ar∕³⁹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³) 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 characterized 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).

2.3 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³ in volume, residing at a depth of 5 to 10 km from the surface (Verma, 1984, 1985a, b, 2000; Verma et al., 1990, 2011; Verma 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, 1983, 2000; Campos-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.

3 Materials and methods

In this work we focus on petrographic investigations including a textural and chemical (mineral chemistry and major-element bulk rock) characterization of the Los Humeros post-caldera-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 LHPCS activity (Figs. 1, 2a–d). 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. 2b–c) 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.

Figure 2 Field photographs of LHPCS volcanic products. (a) E–W panoramic view of Xalapasco Crater; the dashed white line indicates the limit of Cpx-bearing Ol-basalt lavas filling the crater. (b) Intracaldera trachyandesitic lavas outcropping at Los Potreros, south to Los Humeros town. (c) Trachyandesitic lava-dome outcropping inside Los Potreros Caldera, north to Xalapasco Crater. (d) E–W panoramic view from the SE Los Humeros Caldera rim. Dashed lines indicate the two major trachyandesitic lava flows of “El Limón” and “Sarabia”. Pico de Orizaba, Las Derrumbadas, and Cerro Pizarro volcanoes are also indicated. White dots indicate sampling localities.

Concerning the extra-caldera products (Fig. 2d), we selected one sample for each of the four major lava flows: (i) LH15 from El Limón lava flow (Qab), (ii) LH21-2 from Sarabia lava flow (Qta1), (iii) LH17 from Tepeyahualco lava flow (Qtab), and (iv) LH18 from Texcal lava flow (Qb1).

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).

4 Major-element bulk composition

Studied samples show a continuous series from mafic to felsic compositions, with SiO₂ ranging 46.5 wt %–67.6 wt %, and Na₂O+K₂O ranging 3.4 wt %–9.2 wt % (with ${\mathrm{K}}_{\mathrm{2}}\mathrm{O}/{\mathrm{Na}}_{\mathrm{2}}\mathrm{O}<\mathrm{1}$) (Fig. 3a; Table 1). LHPCS mafic rocks (SiO₂<50 wt %; 3 samples) show composition with SiO₂ 46.5 wt %–49.4 wt %, Al₂O₃ 16.2 wt %–17.1 wt %, CaO 9.8–10.7 wt %, MgO 8.0 wt %–8.4 wt % with Mg# (molar $\mathrm{MgO}/[\mathrm{MgO}+{\mathrm{FeO}}_{\mathrm{tot}}]$) equal to 60–61, and Na₂O+K₂O ranging 3.4 wt %–3.5 wt %. LHPCS intermediate rocks ($\mathrm{50}<{\mathrm{SiO}}_{\mathrm{2}}<\mathrm{63}$ wt %; 8 samples) contain 54.4 wt %–62.1 wt % SiO₂, with Al₂O₃ 15.7 wt %–20.7 wt %, Na₂O+K₂O 5.3 wt %–7.1 wt %, MgO 2.2 wt %–3.6 wt % (Mg# 43–51), and low CaO 4.6 wt %–8.5 wt %. LHPCS felsic rocks (SiO₂>63 wt %; 2 samples) show SiO₂ ranging 64.9 wt %–67.6 wt %, associated with Al₂O₃ 15.5 wt %–15.8 wt %, MgO 0.7 wt %–1.2 wt % (Mg#: 26–34), CaO 2.1 wt %–2.8 wt %, and Na₂O+K₂O 8.2 wt %–9.2 wt %.

Table 1 Major-element bulk-rock compositions of studied LHPCS lava samples.

Note: LOI – loss on ignition; Mg# – molar [$\mathrm{Mg}\cdot \mathrm{100}/(\mathrm{Mg}+{\mathrm{Fe}}^{\mathrm{tot}})$].

Figure 3 (a) Satellite image of the LHVC (Image Landsat from Google Earth Pro, 2018, Digital Globe; courtesy of Google) with localization of samples selected for the application of Rayleigh Fractional Crystallization model and for thermobarometry models. (b) Total alkali versus silica (TAS) diagram (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.

On the total alkali versus silica (TAS) diagram (Le Maitre et al., 2002) LHPCS lavas span from basalt to trachyte (Fig. 3b). Los Humeros mafic rocks fall in the “Basalt” field and, following the existing literature (e.g., Barberi et al., 1975; Bellieni et al., 1983; Le Maitre et al., 2002; White et al., 2009; Giordano et al., 2012), can be classified as mildly alkaline (or transitional) basalts and alkali basalts. The high-TiO₂ contents (1.34 wt %–1.5 wt %), together with MgO<12 wt % and low Al₂O₃∕TiO₂ values (average value 11.5), exclude the LH mafic rocks as high-Mg melts (picrites) or komatiites (e.g., Redman and Keays, 1985; Arndt and Jenner, 1986; Le Maitre et al., 2002; Gao and Zhou, 2013; Azizi et al., 2018a, b).

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₂>63 wt %) lava samples fall in the “Trachyte” field. Selected Harker diagrams for major elements are presented in Fig. 3, using SiO₂ wt % as differentiation index. Negative correlations are observed for CaO (Fig. 3c) and Mg# (Fig. 3d), whereas a positive correlation is observed for Na₂O (Fig. 3c).

5 Petrography

5.1 Basalts

LHPCS basalts show vesicle-rich (up to 35 vol %) highly porphyritic (phenocrysts up to 50 vol %) textures (Fig. 4a–d). 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 extra-caldera 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).

Figure 4 Microtextures and fabrics of the LHPCS lavas. (a) Back-scattered electron (BSE) image of the Cpx-free Ol-Basalt fabric, dominated by euhedral unzoned homogeneous Pl + Ol, with major olivine phenocryst characterized by Cr-Spl inclusions. (b) BSE image of Cpx-free Ol-Basalt groundmass highlighting the spinifex to skeletal and dendritic crystallization of olivine, associated with the swallow-tailed morphology of plagioclase. (c) BSE image of Cpx-bearing Ol-Basalt. Normal monotonous zoning at rim is observed for all the main mineral phases (Pl + Ol + Cpx). (d) BSE image of Cpx-bearing Ol-Basalt groundmass characterized by albitic plagioclase, aegirine-pyroxene, Fe-rich olivine, and ilmenite. (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.

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 swallow-tail 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.

Cpx-bearing intracaldera basalts (LH5-2, LH27-1) show euhedral to subhedral plagioclase (ca. 25 vol %), euhedral olivine (ca. 10 vol %–15 vol %), subhedral to anhedral yellow-to-colorless clinopyroxene (ca. 10 vol %–15 vol %), and rare subhedral anorthoclase (<2 vol %) phenocrysts (Fig. 4c) in a holocrystalline groundmass (Fig. 4d). The latter is made up of (in order of abundance) elongated platy plagioclase, olivine, colorless-to-green clinopyroxene, opaque oxides, and rare alkali feldspar. All phenocrysts show core-rim zoning textures (Fig. 4c, d): (i) olivine, plagioclase, and clinopyroxene with homogeneous cores and normal concordant monotonous zoning texture at outer rims; (ii) plagioclase and clinopyroxene with homogeneous cores and low-amplitude euhedral oscillatory zoning texture at rims; (iii) rare plagioclase and clinopyroxene with homogeneous cores and normal concordant step zoning texture at rims; and (iv) very rare plagioclase with patchy cores and normal convolute monotonous zoning texture at rims. Large phenocrysts of olivine (up to 1.5 mm in size) and plagioclase (up to 3 mm in length) are commonly observed. Vesicle size is up to 5 mm in diameter.

5.2 Trachyandesites

LHPCS intermediate volcanic products are lithic free and show low- to medium-porphyritic textures (phenocrysts ranging 10 vol %–40 vol %), with a general fluidal fabric as indicated by orientation of plagioclase and clinopyroxene laths in the groundmass (Fig. 4e, f). Intermediate products vary from poorly vesicular (<10 vol % in LH13) to vesicle-rich (ca. 30 vol % in LH4) lavas. In the highest vesiculated sample (LH4), the size of vesicles (3–5 mm in diameter) is comparable to those of intracaldera basalts (LH5-2, LH27-1). Based on the presence of orthopyroxene (Opx) in the mineral assemblage, trachyandesites can be further subdivided into (i) Opx-free (LH21, LH15) and (ii) Opx-bearing (LH4, LH13, LH17, LH26-1, LH26-2, LH27-2) trachyandesites.

Opx-free trachyandesites show euhedral to subhedral plagioclase (ca. 15 vol %–20 vol %), euhedral to subhedral olivine (ca. 10 vol %), euhedral to anhedral yellow-to-colorless clinopyroxene (ca. 10 vol %–20 vol %), subhedral alkali feldspars (ca. 10 vol %–15 vol %) phenocrysts, in a holocrystalline microcrystalline groundmass composed of elongated platy plagioclase, colorless clinopyroxene, olivine, alkali feldspar, and opaque oxides in order of microlite abundance.

Opx-bearing trachyandesites are generally characterized by euhedral to subhedral plagioclase (ca. 15 vol %–20 vol %), euhedral to subhedral clinopyroxene (ca. 10 vol %–20 vol %), euhedral colorless orthopyroxene (ca. 10 vol %–20 vol %), euhedral to subhedral alkali feldspars (ca. 10 vol %–15 vol %), and euhedral to subhedral olivine (<10 vol %) phenocrysts in a holocrystalline to hypohyaline microcrystalline groundmass made of feldspar (plagioclase and alkali feldspars) microlites, pyroxene (clinopyroxene and orthopyroxene) microlites, olivine microlites, opaque minerals, and glass (Fig. 4f). No olivines are observed in samples collected from the Maxtaloya trachyandesite.

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.

5.3 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 low-SiO₂ (64.93 wt %) LH5-1 trachyte is characterized by euhedral to subhedral phenocrysts of plagioclase (ca. 10 vol %–15 vol %), clinopyroxene (ca. 10 vol %), orthopyroxene (ca. 10 vol %), olivine (ca. 5 vol %–10 vol %), and sanidine (<10 vol %) in a hypohyaline microcrystalline groundmass made of (in order of abundance) sanidine, orthopyroxene, clinopyroxene, rare plagioclase, rare olivine, rare opaque minerals, and very rare glass. All phenocrysts are generally unzoned. Mafic phenocrysts with homogeneous cores surrounded by rims showing normal concordant monotonous zoning textures are also observed. Rare major plagioclase (up to 1.5 mm in length) phenocrysts present patchy cores and normal concordant step zoning texture at rims. Rare major clinopyroxene (up to 1.0 mm in length) and unzoned homogeneous phenocrysts show inclusions of olivine + magnetite.

The high-SiO₂ (67.58 wt %) LH6 trachyte is made up of sanidine (ca. 10 vol %), plagioclase (ca. 5 vol %–10 vol %), and orthopyroxene (ca. 10 vol %) phenocrysts in a fine-grained 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.

6 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).

6.1 Feldspar

In basaltic rocks (Fig. 5a), feldspars are predominantly plagioclase. Plagioclase phenocrysts show anorthitic (X_An = 59 %–81 %, average 67 %) cores and normally zoned (X_An = 42 %–59 %, average 53 %) rims. Orthoclase component (X_Or) is always less than 2 %. Plagioclase microlites in the groundmass show andesine (X_An = 19 %–60 %) composition, with X_Or ranging 1 %–7 %. Alkali feldspars occur as both rare anorthoclase phenocrysts (Ab₆₀Or₃₇) and microlites in groundmass (${\mathrm{Ab}}_{\mathrm{62}-\mathrm{79}}{\mathrm{Or}}_{\mathrm{9}-\mathrm{35}}{\mathrm{An}}_{\mathrm{0}-\mathrm{13}}$).

Figure 5 An-Ab-Or diagrams showing the composition of feldspar in (a) basalts (circles), (b) trachyandesites (diamonds), and (c) trachytes (triangles) of LHPCS lavas.

Plagioclase from trachyandesites (Fig. 5b) have anorthite-rich (X_An = 45 %–87 %, average 67 %) cores and normally zoned (X_An = 27 %–69 %, average 48 %) rims. Cores with X_An in the range 72 %–87 % are observed in all major phenocrysts. Plagioclase core compositions are comparable to those of basalts. The X_Or ranges 1 %–8 %. Plagioclase microlites in groundmass show andesine (X_An = 29 %–63 %) composition with X_Or always less than 10 %. Alkali feldspars occur as (i) anorthoclase (${\mathrm{Ab}}_{\mathrm{59}-\mathrm{68}}{\mathrm{Or}}_{\mathrm{11}-\mathrm{30}}$) and sanidine (${\mathrm{Ab}}_{\mathrm{49}-\mathrm{50}}{\mathrm{Or}}_{\mathrm{43}-\mathrm{48}}$) phenocrysts and (ii) anorthoclase (${\mathrm{Ab}}_{\mathrm{49}-\mathrm{70}}{\mathrm{Or}}_{\mathrm{15}-\mathrm{38}}$) and sanidine (${\mathrm{Ab}}_{\mathrm{38}-\mathrm{48}}{\mathrm{Or}}_{\mathrm{47}-\mathrm{61}}$) microlites in groundmass.

Trachytes (Fig. 5c) show generally unzoned plagioclase phenocrysts with oligoclase-andesine (X_An = 26 %–45 %) composition. Rare An-rich (X_An = 52 %–70 %) cores are reported from major plagioclase phenocrysts in the low-silica trachyte LH5-1. The X_Or is always less than 8 %. Plagioclase microlites in groundmass are rare, with an Ab-rich (${\mathrm{An}}_{\mathrm{21}-\mathrm{30}}{\mathrm{Ab}}_{\mathrm{66}-\mathrm{69}}{\mathrm{Or}}_{\mathrm{4}-\mathrm{10}}$) composition. Alkali feldspar is represented by anorthoclase as phenocrysts (${\mathrm{Ab}}_{\mathrm{65}-\mathrm{66}}{\mathrm{Or}}_{\mathrm{20}-\mathrm{21}}$) and groundmass microlites (${\mathrm{Ab}}_{\mathrm{64}-\mathrm{66}}{\mathrm{Or}}_{\mathrm{21}-\mathrm{24}}$).

6.2 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 low-amplitude 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 clinopyroxene (Cpx2) phenocrysts from few trachyandesites (LH15, LH17, LH26-2).

Cpx1 shows a Mg# of 45–75, Ca of 0.78–0.90 apfu (atom per formula unit), Q+J of 1.84–1.94, and $J/(J+Q)$ of 0.03–0.06, and it can be classified as a Ti-rich augite (${\mathrm{Wo}}_{\mathrm{41}-\mathrm{48}}{\mathrm{En}}_{\mathrm{25}-\mathrm{42}}{\mathrm{Fs}}_{\mathrm{14}-\mathrm{32}}$).

Cpx2 shows a Mg# of 59–84, Ca of 0.20–0.92 apfu, Q+J of 1.77–1.95, and $J/(J+Q)$ of 0.01–0.06, and it can be classified as a diopside-rich augite (${\mathrm{Wo}}_{\mathrm{11}-\mathrm{48}}{\mathrm{En}}_{\mathrm{36}-\mathrm{64}}{\mathrm{Fs}}_{\mathrm{9}-\mathrm{32}}$).

Cpx3 shows a Mg# of 20–86, Ca of 0.27–0.97 apfu, Q+J of 1.57–1.98, and $J/(J+Q)$ of 0.01–0.07, and it can be classified as a diopside-rich augite (${\mathrm{Wo}}_{\mathrm{12}-\mathrm{49}}{\mathrm{En}}_{\mathrm{14}-\mathrm{57}}{\mathrm{Fs}}_{\mathrm{8}-\mathrm{62}}$). The composition of Cpx3 partially overlaps those of Cpx1 and Cpx2 groups, as would be expected for phenocrysts with homogeneous cores (i.e., Cpx1 and Cpx2) and the respective rims (Cpx3) characterized by low-amplitude oscillatory zoning or normal monotonous zoning textures (e.g., Streck, 2008).

Cpx4 shows a Mg# of 31–81, Ca of 0.24–0.87 apfu, Q+J of 1.87–1.97, and $J/(J+Q)$ of 0.01–0.06, corresponding to diopside-rich augite (${\mathrm{Wo}}_{\mathrm{12}-\mathrm{46}}{\mathrm{En}}_{\mathrm{18}-\mathrm{60}}{\mathrm{Fs}}_{\mathrm{11}-\mathrm{38}}$). The composition of Cpx4 partly overlaps that of Cpx3; however, their textural characteristics are completely different.

Cpx5 differs from previous pyroxenes, with a large spread in Mg# ranging 5–73, Ca of 0.03–0.83 apfu, Q+J of 1.51–2.07, and $J/(J+Q)$ of 0.07–0.89. The Cpx5 can be classified as aegirine-augite (Na<0.3 apfu, X_Aeg<0.30; with X_Aeg=Na apfu if $\mathrm{Na}<{\mathrm{Fe}}^{\mathrm{3}+\mathrm{Tot}}$, $X_{\mathrm{Aeg}}={\mathrm{Fe}}^{\mathrm{3}+\mathrm{Tot}}$ apfu if $\mathrm{Na}>{\mathrm{Fe}}^{\mathrm{3}+\mathrm{Tot}}$) to aegirine (Na = 0.68–0.88 apfu, X_Aeg = 0.40–0.88). Cpx5 clinopyroxenes are generally Ti-enriched (TiO₂ up to 2.8 wt %, Ti up to 0.08 apfu) and straddle the $Q+J=\mathrm{2}$ line defining the boundary for “normal” pyroxenes (Morimoto, 1989), thus indicating the presence of a ${{\mathrm{NaR}}^{\mathrm{2}+}}_{\mathrm{0.5}}{{\mathrm{Ti}}^{\mathrm{4}+}}_{\mathrm{0.5}}{\mathrm{Si}}_{\mathrm{2}}{\mathrm{O}}_{\mathrm{6}}$ component (Morimoto, 1988, 1989; Huraiová et al., 2017) (Fig. 6c).

Figure 6 Di-Hd-En-Fs, Q–J, and Ti vs. Na diagrams showing the composition of clinopyroxenes in LHPCS lavas. Symbol shapes follow Fig. 3. (a) Di-Hd-En-Fs diagram for clinopyroxene phenocrysts (Cpx1, Cpx2, Cpx3). (b) Di-Hd-En-Fs diagram for clinopyroxene microlites (Cpx4) and Na-clinopyroxenes (Cpx5). (c) Q–J diagram for pyroxenes with indication of endmembers (Morimoto, 1989). (d) Enlargement of area indicated in (c). (e) Ti vs. Na (apfu) diagram illustrating the compositional differences between clinopyroxenes. (f) Enlargement of area indicated in (e), showing the main augite trend characterizing the evolution from basalts to trachytes and the divergent trend of aegirine-augite and aegirine series.

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 aegirine-augite 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³⁺) content increasing during the magmatic differentiation, whereas the diopside-hedenbergite enrichment could testify to an increase of ferrous iron (Fe²⁺) in magma (e.g., Huraiová et al., 2017).

6.3 Olivine

Olivine is found in all analyzed samples, except for LH27-2 trachyandesite and LH6 trachyte. It consists of idiomorphic (Fig. 4a, c) to skeletal (e.g., Donaldson, 1974; Fowler et al., 2002; Faure et al., 2003; Welsch et al., 2013) (Fig. 4b) phenocrysts and microlites in the groundmass (Fig. 4h). Olivine crystals, both phenocrysts and microcrystals, show homogeneous cores and concordant normal monotonous zoning textures at outer rims. In basalts, olivine shows a continuous compositional range (Fig. 7a) from Fo₈₆Fa₁₄Mtc₀Tep₀ (phenocryst in LH5-2 basalt) to Fo₀₅Fa₉₁Mtc₁Tep₃ (groundmass microlites in LH27-1 basalt). The highest MnO (up to 1.7 wt %) values are systematically found in Fe-rich olivine microlites in basalts. Low monticellite concentration (CaO always <1.0 wt %) in LHPCS samples is typical for magmatic olivine (i.e., Melluso et al., 2014). CaO content positively correlates with the fayalite (FeO) compound (Fig. 5a). Together with the Mg#, the CaO content allows us to discriminate olivine phenocrysts in three coherent compositional clusters: (i) olivine from basalts, with Mg# = 79–87 and CaO = 0.21 wt %–0.73 wt %; (ii) olivine from trachyandesites, with Mg# = 67–80 and CaO = 0.08 wt %–0.43 wt %; and (iii) olivine from trachytes, with Mg# = 58–63 and CaO = 0.16 wt %–0.42 wt %. A minor number of analyzed phenocrysts in basalts show Cr₂O₃ content in the range ca. 0.05 wt %–0.07 wt %. Cr₂O₃ is, instead, below detection limit for almost all analyzed olivine crystals in LHPCS lavas. A minor cluster of peridote Mg-olivine (Fo₉₉Fa₁) xenocrysts, characterized by disequilibrium textures (resorption patterns) at rim, have been identified in LH26-1 trachyandesite lava.

Figure 7 (a) CaO vs. FeO diagram showing the composition of olivine in LHPCS lavas. (b) Di-Hd-En-Fs diagram showing the orthopyroxene chemistry in studied LHPCS lavas. (c) Al^Tot vs. Mg# diagram showing the main compositional differences between orthopyroxene populations from trachytes and trachyandesites. Mn (apfu) contents are also reported for the two populations.

6.4 Orthopyroxene

Orthopyroxene occurs in most of the LHPCS trachyandesite (Fig. 4f) and trachyte samples (Fig. 4g–h). Orthopyroxene phenocrysts are generally unzoned with homogeneous textures. In trachyandesites, they show intermediate (${\mathrm{En}}_{\mathrm{41}-\mathrm{83}}{\mathrm{Fs}}_{\mathrm{14}-\mathrm{55}}{\mathrm{Wo}}_{\mathrm{2}-\mathrm{10}}$) compositions (Fig. 7b), with Mg# of 43–86, Al₂O₃ up to 2.12 wt %, TiO₂ 0.08 wt %–1.33 wt %, and CaO 1.20 wt %–4.72 wt %. Similar compositions (${\mathrm{En}}_{\mathrm{62}-\mathrm{79}}{\mathrm{Fs}}_{\mathrm{18}-\mathrm{33}}{\mathrm{Wo}}_{\mathrm{3}-\mathrm{7}}$) have been obtained for microlites in groundmass (Fig. 7b) with Mg# of 65–81, Al₂O₃ of 0.48 wt %–1.53 wt %, TiO₂ of 0.21 wt %–0.60 wt %, and CaO of 1.35 wt %–3.49 wt %.

In trachytes, orthopyroxene phenocrysts present Mg# ranging 59–65, with low Al₂O₃ (0.18 wt %–0.73 wt %), low TiO₂ (0.11 wt %–0.32 wt %), and CaO (0.81 wt %–1.88 wt %), corresponding to a Fe-rich composition (${\mathrm{En}}_{\mathrm{56}-\mathrm{63}}{\mathrm{Fs}}_{\mathrm{34}-\mathrm{39}}{\mathrm{Wo}}_{\mathrm{2}-\mathrm{4}}$) with a minor Ca-Cpx substitution (Fig. 7b). Orthopyroxene microlites in groundmass (Fig. 7b) show comparable hypersthene (${\mathrm{En}}_{\mathrm{46}-\mathrm{60}}{\mathrm{Fs}}_{\mathrm{35}-\mathrm{45}}{\mathrm{Wo}}_{\mathrm{3}-\mathrm{7}}$) composition with Mg# of 50–63, Al₂O₃ 0.25 wt %–0.82 wt %, TiO₂ 0.19 wt %–0.31 wt % and CaO 1.32 wt %–3.27 wt %.

The compositional variation of orthopyroxene is summarized in the Al^Tot vs. Mg# diagram (Fig. 7c). Orthopyroxene crystals from trachyandesites are characterized by a higher content of Al (apfu) and higher Mg#, whereas those from trachytes are richer in ferrous iron (lower Mg# values) and in manganese (Mn up to 0.04 apfu).

6.5 Spinel and opaque minerals

Basalts show a diversified set of opaque minerals. Phenocrysts are (in order of abundance) (i) Al-spinel (TiO₂ 0.58 wt %–1.00 wt %; Mg# 58–71; Cr# 21–30, with [Cr# = 100 Cr∕(Cr+Al)]), (ii) Ti-magnetite (TiO₂ 1.83 wt %–21.58 wt %; MgO 0.06 wt %–2.19 wt %; MnO 0.44 wt %–0.63 wt %), and (iii) ilmenite (MgO up to 2.18 wt %). Groundmass is characterized by Fe-Ti oxides (ca. 20–30 µm in diameter; Fig. 4f) as ilmenite (MgO 0.27 wt %–1.50 wt %) and Ti-magnetite (MgO 0.18 wt %–1.89 wt %). Cr-spinels (TiO₂ 3.37 wt %–8.55 wt %; Mg# 14–28; Cr# 62–72) are found just as inclusions, up to 200 µm in diameter (Fig. 4a), in larger Mg-rich olivine phenocrysts.

Trachyandesites are characterized by phenocrysts of Ti-magnetite (MgO 0.07 wt %–3.84 wt %), ilmenite (MgO 1.11 wt %–4.79 wt %), and rare rutile (MgO 0.47 wt %). Groundmass microcrystals (ca. 20–30 µm in diameter) show a comparable composition with Ti-magnetite (MgO 0.33 wt %–3.77 wt %), ilmenite (MgO 0.33 wt %–4.79 wt %), and rare rutile (MgO<0.05 wt %). Similar to basalts, Cr-spinels (TiO₂ 6.09–6.47 wt %; Mg# 19–21; Cr# 65–68) are found only as inclusions (100–200 µm in diameter) in major Mg-rich olivine phenocrysts.

In trachytes, Fe-Ti oxides show euhedral to subhedral habit, and, based on chemistry, they are ilmenite (MgO 2.06 wt %–3.31 wt %) and Ti-magnetite (MgO 1.41 wt %–5.47 wt %). Phenocrysts (up to 50–100 µm in diameter; Fig. 4g) and groundmass microcrystals (ca. 15–20 µm in diameter) show the same compositions.

7 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, 2008), alkali feldspar (Putirka, 2008), and clinopyroxene (Putirka et al., 1996, 2003; Putirka, 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, 2008; Mordick 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₂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.

Clinopyroxene-liquid thermometry and barometry, for diopside-augite pyroxenes in basalts and trachyandesites (groups Cpx1, Cpx2, Cpx3, Cpx4), were calculated by the application of the Jd-DiHd exchange thermometer (Putirka et al., 1996, 2003) using Eq. (33) in Putirka (2008) and the Al-partitioning barometric model Eq. (32c) in Putirka (2008). Clinopyroxene-liquid thermometry and barometry, for diopside-augite pyroxenes in trachytes (groups Cpx3 and Cpx4), were calculated by the application of the Jd-DiHd exchange thermometer (Putirka et al., 1996, 2003; Putirka, 2008) recalibrated for alkaline differentiated magmas using Eqs. (Talk33) and (Palk2012), respectively, in Masotta et al. (2013). Clinopyroxene-liquid thermometry and barometry, for augite-aegirine pyroxenes (Cpx5), were calculated integrating Eqs. (33) and (32c) in Putirka (2008) with equations Eqs. (Talk2012) and (Palk2012) in Masotta et al. (2013). Results of mineral-melt equilibrium tests (Figs. 8, 9, 10), hygrometry calculations (Fig. 10), and geothermobarometric estimates are presented contextually in supplementary mineral chemistry tables. A summary of the thermobarometry estimates is reported in a pressure–temperature diagram (Fig. 11).

Figure 8 Rhodes diagrams showing the results of the test of equilibrium liquid and olivine (a) and orthopyroxene (b). The partitioning of Fe-Mg between mineral and liquid (Fe-Mg exchange coefficient) or $K_{\mathrm{D}}^{\mathrm{Min}-\mathrm{Liq}}$(Fe-Mg) is shown (black lines). The accepted range of equilibrium constant values for both panels (a) and (b) is indicated by dashed lines. $K_{\mathrm{D}}^{\mathrm{Min}-\mathrm{Liq}}$(Fe-Mg) values are from Putirka (2008). Nominal melt compositions are selected from whole-rock analyses. Vectors of olivine and orthopyroxene removal from melt and closed system differentiation are redrawn after Putirka (2008) and references therein. Fields of “Xenocrystal/antecrystal” and “Late crystallization” are also indicated. Symbols and colors refer to Fig. 7.

Figure 9 Clinopyroxene-melt equilibrium tests; (a) DiHd: diopside-hedenbergite, (b) EnFs: enstatite-ferrosilite, and (c) CaTs: Ca-Tschermak components. Equilibrium associated with observed components in pyroxenes are paired with predicted components in respective hosting melts. The accepted range of equilibrium is indicated in each figure by dashed lines. Nominal melt compositions for clinopyroxene are selected from whole-rock analyses. Symbols and colors refer to Fig. 6.

7.1 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 Fe-Mg between mineral and liquid is known as Fe-Mg exchange coefficient, or $K_{\mathrm{D}}^{\mathrm{Min}-\mathrm{Liq}}$(Fe-Mg) (defined as $K_{\mathrm{D}}^{\mathrm{Min}-\mathrm{Liq}}\left(\text{Fe-Mg}\right)=\left[{\mathrm{MgO}}^{\mathrm{Liq}}{\mathrm{FeO}}^{\mathrm{Min}}\right]/\left[{\mathrm{MgO}}^{\mathrm{Min}}{\mathrm{FeO}}^{\mathrm{Liq}}\right]$, where Liq is the liquid composition, Min is the mineral composition, and MgO and FeO are molar fractions; Roeder and Emslie, 1970; Langmuir and Hanson, 1980; Putirka, 2005a, 2008). It is used here to test the equilibrium between mafic minerals (olivine, orthopyroxene and clinopyroxene) and liquid (e.g., Maclennan et al., 2001; Putirka, 2008; Stroncik et al., 2009; Aulinas et al., 2010; Keiding and Sigmarsson, 2012; Melluso et al., 2014; Feng and Zhu, 2018).

We calculated $K_{\mathrm{D}}^{\mathrm{Min}-\mathrm{Liq}}\left(\mathrm{Fe}-\mathrm{Mg}\right)$ values using (i) Eq. (17) in Putirka (2008) for Ol and Opx, (ii) temperature-dependent Eq. (35) in Putirka (2008) for diopsidic-augitic Cpx in basalts and trachyandesites, and (iii) the Na-corrected Eq. (35a) in Masotta et al. (2013) for Na-rich Cpx5 group and for all Cpx from LH5-1 and LH6-1 trachytes. The calculated $K_{\mathrm{D}}^{\mathrm{Min}-\mathrm{Liq}}(\mathrm{Fe}-\mathrm{Mg})$ values for olivine and orthopyroxene are plotted in a Rhodes diagram (Dungan et al., 1978; Rhodes et al., 1979; Putirka, 2005a, 2008) to graphically test the equilibrium between Ol (Fig. 8a) or Opx (Fig. 8b) and the respective hosting melts (Liq). Furthermore, the Rhodes diagram is useful to recognize (i) presence of xenocrystals and/or antecrystals, (ii) late or groundmass crystallization, (iii) crystal removal (decrease of Mg#^Liq only), and (iv) closed system crystallization (decrease of Mg#^Min only) by deviations of the measured compositions from the expected ones (Rhodes et al., 1979; Putirka, 2008; Melluso et al., 2014).

The calculation of $K_{\mathrm{D}}^{\mathrm{Cpx}-\mathrm{Liq}}\left(\mathrm{Fe}-\mathrm{Mg}\right)$ does not consider variations of Ca and Al contents in Cpx (Rhodes et al., 1979; Putirka, 1999, 2005b, 2008). Therefore, a further equilibrium test was achieved through the comparison of analyzed Cpx compositions (as expressed by the components EnFs, DiHd, and CaTs, where CaTs is Ca-Tschermak) with component contents predicted from melt composition (e.g., Putirka, 2008; Mollo et al., 2010; Jeffery et al., 2013; Barker et al., 2015; Ellis et al., 2017). Normative components of Cpx were calculated following the scheme proposed in Putirka et al. (1996) and Putirka (2008). Calculation of Cpx components based on melt composition was performed using equations (Eq. 3.1a) for DiHd, (Eq. 3.2) for EnFs and (Eq. 3.4) for CaTs in Putirka (1999). A graphical presentation (e.g., Jeffery et al., 2013; Barker et al., 2015) of this test is shown in Fig. 9.

The partitioning of An-Ab between mineral and liquid is known as An-Ab exchange coefficient, or $K_D^{\mathrm{Pl}-\mathrm{Liq}}\left(\mathrm{An}-\mathrm{Ab}\right)$ (defined as $K_{\mathrm{D}}^{\mathrm{Pl}-\mathrm{Liq}}\left(\mathrm{An}-\mathrm{Ab}\right)=\left[X{\mathrm{Ab}}^{\mathrm{Pl}}X{\mathrm{AlO}}_{\mathrm{1.5}}^{\mathrm{Liq}}X{\mathrm{CaO}}^{\mathrm{Liq}}\right]/\left[X{\mathrm{An}}^{\mathrm{Pl}}X{\mathrm{NaO}}_{\mathrm{0.5}}^{\mathrm{Liq}}X{\mathrm{SiO}}_{\mathrm{2}}^{\mathrm{Liq}}\right]$, 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).

Figure 10 Plagioclase-melt equilibrium test. Equilibrium associated with anorthite (X_An) component in plagioclase are paired with predicted anorthite in melt. Nominal melt compositions for plagioclase are selected from whole-rock analyses. Calculated water concentrations using plagioclase-melt hygrometer (Putirka, 2008) are reported in diagrams with isolines (graded blue lines). Symbols and colors refer to Fig. 5.

In summary, we accept (i) Ol with $K_{\mathrm{D}}^{\mathrm{Ol}-\mathrm{Liq}}\left(\mathrm{Fe}-\mathrm{Mg}\right)=\mathrm{0.30}\pm \mathrm{0.06}$ (Roeder and Emslie, 1970; Putirka, 2005a, 2008, and references therein) (Fig. 8a), (ii) Opx with $K_{\mathrm{D}}^{\mathrm{Opx}-\mathrm{Liq}}\left(\mathrm{Fe}-\mathrm{Mg}\right)=\mathrm{0.29}\pm \mathrm{0.06}$ (Putirka, 2008, and references therein) (Fig. 8b), (iii) Cpx with $K_{\mathrm{D}}^{\mathrm{Cpx}-\mathrm{Liq}}\left(\mathrm{Fe}-\mathrm{Mg}\right)=\mathrm{0.28}\pm \mathrm{0.08}$ (Putirka, 2008) and verifying the one-to-one (±0.1) relationship between predicted vs. observed normative components (EnFs, DiHd and CaTs) for at least two of the monitored components (Fig. 9), (iv) Pl with $K_{\mathrm{D}}^{\mathrm{Pl}-\mathrm{Liq}}\left(\mathrm{An}-\mathrm{Ab}\right)=\mathrm{0.27}\pm \mathrm{0.11}$ for T>1050 ^∘C and $K_{\mathrm{D}}^{\mathrm{Pl}-\mathrm{Liq}}\left(\mathrm{An}-\mathrm{Ab}\right)=\mathrm{0.10}\pm \mathrm{0.05}$ for T<1050 ^∘C (Putirka, 2008) or falling within ±0.1 of the one-to-one relationship between predicted vs. observed An components (Fig. 10), and (v) Afs with $K_{\mathrm{D}}^{\mathrm{Afs}-\mathrm{Liq}}\left(\mathrm{An}-\mathrm{Ab}\right)=\mathrm{0.27}\pm \mathrm{0.18}$ (Putirka, 2008). All mineral-liquid pairs exceeding the accepted exchange coefficient values for Ol, Cpx, Opx, and Fsp were discarded for thermobarometric analyses.

7.2 Pre-eruptive H₂O^Liq content estimates

Thermobarometric models for volcanic systems require an initial estimate of the pre-eruptive water concentration (wt %) in melt (H₂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₂O^Liq wt % values (±1σ standard deviation of the weighted mean) are plotted as isolines in Fig. 10. The hygrometer of Putirka (2008) indicates (Fig. 10) (i) H₂O^Liq negative values in basalts, from −0.20 wt % to −0.40 wt %, with a weighted mean of $-\mathrm{0.37}\pm \mathrm{0.20}$ wt % (mean square weighted deviation, MSWD = 0.0026; n=95); (ii) trachyandesites pre-eruptive water content in the H₂O^Liq range 0 wt %–1.40 wt % (weighted mean of 0.57±0.13 wt %, MSWD = 0.13, n=245); and (iii) trachytes with the highest water concentration (H₂O^Liq: 1.40 wt %–1.90 wt %; weighted mean of 1.46±0.32 wt %, MSWD = 0.059, n=37). Following the approach of Keiding and Sigmarsson (2012), negative values in basalts are interpreted as anhydrous melt compositions. Coherently with the existing literature (e.g., Webster et al., 1999), the anhydrous character is then assumed as a H₂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₂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).

On the contrary, existing studies (e.g., Kushiro, 1969; Sisson and Grove; 1993; Yang et al., 1996; Putirka, 2005a, b, 2008; Kelley and Barton, 2008; Keiding and Sigmarsson, 2012) demonstrated a negligible effect of water for basaltic and intermediate melts showing H₂O^Liq ranging 0 wt %–1 wt %.

7.3 Thermobarometry results

7.3.1 Basalts

When applied to phenocryst cores, the Pl-Liq thermobarometry (Fig. 11a–c) shows that all LHPCS basaltic materials have magmatic anhydrous temperature in the range 1230–1266 ^∘C (weighted mean of 1250±5 ^∘C, ±1σ standard deviation of the weighted mean, MSWD = 0.112, n=95). Pressure estimates are in the range 6.5–8.7 kbar (weighted mean of 7.9±1.1 kbar, ±1σ standard deviation of the weighted mean, MSWD = 0.024, n=28) for LH18 Ol-basalts and 7.2–10.3 kbar (weighted mean of 9.2±0.7 kbar (±1σ), MSWD = 0.064, n=67) for LH5-2 and LH27-1 Ol-Cpx-basalts. Olivine-melt equilibrium (Fig. 11a–c), for the olivine compositional range of Fo 80 %–85 %, yields a temperature window of 1240–1297±27 ^∘C (±1σ), consistent with the results obtained with Pl-Liq thermometry. The Cpx thermobarometry (Fig. 11a, c), for both Cpx1 (phenocryst cores) and Cpx3 (phenocryst rims and unzoned phenocrysts), provides temperatures of 1006–1209 ^∘C (weighted mean of 1124±12 ^∘C (±1σ), MSWD = 3.4, n=82). Pressure estimates are in the range 3.1–11.5 kbar (weighted mean of 7.6±0.8 kbar (±1σ), MSWD = 2.7, n=36) for Cpx1 and 2.5–7.7 kbar (weighted mean of 4.0±0.8 kbar (±1σ), MSWD = 0.63, n=14) for Cpx3. Thermobaric estimates for Cpx4 (microlites in groundmass) indicate shallow conditions (0.3–3.0; weighted mean of 1.6±1.2 kbar (±1σ), MSWD = 0.38, n=6) for temperatures (1006–1123 ^∘C; weighted mean of 1060±54 ^∘C (±1σ), MSWD = 2.9, n=6) comparable to those obtained for Cpx1 and Cpx3. Higher temperature estimates (1067–1221 ^∘C; weighted mean of 1157±53 ^∘C (±1σ), MSWD = 2.4, n=7) at low-pressure conditions (0.4-4.7; weighted mean of 2.9±1.1 kbar (±1σ), MSWD = 0.83, n=7) are instead obtained for a limited number of Cpx5 (aegirine-rich) compositions (Fig. 11a, c).

Figure 11 Thermobarometric estimates; (a) a summary of the results obtained from thermobarometry models applied to Los Humeros post-caldera-stage lavas. Symbols refer to whole-rock chemistry compositions, whereas colors of different phases refer to mineral chemistry diagrams. Green-shaded fields show the results of olivine-liquid thermometry. Blue dashed isolines represent the results of plagioclase-melt hygrometer. Yellow-shaded fields indicate the pressure–temperature domain of crystallization of feldspars in groundmass. Basalt liquidus curve, Ol + Cpx + Pl + Spl + Liq stability fields and point “A” (basalt liquidus in equilibrium with mantle peridotite mineral assemblage of Ol + Cpx) are redrawn after Grove (2000). Opx-in stability curve is redrawn after Wallace and Anderson Jr. (2000). Schematized results are presented separately for (b) Cpx-free Ol-basalt, (c) Cpx-bearing Ol-basalts, (d) trachyandesites, and (e) trachytes.

7.3.2 Trachyandesites

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

Opx-free trachyandesites LH15 and LH21-2 (El Limón and Sarabia lava flows, respectively) are characterized by (i) plagioclase phenocryst cores crystallized at temperatures of 1190–1263 ^∘C (weighted mean of 1248±7 ^∘C (±1σ), MSWD = 1.09, n=39) and pressures of 4.8–9.4 kbar (weighted mean of 7.7±0.9 kbar (±1σ), MSWD = 0.14, n=39); (ii) comparable temperature (1193–1263 ^∘C; weighted mean of 1227±37 ^∘C (±1σ), MSWD = 2.3, n=6) and pressure (6.7–9.6 kbar, mean value of 7.8±2.4 kbar (±1σ), MSWD = 0.101, n=6) obtained for rare phenocryst rims and microlites at equilibrium; (iii) olivine-melt equilibrium (with Fo: 75 %–80 %) showing a temperature window of 1030–1055 ± 27 ^∘C (±1σ); (iv) rare Cpx2 (clinopyroxene phenocryst cores) showing equilibrium with melt and yielding temperature of 1061–1239 ^∘C (weighted mean of 1116±29 ^∘C (±1σ), MSWD = 2.3, n=12) and pressures of ca. 2.9–8.3 kbar (weighted mean of 5.2±1.2 kbar (±1σ), MSWD = 1.5, n=12); (v) Cpx3 (rims of and unzoned phenocrysts) showing equilibrium with melt and yielding thermobarometric results (temperatures of 938–1139 ^∘C, with weighted mean of 1074±15 ^∘C (±1σ), MSWD = 1.9, n=32; and pressures of 1.0-4.4 kbar, with weighted mean of 2.8±0.5 kbar (±1σ), MSWD = 0.22, n=32); (vi) Cpx4 (groundmass microcrystals) compositions indicating, with respect to Cpx3, comparable temperatures (1026–1127 ^∘C, with weighted mean of 1059±16 ^∘C (±1σ), MSWD = 0.71, n=14) at lower pressure conditions (0.3–3.6 kbar with weighted mean of 1.4±0.8 kbar (±1σ), MSWD = 0.35, n=14). The unique Cpx5-liquid pair at equilibrium yielded P–T conditions of 5.6±1.5 kbar and 1122±30 ^∘C.

Thermobarometric estimates (Fig. 11a, d) for Opx-bearing trachyandesites (LH4, LH13, LH17; LH26-1; LH26-2; LH27-2) show overlapping P–T conditions for plagioclase populations with (i) phenocryst cores crystallizing at temperatures of 1145–1228 ^∘C (weighted mean of 1187±4 ^∘C (±1σ), MSWD = 1.17, n=166) and pressures of 4.1–7.7 kbar (weighted mean of 5.8±0.5 kbar (±1σ), MSWD = 0.059, n=166) and (ii) phenocryst rims and microcrystals forming at temperatures of 1140–1224 ^∘C (weighted mean of 1168±8 ^∘C (±1σ), MSWD = 0.92, n=34) and pressures of 4.4–8.5 kbar (weighted mean of 6.4±1.0 kbar (±1σ), MSWD = 0.14, n=34). Lower temperatures (1050–1090 ± 27 ^∘C (±1σ)) are obtained using olivine-liquid (Fo 70 %–80 %) equilibrium model.

Thermobarometers applied to pyroxenes indicate (i) Cpx2 (phenocryst cores) crystallizing at temperatures of 979–1204 ^∘C (weighted mean of 1060±8 ^∘C (±1σ), MSWD = 1.8, n=101) and pressures of 3.4–11.5 kbar (weighted mean of 7.0±0.3 kbar (±1σ), MSWD = 0.94, n=101), (ii) Cpx3 crystallizing at temperatures of 959-1106 ^∘C (weighted mean of 1026±6 ^∘C (±1σ), MSWD = 1.3, n=145) and pressures of 1.2–6.9 kbar (weighted mean of 4.3±0.2 kbar (±1σ), MSWD = 0.72, n=145), (iii) rare Cpx4 showing general equilibrium with melt and forming at P–T conditions of temperatures of 920–1123 ^∘C (weighted mean of 1020±21 ^∘C (±1σ), MSWD = 2.7, n=24) and pressures of 0.1–3.4 kbar (weighted mean of 1.8±0.6 kbar (±1σ), MSWD = 0.56, n=24), and (iv) Opx yielding crystallization conditions, for both phenocrysts and microlites, of temperatures of 1048–1123 ^∘C (weighted mean of 1078±5 ^∘C (±1σ), MSWD = 0.24, n=129) and pressures of 0–2.8 kbar (weighted mean of 1.1±0.6 kbar (±1σ), MSWD = 0.057, n=84). In all trachyandesite samples, temperature estimates obtained through the Ol-Liq model and the Cpx-Liq model are comparable (Fig. 11a, d), whereas the Pl-Liq model shows higher temperature values. These results can be interpreted as an earlier plagioclase crystallization with respect to olivine and clinopyroxene. Furthermore, orthopyroxene (Opx) can be considered a tracer of trachyandesitic magma stagnation at shallow depths, as indicated by the lower pressure estimates obtained through the Opx-liquid barometer.

7.3.3 Trachytes

Magmatic P–T conditions (Fig. 11a, e) of trachytic (LH5-1 and LH6) melts are defined by (i) plagioclase crystallization at temperature of 1050–1094 ^∘C (weighted mean of 1069±6 ^∘C (±1σ), MSWD = 0.39, n=37) and pressures of 4.7–9.0 kbar (weighted mean of 6.5±1.0 kbar (±1σ), MSWD = 0.20, n=37), (ii) olivine-liquid (Fo 55 %–65 %) regression indicating olivine crystallization at 900–920 ± 27 ^∘C (±1σ), and (iii) clinopyroxene crystallization, both phenocrysts (Cpx3) and groundmass (Cpx4), at temperature of ca. 955 ^∘C (weighted mean of 956±14 ^∘C (±1σ), MSWD = 0.00056, n=17) and very shallow depth conditions (pressure weighted means of 2.3±0.9 kbar (±1σ), MSWD = 0.047, n=10, and 1.6±1.1 kbar (±1σ), MSWD = 0.04, n=7, for Cpx3 and Cpx4, respectively). Shallow-depth conditions are also obtained for orthopyroxene crystallization with temperatures of 960–1006 ^∘C (weighted mean of 990±7 ^∘C (±1σ), MSWD = 0.28, n=49) and pressures of 0.2–3.6 kbar (weighted mean of 1.6±0.9 kbar (±1σ), MSWD = 0.101, n=35). The alkali-feldspar-liquid thermometer provided temperature estimates always <500 ^∘C, here interpreted as feldspar re-equilibration in subsolvus and/or subsolidus post-eruptive conditions (Nekvasil, 1992; Brown and Parsons, 1994; Plümper and Putnis, 2009; Kontonikas-Charos et al., 2017; Latutrie et al., 2017). Interestingly, temperatures obtained through the Pl-Liq model are higher than those obtained with Ol-Liq, Cpx-Liq, or Opx-Liq, suggesting an earlier crystallization of plagioclase with respect to mafic minerals. Moreover, the Pl-Liq models indicate thermobaric estimates comparable to those obtained for trachyandesitic rocks.

8 Discussion

8.1 Major-element mass balance modeling

Based on the documenting of textural evidence – (i) Cpx-bearing 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; Lucci et al., 2016). The FC modeling is focused on the following hypotheses: (i) a direct cogenetic relationship between all LHPCS basalts and (ii) a common genesis for all LHPCS trachyandesites and trachytes through differentiation via fractional crystallization starting from the same basaltic parental melt.

Major-element mass balance models (e.g., Bryan et al., 1969) can be used to test and define relative proportions of phases involved in Rayleigh fractional crystallization (RFC, Daughter = Parent – fractionating assemblage) and crystal accumulation (Cumulate = Melt + accumulated assemblage) hypotheses (e.g., White et al., 2009; Moghadam et al., 2016; Lucci et al., 2016).

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²) as

$$\begin{array}{}\text{(1)}& \sum r^{\mathrm{2}}=\sum _{i-\mathrm{1}}^n{\left(b_i^{\prime }-b_i\right)}^{\mathrm{2}}.\end{array}$$

RFC and Cumulate model results are considered acceptable when Σr²<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₂-TiO₂-Al₂O₃-FeO^*-MnO-MgO-CaO-Na₂O-K₂O. The LH5-2 Cpx-bearing basalt, with the lowest SiO₂ 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 RFC modeling has been applied to all studied trachyandesites and trachytes. It was verified that a fractionation of the Pl + Cpx + Ol + Spl assemblage in the range of (i) 45 wt %–63 wt % (Σr² 0.37–0.92) is necessary to produce Opx-free trachyandesites, (ii) 59 wt %–69 wt % (Σr² 0.38–0.92) is capable to produce Opx-bearing trachyandesites, and (iii) 73 wt %–74 wt % (Σr² 0.88–0.91) is requested to produce trachytes. The crystal accumulation has been tested to verify the genetic linkage between Cpx-free basalt (LH18-1) and Cpx-bearing basalts (LH5-2, LH27-1). It was verified that a crystal accumulation of the Pl + Cpx + Ol + Spl assemblage in the range of 16 wt %–17 wt % (Σr² 0.05–0.15), with Cpx ranging 5 wt %–7 wt %, can produce the LHPCS Cpx-bearing basalts.

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₂O < 1 wt %; e.g., Webster et al., 1999) mafic parental melt. Integrating FC model and hygrometer (Putirka, 2008) results, LHPCS trachytes show H₂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₂O in the range 0.3 wt %–0.5 wt %.

8.2 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 high-temperature (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 olivine-supersaturated 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).

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.

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.

Where these magmas erupted as intracaldera basalts (LH5-2, LH27-1), they are characterized by the further crystallization of (i) progressively Fe-rich olivine (up to X_Fo = 17 %–20 %), (ii) Ab-rich plagioclase (X_An = 25 %–30 %), (iii) Cpx3 unzoned homogeneous phenocrysts and overgrowth (normal monotonous and normal low-amplitude oscillatory zoning textures) on Cpx1-cores, (iv) Cpx4 (Di-rich) microcrystals and microlites, and (v) Cpx5 (Aeg-Aug) Na clinopyroxenes. This mineral assemblage (mineral chemistry and textures), together with the obtained thermobarometric results, describes a near-isothermal magma uprising within a narrow temperature window of ca. 1070–1150 ^∘C. Such crystal-bearing magmas ascend from the deeper reservoir to intermediate and shallower stagnation levels, where different phases would crystallize, before the eruption (e.g., Feng and Zhu, 2018). In particular, (i) the homogeneous unzoned cores of phenocrysts represent the early crystallization at equilibrium with the melt; (ii) the normal low-amplitude oscillatory zoning textures at rims of Pl and Cpx phenocrystals indicate a kinetically driven crystallization (e.g., Ginibre et al., 2002; Streck, 2008; Renjith, 2014), whereas the normal monotonous zoning textures at rims observed in many Pl, Cpx, and Ol phenocrysts indicate a fast growth during ascent of the magma (e.g., Streck, 2008); (iii) microlite formation indicates water-exsolution-driven crystallization (e.g., Rutherford, 2008; Renjith, 2014) during a relatively rapid ascent or eruption process (e.g., Renjith, 2014); and (iv) the similarity of compositions between Pl and Cpx phenocrysts rims and microlites confirms that there were essentially no major changes in the temperature of any of these basaltic magmas during the ascent (e.g., Rutherford, 2008). This scenario of rapid ascent of LHPCS basaltic magmas is also supported by the observed high-vesicularity textures, interpreted as bubble-growth processes during a relatively fast magma rise precluding exsolved volatiles to escape (e.g., Sparks, 1978; Sparks et al., 1998; Rutherford and Gardner, 2000; Rutherford, 2008; Costa et al., 2013; Feng and Zhu, 2018).

Fractional crystallization of An-rich plagioclase, Fo-rich olivine, Ti-rich augite, and spinel (Pl + Ol + Cpx + Spl in RFC models) in the primary Cpx-bearing basaltic magmas produces residual melts (ca. 30 wt %–55 wt %) of trachyandesitic compositions. These evolved buoyant melts will be prone to leave the basaltic reservoir to produce shallower intrusions in a vertically extensive magmatic system (e.g., Jackson et al., 2018), carrying early-formed phenocrysts (i.e., anorthitic plagioclase antecrystals) to the intermediate reservoir and stall. Within this intermediate, vertically distributed, layered storage system in the middle crust, Cpx2 clinopyroxene and all the rest of the plagioclase phenocrysts start to crystallize, producing progressively evolved felsic residual melts able to migrate upward in the feeding system or erupt (e.g., Freundt and Schminke, 1995; Patanè et al., 2003; Klügel et al., 2005; Stroncik et al., 2009; Aulinas et al., 2010; Dahren et al., 2012; Keiding and Sigmarsson, 2012; Scott et al., 2012; Jeffery et al., 2013; Coombs and Gardner, 2001; Barker et al., 2015; Feng and Zhu, 2018). Similarly to LHPCS basalts, the phenocryst morphologies and textures, together with the microlites compositions and the vesicle-rich textures described in trachyandesitic melts, suggest a nearly isothermal rapid ascent precluding exsolved volatiles to escape and producing water-exsolution-driven crystallization (e.g., Rutherford, 2008; Renjith, 2014).

The shallowest magma stagnation level (<3 kbar; mean 1.5 kbar) has been here interpreted as a complex magma 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 LHPCS 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 assimilation 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 mineral-melt 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 An-rich cores can be crystallized in two possible scenarios. The first one is related to the H₂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 An-rich 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 An-rich 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 An-rich 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) high-temperature 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.

8.3 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⁻³ for unconsolidated sediments via about 3000 kg m⁻³ for the lower crust to 3300 kg m⁻³ 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⁻³), (ii) middle sediments (thickness: 0.5 km, density 2370 kg m⁻³), (iii) upper crust (thickness: 13.6 km, density 2740 kg m⁻³), (iv) middle crust (thickness: 15.3 km, density 2830 kg m⁻³), and (v) lower crust (thickness: 13.6 km, density 2920 kg m⁻³). Inferred (seismic) Moho depth is reported at −41.7 km with an upper mantle density of 3310 kg m⁻³ (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 deep-seated anhydrous Ol-basalt reservoir at depths of ca. 28–33 km (7.6–9.2 kbar), at the boundary between lower and middle crust, below the caldera is recorded by (i) An-rich Pl cores (X_An = 50 %–70 %) and (ii) Ti-rich augitic Cpx1 cores (Mg# up to 75, TiO₂ up to 4.57 wt %). For this mafic reservoir, the overlapping of the calculated anhydrous temperature estimates as derived from Pl-Liq, Cpx1-Liq, and Ol-Liq pairs spans ca. 1000–1300 ^∘C. The highest anhydrous temperature values are derived from the Cpx-free Ol-basalt Texcal lava flow (LH18), where the convergence of Pl-Liq thermobarometry and Ol-Liq thermometry models indicate conditions of temperature of ca. 1230–1270 ^∘C at a pressure of ca. 8 kbar. Lower anhydrous temperatures of ca. 1000–1210 ^∘C are obtained at comparable average pressure values for Cpx-bearing intracaldera Ol-basalts (LH5-2; LH 27-1). These results are in agreement with existing literature on the near-liquidus melting behavior of high-Al basaltic magmas (Mg# ca. 60–70 and Al₂O₃: 17 wt %–19 wt %) under dry conditions (e.g., Thompson, 1974; Grove et al., 1982; Crawford et al., 1987; Bartels et al., 1991; Grove, 2000). At 1250–1300 ^∘C and ca. 10 kbar (point A in Fig. 11) the basaltic melt is in equilibrium with a mantle peridotite mineral assemblage of olivine + clinopyroxene (Kushiro and Yoder Jr., 1966; Presnall et al., 1978; Grove et al., 1982; Fuji and Scarfe, 1985; Takahashi, 1986; Falloon and Green, 1987; Bartels et al., 1991; Sisson and Layne, 1993; Wagner et al., 1995; Grove et al., 1997; Grove, 2000; Kinzler et al., 2000). Following the models proposed by Thompson (1974), Bartels et al. (1991), and Grove (2000), a temperature decrease would lead primary melts to pass the “dry basaltic liquidus” and start the crystallization of Ol + Pl (higher temperatures) or Ol + Cpx + Pl (lower temperatures) assemblages (see stability fields in Fig. 11). Given the ubiquitous presence in all LHPCS basalts of well-developed euhedral to subhedral olivine crystals (both as phenocrysts and microlites) at equilibrium with anorthitic plagioclase, it is possible to exclude that crystallization history started at depths > ca. 10–12 kbar where olivine is not a stable phase and the primary assemblage would be characterized only by Cpx + Pl + Spl in equilibrium with melt (Kushiro and Yoder Jr., 1966; Thompson, 1974; Presnall et al., 1978; Bartels et al., 1991; Grove, 2000).

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₂ 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 nucleation-dominated 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 earlier-formed 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).

The fourth shallowest storage zone located at depths of ca. 3–7 km (weighted mean value of 1.5±0.2 kbar (±1σ), MSWD = 0.24, n=177; P ranging ca. 0.1–4.5 kbar) is required to explain the presence of (i) Cpx4 clinopyroxene (microcrystals and microlites) in all LHPCS lavas, (ii) Aeg-rich Cpx5 clinopyroxene in basalts, and (iii) Fe olivine (Fo = 55 %–65 %) and orthopyroxene in Opx-bearing evolved LHPCS lavas. Magmas in this shallow storage system show a wide range of temperature values calculated for hydrous melts: (i) ca. 1060 ^∘C (weighted mean value, MSWD = 2.4, n=7) for Aeg-rich Cpx5 crystallization in basalts; (ii) ca. 1070 ^∘C (weighted mean value, MSWD = 1.09, n=168) for Cpx4, and Opx crystallization in trachyandesites; and (iii) ca. 965 ^∘C (weighted mean value, MSWD = 2.2, n=78) for olivine, Cpx4 and Opx crystallization in trachytes. Thermobaric estimates obtained for Aeg-rich Cpx5 agree with those calculated for transitional basalts at Pantelleria (White et al., 2009, and references therein), whereas orthopyroxene crystallization conditions overlap with the existing literature for intermediate rocks (e.g., Rutherford et al., 1985; Wallace and Anderson Jr., 2000; Reubi and Nicholls, 2004; Allan et al., 2013; Jeffery et al., 2013). The broad distribution of melt chemistry from basalt to trachyte, together with the obtained thermobaric estimates, defines a shallow magma storage environment characterized by progressive accumulation of small, locally interconnected magma pockets and batches (e.g., Reubi and Nicholls, 2004; Jeffery et al., 2013) dispersed in the upper crust (<10 km) with a possible magma plexus at a depth of 2–4 km under the caldera (e.g., Armienti et al., 1989; Freundt and Schminke, 1995; Pietruszka and Garcia, 1999; Patanè et al., 2003; Klügel et al., 2005; Stroncik et al., 2009; Dahren et al., 2012; Jeffery et al., 2013; Coombs and Gardner, 2001, 2004).

8.4 Standard Model versus multilayered magmatic plumbing system

Existing conceptual models for LHVC are based on the Standard Model (sensu Gualda and Ghiorso, 2013), considering a single bowl-shaped long-lived melt-dominated magma chamber of 1000–1500 km³, at depths of 5 to 10 km (Verma, 1983, 1984, 1985a, b; Verma and López, 1982; Verma et al., 1990, 2011; Verma and Gomez-Arias, 2013; Verma and Andaverde, 1995; Carrasco-Núñez et al., 2018). However, these models mainly refer to the Los Humeros caldera-stage activity (Carrasco-Núñez et al., 2018, and references therein), lasting ca. 130 kyr, where the major caldera-forming events (Xaltipan and Zaragoza ignimbrites, 115 and 15 km³ DRE, respectively) and the large Plinian eruptive episode (Faby Tuff, 10 km³ DRE) necessitated feeding from a huge, voluminous magma chamber (Carrasco-Núñez and Branney, 2005; Carrasco-Núñez et al., 2018).

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³ of mafic lavas, 10 km³ of intermediate and felsic lava, and 1 km³ 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 post-caldera 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 Member 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., Nairn et al., 1998; 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 magma-ponding 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.

8.5 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³) 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³) 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 background-positive 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 LHPCS 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.

9 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.

Appendix A: Analytical details

A1 Petrography of volcanic samples

Rock magmatic fabrics, textures, and mineral assemblages were studied on polished thin sections, using a Nikon Eclipse 50iPol polarized light microscope (PLM) equipped with a Nikon Ds-Fi2 CCD camera (Nikon, Tokyo, Japan) and Nikon Nis-Elements software (Ver4.30.01), at Laboratorio di Microtettonica, Dip. Science, Università Roma Tre (Rome, Italy). Mineral abbreviations follow 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. (2015, 2018a, b). 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 Developments). 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.

Validation of mineral chemistry results was also achieved through opportune comparisons with the existing literature for (i) Mg-olivine (e.g., Hirano et al., 2004; White et al., 2009; Giordano et al., 2012; Melluso et al., 2014); (ii) Fe-olivine (e.g., Aldanmaz, 2006; White et al., 2009; Melluso et al., 2010; Giordano et al., 2012); (iii) aegirine-augite clinopyroxene (Cpx5 group; e.g., Piilonen et al., 1998; White et al., 2009; Njonfang et al., 2013); (iv) augite-diopside clinopyroxene (Cpx1-4 groups; e.g., Dawson and Hill, 1998; Aldanmaz, 2006; Melluso et al., 2010, 2014); (v) orthopyroxene (e.g., Papike et al., 1995; Aldanmaz, 2006; Carvalho and de Assis Janasi, 2012; Hu et al., 2018); (vi) feldspar (e.g., Keil et al., 1972; Giordano et al., 2012; Innocenti et al., 2013; Njonfang et al., 2013); and (vii) spinel and opaque minerals (e.g., Melluso et al., 2014). Back-scattered electron (BSE) imaging was obtained by using the same electron microprobe with operating conditions of 15 kV and 50 nA. Mineral structural formulae of feldspar, olivine, and spinel were calculated with the software CalcMin_32 (Brandelik, 2009). Mineral structural formulae of orthopyroxene were calculated following Putirka et al. (1996) and Putirka (2008). Clinopyroxene formula has been calculated following procedures reported in Putirka et al. (1996), Putirka (2008), and Masotta et al. (2013). Clinopyroxenes were then classified integrating the Wo-En-Fs scheme (Morimoto, 1989) and J vs. Q scheme (Morimoto, 1988, 1989) with J = 2 Na apfu and $Q=(\mathrm{Ca}+\mathrm{Mg}+{\mathrm{Fe}}^{\mathrm{2}+})$ apfu. Aegirine (X_Aeg) component correction, for Na-rich Cpx (aegirine-augite series), followed the scheme (X_Aeg = Na apfu if $\mathrm{Na}<{\mathrm{Fe}}^{\mathrm{3}+\mathrm{Tot}}$, $X_{\mathrm{Aeg}}={\mathrm{Fe}}^{\mathrm{3}+\mathrm{Tot}}$ apfu if $\mathrm{Na}>{\mathrm{Fe}}^{\mathrm{3}+\mathrm{Tot}}$) proposed by Putirka et al. (1996), Putirka (2008), and based on ${\mathrm{Fe}}^{\mathrm{2}+}-{\mathrm{Fe}}^{\mathrm{3}+}$ correction of Lindsley (1983).