A petrological and conceptual model of Mayon volcano (Philippines) as an example of an open-vent volcano

Abstract

Mayon is a basaltic andesitic, open-vent volcano characterized by persistent passive degassing from the summit at 2463 m above sea level. Mid-size (< 0.1 km³) and mildly explosive eruptions and occasional phreatic eruptions have occurred approximately every 10 years for over a hundred years. Mayon’s plumbing system structure, processes, and time scales driving its eruptions are still not well-known, despite being the most active volcano in the Philippines. We investigated the petrology and geochemistry of its crystal-rich lavas (~ 50 vol% phenocrysts) from nine historical eruptions between 1928 and 2009 and propose a conceptual model of the processes and magmatic architecture that led to the eruptions. The whole-rock geochemistry and mineral assemblage (plagioclase + orthopyroxene + clinopyroxene + Fe-Ti oxide ± olivine) of the lavas have remained remarkably homogenous (54 wt% SiO₂, ~ 4 wt% MgO) from 1928 to 2009. However, electron microscope images and microprobe analyses of the phenocrysts and the existence of three types of glomerocrysts testify to a range of magmatic processes, including long-term magma residence, magma mixing, crystallization, volatile fluxing, and degassing. Multiple mineral-melt geothermobarometers suggest a relatively thermally buffered system at 1050 ± 25 °C, with several magma residence zones, ranging from close to the surface, through reservoirs at ~ 4–5 km, and as deep as ~ 20 km. Diffusion chronometry on > 200 orthopyroxene crystals reveal magma mixing timescales that range from a few days to about 65 years, but the majority are shorter than the decadal inter-eruptive repose period. This implies that magma intrusion at Mayon has been nearly continuous over the studied time period, with limited crystal recycling from one eruption to the next. The variety of plagioclase textures and zoning patterns reflect fluxing of volatiles from depth to shallower melts through which they eventually reach the atmosphere through an open conduit. The crystal-rich nature of the erupted magmas may have developed during each inter-eruptive period. We propose that Mayon has behaved over almost 100 years as a steady state system, with limited variations in eruption frequency, degassing flux, magma composition, and crystal content that are mainly determined by the amount and composition of deep magma and volatile input in the system. We explore how Mayon volcano’s processes and working model can be related to other open-vent mafic and water-rich systems such as Etna, Stromboli, Villarrica, or Llaima. Finally, our understanding of open-vent, persistently active volcanoes is rooted in historical observations, but volcano behavior can evolve over longer time frames. We speculate that these volcanoes produce specific plagioclase textures that can be used to identify similar volcanic behavior in the geologic record.

Introduction

Open-vent, mafic volcanoes are characterized by prominent volcanic gas plumes and persistent passive degassing from summit vents and/or lava lakes during quiescence, often with limited or no magma erupted (Rose et al. 2013). This persistent degassing is punctuated by a diverse range of eruptive behavior, from effusive and almost dome-like extrusion of viscous lavas to explosive eruptions with a wide range of compositions (Rose et al. 2013). Open-vent volcanoes are found in all tectonic settings and include some of the most active terrestrial volcanoes, such as Stromboli (Italy), Villarrica (Chile), Merapi (Indonesia), Erta Ale (Ethiopia), and Kilauea (USA), among others (e.g., Rose et al. 2013; Girona et al. 2014; Vergniolle and Métrich 2021). Many subduction zone open-vent volcanoes typically produce small to mid-sized eruptions (0.1 km³; Volcanic Explosivity Index, VEI, 1 to 3), but they are also capable of larger eruptions (> VEI 3) that occur less frequently (Rose et al. 2013). Most magmas are basalt to basaltic andesite, and crystal-rich lavas with high phenocryst contents (20–45 vol%) are the volumetrically dominant deposits.

The near continuous degassing of these volcanoes implies that there is a more or less permanent physical connection between the surface and the magma reservoir(s) that supply the volcanic gas. This condition allows the use of volcanic gas and fluxes (e.g., Allard 1997; Allard et al. 2008; Girona et al. 2014) for an increased understanding of the processes occurring in the plumbing system. The gas emissions during quiescence have been explained by magma convection in the conduit and degassing near the earth’s surface (Kazahaya et al. 1994; Stevenson & Blake 1998; Burton et al. 2007; Allard et al. 2016) and/or by gas percolation through the melt in viscous crystal-rich conduits (Boudreau 2016; Ruth et al. 2016; Edmonds et al. 2019). The daily gas flux that is produced by these volcanoes can be very large, on the order of several thousand or tens of thousands of tonnes per day (mainly H₂O, but also CO₂ and SO₂). Such mass release may have a large influence on the plumbing system and therefore be one of the drivers of renewed intrusions that lead to eruptions (Girona et al. 2014, 2015). There are many unanswered questions about open-vent volcanic systems including: what processes keep these volcanic systems frequently erupting? How does the magma composition remain so homogenous over multiple eruptions? How do we identify evidence of similar open-vent behavior in the geologic past? Is it possible to propose a robust generalized model for how these volcanoes work that would be useful toward identifying a suite of analogue volcanoes (Tierz et al. 2019)?

Mayon volcano (13.257°N, 123.685°E, 2463 m above sea level) is located on Luzon Island in the Philippines, ~ 300 km southeast of the capital, Manila (Fig. 1). It is a famously steep-sided and nearly symmetrical conical volcano, and the most active in the Philippines. Nearly 50 eruptions have occurred since written records began in 1616 (Catane et al. 2003), with the most recent eruption occurring in February 2018. Eruption deposits from ~ 1250 A.D. to 1978 are dominated by andesite to basaltic andesite lavas (Newhall 1977, 1979). Using whole-rock geochemical data, Newhall (1979) highlighted the occurrence of cyclical chemical variations with more primitive magma erupted after higher intensity eruptions. Such variations were explained by periodic basaltic magma injection and subsequent fractional crystallization at generally shallow depths. Mayon lavas erupted after 1979 have either been mainly investigated within the context of Philippine and Bicol arc-scale processes (Castillo & Newhall 2004; McDermott et al. 2005; DuFrane et al. 2006), or are cursory geochemical and petrological reports (Listanco et al. 1999; Catane and Mirabueno 2001; Arpa et al. 2006).

Fig. 1

a Regional location map of Mayon volcano (red circle) and other volcanoes (gray circles) within the Philippine Archipelago. b The hillshade map shows the location of main edifice and the two satellite scoria cones. Ligñon Hill is to the south of the main edifice, whereas Paulog Cinder Cone is to the west of the summit (shown with stars, respectively). inset) The photo shows typical open-vent degassing and symmetrical morphology of the main edifice. Regional location map made with GeoMapApp (www.geomapapp.org)/ Creative Commons license by (Ryan et al. 2009)

Here we present a detailed geochemical and petrological dataset for multiple Mayon eruptions between 1928 and 2009, that mainly produced lava flows. We use these data to reveal the timing of crustal magmatic processes that may drive its eruptions. Deciphering the plumbing system structure, location of magmatic environments, and the timing of magmatic processes at Mayon provides spatial and temporal constraints to help monitoring efforts for better anticipating its eruptions. Furthermore, by investigating plumbing systems over time, and comparing multiple volcanoes, we aim to develop a generalized conceptual model of how persistently degassing and open-vent volcanoes work.

Geologic background

Mayon volcano is situated within the Bicol Volcanic Arc, which forms from the westward oblique subduction of the Philippine Sea Plate into the Philippine Trench (McDermott et al. 2005; Lagmay et al. 2005; Andal et al. 2005). Ultramafic and metamorphic rocks constitute the deep basement, which is overlain by limestone, other marine sediments, and volcaniclastic sediments (Andal et al. 2005; McDermott et al. 2005). The mantle source has been modified by sediment and fluid inputs from the subducting plate as indicated by trace element and isotopic studies (e.g., Castillo & Newhall 2004; McDermott et al. 2005; DuFrane et al. 2006). The oblique convergence also results in the Philippine Fault, a large left-lateral strike slip system with splay faults like the San Vicente-Linao Fault that cuts across the southwestern side of Mayon (Newhall 1979; Lagmay et al. 2005). Volcanism in the Bicol Arc began ~ 70 Ma and continued to 1 Ma (Andal et al. 2005). The current stage of active volcanism resumed less than 0.1 Ma and continues to this day (Andal et al. 2005; Mirabueno et al. 2006). In addition to the Mayon stratocone, there are several scoria cones 11–14 km from the summit vent, two of which have been previously investigated (Maximo, 2011). Ligñon Hill (156 m a.s.l.) is located south of the main edifice and within Legaspi City limits. Mayon Volcano Observatory, managed by the Philippine Institute of Volcanology and Seismology (Philippine Institute of Volcanology and Seismology 2019b), is located at the summit of Ligñon Hill. Paulog cinder cone (120 m a.s.l.) is located east of the main edifice (Maximo 2011). The ages for these scoria cones are not well constrained, but paleosol radiometric ages adjacent to Paulog cinder cone provide a maximum age of ~ 20 ka (Mirabueno et al. 2006).

Eruptive activity in the last 200 years has ranged from minor ash explosions to Strombolian and Vulcanian eruptions (VEI 2–3) (Global Volcanism Program 2013; Philippine Institute of Volcanology and Seismology 2020). In addition to the more frequent mildly explosive eruptions, Mayon has also produced much more explosive events in 1814 and 1897 (Philippine Institute of Volcanology and Seismology 2020). Here we focus on the more common eruptions that have produced many lava flows, although the more explosive and larger events are also an important component of the Mayon dynamics. Common hazards include lava flows, ash fall, pyroclastic flows, and lahars (Catane et al. 1991). Monsoonal rains make lahars a major risk to the nearly 1 million people who live within a 30 km radius (Paguican et al. 2009). The volcano is currently monitored with geophysical (broadband, short period, and infrasound), geodetic (electronic tilt, and global positioning system, GPS), and geochemical (mini-DOAS, multi-gas) instrumentation.

Inter-eruption activity is characterized by persistent outgassing from the summit region, between 400 and 600 tonnes of SO₂ per day, measured via mini-DOAS (Philippine Institute of Volcanology and Seismology 2019a, b; Arelleno et al. 2021). Mayon crater often shows incandescence, indicating that partially crystallized magma is very close to the surface, and probably filling the upper parts of the conduit (Girona et al. 2015). Sustained outgassing requires a relatively high flux of magma and volatiles from deeper parts of the system, a process which has been proposed at other mafic open-vent volcanoes (Etna; Allard 1997). Gas transport through the magmatic system and upper parts of the conduit may occur via magma convection (Shinohara 2008; Witter et al. 2004; Beckett et al. 2014; Sheehan & Barclay 2016), although the high viscosity of the magma implied by the crystal-rich nature of Mayon lavas may require gas percolation through a high permeability crystal-network (Edmonds et al. 2003; Holland et al. 2011) or some combination of magma convection and gas percolation (Burton et al. 2007).

Our understanding of the Mayon plumbing system is partly based on geophysical datasets. Gravity and tilt data suggest magma is stored at 4–5 km below sea level, which corresponds to ~ 6–8 km below the summit (Camacho et al. 2007; Vajda et al. 2012; D. Hidayat personal communication). Volcano-tectonic earthquakes between 2005 and 2016 occurred throughout the crustal column, but cluster between 6 and 12 km (Philippine Institute of Volcanology and Seismology 2020). Finally, preliminary receiver function analysis also suggests a large magma body around 20 km depth (D. Hidayat personal communication; S. Wei personal communication).

Sampling and analytical methods

The limited geographic extent of the tephra deposits from the low VEI eruptions, coupled with significant rain-induced erosion severely limited tephra sampling at Mayon. Moreover, the compositional similarity adds complexity in distinguishing tephras deposited by different eruptions. Thus, we focused mainly on lava flows and bomb samples from nine main edifice eruptions between 1928 and 2009 (Table 1). Whole-rock geochemical data from the 1814 and 1897 eruptions, as well as from the literature are also included for context. For eruptions before and including the one in 1978, we analysed the same samples described in Newhall (1979). Samples from recent eruptions (1984, 1993, 2000, 2006, and 2009) were collected from exposed outcrops on Mayon over several field excursions. Multiple thin sections from each eruption were prepared, and high-resolution images (3200 dpi) of the entire sections were obtained using a flatbed scanner. We used JMicroVision software (Roduit 2007; http://jmicrovision.com/) and counted a minimum of 500 points per thin section to determine phenocryst modal abundances as vesicle-free vol%. Backscattered electron (BSE) images of thin sections and quantitative mineral chemistry were collected on samples from the 1814, and 1928–onward.

Table 1 Summary table of key parameters for Mayon eruptions studied in this paper

Textural definitions, abbreviations, and mineral end members

We classified minerals into phenocrysts (> 300 μm), microphenocrysts (100–300 μm), and microlites (< 100 μm). Mineral chemistry is discussed in terms of end-members. For plagioclase we used Anorthite, An mol% = 100 × Ca/(Ca + Na + K). We discuss orthopyroxene and clinopyroxene in terms of their Mg#, where Mg# (mol%) = 100 × Mg/(Mg + Fe^t) where we computed all iron (Fe^t) as Fe²⁺. We also provide the molar end-members of wollastonite Wo = 100 × Ca/(Ca + Mg + Fe^t), enstatite En = 100 × Mg/(Ca + Mg + Fe^t) and ferrosilite Fs = 100 × Fe/(Ca + Mg + Fe^t) after Morimoto (1989). For olivine we used the Forsterite content, Fo (mol%) = 100 × Mg/(Mg + Fe^t). For Fe-Ti oxides group we calculated the end members between spinel (X_sp MgAl₂O₄), magnetite (X_mag Fe₃O₄), and ulvöspinel (X_ulv Fe²⁺TiO₄).

Whole-rock geochemistry

Whole-rock geochemical analyses of major and trace elements were obtained at ActLabs (Canada) with the 4Lithores analytical package (see www.actlabs.com for details). Samples were powdered and fused using lithium metaborate/tetraborate as a flux and were analyzed with the Perkin Elmer Sciex ELAN Inductively Coupled Plasma-Mass Spectrometer (ICP-MS) system. Major elements were measured with Inductively Coupled Plasma-Optical Emission Spectrometer (ICP-OES) whereas trace elements were measured with ICP-MS. The complete geochemical dataset and details on reproducibility and precision of the analysis can be found in the supplementary Table S1.

Electron images and microprobe analysis

We used the BSE images to identify compositional differences throughout a given crystal and collected compositional data to capture the observed variability. High-resolution BSE images (at least 2048 pixels) were collected on a JEOL-JXA-8530F electron probe microanalyser (EPMA) and a JEOL JSM-7800F scanning electron microscope (SEM), both hosted at Nanyang Technological University, Singapore. Phenocryst compositions from each eruption were analysed via EPMA with a 15 kV accelerating voltage, beam current from 10–30 nA, and beam diameter from point to 5 µm, depending on the analyzed mineral. Core and rim analyses were collected in all minerals if grey scale heterogeneities were observed in BSE images. Core to rim traverses of orthopyroxene and plagioclase crystals were acquired perpendicular to the zoning boundary and with a spacing between individual points ranging between 2.5 and 4 µm. In addition, we also used calibrated BSE images of orthopyroxene to model Fe–Mg compositional profiles (e.g., Ginibre et al. 2002; Martin et al. 2008; Fabbro 2014). Full details for EPMA analytical procedures are provided in the supplementary Table S2.

Secondary ion mass spectrometry

Selected orthopyroxene and clinopyroxene crystals from the 2000 eruption were analyzed for ¹H concentrations. Gold-coated thin sections were loaded into the Cameca ims4f secondary ion mass spectrometer (University of Edinburgh, UK). ¹H, ²⁵ Mg, ³⁰Si, and ²⁷Al were measured using a 15 kV O⁻ ion beam and accelerating voltage of 4500 eV with a -75 eV offset. Beam diameter was about 50 μm for a 4-min pre-sputtering, after which each element peak position was confirmed. The beam diameter was reduced to 8–10 μm for the analysis. Analyses were collected from rims and cores. Each analysis had a minimum of 10 cycles where on peak count times ranged from 2 to 5 s. We normalized the counts based on ³⁰Si and used glass standard-based calibration curves to obtain H⁺ concentrations (for calibration curves and standards, see supplementary Table S3).

Magmatic H₂O content

We measured the H⁺ concentrations in orthopyroxene and clinopyroxene and used equilibrium partitioning relations to calculate the H₂O contents of the melt (Wade et al. 2008; O’Leary et al. 2010). We analysed seven clinopyroxene and 34 orthopyroxene crystals from the 2000 eruption, including crystal cores and rims. To determine the H⁺ partition coefficients, we used Eq. 11 for clinopyroxene, and Eq. 12 for orthopyroxene from O’Leary et al. (2010) and a temperature of 1050 °C. H⁺ diffusion in pyroxenes is relatively fast where proton-polaron diffusion is the dominant mechanism and is affected by the presence of Al³⁺ (Wade et al. 2008; O’Leary et al. 2010; Ferriss et al. 2016; Lloyd et al. 2016). It is likely that the H⁺ in pyroxenes partially re-equilibrated during ascent and cooling of the samples on the earth’s surface (Lloyd et al. 2016) and thus the values we report are minimum melt H₂O concentrations.

Diffusion chronometry modeling

Diffusion modeling of chemical gradients in phenocrysts is becoming a routine method to determine magmatic timescales (Costa et al. 2008). Zoned plagioclase in Mayon lavas are an ideal target for such studies due to their abundance in the lava flows. However, plagioclase have complex textures, with many dissolution zones and small mineral inclusions of pyroxenes, oxides and recrystallized melt inclusions. Such features make it impossible for us to model the diffusion of elements such as Mg in plagioclase (Costa et al. 2003) because the individual analysis from electron microprobe traverses may include multiple phases. As a result, the measured concentration changes cannot be simply approximated by volume diffusion in plagioclase. We have thus focused instead on the Fe–Mg zoning in orthopyroxene, with a total of 239 traverses (187 reverse zoned and 52 normal zoned) which we have modeling using the diffusion equation and standard methodology (e.g., Allan et al. 2013; Fabbro et al. 2017; Flaherty et al. 2018; see also methods). The input parameters included (their specific values are in the results section): temperature, oxygen fugacity to calculate the diffusion coefficient, and the model initial and boundary conditions, where temperature is typically the largest source of uncertainty (Costa et al. 2008). The initial shape of the model was based on the observed chemical profile in the orthopyroxene; these were typically compositionally homogenous cores and rims separated by a step function (Fig. S1a-c). Finally, we assumed an oxygen fugacity at the nickel—nickel oxide buffer which is consistent with arc magmas (Toplis & Carroll 1995). Details about the model assumptions are provided in the supplementary materials.

Results

Whole-rock petrography

Throughout its last 100 years of eruptive history, Mayon has produced remarkably homogenous lavas in terms of mineral assemblage and phenocryst content. All lavas are porphyritic, with 45 to 52 vol% phenocrysts of plagioclase + orthopyroxene + clinopyroxene + Fe-Ti oxide ± olivine (Table 1, Fig. 2a). Plagioclase is the most abundant mineral (75–84 vol% of the total) and is found as phenocrysts, microphenocrysts, microlites, and in glomerocrysts. Pyroxenes (orthopyroxene + clinopyroxene) make up 16–24 vol% of the mineral assemblage. Orthopyroxene (~ 90–95 vol% of the pyroxene population) is found as phenocrysts, microphenocrysts, and microlites, and in glomerocrysts. Clinopyroxene occurs primarily as phenocrysts, microphenocrysts, and in glomerocrysts. Fe-Ti oxides are < 5 vol% of the mineral assemblage and is present as phenocrysts, microphenocrysts, microlites, as inclusions in other minerals, and in the glomerocrysts. Olivine is rare and found as microphenocrysts and included in orthopyroxene. Plagioclase, orthopyroxene, and clinopyroxene show a variety of textures and chemical zoning which are discussed below. Finally, previous works report xenocrystic pargasitic hornblende poikilitically enclosing high-Ca plagioclase (An_93–95) (Newhall 1977) and amphibole phenocrysts in lavas from the 2006 eruption (Maximo 2011).

Fig. 2

a False color X-ray intensity maps combining Ca (teal), Mg (yellow), and Fe (magenta) showing the representative mineralogy and textures at Mayon. Scale bar is 5 mm. b–e Representative backscatter electron (BSE) images of fracture textures in different minerals. Specific fractures highlighted with arrows. All scale bars are 100 μm.  f–h Three different types of glomerocrysts commonly observed in lavas. All scale bars are 100 μm. In all panels—Cpx: clinopyroxene, Opx: orthopyroxene, Plag: plagioclase, Ox: Fe-Ti oxides, Ol: olivine

Healed microfractures are ubiquitous within plagioclase phenocrysts, and commonly observed in orthopyroxene and clinopyroxene phenocrysts (Fig. 2b–e). Many are concentrated near crystal corners, especially when the interfacial angle is acute, where they may intersect crystal twinning or crystal boundaries (Fig. 2b, d). In multi-pyroxene crystal clusters, fractures occur along crystal boundaries, often showing chemical zoning (Fig. 2c). Finally, fractures are observed linking various melt inclusions and the matrix glass, especially in plagioclase crystals (Fig. 2b, e). In many cases, fractures are “in-filled” which provides evidence for post-fracture “healing” at magmatic conditions. These fractures may result from a range of processes including melt inclusion decrepitation, glomerocryst breakup, flow shearing (Bindeman 2005), and brittle fracturing (Taddeucci et al. 2021).

The samples also contain three types of glomerocrysts, or crystal clusters (Fig. 2f–h). The most abundant glomerocryst, Type 1, is composed of orthopyroxene + clinopyroxene + Fe-Ti oxides ± plagioclase ± olivine. The equigranular pyroxenes often host large Fe-Ti oxides inclusions. Type 2 glomerocryst lacks clinopyroxene and is composed of orthopyroxene + olivine + Fe-Ti oxides ± Plagioclase. Textures in Type 2 glomerocrysts have orthopyroxene mantling olivine, commonly with associated symplectite-textured Fe-Ti oxides. Plagioclase may be found as an interstitial phase, and/or a large component of the glomerocryst. Type 3 glomerocrysts are the least abundant and are composed of plagioclase + orthopyroxene + Fe-Ti oxides ± olivine. Plagioclase and orthopyroxene are generally equigranular. Fe-Ti oxides in Type 3 glomerocrysts are hosted in orthopyroxene (Fig. 2h). The abundance of fractures and glomerocrysts is likely a reflection of the high crystal content of the subvolcanic plumbing system of Mayon volcano which behaves as a rigid crystal mush (Holness et al. 2019).

Whole-rock, matrix glass, and melt inclusion compositions

The whole-rock compositions erupted from the main edifice of Mayon from the eruptions that we have studied are high-alumina basaltic andesites (Figs. 3a, 4a) with remarkably homogenous major element contents (see supplementary Table S1 for the complete geochemical dataset). Whole-rock SiO₂ and Mg# show limited variations from 1928–2009, from about 54–55 wt% and 47–49 mol%, respectively (Fig. 3b). These compositions overlap with those reported in previous studies (Castillo and Newhall 2004; McDermott et al. 2005; Andal et al. 2005; DuFrane et al. 2006; Maximo 2011). Previously reported plagioclase-hosted melt inclusions in lavas from the main edifice and matrix glass range from basaltic andesite to dacite with SiO₂ of 53–65 wt% and MgO from 0.3 to 4.6 wt% (Newhall 1977; Mirabueno 2001). The scoria cones have the most primitive whole-rock compositions and olivine-hosted melt inclusions, with SiO₂ ranging from 48 to 52 wt%, and MgO ranging from 5.1 to 9.3 wt% (Maximo 2011).

Fig. 3

a TAS diagram showing previous work (main edifice: light gray symbols, scoria cones: dark gray symbols) and this study (orange symbols). Whole rocks are basalt to basaltic andesite. Melt inclusion and matrix glass are generally more evolved. b Time-series showing that SiO₂ and Mg# have not changed from 1928 to 2009. Previous work from Moore and Melson 1969; Newhall 1977; Newhall 1979; Magalit and Ruelo 1985; Listanco et al. 1999; Jentzsch et al. 2001; Mirabueno 2001; Castillo and Newhall 2004; McDermott et al. 2005; Arpa et al. 2006; DuFrane et al. 2006; Maximo 2011

Fig. 4

Symbols the same as in Fig. 3a–c) Harker diagrams showing Al₂O₃, MgO, CaO. Fractional crystallization MELTS modeling completed with rhyolite-MELTS v1.2 (Ghiorso and Gualda 2015) using a parent composition from the scoria cones (PAC4; Maximo 2011) and varied pressures and with 2 wt% of water and oxygen fugacity at the Ni-NiO reaction buffer (NNO). MELTS modeling results follow the general whole-rock trend for Al₂O₃, MgO, and CaO (with offsets). However, the K₂O model shows that simple fractionation cannot reproduce the whole-rock data of the Mayon edifice, suggesting that the scoria cones are not the dominant parental melts

Major element Harker diagrams of the scoria cones, Mayon whole-rock, matrix glass, and melt inclusions show negative linear trends between some oxides, such SiO₂ vs. MgO, and CaO, (Fig. 4b, c). MELTS (Ghiorso and Gualda 2015) fractional crystallization modeling using the scoria cone basalt as an initial composition shows that MgO, CaO and Al₂O₃ results generally follow the trends outlined by the main edifice lavas. However, other elements such as K₂O have significantly higher contents that do not match the main edifice trends.

Whole-rock trace element compositions show significant scatter (Fig. 5a–d); the scoria cone samples have comparable to much higher Rb, Sr, La, and Zr contents than the basaltic andesite lavas of the main edifice. Such differences further indicate that the scoria cones are not appropriate parental magmas as incompatible elements should increase with crystal fractionation of anhydrous minerals. Given the high crystal content of the Mayon lavas (especially plagioclase) it seems possible that their incompatible trace elements would have been diluted by the addition of crystals (Maximo, 2011). However, the cone samples have comparable to higher Sr and Eu/Eu* ($Eu/{Eu}^{*}= {Eu}_{N}/{\left({Sm}_{N}\times {Gd}_{N}\right)}^{0.5})$, as well as different ratios of incompatible elements (e.g., La vs Zr) than the Mayon lavas we have analyzed, and thus these relations cannot simply be explained by plagioclase accumulation. Therefore, it appears that we have not sampled the dominant mafic parental magmas of the main Mayon basaltic andesite lavas, and that the basaltic magmas that produced the scoria cones originated from different degrees of partial melting, sediment contribution, and fluid contributions (Castillo and Newhall 2004).

Fig. 5

Symbols the same as in Fig. 3a–d) Trace element variation diagrams vs Zr. The scoria cones (symbols circled in b, d) have different Rb, and Sr trends when compared with the whole-rock trace element data; this suggests that scoria cones cannot be the dominant parental magma for the main edifice. Sr and Eu/Eu* data indicate varying degrees of plagioclase fractionation and accumulation. The two trends in La vs Zr indicate varying degrees of fluid input and partial melting produce the magmas erupted from the main edifice

Mineral compositions, zoning, and textures

Although the whole-rock chemistry and mineral assemblage of Mayon lavas erupted in the last century do not vary through time, mineral compositions and textures vary significantly within a single sample. The crystal compositions and chemical zoning are key to unraveling the structure of the plumbing system, processes and time scales that lead to eruptions at Mayon, and thus are described in some detail in the following sections.

Plagioclase

Plagioclase phenocrysts are commonly euhedral to subhedral and exhibit a variety of textures, often within individual crystals, including multiple compositionally distinct zones, patchy and sieve textures, and recrystallized melt inclusions (Fig. 6). Crystal cores can be rounded to angular, but most are subrounded. Plagioclase shows a wide compositional range, from An₂₉ to An₉₄ (Fig. 6), and the compositional histograms show three main modes at about An₆₂, An₇₈, and An₉₀. Crystal rims are consistently An_60–65, whereas the crystal interiors tend to be more calcic (> An₆₅). Most crystal cores are at about An₇₈ rather than An₉₀, with very few crystals of core compositions < An₆₅. The lowest An (< An₅₅) are interstitial plagioclase crystals found within the multi-phase glomerocrysts. An contents do not monotonically increase or decrease from crystals’ core to rim. Rather, many crystals show multiple compositional bands of high An that are 10 s of micrometers thick within lower An areas, with the compositional bands parallel to the crystal margins. Reverse zoned plagioclase crystals (An increasing from core to rim) are rare; we only observed one crystal with a core composition of An₅₅, intermediate zones ranging from An_52–70, and the rim composition of An₆₀ (see Fig. S2).

Fig. 6

Summary of compositions and textures of plagioclase including textures, chemistry and chemical distribution. a, b Compositional distribution by year (a) and then combined into a master dataset (b). Bin width for both is 2 mol% An. Note three peaks in the master set at An₉₀, An_77–78, and An_60–62. c-e BSE images showing representative plagioclase textures. All scale bars are 100 μm. f–h An content from electron microprobe traverses in c–e

Combining textural and compositional observations we further classified the plagioclase into three groups. Group I crystals have a high An (An_80–90) interior with limited to no textural features (Fig. 6c, f) that is mantled by a low An (< 65) rim with a sharp boundary between them. Group II crystals are similar to Group I, except that the high An interior has a boxy or patchy texture filled with low An plagioclase (Fig. 6d, g). Group III crystals are characterized by multiple, alternating An-rich and An-poor bands (Fig. 6e, h). The high-to-low An transitions are often sharp. In contrast, low-to-high An transitions are more progressive, often with uneven dissolution-like textures. These lower An regions with dissolution features also contain melt pockets (recrystallized to various mineral assemblages), multi-mineral assemblages of small (< 20 µm) pyroxenes, and Fe-Ti oxides (e.g., Fig. 6e), and vesicles. Fine-scale oscillatory zoning displaying small variations in An content (< 5 mol%) is often observed in the low An zones (e.g., Fig. 6c). Finally, other minerals such as orthopyroxene, olivine, and Fe-Ti oxides are commonly included in the interior of the plagioclase phenocrysts.

Orthopyroxene

Orthopyroxene phenocrysts are predominantly euhedral to subhedral whereas anhedral orthopyroxene inclusions are often found in plagioclase phenocrysts. Compositions range from Mg# 65 to 83 (with end-members of Wo_0–14, En_56–78, Fs_19–39) and show a normal distribution, when disregarding zoning (Fig. 7a–d). We note that some crystals are pigeonite and likely out of equilibrium with the melt. Reverse and normal zoned orthopyroxene coexist in each sample (Fig. 7c–h); reverse zoned crystals are more common. Based on the population statistics, Mg# varies by 2 mol% between core and rims in the reverse zoned orthopyroxene, whereas Mg# varies by 8 mol% in normal zoned orthopyroxene (Fig. 7c, d). Fine-scale oscillatory zoning and sector zoning in Mg# are also observed (Fig. S3a, b). Orthopyroxene crystals rarely show multiple zones of alternating Mg# (Fig. S3c, d). Textures typically associated with rapid growth (e.g., dendritic/skeletal features at the apices of crystal faces and aligned melt inclusions) are found in normal zoned and smaller crystals (Fig. S3c). Some orthopyroxene exhibit exsolution lamellae in their cores (Fig. S3e, f); these may be indicative of long-term storage at near solidus temperature, and imply the existence of crystal-rich environments. Finally, Fe-Ti oxides and olivine are common mineral inclusions in orthopyroxene. The majority of orthopyroxene cores are subrounded (often preferentially along the longest visible axis) with an oblate morphology, indicative of dissolution. Some reversely zoned orthopyroxene are more complexly rounded and exhibit non-uniform dissolution textures (Fig. S3g-j).

Fig. 7

Summary of compositions and textures of orthopyroxene including textures, chemistry, and chemical distribution. a, b Compositional distribution by year (a) and then combined into a master dataset (b). Bin width for both 2 mol% Mg#. c, d Compositional distribution separated by zoning. Three peaks emerge at Mg# 67–68, 73–74, and 75–76. Note the difference in y-axis between c) and d). The intermediate blue denotes the overlapping distributions between the core and rims. e, g BSE images of representative orthopyroxene textures. f, h Mg# distribution for a reverse and normal zoned crystal, respectively

Clinopyroxene

Clinopyroxene phenocrysts are subhedral to euhedral and often have simple basal twins (Fig. 8a–c). The clinopyroxene composition ranges from Mg# 65 to 89 (Wo_36–45, En_39–69, Fs_13–31) (Fig. 8d). Reverse zoning in Mg# is more common than normal zoning (Fig. 8a, b, e, f). Clinopyroxene crystals with multiple changes in Mg# are rare and have subrounded or resorbed cores. Normal zoned clinopyroxene crystals show changes of about 8 mol% Mg# between cores (at about 82) and rims (about 74; Fig. 8e). Reverse zoned clinopyroxene show a much narrower compositional range of 2 mol% Mg# between rims (about 82) and cores (about 80; Fig. 8f).

Fig. 8

Summary of compositions and textures of clinopyroxene including textures and chemical distributions. a, b BSE images of representative clinopyroxene textures. c Mg# distribution combined from nine eruptions with 2 mol% Mg# bin widths. d, e Distribution separated by zoning. Three peaks emerge at Mg# 81–82, Mg# 79–80, and Mg# 73–74. Note the difference in y-axis between d) and e). The intermediate blue denotes the overlapping distributions between the core and rims

Olivine and Fe-Ti oxide

Rare olivine microphenocrysts are subhedral to anhedral, with a compositional range from Fo₆₄ to Fo₇₉ (Fig. 9a, b, d). They are normal zoned (cores about 78, rims about 68), with zoning lengths ranging from 8 to 40 µm (average of ~ 15 µm). We found one reverse zone olivine crystal with a 2 mol% difference between core (Fo₇₇) and rim (Fo₇₉). Olivine crystals in the Type 2 glomerocrysts are anhedral and unzoned and are often associated with Fe-Ti oxide and mantled by orthopyroxene. We did not observe deformed olivine. However, it has been reported in deposits from the 1814 eruption (Mirabueno 2001) and is often associated with crystal compaction in a magma mush (Holness et al. 2019).

Fig. 9

Summary figure for olivine and Fe-Ti oxides including textures and chemical distribution. a BSE images showing rare olivine microphenocrysts with normal zoning. Although we did not measure diffusion in olivine because of their paucity, the diffusion lengths of the observed profiles ranged from 8 to 19 μm. b Fo distribution from all analysed olivine. Distribution separated by olivine type (zoned—rim and core, and inclusion). Three peaks emerge at Fo_75–76, Fo_73–74, and Fo_67–68. c BSE images of representative Fe-Ti oxides textures including their presence as a phenocryst phase, inclusions, and exsolution lamellae. d BSE image showing the peritectic, symplectite texture of orthopyroxene-olivine-Fe-Ti oxides

Fe-Ti oxide is subhedral or anhedral and commonly included within orthopyroxene and clinopyroxene (Figs. 2a, 7g, 8b, c, 9c), or sometimes in plagioclase. Fe-Ti oxides phenocrysts may exhibit exsolution lamellae (Fig. 9d). We observe abundant Fe-Ti oxides microlites in the groundmass and as inclusions within the glomerocrysts. Those in type 2 glomerocrysts show symplectite textures and are associated with olivine inclusions (see Figs. 2f, 9d). The majority of Fe-Ti oxides are magnetite, ranging in composition from spinel (X_sp 0.08–0.22) to magnetite (X_mag 0.53–0.70) to ulvöspinel (X_ulv 0.13–0.32). Ti–rich magnetite is occasionally observed as exsolution lamellae with compositions of spinel (X_sp = 0.09–0.13), magnetite (X_mag = 0.11–0.27) and ulvöspinel (X_ulv = 0.61–0.74).

Assessing the intrinsic variables of the magmas and plumbing system

Calculating equilibrium melt compositions

We estimated the composition of the melt from which the dominant plagioclase compositions might have grown (An₉₀, An₇₈, An₆₂), using their K concentrations, the plagioclase-melt partition coefficient from Bindeman et al. (1998), and making use of the mainly linear relationships between K and Si concentration shown by the bulk-rocks and various glass compositions (method after Fabbro et al. 2017; Fig. 10a, b). The K crystal-melt partitioning of Bindeman et al (1998) is weakly dependent on temperature, therefore for simplicity we report the results at 1050 °C (see Fig. S5 for all results including the limited effect of a temperature variation from 1000 to 1100 °C). We found that the prominent peak of An₆₂ plagioclase likely formed in andesitic to dacitic melts (Fig. 10), whereas the An₇₈ plagioclase formed in basaltic andesitic to andesitic melts that overlap with the observed whole-rock compositions. Melts in equilibrium with the most calcic plagioclase grew from liquids with lower K and Si than the observed whole-rock compositions, and may originate from basaltic melts that reside below Mayon but which have not erupted, as has been found elsewhere (Izbekov et al. 2002). Alternatively, the high An plagioclase may have originated from a hydrous basaltic andesitic melt.

Fig. 10

a–b Equilibrium melt composition calculated from plagioclase partition coefficients (after Bindeman et al. 1998). a The equilibrium melt composition for K₂O decreases with increasing An content. b Whole rock, matrix glass, and melt inclusion K₂O content vs SiO₂; symbols as in Fig. 3, 4. Note the same y-axes for a) and b) to directly compare the equilibrium melt K₂O and sample K₂O compositions. The dashed line in b) projects the sample trend to basalt compositions to estimate the equilibrium melt associated with An₉₀. Bold vertical lines show the estimated equilibrium melt based on the intercept with the regression line

Calculated melt H₂O contents

H₂O concentration in pyroxenes ranges from 18 to 340 ppm, with lower concentrations in clinopyroxene than in orthopyroxene (Fig. 11; see supplementary Table S8, and Fig. S6). The calculated H₂O in the melt ranges from about 0.15 wt% H₂O to 3.5 wt% H₂O, but most data are < 1 wt% H₂O. These values are consistent with H₂O contents reported for mafic open-vent volcanoes from subduction zones (e.g., Bouvet de Maisonneuve et al. 2012a; Plank et al. 2013; Ruth et al. 2016). Such relatively low values may partially reflect H⁺ loss via diffusion out of the crystal during magma ascent or cooling at the earth’s surface. Low H⁺ contents may also indicate re-equilibration due to prolonged residence in degassed magma within the shallow plumbing system.

Fig. 11

H₂O content in the magma calculated from measured H₂O contents in clinopyroxene and orthopyroxene. Melt H₂O contents calculated from orthopyroxene range from < 0.5 to ~ 3.5 wt%. clinopyroxene-derived melt H₂O contents are more constant at less than 1 wt%. We would expect higher H₂O contents derived from the clinopyroxene (Wade et al. 2008); therefore, we posit these have been re-equilibrated and the H₂O contents are minimum values

The relatively low H₂O contents are consistent with the virtual absence of amphibole phenocrysts in Mayon lavas, although it could still be an important mineral growing at the deeper part of the system (Blatter & Carmichael 1998; Moore & Carmichael 1998).

We also calculated melt H₂O contents using the plagioclase-melt hygrometer (Waters and Lange 2015). For An₆₂ and andesitic melt composition (Newhall 1977, sample 6–1), H₂O contents vary from 1.8 to 3.1 wt% for temperatures from 1100 to 1000 °C, respectively. For An₇₆ and the average whole rock composition (basaltic andesite), H₂O contents vary from 2.5 to 3.8 wt% for temperatures from 1100 to 1000 °C. Finally, for An₉₀ and basaltic melt composition (Newhall 1977, sample 91-LH), H₂O contents vary from 2.9 to 4.2 wt% for temperatures from 1100 to 1000 °C.

Calculated temperatures-geothermometry

The varied mineral assemblage allowed us to take advantage of multiple geothermometers to calculate temperature. Core and rim mineral compositions were combined with whole-rock and matrix glass compositions, respectively; we acknowledge that some ambiguity exists on the liquid compositions, and the effects of water in the melt. We calculated temperatures assuming 1 and 3 wt% H₂O contents in the melt (based on the calculated melt H₂O content from the H⁺ in orthopyroxene and clinopyroxene, and plagioclase-melt hygrometer results) and discarded all results outside the published acceptable equilibrium crystal-melt partitioning (K_D) values (Table 2). Finally, for the clinopyroxene-orthopyroxene thermometer, we measured touching crystal pairs.

Table 2 Summary of the textural relations and crystal-melt pairs we used for various geothermometers

The clinopyroxene-orthopyroxene temperatures of all studied eruptions varied from 940 to 1060 °C, with an average temperature of 1000 ± 20 °C (Fig. 12a). Clinopyroxene rim-matrix melt results (at 3 wt% H₂O) partly overlap with the clinopyroxene-orthopyroxene results, with a peak around 1025 °C. The clinopyroxene core-whole rock melt, orthopyroxene rim-matrix melt, and orthopyroxene core-whole rock thermometers largely overlap, regardless of H₂O content, with peak temperatures of about 1050 to 1100 °C. The plagioclase thermometers have the widest range in temperatures with respect to cores vs. rims, and H₂O content. The An_60–62 plagioclase rim-matrix melt temperature (at 3 wt% H₂O) is slightly less than 1050 °C, whereas An_74–77 and An_>88 plagioclase core-whole rock melt temperatures (at 3 wt% H₂O) are around 1075–1080 °C. Generally, in all minerals the core-whole rock melt temperatures are slightly higher than the rim-matrix melt temperatures. We conclude that Mayon magmas erupted in the last 100 years record a relatively limited range of temperatures, most are 1050 ± 25 °C (Fig. 12b).

Fig. 12

Summary of geothermometry results. a Geothermometry results from individual eruptions calculated with the clinopyroxene-orthopyroxene thermometer in Putirka (2008). The average temperature over the twentieth century was 1002 ± 22 °C. b The temperature range at Mayon calculated using multiple mineral thermometers from Putirka (2008) and separated in core + whole rock (WR), and rim + matrix glass (MG). Equations given in the image. The orthopyroxene-melt, clinopyroxene-melt thermometers are generally within 50 °C. The plagioclase-melt thermometers show a wide range based on water content of the melt

Calculated pressure-geobarometry

We used various approaches to constrain the pressure of crystallization and likely range of magma storage depths below Mayon. We first calculated the volatile melt saturation pressure using the MagmaSatApp and the solubility model of Ghiorso and Gualda (2015). With the average whole rock composition and using only the calculated melt H₂O contents from the H⁺ concentrations in clinopyroxene and orthopyroxene, we obtained minimum H₂O saturation pressures from 0.15 MPa to 154 MPa. We also applied the clinopyroxene-orthopyroxene barometer, clinopyroxene-melt barometer, and Orthopyroxene-melt barometer (using crystal cores only and 3 wt% H₂O) based on the mineral analyses that fulfilled the accepted K_D values by Putirka (2008). We found that these mineral barometers provide information about the deeper plumbing system (Fig. 13). The clinopyroxene-melt and orthopyroxene-melt barometer results range from 150 to > 600 MPa (about 7.5–25 km). The clinopyroxene-melt barometer results show a main peak at 416 MPa (18 km), very similar to those of orthopyroxene-melt (422 MPa, 18 km). The clinopyroxene-orthopyroxene barometer results also ranged from ~ 200 to > 600 MPa (9.5–25 km), with the main peak at 457 MPa (20 km) and a smaller peak at ~ 330 MPa (15 km). The difference between the major element barometers of the pyroxenes, and that calculated from the H⁺ content and water saturation pressure is likely due to H⁺ loss from the pyroxenes during magma ascent or cooling at the earth’s surface.

Fig. 13

Histogram results compiled from multiple geobarometers (Putirka 2008, 2016; Li et al. 2020) and the minimum H₂O saturation pressure (from melt H₂O contents calculated from the orthopyroxene H₂O contents). The dashed gray line shows the depth of the amphibole barometry result (430 MPa). These data show that minerals are sourced from nearly the entire crustal column beneath Mayon

Hornblende stability requires significant H₂O contents in the melt (typically > 3 wt%; Krawczynski et al. 2012), which can be easily achieved at deeper storage zones (Melekhova et al. 2017; Goltz et al. 2020). We used the hornblende composition reported in Newhall (1977) and the geobarometer (Eq. 7a) of Putirka (2016) to obtain a first-order estimate of the depth of hornblende crystallization. We tested several melt compositions from basalt to andesite to determine which was in equilibrium with the amphibole using the K_D (Si-Al) and K_D (Fe–Mg) relationships following Li et al. (2020). We found that amphibole is in equilibrium with a basaltic andesite composition, with 54–57 wt% SiO₂, ~ 18 wt% Al₂O₃, and FeO* from 7 to 16 wt%; these values overlap with reported whole rock composition (Fig. 4a). The geobarometer indicates that hornblende grew at 430 MPa, corresponding to a depth of about 19 km overlapping with the main mode of depths from pyroxene-based barometers.

Assessing magmatic timescales with diffusion chronometry

Magmatic temperature is the most important parameter that affects the values of the diffusion coefficient. We modeled diffusion isothermally with 1030 °C for rims with Mg# < 70, and 1070 °C for rims > Mg# 70 with ± 30 °C for the uncertainty; these values encompass the natural temperature variability in the magmas from the 9 eruptions. Orthopyroxene crystals from all eruptions gave a range of timescales (Fig. 14a; supplementary data) from a few days to about 65 years, although the majority (ca. 88%) of the timescales are less than about 10 years, which is within the average inter-eruptive repose time of Mayon. Timescales from normal and reverse zoned crystals generally overlap, but the shortest timescales (< 60 days) are exclusively from reverse zoned crystals (Fig. 14a). We also observe an almost exponential increase of the number of timescales about three years prior to eruption (Fig. 14a inset). We calculated the injection/mixing date for each crystal by subtracting the diffusion times from the eruption date (Fig. 14b), placing the diffusion timescales in an absolute-time reference frame. For this we assume that all crystals were ejected during the first day of the eruption and that cooling time of the lava flow did not affect the calculated timescale. We find that individual eruptions show the same pattern as the aggregate histogram, with many more timescales as time gets closer to the eruption date. These data also suggest limited crystal recycling between eruptive events and that the processes that create the zoning occur between eruptive episodes. We find that only four eruptions had significant crystal recycling from previous events (Fig. 14c). For these eruptions, 18–23% orthopyroxene crystals were recycled from previous events, with the 2009 eruption having up 45% of recycled crystals.

Fig. 14

a Master dataset compiled from diffusion timescales from all eruptions; note the log scale on x-axis. Uncertainty in the diffusion timescales are shown with the horizontal bars; this is associated mostly with the ± 30 °C uncertainty from the temperature data. The majority of the timescales are from reversely zoned crystals. The shortest timescales are always reversely zoned. The longest may be either reversely or normal zoned. Inset: Histogram of timescale results, 1-year bins. The majority of timescales recorded a mixing event within the average inter-eruption repose period. Note the general uptick in diffusion-inducing events within two years of an eruption at time 0. b Injection/mixing date determined from diffusion chronometry. Increased number of injections occurred within 1–2 years prior to each eruption. c A subset of the timescale dataset showing the mixing year for the recycled crystals from each eruption. ** Mayon erupted in 1938. Although we show it here for reference, we were unable to sample the lava flow

Discussion

The magmatic processes within Mayon’s plumbing system have produced the same whole-rock composition and mineral assemblage over the last hundred years, suggesting a stable structure and a steady-state system. Some of the key processes and observations include the flux of gases through the magmatic system to sustain its permanent volcanic gas plume, the extensive crystallization of mainly plagioclase that is revealed by the crystal-rich content of its lavas, and the repetitive magma intrusions that are required for the frequent eruptions, and which we propose are recorded in the plagioclase and pyroxene crystals. Our understanding of subvolcanic plumbing systems has significantly changed over the last two decades. It is now recognized that magmatic systems are mainly made of magma mushes dominated by crystal networks rather than large pools of melt (Marsh 2006; Holness et al. 2019). The magma mushes may be vertically extensive, forming a network distributed throughout most, if not the entire crust (i.e., the so-called “transcrustal magmatic system”, Cashman et al. 2017). The presence of such crystal-rich zones has been proposed using textural observations of crystals and aggregates (Holness et al. 2019). Moreover, it has also become apparent that the crystal cargo has components recycled from previous eruptions/magmas from a variety of sources and ages, and not necessarily those related to the melt that carried the crystals to the surface (Davidson et al. 2007). The recognition that the plumbing system may be very crystal-rich has sparked the need for new numerical simulations that are capable of exploring the dynamics of intrusions in crystal-network dominated systems (Schleicher et al. 2016; Carrara et al. 2020). The high crystallinity of the system can be due to simply long sustained cooling, although another major conceptual shift has been the recognition that fluxing of volatiles from deeper to shallower parts of the system may be also a very important process in creation of crystal mushes. More specifically, volatile fluxing and percolation of CO₂ from the deeper to the shallower parts of the system can drive crystallization by decreasing the H₂O activity, and has been documented at various volcanoes (e.g., Métrich et al. 2004; Yoshimura and Nakamura 2011; Caricchi et al. 2018).

In this section, we bring together the textural and geochemical data of Mayon’s crystal cargo to reveal the main components of the structure below the volcano, including the processes and time scales that are active during repose and those that may lead to eruption. We integrate our observations from Mayon to those from those of other open-vent volcanoes from the literature and discuss the importance of developing an interpretative framework of mafic, open-vent systems.

The plumbing system and associated magmatic processes

Several observations suggest that the Mayon plumbing system consists of magma mush regions, especially in the shallower parts of the system where most of plagioclase crystallizes (Fig. 15). One is the high phenocryst crystal content of the lavas, which is dominated by plagioclase (> 30 wt%) and the common occurrence and textural variety of glomerocrysts as has also been reported in other systems (e.g., Arculus and Wills 1980; Stamper et al. 2014; Holness et al. 2019; McIntire et al. 2019). Another is the presence of phenocrysts with healed microfractures, likely due to deformation and brittle failure of the magma mush (Bindeman 2005; Holness et al. 2019; Taddeucci et al. 2021), which would not occur in a liquid-dominated system. The orthopyroxene exsolution lamellae and orthopyroxene-olivine-Fe-Ti oxide peritectic textures suggest prolonged storage at near solidus conditions as can be expected in a crystal-rich system. Finally, despite the remarkably uniform whole rock geochemistry, the crystal-rich mushes likely host a wide range of melt compositions, from basaltic andesite to dacite as we have found in the calculated melt (Fig. 10), and the melt inclusion data from the literature (Fig. 3). Similar observations have been reported in other systems (Kilgour et al. 2013; Ruth et al. 2016; Stock et al. 2020).

Fig. 15

Plumbing system schematic based on the results in this paper. Bold black lines show potential magma pathways for a given eruption. Note that pathways may or may not be repeated, thus limiting crystal recycling between eruptions. Mineral chemistry suggests a trans-crustal plumbing system with at least two magma mush zones with a larger proportion of melt. The diffusion/injection timescales and limited amount of crystal recycling suggests that each eruption is derived from a small portion of the system that is activated and mixed during the inter-eruption repose interval. The remaining parts of the system is then in varying states of evolution and cooling, as evidenced by the glomerocrysts, evolved melt inclusions, and near solidus textures in orthopyroxene and Fe-Ti oxides

Moving upward from depth and towards the earth’s surface, we infer the existence of one or two mid to deep crustal regions of magma ponding and crystallization by the presence/absence of some minerals and geothermobarometry calculations (Fig. 13). The available amphibole barometry places this region at ~ 19 km below the crater. This region could be where basalt or hydrous basaltic andesite melts enter the system and crystallize An₉₀ plagioclase to evolve towards high-Al basaltic andesite that Mayon erupts. The high water contents and pressure inhibits significant plagioclase crystallization (Gust and Perfit 1987).

Clinopyroxene, orthopyroxene, and olivine crystallization like occurred at these depths as well. Although the main edifice magmas fractionated from a basalt, the alkalic compositions of the scoria cone basalts indicates that they are not the parental magmas. Our inferences agree with experimental phase equilibria results from mafic subduction zone magmas (Moore and Carmichael 1998; Blatter and Carmichael 1998; Blatter et al. 2013; Melekhova et al. 2017) and could also explain the absence of plagioclase within the Type 1 glomerocrysts. The main orthopyroxene-melt, clinopyroxene-melt, and clinopyroxene-orthopyroxene barometry results indicate magma storage depths of 18–20 km and support the presence of a deep magma mush region at Mayon. The geobarometry results for samples from the deeper plumbing system (> 15 km) generally overlap with the preliminary receiver function dataset (S. Wei, personal communication).

The shallowest melt storage below Mayon extends from about near surface conditions to about 4–8 km depth (Figs. 13, 15) as indicated by our geothermobarometry and volatile saturation pressure calculations. Magma/melt mixing/mingling between the injecting basaltic andesite and andesitic/dacitic magma potentially produces the observed linear whole rock geochemistry trends. An_60–62 plagioclase crystallization from andesitic/dacitic melts and re-equilibration of deeply derived pyroxene likely occurs in this region. Such depths are consistent with textures observed within Type 2 glomerocrysts (orthopyroxene + olivine + Fe-Ti oxides ± plagioclase) which show a peritectic texture, suggesting the reaction where orthopyroxene and Fe-Ti oxides begin to crystallize at the expense of olivine. Such reactions occur within more evolved melts which often develop at shallow depths (Blatter and Carmichael 1998). Moreover, the equilibrium between plagioclase and orthopyroxene in Type 3 glomerocrysts also suggests mid-crustal to shallower conditions (< 6–8 km), where H₂O contents are low enough (< 4 wt%) to allow significant plagioclase crystallization (Moore and Carmichael 1998). Our geobarometry results do not match the locations of inferred gravity anomalies in the shallow system (6–8 km depth) (Camacho et al. 2007; Vajda et al. 2012). This may reflect the uncertainties in the geobarometry method and/or the geophysical approaches. Additional petrological sampling and analysis to refine the shallow system plumbing system may lead to better depth correlations.

Magma injection and mixing are dominant processes within the Mayon system as evidenced by the large textural variety and compositional zoning of the phenocrysts, especially the plagioclase. The presence of normal and reverse zoned clinopyroxene and orthopyroxene crystals in a single hand specimen can be simply explained by recharge of relatively hotter and more Mg-rich magma into a more evolved and cooler reservoir (Streck 2008). The reverse Mg# zoning in orthopyroxene crystals spans about 2 mol%, and mainly consists of a homogenous core surrounded by a rim of a different composition (Fig. 7). However, the cores of many crystals are subrounded and suggest that crystal dissolution may have erased the past evidence for multiple mafic recharge events (Fig. 7e, and Fig. S3). Thus, the orthopyroxene crystals are likely only recording the last magma replenishment event prior to eruption. The replenishment melt was probably hotter and more volatile-rich than the resident melts, although mixing between resident evolved melt and incoming mafic melt means that the crystal rims may record only the mixed melts, rather than the intruding melt. We may infer some information of the resident and intruding magmas using the temperature difference between the different core and rim thermometer results (at the same H₂O content) (Fig. 11) which varies ~ 50 °C from approximately 1050 to 1100 °C. This shows the injecting magma may not have been significantly hotter than the resident magma, and that the difference is mainly related to storage, crystallization, and degassing at shallower conditions. Plagioclase composition and textures also record replenishment of more primitive magma. Using the K concentrations in plagioclase and its melt partition coefficient we found that the An₉₀ plagioclase likely grew from a basalt with 48–52 wt% SiO₂ (or a hydrous basaltic andesite), whereas An₇₈ plagioclase is a proxy for the main injecting magma of basaltic andesite to andesite composition (Fig. 10b). Similarly, we view the An₆₂ plagioclase as a proxy for the shallow resident magma with an andesite to dacite composition. However, the interpretation of An variations is complicated by the major influence of H₂O activity on An and the role of fluid fluxing as we discuss below.

The relatively simple Fe–Mg zoning of the mafic minerals contrasts with complex zoning and textures of many plagioclase phenocrysts (Fig. 6). Plagioclase stability and composition vary as a function of temperature and melt composition, but the H₂O content, and more specifically the H₂O activity has also a very strong effect, in contrast with the much lesser role in the mafic minerals (Sisson and Grove 1993; Waters and Lange 2015). The high abundance of plagioclase phenocrysts and the variety of plagioclase textures and chemistry (Fig. 6) reflect a major role for changes in the abundance and composition of the volatiles in the Mayon plumbing system. Addition of more primitive (higher Ca/Na) and relatively hotter, mafic melt would tend to first partially melt any pre-existing lower An plagioclase crystals producing boxy, sieve textures and patchy zoning followed by relatively higher An plagioclase crystallization (see Fig. 6). We favor this process to explain the large-scale compositional changes in Group II and large-scale oscillatory zoning in Group III plagioclase. Such textures and observed void spaces (an inferred exsolved vapor phase) have been linked to dissolution and rapid regrowth owing to volatile oversaturation in shallow magma storage regions (Izbekov et al. 2002, 2004; Landi et al. 2004; Streck 2008). Moreover, crystals from the deeper system (An_80–90) may be entrained and mixed into the resident magma (Group I plagioclase). The repetition of these large-scale zoning features suggests repeated injection events, perhaps over many years, which is supported by the diffusion timescale results (Fig. 14b).

In contrast, lower An contents may be derived from shallow H₂O degassing, and/or CO₂ flushing. Mayon emits ~ 700 tonnes/day CO₂ (Allard et al. 2017; Aiuppa et al. 2019) and flushing from depth would decrease the H₂O activity in the melt leading to additional crystallization of low An plagioclase (Caricchi et al. 2018). This process could explain the finer-scale oscillatory zoning associated with small An fluctuations; the overall low An content suggests that this type of degassing occurs at shallow pressures. The range of plagioclase compositions and textures observed in Mayon and also recorded at other open-vent volcanoes (Stromboli – Landi et al. 2004; Etna – Viccaro et al. 2010; Llaima – Bouvet de Maisonneuve et al. 2012b) may be the result of a complex interplay between arrivals of fluids and melts from depth with cooling and melt degassing that is needed to sustain the vigorous volcanic gas plumes (Giuffrida et al. 2017). Distinguishing the finer details of how the different processes have affected the An and textural patterns of plagioclase would necessitate sophisticated numerical models of plagioclase growth and dissolution due to changing variables (Wallace and Bergantz 2002), which are beyond the scope of this paper.

Timing of magma intrusions and crystal recycling

Modeling the pyroxene Fe–Mg zoning patterns using diffusion chronometry from multiple eruptions provides unique insights into the timing of sustained and/or repeated magma injection and mixing events. We found that almost 90% of the diffusion timescales show evidence for magma replenishment/mixing of less than about 11 years (Fig. 14), which is the average repose interval of Mayon. Mixing and injection generally occurred in the deeper plumbing system (~ 15–20 km), based on the orthopyroxene-melt barometric results. This suggests that plumbing system is almost continuously replenished between eruptions by new melt from depths greater than 18–20 km. Such additions of mass, heat, and fluids likely allows for the relatively narrow range of temperature and continues degassing which in turn maintains the open plumbing system. Moreover, we also found that all eruptions were preceded by increased injection frequency and/or mixing on the order of 1–2 years before eruption onset, which suggest that eruptions were preceded by an increased amount of melt replenishment that could have been driven by the accumulated amount of gas lost (Girona et al. 2015).

Recalculation of the time scales into absolute magma intrusion times shows that between 10 and 20% of the pyroxene crystals found in each eruption are recycled from previous events. For example, one crystal erupted in 2009 is inherited from processes that started in 1989 (Fig. 14b, c, see supplementary Table S9 for exact mixing dates). Given the crystal-rich nature of the Mayon plumbing system one could expect a larger amount of recycled pyroxene crystals, although most of the crystal mush is likely made by plagioclase. There are several possible scenarios to explain the small number of crystals with diffusion time scales that are longer than the repose times. One is that older pyroxenes may be dissolved by each renewed intrusion, as orthopyroxene dissolution is relatively fast if we assume that it is mainly dependent on the Mg diffusion in the liquid as it is for olivine (Chen and Zhang 2008). This scenario would require significantly hotter mafic magma input, which is not supported by the thermometry data. Another possibility is that the relationship between the crystal-mush and the intruding melt is very heterogenous and that each intrusion-eruption “couple” only samples a localized, small volume part of the system (Fig. 15). Recent numerical modeling simulations of melt intrusion into crystal-rich mafic mushes highlights that crystals entrained during an injection event may be redistributed to distal regions of the active crystal mush and these crystals may only be entrained during a later injection event (Schleicher et al. 2016 and their supplementary data). Thus, the long-lived magma mush “traps” these crystals, allowing for decades of diffusion before being entrained during a later injection event and erupted. Finally, Mayon eruptions between 1928 and 2009 have been relatively small to moderate (VEI 2–3) and it is most likely that the majority of the intruded magma volume is erupted, so that only a limited number of crystals are left over in the plumbing system. This implies that the magma mush could have returned to the pre-eruptive state relatively quickly, in the few years during the repose period during which addition of CO₂ rich fluids to the shallower parts of the system followed by vigorous degassing towards the surface promotes extensive plagioclase crystallization.

Comparing the plumbing system of Mayon to other persistently active volcanoes

Our conceptual model for the Mayon plumbing system is that magma extends from near surface conditions to depths of about 30 km, although there is not necessarily significant melt at all storage depths (Figs. 13, 15). Given that Mayon is an open-volcano with a permanent volcanic gas plume, the various parts of the magmatic system are well connected (Girona et al. 2015), and allow CO₂-rich fluids to potentially percolate from the deeper to the shallower parts of the system. Such volatile influx could lead to additional plagioclase crystallization as it decreases the activity of H₂O in the melt and promotes additional degassing of H₂O, CO₂ and SO₂ towards the atmosphere (Giuffrida et al. 2017; Caricchi et al. 2018). During repose periods such volatile fluxing and some melt additions leads to the complex plagioclase textures, with numerous dissolution zones followed by crystallization. Pyroxenes may also dissolve and regrow more mafic rims, producing reverse zoning and at some point, due to an increasing amount of melt replenishment the pressure is large enough that the volcano erupts. The eruptions are low-moderate intensity (VEI 2–3), and predominantly involve the conduit and upper parts of the system and bring to the surface a range of crystal glomerocrysts. The eruption frequency and crystal content could be controlled by the amount of magma (and volatile) transfer through the plumbing system (Locke et al. 2003; Girona et al. 2015; Caricchi et al. 2018; McCarthy et al. 2020). The larger and more explosive eruptions that also occur in Mayon such as those in 1814 and 1897 are likely affected by additional processes and thus would necessitate a different conceptual model (for example Merapi 2006 vs 2010, Costa et al. 2013; Rose et al. 2013).

Other studies of open-vent volcanoes such as Stromboli, Etna or Llaima, among others found that their plumbing system often consists of multiple, vertically distributed magmatic environments or reservoirs, often extending to mid-crustal depths, if not deeper (Métrich et al. 2010; Bouvet de Maisonneuve et al. 2012a; Ruth et al. 2016). Smaller, lower intensity eruptions mostly involve magma from the shallow reservoir, but more primitive, deeply derived magmas may also be found in the most powerful eruptions (Métrich et al. 2021). At Mayon, we have not found evidence for these more primitive magmas, but they might be present in the more powerful sub-Plinian to Plinian eruptions such as in 1814. Many of these open-vent systems have high eruption frequencies on the order of days to weeks (Stromboli), once every 1–2 years (Etna), every 10 years (Mayon), and 10–50 years (Llaima) (Table 3). Moreover, these persistently active volcanoes erupt often with highly crystalline eruptive products (plagioclase dominant, 10–60 vol%) and generally homogenous whole rock major element composition over time; although some evolutionary trends have also been documented (e.g., Stromboli—Pompilio et al. 2012). Lavas from these open-vent volcanoes exhibit mineral textures (especially plagioclase) that are like those observed at Mayon (Landi et al. 2004; Viccaro et al. 2010). Diffusion chronometry results in orthopyroxene, clinopyroxene, olivine and plagioclase from many of these systems show that magma injection and mixing consistently occurs months to a few years before a given eruption, within the repose interval for these volcanoes, and often with increasing mixing within the last year to days before a given eruption (Kilgour et al. 2014; Kahl et al. 2015; Bouvet de Maisonneuve et al. 2016; Ruth et al. 2018; Oeser et al. 2018; Métrich et al. 2021). Moreover, for the systems discussed here, crystals with the longest residence times are erupted from volcanoes with the longest repose interval, which again illustrates that persistently active volcanoes are steady state systems, where volume inputs and outputs into a magma reservoir are largely balanced (e.g.,Armienti et al. 2007; Kelly et al. 2008; Firth et al. 2014).

Table 3 Comparison of some volcanological and petrological attributes between several persistently active mafic volcanoes

We propose that the variety of attributes from open-vent volcanic systems such as eruption frequency, degree of magmatic differentiation, amount of degassing, and crystallinity are mainly controlled by the amount of magma and volatile added from depth through the whole plumbing system. In addition to these simple range of variables, the regional tectonic or the stresses created by the volcanic edifice itself also modulate the amount of melt or magma intrusion that is required for eruption (Fig. 16, Table 3). Thus, volcanoes such as Etna and Stromboli with higher eruptive frequencies are perhaps driven by higher melt and volatile fluxes from depth, which lead to more primitive magma compositions, shorter residence times, and variable lava crystallinities. Decreasing the amount of melt and volatile transfer leads to lower eruption frequencies, more evolved rock compositions (high-Al basaltic andesites), larger amount of shallow crystallization, and more effusive eruption of crystal-rich lavas.

Fig. 16

Schematic cartoon compiling the key characteristics observed at several open-vent, persistently active volcanoes. The x-axis in the mineral residence time vs number of crystals is not linked to specific data, but rather shows the relative abundance of long vs short time scales between, in this case, Stromboli and Mayon. E = t0, or the eruption day. Table 3 provides the data and references supporting this figure

Can we identify persistently active volcanic systems in the geologic record?

The definition of open-vent volcanoes depends on our direct observations of their prominent degassing plumes and thus is limited to the historical record. It seems unlikely that this behavior persists for their entire volcanic lifecycle, and it could be important to identify whether a given volcano that is now not active or extinct used to be an open-vent. It appears that some petrological and textural features of the minerals could be the “hall-mark” of an open-vent system. The textures of Group III plagioclase (e.g., Fig. 6e), in particular, the combination of oscillatory zoning with large An jumps and boxy/sieve textures, show some potential as a fingerprint for persistently degassing behavior; a similar concept was suggested by Giacomoni et al. (2014). Plagioclase textures similar to those discussed here have been observed at other persistently degassing volcanoes in arcs such as Llaima (Bouvet de Maisonneuve et al. 2012b their Fig. 7c–d), Stromboli (Landi et al. 2004 their Fig. 3a–f; Armienti et al. 2007 their Fig. 2a–c; Viccaro et al. 2021 their Fig. 3), Bezymianny (Shcherbakov et al. 2011 their Fig. 4b), Etna (Viccaro et al. 2010 their Fig. 3e, 5a, 6d; Giacomoni et al. 2014 their Fig. 9c). With the advent of machine learning approaches for textural analysis (Cheng et al. 2017), a quantitative textural study comparing plagioclase textures between open and closed vent systems would shed more light. Nevertheless, we speculate that the combination of multiple, large-scale oscillatory zoning and boxy/sieve textures that seems unique to persistently active volcanoes in arc settings could be useful to identify volcanoes that experienced open-vent behavior in the geologic past.

Conclusions

The frequent eruptive activity at Mayon volcano serves as a natural laboratory that allows us to investigate the plumbing system dynamics of open-vent, persistently active volcanoes. For the past 100 years Mayon has produced similar mild-intensity eruptions, and remarkably homogenous high-crystallinity basaltic andesite lavas. However, textural observations, phenocryst mineral chemistry, and Fe–Mg diffusion chronometry in orthopyroxene all reveal complex magmatic processes (injection, mixing, H₂O and likely CO₂ fluxing and degassing) throughout the plumbing system and over different time scales. Our detailed observations suggest a vertically extensive plumbing system dominated by magma mushes in the shallow and intermediate crust. Steady state magma injection and mixing provides the heat and gases that maintain the persistently active behavior observed at the summit.

A comparison of persistently active mafic volcanoes reveal they share some key characteristics: a plumbing system that extends to mid-crustal depths and often deeper, whole rock chemical homogeneity over decades, highly crystalline lavas, and a timing of magma intrusion that is consistent with their eruption frequencies. We propose that the range of eruption frequencies, crystallinities, and volatile fluxes are driven by differences in the amount of magma and volatiles supplied at depth and percolating through the plumbing system. Given that a volcano’s eruptive behavior may evolve over its lifespan, an important question is whether persistently active volcanoes could be identified in the geologic record. We propose that the combination of multiply zoned and sieve-textured plagioclase crystals observed in many persistently active systems today could be an indicator of persistent activity at volcanoes in the geologic past.