Petrographic and Geochemical Evidence for a Complex Magmatic Plumbing System beneath Bagana Volcano, Papua New Guinea

Abstract

Bagana is a persistently active stratovolcano located on Bougainville Island, Papua New Guinea. Characteristic activity consists of prolonged lava effusion over months to years, with occasional shifts to explosive vulcanian or sub Plinian eruptions that threaten surrounding communities. Satellite observations have shown that Bagana is a major SO₂ emitter, particularly during eruptive intervals. Despite persistent and potentially hazardous activity, no previous geophysical, petrological, or geochemical studies have constrained the magma storage conditions and reservoir processes at Bagana. To address this knowledge gap, we present new bulk rock major, trace element, and radiogenic isotope data, plus mineral phase major element compositions, for Bagana lavas erupted in 2005 and 2012 and ash erupted in 2016. We use our new data to understand the magmatic processes controlling the typical effusive activity and provide the first estimates of magma storage conditions beneath Bagana. The basaltic andesite bulk rock compositions (56–58 wt% SiO₂) of our Bagana lavas reflect accumulation of a plagioclase + clinopyroxene + amphibole + magnetite + orthopyroxene crystal cargo by andesitic-dacitic (57–66 wt% SiO₂) carrier melts. Constraints from clinopyroxene and amphibole thermobarometry, amphibole hygrometry, and experimental petrology suggest that the high-An plagioclase + clinopyroxene + amphibole + magnetite assemblage crystallizes from basaltic-basaltic andesite parental magmas with >4 wt% H₂O, over a temperature interval of ~1100–900°C, at pressures of ~130–570 MPa, corresponding to ~5–21 km depth. Continued crystallization in the magma storage region at ~5–21 km depth produces andesitic to dacitic residual melts, which segregate and ascend towards the surface. These ascending melts entrain a diverse crystal cargo through interaction with melt-rich and mushy magma bodies. Degassing of carrier melts during ascent results in crystallization of low-An plagioclase and the formation of amphibole breakdown rims. The radiogenic isotope and trace element compositions of Bagana lavas suggest that parental magmas feeding the system derive from an enriched mantle source modified by both slab fluids and subducted sediments. Our findings suggest that the prolonged lava effusion and persistently high gas emissions that characterise Bagana’s activity in recent decades are sustained by a steady state regime of near-continuous ascent and degassing of magmas from the crustal plumbing system. Our characterisation of the Bagana magmatic plumbing system during effusive activity provides a valuable framework for interpreting ongoing monitoring data, and for identifying any differences in magmatic processes during any future shift to explosive activity.

Bagana volcano, geochemistry, igneous petrology, magmatic plumbing systems, subduction zone magmatism

INTRODUCTION

Bagana is a stratovolcano located on Bougainville Island in southeastern Papua New Guinea and is among the youngest and most active volcanoes in Melanesia (Bultitude, 1976; Global Volcanism Program, 2023a). The volcano occupies a remote location in the densely forested interior of Bougainville and recent studies have drawn heavily on satellite or airborne observations (McGonigle et al., 2004; McCormick et al., 2012; Wadge et al., 2012, 2018; McCormick Kilbride et al., 2019, 2023). Bagana’s characteristic activity is the persistent, pulsatory effusion of blocky lava flows that spill from a central summit crater to form compound flow fields on the volcano’s flanks (Wadge et al., 2012). Inter-comparison of digital elevation models constructed over various intervals within the past 70 years suggest that the Bagana edifice is both young and fast-growing. Based on a mean extrusion rate of ~1.0 m³ s⁻¹, which is inferred to be representative of the long-term extrusion rate over the lifetime of Bagana, the whole edifice could have grown in only 300–500 years (Wadge et al., 2018). Spectroscopic remote sensing measurements have shown that Bagana is a major source of sulfur dioxide (SO₂) emissions to the atmosphere, globally the third largest volcanic source in 2005–15 (McGonigle et al., 2004; McCormick et al., 2012; Carn et al., 2017). Outgassing is highest (≥10³–10⁴ t d⁻¹ SO₂) during eruptive intervals and falls to more modest levels (≥10² t d⁻¹ SO₂) during inter-eruptive pauses (McCormick Kilbride et al., 2019, 2023). Bagana also exhibits explosive behaviour, from minor ash venting to infrequent vulcanian or sub Plinian eruptions. In July 2023, a sequence of explosive eruptions, accompanied by pyroclastic flows and generating ash clouds reaching 18 km altitude, resulted in widespread evacuations from several nearby villages (Global Volcanism Program, 2023b, 2023c).

Consistent seismic monitoring was conducted for Bagana volcano between 1980 and 1990 by Rabaul Volcanological Observatory (RVO) staff, Papua New Guinea, with a single seismic station located several kilometres from the summit crater on the southwest side of the volcano. Monitoring ceased in 1990 during civil unrest on Bougainville and was restored later for a few years until around 2013. Seismic monitoring has now been re-established following the July 2023 explosions. During the ~14 years of monitoring, seismicity consisted of high frequency volcano-tectonic earthquakes, low frequency volcanic earthquakes, volcanic tremors, and rockfall events. Rockfall events, generated by lava fragments detaching from the front of active lava flows, were the dominant source of seismicity and occurred more frequently than high frequency volcano-tectonic earthquakes.

These growing observational datasets are not yet supported by a conceptual understanding of the processes controlling Bagana’s activity, specifically: its putative youth; the multi-decadal persistence and dominantly effusive mode of eruption; the causes and timing of transitions to explosive behaviour; or the origin and abundance of magmatic volatiles and how these sustain the long-lived summit outgassing. The Bagana ‘lava cone’ edifice has been suggested to represent the first stages of growth of a new polygenetic arc volcano that may or may not share a plumbing system with the nearby pyroclastic shield volcano Billy Mitchell (Rogerson, 1989; Wadge et al., 2018). However, geophysical evidence for connectivity between these two volcanic systems has yet to be established and satellite-based geodetic studies are hampered by high atmospheric water vapour and dense forests (Wadge et al., 2018). Previous petrographic and geochemical analyses of Bagana lavas are limited to whole-rock major oxide and trace element data, modal mineral abundances, and a single ⁸⁷Sr/⁸⁶Sr ratio (Blake & Miezitis, 1967; Taylor et al., 1969; Page & Johnson, 1974; Bultitude et al., 1978; Rogerson, 1989). To our knowledge, no petrographic or geochemical studies of Bagana lavas have been published since Rogerson (1989). Since then, the accuracy, precision, and detection limits of analytical techniques have improved (e.g. Weyer et al., 2002; Rinaldi & Llovet, 2015), new thermodynamic models have been developed to investigate magma storage conditions (e.g. Putirka, 2008), and conceptual models of magma plumbing systems have evolved (e.g. Cashman et al., 2017; Edmonds et al., 2019). A study using modern analytical and modelling techniques offers an enticing opportunity to advance our understanding of the Bagana magmatic system.

In this contribution, we present a petrographic and geochemical study of recent Bagana lavas to identify magmatic processes and determine magma storage conditions within the subvolcanic plumbing system. We present new bulk-rock major, trace element and radiogenic isotope data for Bagana lavas erupted in 2005 and 2012, and ash erupted in 2016. We present the first major element compositions of mineral phases in Bagana lavas. We use mineral chemistry, textural observations and thermobarometric calculations to investigate sub-volcanic processes and provide the first estimates of magma storage conditions (pressure, temperature, melt H₂O content) for the Bagana plumbing system. We find that Bagana lava bulk compositions reflect the accumulation of a mafic crystal cargo by andesitic-dacitic carrier melts. We suggest that persistent effusive activity at Bagana is fed by near-continuous ascent of magmas from a crustal storage region at ~5–21 km depth, accompanied by degassing. Finally, we use our radiogenic isotope and trace element data to suggest that recycled subducting slab components contribute to the chemical signature of parental magmas that supply the Bagana magmatic system. Our characterization of the Bagana plumbing system provides a valuable basis for interpreting ongoing and future monitoring data, investigating the processes controlling transitions from long-lived effusive behaviour to intense explosive eruptions, and exploring the potential link between Bagana and Billy Mitchell.

GEOLOGICAL SETTING

Bagana volcano is located on Bougainville, the largest island in the Solomon archipelago and an autonomous region of Papua New Guinea (Fig. 1). Bougainville is aligned northeast-southwest along the southwestern margin of the Pacific Ocean and is composed mostly of Oligocene to Quaternary volcanic rocks, derived sediments, and minor Miocene and early Pleistocene limestones (Fig. 1, Blake & Miezitis, 1967). The axial spine of Bougainville is a chain of 17 post-Miocene volcanoes, grouped in 10 main centres (Fig. 1, Blake & Miezitis, 1967; Blake, 1968). Bagana is the only volcano of this chain with a confirmed eruption in the past century (Bultitude et al., 1978).

Fig. 1

a) Regional tectonic context of Bougainville, redrawn from Holm et al. (2016). Bagana is shown by the red star. AP = Adelbert Plate, NB = New Britain, NI = New Ireland, Bo. = Bougainville; PP = Papuan Peninsula. Active convergent margins are marked by black filled triangles; inactive convergent margins by open triangles. Topography and bathymetry are from the ETOPO1 Global Relief Model (https://www.ngdc.noaa.gov/mgg/global/global.html). b) Geological map of Bougainville, redrawn from Blake & Miezitis (1967); Blake (1968). Bagana is shown by the red triangle (location 6).

Quaternary volcanism on Bougainville is the consequence of north-eastward subduction of the Solomon Sea microplate beneath the Pacific plate. The older volcanic rocks of Bougainville form part of the Melanesian arc related to the southwestward subduction of the Pacific plate beneath the Australian plate, which ceased at ~26–20 Ma due to docking of the Ontong Java plateau (Petterson et al., 1999; Holm et al., 2016). The ocean trench southwest of Bougainville reaches 8000 m in depth. Total convergence between the Australian and Pacific plates, variously apportioned across the dividing microplate boundaries, is 10–11 cm yr⁻¹ (Holm et al., 2016). The Bougainville arc front volcanoes are roughly 130 km from the trench, and the downgoing slab dips steeply, up to 76°, beneath the island (Syracuse & Abers, 2006).

The Bagana edifice (summit elevation 1850 m) is dominantly built from overlapping lava flows with minor intercalated pyroclastic deposits. Lava flows from Bagana are typically blocky, slow moving, and can be continuously extruded over periods of months to years (Bultitude, 1976; Bultitude & Cooke, 1981; Wadge et al., 2018). Bagana has been near-continuously active since at least 1947, before which activity is poorly documented (Bultitude, 1976; Wadge et al., 2018). The last major eruption lasted from 2000 to 2014, involving the emplacement of three major lava flow fields, initially to the southwest, then the west, and then to the northeast (Wadge et al., 2018). Explosive eruptions generating widespread ashfall occurred in 2007, 2014, and 2023. The July 2023 eruptions generated pyroclastic density currents and had a substantial impact on surrounding communities, including evacuations, ashfall on local agricultural land, and the contamination of water supplies (Global Volcanism Program, 2023b, 2023c).

SAMPLES

We sampled Bagana lavas in September 2016. Our sample sites (Fig. 2) are at the distal ends of two major lava flow fields on the lower western slopes of the edifice whose ages are given as 2000–08 and 2010–11 by Wadge et al. (2018). The lava flows are blocky and unstable, making access to higher points on the edifice dangerous. We assess the most likely eruption dates of our samples to be 2005 and 2012, based on sampling locations, satellite observations (Wadge et al., 2018), and discussions with residents of the villages nearest to Bagana’s western flanks. We have four lava samples from each of 2005 and 2012. During fieldwork in 2016, we witnessed small ash venting eruptions from a flank vent on the volcano, and we collected samples by brushing fresh deposits from sago palm leaves in the downwind direction beneath the plume (distance to summit ∼2 km).

Fig. 2

a) Aerial photograph of Bagana, viewed from the northwest, in September 2019 (McCormick Kilbride et al., 2023). The 2010–12 and 2000–08 lava flow fields are indicated. b) Field photograph showing the tip of the southern lobe of the 2010–12 lava flow field, author L.C. Salem for scale. c) Satellite image courtesy of Bing Images labelled to highlight three adjacent volcanic edifices and the 2010–12 and 2000–08 lava flow fields. Stars indicate our lava sampling locations.

ANALYTICAL METHODS

Whole-rock major and trace elements

We measured whole-rock major and trace element compositions at ETH Zürich. For XRF analyses, we crushed and powdered 1 g material per sample, then fused rock powders into glass beads with lithium tetraborate-metaborate flux, following the setup of Troch et al. (2018). We analysed major elements using a PANalytical AXIOS wavelength dispersive X-ray fluorescence spectrometer. Relative errors are <1% for the majority of major elements and < 2% for Na₂O. We measured trace element compositions on broken internal surfaces of the same glass beads analysed by XRF. We analysed trace elements using an Excimer 193 nm (ArF) GeoLas laser system coupled to a Perkin Elmer Nexion 2000 DRC quadrupole ICP-MS. Measurements were performed using a laser repetition rate of 10 Hz, laser fluence of 10 J cm⁻², and a He-Ar carrier gas (flow rate ~1 L min⁻¹). We used a scheme of 30 s gas blank recording followed by 40 s target ablation, following the primary standard-flux blank-secondary standard-sample bracketing method of Cortes-Calderon et al. (2022). We analysed primary standard NIST SRM610 glass (Jochum et al., 2011) to test for calibration and instrumental drift, using a 40 μm spot size. We analysed flux blank using a 100 μm spot size to correct for matrix effects on the glass beads. We assessed the precision and accuracy of trace element analyses via repeat measurements of USGS standard reference material BCR-2G (Supplementary File S1), using a 100 μm spot size. We measured the majority of elements to within 15% of preferred values (from the GeoREM database, Jochum et al., 2005) excluding Pb (21%). Reproducibility of the analyses varied from <1–11% relative depending on the element. We made three individual measurements per sample, using a 100 μm spot size and averaged the results. We carried out data reduction using the MATLAB-based SILLS program (Guillong et al., 2008). Raw counts were converted to concentrations using SiO₂ concentrations measured by XRF as an internal standard.

Whole-rock radiogenic isotopes

We measured radiogenic isotope compositions at New Mexico State University, USA. Whole-rock powders (200–400 mg of material, subsampled from bulk rock) were dissolved using HF, HNO₃, and HCl. Strontium was purified using cation exchange resin and 2.5 M HCl. Strontium was loaded onto pre-outgassed and clean rhenium filaments with phosphoric acid and tantalum oxide. We analysed Sr isotopes using Thermal Ionization Mass Spectrometry (TIMS) and five Faraday collectors in dynamic mode and ⁸⁸Sr = 3.0 V. Sr isotopes were normalized to ⁸⁶Sr/⁸⁸Sr = 0.1194 and corrected for any Rb present during the analysis.

Rare earth elements (REEs) were purified using cation exchange resin and 6.0 M HCl. Neodymium was then separated from the remaining REEs using HDEHP resin and 0.25 M HCl. Neodymium was loaded onto pre-outgassed and clean rhenium triple filaments. We analysed Nd isotopes using a Thermo Scientific Neptune MC-ICP-MS using seven Faraday collectors in dynamic mode with ¹⁴⁴Nd = 0.5 V to 3.0 V. Nd isotopes were normalized to ¹⁴⁶Nd/¹⁴⁴Nd = 0.7219 and corrected for any Sm present during the analysis. Lead was purified using anion exchange resin and 1.0 M HBr. Lead was dissolved in 1 ml of 2% HNO₃ and Tl was introduced into the sample to obtain a Pb/Tl ratio between 2–5 (Wolff & Ramos, 2003). We analysed Pb isotopes using a Thermo Scientific Neptune MC-ICP-MS in static mode using five Faraday collectors. Pb isotope results were normalized to ²⁰³Tl/²⁰⁵Tl = 0.41892 and corrected for Hg during analysis. We obtained the following values for secondary standards measured alongside samples: ⁸⁷Sr/⁸⁶Sr = 0.710271 ± 0.000035 (2σ, n = 2) for NBS987, ¹⁴³Nd/¹⁴⁴Nd = 0.511872 ± 0.000011 (2σ, n = 2) for La Jolla, and ²⁰⁶Pb/²⁰⁴Pb = 16.931 ± 0.001, ²⁰⁷Pb/²⁰⁴Pb = 15.484 ± 0.003, and ²⁰⁸Pb/²⁰⁴Pb = 36.674 ± 0.006 (2σ, n = 2) for NBS 981.

Whole-rock oxygen isotopes

We measured oxygen isotopic compositions at the University of Cape Town, South Africa, using a Thermo DeltaXP mass spectrometer. We followed the sample preparation and analytical procedures of Heap et al. (2022). We dried 10 mg aliquots of bulk rock powders overnight at 50 °C, and then in a nickel reaction vessel under vacuum. The dried powders were reacted with CIF₃ to extract the oxygen, which was converted to CO₂ by passing it over a high-temperature platinised carbon rod. Duplicates of the internal quartz standard MQ were run with unknowns and used to calibrate the raw data to the Standard Mean Ocean Water scale, using a δ¹⁸O value of 10.1‰ for MQ (calibrated against the quartz standard NBS-28 (δ¹⁸O = 9.64‰)). Long-term repeated analyses of MQ indicate that analytical uncertainties are approximately ±0.20‰ (2σ).

QEMSCAN images and mineral modal abundance estimates

We produced QEMSCAN images using a Quanta 650F field emission gun scanning electron microscope, equipped with two Bruker XFlash 6130 energy-dispersive X-ray spectrometers (EDS) at the Department of Earth Sciences, University of Cambridge. Using a predefined backscattered electron (BSE) threshold (calibrated against quartz, gold, and copper standards) together with EDS spectra generated per pixel, each pixel is then binned according to the species identification protocol files in the iDiscover software package. The analysis was performed at 25 kV and 10 nA, with a resolution (i.e. pixel size) of 5 μm². The system is set to count 2000 X-ray counts per pixel to aid fast analysis. We estimated mineral modal abundances via point counting (results in Supplementary File S1). We performed point counting by overlaying a grid of 1000 points over QEMSCAN images, which cover a ~13 × 13 mm area of each analysed thin section.

Electron probe microanalysis

We obtained mineral phase major element compositions using electron probe microanalysis (EPMA) on the same polished thin sections imaged using QEMSCAN. We carried out EPMA analyses using a Cameca SX-100 electron microprobe at the Department of Earth Sciences, University of Cambridge. We analysed all mineral phases using an accelerating voltage of 15 kV. We used a spot size of 5 μm for all phases except oxides (1 μm). We used a beam current of 10 nA for plagioclase and amphibole, 20 nA for pyroxenes and 40 nA for oxides. On-peak and background counting times for each element analysed were varied by mineral phase and are listed in Supplementary File S1. The following standards were used for calibration; jadeite (Na), diopside (Si, Ca), orthoclase (K), rutile (Ti), fayalite (Fe), Mn (Mn), St Johns olivine (Mg), corundum (Al), Cr₂O₃ (Cr).

Complete data for unknowns and secondary standards are provided in Supplementary File S1, along with 1σ precisions from counting statistics. We monitored accuracy by analysing the following secondary standards from the Smithsonian National Museum of Natural History: anorthite, Kakanui augite, Kakanui hornblende, plagioclase, chromium augite, diopside (Jarosewich et al., 1980, 1987). We typically measured major element (SiO₂, Al₂O₃, MgO, CaO, FeO in pyroxene and amphibole) concentrations within 4% of preferred values (Jarosewich et al., 1980, 1987) and minor elements (TiO₂, Na₂O, MnO, K₂O in amphibole, FeO in plagioclase) within 1–23%. Uncertainty on measurements of unknowns shows that we typically measured major element oxides with 1σ precisions better than 4%. Minor element oxides were measured with more variable precision, dependent on concentration. For example, Na₂O was measured with 1σ precision better than 7% in plagioclase and amphibole, and 15% in clinopyroxene.

We carried out additional EPMA analyses on different crystals of plagioclase, clinopyroxene and amphibole using a Cameca SX-100 electron microprobe in the Department of Earth and Environmental Sciences, University of Manchester. These analyses primarily focused on obtaining high-precision data for minor elements in clinopyroxene that have a strong influence on the precision of thermobarometric calculations (Wieser et al., 2023b). We analysed all mineral phases using an accelerating voltage of 15 kV. We used a spot size of 1 μm for clinopyroxene, 5 μm for amphibole and a defocussed 10 μm spot for plagioclase. On-peak and background counting times for each element analysed were varied by mineral phase and listed in Supplementary File S1. We used a beam current of 10 nA for plagioclase and amphibole. For clinopyroxene, we analysed Si, Ti, Al, Fe, Mg, Ca, and K at 10 nA and Na, Mn, Cr, Ni, and P at 40 nA. The following standards were used for calibration; jadeite (Na), wollastonite (Si), anorthite (Ca, Al), orthoclase (K), rutile (Ti), fayalite (Fe), tephroite (Mn), periclase (Mg), Cr₂O₃ (Cr), apatite (P), NiO (Ni).

We monitored accuracy by analysing the secondary standards anorthite, Kakanui augite and Kakanui hornblende (Jarosewich et al., 1980). We typically measured major element concentrations within 4% of preferred values and minor elements within 1–23%. For unknowns, we measured major elements with 1σ precision better than 2% and minor elements with more variable precision, dependent on concentration. For example, Na₂O was measured with 1σ precisions better than 2% in plagioclase and amphibole, and 4% in clinopyroxene. The compositions of clinopyroxene and amphibole measured in the two different laboratories closely overlap (Supplementary Fig. S1).

The reproducibility of plagioclase An (An = Ca/(Ca + Na + K)) and clinopyroxene and amphibole Mg# (Mg# = Mg/(Mg + Fe), all Fe as Fe²⁺) measurements was assessed based on the repeat analyses of secondary standards plagioclase, Kakanui augite, and Kakanui hornblende. Reproducibility (expressed as (2SD/mean) × 100) of plagioclase An was less than 2.5% relative (Supplementary File S1). Reproducibility of clinopyroxene Mg# was less than 0.5% relative and amphibole Mg# less than 6% (Supplementary File S1).

THERMOBAROMETRY, HYGROMETRY, AND CHEMOMETRY CALCULATIONS

We performed all thermobarometry and hygrometry calculations using the Python3 tool Thermobar (version 1.0.13, Wieser et al., 2022). We did not find groundmass glass in our Bagana lavas and melt inclusions in our samples were too small for analysis. To our knowledge, melt inclusions from Bagana have not been analysed in previous studies, which makes identifying liquid compositions in equilibrium with mineral phases difficult. Therefore, we used thermobarometric and hygrometric formulations that require only the composition of mineral phases as input data.

To estimate pressures and temperatures of clinopyroxene crystallization, we used the machine learning-based clinopyroxene-only thermobarometer of Jorgenson et al. (2022). The quoted standard error of estimate (SEE) is ±320 MPa, 72.5°C. However, this thermobarometric model has the advantage of allowing the user to calculate the uncertainty on individual pressure–temperature estimates, which can be significantly smaller than the quoted model SEE (Jorgenson et al., 2022). When applied to a test dataset of clinopyroxene from experiments performed at pressure–temperature-melt H₂O conditions relevant to arc settings, the Jorgenson et al. (2022) thermobarometer produced the smallest mismatch between calculated and experimental pressures, compared with other clinopyroxene-only thermobarometers (Wieser et al., 2023a). In addition, the Jorgenson et al. (2022) thermobarometer requires only clinopyroxene compositions as input data, unlike models such as Wang et al. (2021) that require an additional assumption of melt H₂O content to calculate temperature. Our reported pressure–temperature estimates represent the mean value of the machine learning voting distribution for each input clinopyroxene composition, with uncertainty shown as the interquartile range of the voting distribution (as recommended by Jorgenson et al., 2022).

To estimate pressures and temperatures of amphibole crystallization, we used the amphibole-only barometer of Ridolfi (2021) and thermometer of Ridolfi & Renzulli (2012), with pressure calculated first and subsequently used in the equation to calculate temperature. Using these equations, uncertainties on pressure and temperature estimates are quoted as ±12% and ± 23.5°C respectively (Ridolfi & Renzulli, 2012; Ridolfi, 2021). When performing calculations in Thermobar, we applied the filter recommended by Ridolfi (2021) to remove any amphibole analyses with low oxide totals or which fall outside the compositional range of the calibration dataset. Our clinopyroxene and amphibole data fall within the compositional range of the calibration datasets for our chosen thermobarometers (Supplementary Fig. S2). We estimated H₂O contents and major element compositions of melts in equilibrium with amphibole using the formulations of Ridolfi (2021) (uncertainty of ±14%) and Humphreys et al. (2019), respectively. Uncertainties on calculated melt major element compositions are: SiO₂ (±3.73 wt%), Al₂O₃ (±1.31 wt%), TiO₂ (± 0.64 wt%), MgO (±1.12 wt%), CaO (±1.45 wt%), FeO (±2.02 wt%), K₂O (±0.76 wt%).

RESULTS

Petrography

Bagana lavas erupted in 2005 and 2012 are highly porphyritic with 38–55 vol% crystals, which are divided into macrocrysts (>500 μm length) and mesocrysts (100–500 μm length), after Zellmer (2021),. The mineral assemblage consists of plagioclase (61–76 vol% of total crystal assemblage, macrocrysts + mesocrysts), clinopyroxene (12–27 vol%), amphibole (1–14 vol%), orthopyroxene (3–6 vol%) and magnetite (5–6 vol%), plus very minor amounts of apatite. We found rare individual olivine crystals in two samples. The groundmass is largely microcrystalline and dominated by plagioclase microlites, along with minor oxides, pyroxenes, and sparse apatite. These characteristics are consistent with lavas from Bagana erupted pre-1989 (Blake, 1968; Bultitude et al., 1978; Rogerson, 1989).

We distinguish three populations of plagioclase crystals in Bagana lavas based on zoning patterns (Figs. 3 and 4). Plagioclase crystals of population 1 are subhedral to euhedral and characterised by high Ca bytownite (An_72–86) cores with sharp transitions to more sodic labradorite to andesine (An_47–56) rims (Figs. 3, 4, and5). Larger population 1 plagioclase crystals (>1 mm) often show elongate, sieve-textured cores with inclusions of glass and rarely clinopyroxene, with narrow rims. Smaller population 1 plagioclase crystals (0.1–0.8 mm) have more variable rim thicknesses and we interpret these as cross-sections of larger crystals due to their identical chemical characteristics. Population 2 plagioclase crystals (0.5–2 mm, most >1 mm) are subhedral to euhedral and characterised by cores (which may be homogenous or show internal zonation) displaying a range (bytownite or labradorite, An_56–78) of compositions and fine-scale (typically <10–20 μm bands) oscillatory zoned mantles and rims (Figs. 3, 4, and5). Some population 2 plagioclase crystal mantles have thick bands with distinctly higher Ca compared with adjacent zones, some of which show sieve textures (Fig. 5c). Some population 2 plagioclase crystals have high-Ca bytownite cores that may be sieved or partially resorbed (Fig. 4). Population 3 plagioclase crystals (0.1–0.8 mm) are subhedral to euhedral and unzoned (Fig. 4), with compositions overlapping with population 1 and 2 plagioclase rims.

Fig. 3

QEMSCAN images (5 μm resolution) showing typical textures and mineral assemblage of Bagana lavas, from 2005 (upper) and 2012 (lower). Typical examples of plagioclase population 1, mainly found in 2005 lavas, and plagioclase population 2, mainly found in 2012 lavas, are labelled.

Fig. 4

QEMSCAN Ca maps (5 μm resolution) showing typical examples of plagioclase from populations 1, 2 and 3. A brighter colour corresponds to a greater concentration of Ca (higher plagioclase An content). Population 1 plagioclase are characterised by high Ca cores and low Ca rims. Population 2 plagioclase have cores with variable Ca content and oscillatory zoned mantles. Note high Ca sieved cores in population 2 plagioclase from sample 2012–1. Population 3 plagioclase are unzoned.

Fig. 5

a) BSE image of typical population 1 plagioclase with sieved core, from sample 2005–1. b) BSE image of typical population 2 plagioclase (centre), with a high-An core (paler grey) surrounded by an oscillatory zoned mantle (darker grey), from sample 2012–2. c) BSE image of population 2 plagioclase with a mantle containing a thick high-An band with sieve texture (paler grey), between lower-An zones (darker grey), from sample 2012–2. d) XPL image of a typical clinopyroxene, displaying both concentric and sector zoning, with magnetite and plagioclase inclusions, from sample 2005–1. e) BSE image of typical clinopyroxene, displaying concentric zoning, with magnetite and plagioclase inclusions, from sample 2012–2. f) BSE image of typical amphibole, with unaltered core and breakdown rim composed of plagioclase, clinopyroxene, orthopyroxene, and magnetite, from sample 2005–2. g) BSE image of amphibole displaying internal alteration, with the interior partially replaced by microcrystalline plagioclase, pyroxenes, and magnetite, from sample 2012–2. h) PPL image of typical glomerocryst, composed of clinopyroxene (CPX), plagioclase, and magnetite (Mt). Note the plagioclase filling interstices between clinopyroxene. Red scale bars (b,c,d,e,g,h) = 500 μm. White scale bars(a,f) = 100 μm.

Plagioclase crystals from 2005 lavas are mostly population 1 (61 vol% of plagioclase macrocrysts and mesocrysts) with a smaller fraction of population 3 (39 vol%). The majority of plagioclase crystals from 2012 lavas are from population 2 (87–96 vol%) with a minor contribution from population 3 (4–12 vol%) and sparse population 1 (≤1 vol%). Sample 2012–3 is an exception, with 66 vol% population 1, 5 vol% population 2 and 29 vol% population 3 plagioclase crystals.

Clinopyroxene crystals are predominantly macrocrysts (0.2–2.3 mm, most >1 mm) and mostly euhedral with planar crystal faces in contact with the groundmass (Figs. 3 and 5). Clinopyroxene crystals contain abundant inclusions of magnetite, glass, and occasionally plagioclase, while sparse orthopyroxene inclusions are found in clinopyroxene crystals from 2012 lavas only. The majority of clinopyroxene crystals show concentric oscillatory zoning (typically ~30–150 μm-thick zones), which is sometimes superimposed upon sector zoning (Fig. 5d, e).

Amphibole crystals are predominantly macrocrysts (0.2–2.7 mm, most >0.5 mm), subhedral to euhedral and surrounded by breakdown rims of variable thickness (54–105 μm) composed of plagioclase, clinopyroxene, oxides ± orthopyroxene (Fig. 5f, g). Some crystals also exhibit extensive internal alteration in the form of embayment or complete replacement of the interior by microcrystalline plagioclase, pyroxenes, oxides, and, in rare cases, olivine (Fig. 5g).

Orthopyroxene mesocrysts (maximum length 0.5 mm) are generally euhedral and occur as free groundmass crystals (Fig. 3). Magnetite mesocrysts (maximum length 0.5 mm) are subhedral to euhedral and occur as free groundmass crystals or commonly as inclusions in clinopyroxene.

Glomerocrysts are present in all our lava samples (Figs. 3 and 6). The majority of glomerocrysts consist of the following assemblages: clinopyroxene + magnetite (31% of all observed glomerocrysts in examined thin sections), clinopyroxene + plagioclase ± magnetite (20%), clinopyroxene + plagioclase + amphibole ± magnetite (9%), clinopyroxene + amphibole ± magnetite (10%), and plagioclase only (29%) (Fig. 6). In 2012 lavas, we found individual glomerocrysts composed of amphibole + plagioclase and clinopyroxene + magnetite + olivine (Fig. 6f). In some glomerocrysts, crystals are joined along near planar boundaries, whereas others contain interlocking crystals with random orientations (Fig. 6b,c,d). Clinopyroxene in glomerocrysts are texturally similar to individual crystals, with many displaying concentric and/or sector zoning (Fig. 6). Different zoning patterns are present in adjacent clinopyroxene crystals within the same glomerocryst (Fig. 6d). Amphibole do not have breakdown rims at shared boundaries with clinopyroxene crystals in glomerocrysts (Fig. 6e). Plagioclase crystals in polymineralic glomerocrysts are typically bytownite and often display more sodic labradorite to andesine rims in contact with the groundmass (Fig. 6a). The majority of plagioclase-only glomerocrysts are formed of population 2 plagioclase crystals. In some glomerocrysts, plagioclase occupies interstices between clinopyroxene crystals (Figs. 5h and 6a,b).

Fig. 6

a) Extract of QEMSCAN map (5 μm resolution) showing a glomerocryst composed of clinopyroxene + plagioclase + magnetite + minor amphibole. Note the low-An rims on high-An plagioclase in contact with the groundmass (marked by the red arrows) and interstitial plagioclase. From sample 2005–2. b) XPL image of the glomerocryst in a). The interlocking texture of plagioclase and clinopyroxene is marked by the small red arrows. c) XPL image of glomerocryst composed of clinopyroxene + plagioclase + magnetite. An example of a near planar crystal boundary between clinopyroxene and plagioclase is marked by the small red arrow. From sample 2005–2. d) XPL image of glomerocryst composed of clinopyroxene + plagioclase + magnetite. Examples of near planar crystal boundaries are marked by the small red arrows. Note the presence of adjacent clinopyroxene crystals with different zoning patterns. From sample 2012–3. e) XPL image of glomerocryst composed of clinopyroxene + amphibole + magnetite. From sample 2005–1. f) XPL image of glomerocryst composed of clinopyroxene + magnetite + olivine. The labelled clinopyroxene (CPX) crystals adjacent to the olivine (Ol) have high Mg# (>80) compositions. From sample 2012–2. g) XPL image of glomerocryst composed of clinopyroxene + magnetite from sample 2005–2. h) XPL image of plagioclase-only glomerocryst from sample 2012–2. Red scale bars (a,b,c,e) = 500 μm. White scale bars (d,f,g,h) = 100 μm.

Whole-rock geochemistry

We report the major, trace element and radiogenic isotopic compositions of 2005–12 Bagana lavas and 2016 ash in Supplementary File S1. Figure 7 shows the major element compositions of 2005–12 Bagana lavas compared with historical eruptions (pre-1989, Taylor et al., 1969; Bultitude et al., 1978; Rogerson, 1989). Our Bagana lavas occupy a narrow compositional range of ~56–58 wt% SiO₂, and straddle the compositional boundary between basaltic andesites and andesites (LeMaitre, 2002) with lavas erupted in 2012 being slightly more evolved. The 2016 ash samples have an andesitic composition with 60 wt% SiO₂. All Bagana lavas plot in the medium-K field of Gill (1981) (Fig. 7g). Our Bagana lavas overlap the compositions of historical eruptions, which form linear arrays on major element oxide vs SiO₂ variation diagrams (Fig. 7), with MgO, FeO, CaO, and TiO₂ correlated negatively with SiO₂ and Al₂O_3, Na₂O and K₂O showing weak positive correlations, while P₂O₅ is constant (Supplementary File S1).

Fig. 7

Major oxides vs SiO₂ for Bagana lavas and ash from this study and historic eruptions (triangle symbols, data from Taylor et al., 1969; Bultitude et al., 1978; Rogerson, 1989). All data have been normalized to 100% volatile-free with all Fe as FeO. Uncertainties for our Bagana data are smaller than symbol size. The low-K, medium-K, and high-K fields in (g) are taken from Gill (1981). Also shown are calculated groundmass compositions for Bagana lavas and calculated compositions of melts in equilibrium with amphibole (AEMs). Uncertainties on AEMs are listed within the text.

Large ion lithophile element (LILE) e.g. Ba, Rb, and some high field strength element (HFSE) e.g. Zr concentrations are higher in the 2016 ash (Ba = 305 ppm) than 2005–12 lavas (Ba = 226–260 ppm), while Nb (4–6 ppm), Y (15–19 ppm), and REEs are similar (Fig. 8, Supplementary Fig. S3). The 2005–12 lava trace element compositions overlap closely with historical eruptions, except for lower Rb (Fig. 8, Supplementary Fig. S3). Lavas from 2005 and 2012 have very similar trace element compositions that are strongly enriched in LILEs relative to N-MORB, typical for island arc lavas (e.g. Pearce & Peate, 1995) (Fig. 8). All lavas and ash show small negative Eu anomalies and HFSE contents are slightly enriched relative to N-MORB (Fig. 8). We analysed Sr, Nd, and Pb radiogenic isotopes in one lava sample from both 2005 and 2012 and one 2016 ash sample. The three samples have near identical compositions, with ⁸⁷Sr/⁸⁶Sr = 0.703662–0.703687, ¹⁴³Nd/¹⁴⁴Nd = ~0.512906–0.512924, ²⁰⁶Pb/²⁰⁴Pb = 18.635–18.647, ²⁰⁷Pb/²⁰⁴Pb = 15.522–15.524, and ²⁰⁸Pb/²⁰⁴Pb = 38.317–38.328 (Fig. 9). These samples were also analysed for bulk rock oxygen isotopic compositions, returning δ¹⁸O values of 8.3‰ (2005 lava), 6.7‰ (2012 lava), and 7.6‰ (2016 ash).

Fig. 8

a) N-MORB normalized trace element patterns for Bagana lavas and ash from this study. N-MORB normalizing values from Sun & McDonough (1989). (b) Ba and (c) Rb vs SiO₂ for Bagana lavas and ash from this study and historic eruptions (data from Taylor et al., 1969; Bultitude et al., 1978; Rogerson, 1989). Symbols are the same as Fig. 7. d) Th/Yb vs Ba/Th and e) Th/Yb vs Nb/Yb of Bagana lavas and ash (yellow triangles, this study), compared with basalts from nearby arc systems. Uncertainties for our Bagana data (based on reproducibility of secondary standard BCR-2G, expressed as (2SD/mean) × 100, Supplementary File S1) are either shown by the representative error bar or are smaller than symbol size if an error bar is not shown. The MORB array in (e) is taken from Pearce & Peate (1995). Solomon Island arc basalt data from Schuth et al. (2004, 2009), Petterson et al. (2011). Tabar-Lihir-Tanga-Feni (TLTF) basalt data from Stracke & Hegner (1998) and Müller et al. (2001). New Britain arc basalt data from Woodhead et al. (1998).

Fig. 9

a)  ¹⁴³Nd/¹⁴⁴Nd vs ⁸⁷Sr/⁸⁶Sr and b)  ²⁰⁸Pb/²⁰⁴Pb vs ²⁰⁶Pb/²⁰⁴Pb for Bagana lavas and ash (yellow triangles) compared with lavas from nearby arcs and potential representatives of mantle source and subduction components. Uncertainties for our Bagana data are smaller than symbol size. Data for Manus Basin MORB, Solomon Sea plate (SSP) basalts and sediments and New Britain arc lavas from Woodhead et al. (1998). Woodlark basin N-MORB data from Park et al. (2018). Solomon Island arc data from Schuth et al. (2004, 2009, 2011). Tabar-Lihir-Tanga-Feni (TLTF) lava data from Stracke & Hegner (1998) and Kamenov et al. (2005). TLTF mantle xenolith data from Kamenov et al. (2008). The thick black line in (b) shows the boundary between the Pacific and Indian mantle domains, from Kempton et al. (2002).

Mineral chemistry

Plagioclase

Plagioclase crystals from population 1 show different An contents between cores (An_72–86) and rims (An_47–56, Fig. 10a). Plagioclase crystals from population 2 show a wider range of core compositions (An_56–78) and generally slightly less evolved (An_50–62) rims compared with population 1 (Fig. 10a). Unzoned population 3 plagioclase crystals can be divided into two groups, An_46–52 and An_57–62 (excluding one crystal with An₆₆), corresponding to the rim compositions of population 1 and population 2 plagioclase crystals respectively (Fig. 10a). Minor elements (FeO, TiO₂) show a similar compositional range for cores and rims in all plagioclase populations (Supplementary File S1).

Fig. 10

a) K₂O vs An (An = Ca/(Ca + Na + K)) for Bagana plagioclase crystals. The smaller circle symbols with green outline show the compositions of population 3 plagioclase. b) Results of Ab-An exchange test for equilibrium between plagioclase and bulk-rock compositions (Putirka, 2008). Plagioclase crystals are considered in equilibrium with bulk-rocks if the calculated K_D(An-Ab)^Pl-liq value lies within the range 0.28 ± 0.11, shown by the red shaded region. A pressure of 400 MPa was assumed for the calculations, though no difference in results was observed when a pressure of 200 MPa was tested. Plagioclase crystal compositions from 2005 and 2012 were compared with average bulk-rock compositions from the respective eruptions. c), d) and e) Compositional profiles showing variation in An across plagioclase crystals, highlighting typical zoning patterns of the different plagioclase populations. The red lines across crystals mark the locations of profiles. Small red numbers in some BSE images mark locations of EPMA point analyses. Red scale bars = 200 μm. Uncertainty on plagioclase K₂O and An is shown by the representative error bars. Uncertainty on K₂O is based on 1σ precisions from counting statistics (Supplementary File S1). Uncertainty on An is based on reproducibility of the plagioclase secondary standard An composition (Supplementary File S1).

Figure 10c illustrates the sharp normal zoning (~30 mol% An decrease between core and rim) typical of population 1 plagioclase. Figure 10d,e shows examples of compositional profiles from population 2 plagioclase, highlighting the variable core compositions. Although the spot size used for EPMA analyses (5–10 μm) may not always capture individual growth bands, profiles show fluctuations in An content of typically <10 mol% across oscillatory zoned regions towards rims. The occasional thick, high-An bands in some population 2 plagioclase mantles have approximately 10 mol% greater An (~An₇₀) than adjacent zones (~An₆₀) (Supplementary Fig. S4).

Clinopyroxene and orthopyroxene

Clinopyroxene from 2005–12 lavas are predominantly classified as augite, with occasional diopsides (Wo_40–48, En_38–47, Fs_10–17) and are chemically similar (Fig. 11b, c, d). Most clinopyroxene show concentric zoning (Fig. 5d, e) and core refers to the composition of the innermost concentric zone, measured at a point in the centre of a crystal, while rim refers to the outermost concentric zone, measured close to the edge of a crystal. Al₂O₃, Na₂O, and TiO₂ (not shown) broadly increase with decreasing Mg# (Fig. 11b, d). Cr₂O₃ contents are very low (<0.1 wt%) except for Mg# >80 cores (up to 0.64 wt%, Supplementary File S1). The range of chemical compositions within clinopyroxene from glomerocrysts and free crystals is similar (Fig. 12, Supplementary Fig. S4), except for the rare high Mg# (>80), high Cr₂O₃ cores that are only found in glomerocrysts (Fig. 6).

Fig. 11

a) Test for Fe-Mg equilibrium between clinopyroxene and bulk-rock compositions. Clinopyroxene are considered in equilibrium with bulk-rocks if the calculated K_D(Fe-Mg)^cpx-liq value lies within the range 0.28 ± 0.08, marked by the dashed lines. Clinopyroxene compositions from 2005 and 2012 were compared with average bulk-rock compositions from the respective eruptions. b) Al₂O₃  c) CaO and d) Na₂O vs Mg# (Mg# = Mg/(Mg + Fe), all Fe as Fe²⁺) for Bagana clinopyroxene. e) Al₂O₃ and f) TiO₂ vs Mg# for Bagana amphibole. Uncertainties are shown by the representative error bars. Uncertainties on major element oxides are based on 1σ precisions from counting statistics (Supplementary File S1). Uncertainties on Mg# are based on reproducibility of the Kakanui augite and Kakanui hornblende secondary standard Mg# compositions (Supplementary File S1). Where no error bar is shown, uncertainties are smaller than symbol size.

Fig. 12

a), b), c) and d) compositional profiles across clinopyroxene crystals, displaying a variety of zoning patterns typical of the wider clinopyroxene population. Black circles and line show clinopyroxene Mg# (Mg# = Mg/(Mg + Fe), all Fe as Fe²⁺). Red triangles and line show clinopyroxene Al₂O_3, orange diamonds and line clinopyroxene TiO₂. The red lines across crystals mark the locations of profiles. Small red numbers in some BSE images mark locations of EPMA point analyses. Red scale bars = 500 μm. Uncertainties are smaller than symbol size.

Clinopyroxene displays a variety of chemical zoning patterns and consistent trends are not observed between crystals from the same sample or individual eruptions (Fig. 12). Concentric zoning is most clearly defined by fluctuations in Mg# (up to 8 mol% variation within a single crystal) and Al₂O₃ (up to 4 wt% variation), with zones of higher Al₂O₃ and TiO₂ broadly corresponding with lower Mg# (Fig. 12). Some clinopyroxene profiles capture a final, very narrow low Mg# (~72) outermost rim (Fig. 12c and Supplementary Fig. S4). These cannot be resolved on all crystals due to the large analytical spot size used (1 or 5 μm). Some glomerocryst clinopyroxene crystals have relatively homogenous core regions with contrasting compositions to concentric zoned mantles (Fig. 12b, Supplementary Fig. S4). Orthopyroxene mesocrysts display a narrow compositional range (Wo_3–6, En_69–75, Fs_22–28, Mg# 72–77, Supplementary File S1).

Amphibole

Due to the presence of breakdown rims, most amphibole measurements were taken from unaltered interiors. The few low-Mg# (<70), Fe-rich amphibole represent points directly adjacent to breakdown rims (Fig. 11e, f). The majority of amphibole are classified as magnesio-hastingsites with a few tschermakitic pargasites, according to the scheme of Ridolfi et al. (2010). Amphibole compositions (excluding points adjacent to breakdown rims) form a broad cluster (Mg# 70–79), with slightly greater compositional variation in 2012 amphiboles compared with 2005, particularly in Al₂O₃ (Fig. 11e, f). Amphibole generally show minimal chemical zonation with less than ~1 wt% fluctuations in major oxides (Al₂O₃, FeO, and MgO) and < 0.5 wt% fluctuations in minor oxides (TiO₂, Na₂O) across profiles (Supplementary Fig. S5).

DISCUSSION

Bulk-rock compositions and crystal accumulation

The high crystallinity of the groundmass in our Bagana lavas hinders direct measurement of the erupted melt composition. Therefore, we calculated approximate groundmass/erupted melt compositions for 2005–12 Bagana lavas by mass balance using bulk-rock compositions and representative mineral compositions measured by EPMA. Groundmass compositions were calculated using the following equation for each major element oxide, following the approach of D’Mello et al. (2023):

$$ \begin{align*} \text{C}_{\text{GM}}=(\text{C}_{\text{WR}}-\text{XC}_{\text{MIN}})/(1-\text{X}), \end{align*} $$

where C_GM = groundmass composition (wt% oxide), C_WR = bulk-rock composition (wt% oxide), C_MIN = average mineral assemblage composition (wt% oxide), and X = proportion of crystals. The value of X was set at 0.46 based on the average proportion of macrocrysts + mesocrysts, which varies from 38–55 vol% (Supplementary File S1). Details of the mineral phase compositions and proportions used to calculate C_MIN are listed in Table 1. Since the compositions and proportions of the mineral phases in 2005 and 2012 lavas are broadly similar (Figs. 10 and 11), our calculations assume the same average mineral assemblage composition for all samples. Without undertaking time-consuming detailed chemical mapping of thin sections, it is difficult to accurately constrain the mineral compositional variation within individual samples and, hence, the bulk composition of the mineral assemblage. Therefore, our use of an average mineral assemblage composition and proportion means that our groundmass calculations are subject to significant uncertainty and are considered an approximation of the erupted melt composition. The calculated groundmass compositions are dacitic (~63–66 wt% SiO₂, after LeMaitre, 2002) and considerably more evolved than bulk-rock compositions (Fig. 7).

Table 1

Representative mineral phase compositions and proportions used to calculate the average mineral assemblage composition (C_MIN), used for calculation of Bagana lava groundmass compositions

Mineral phase	Proportion	Composition
		SiO2	TiO2	Al2O3	FeO	MnO	MgO	CaO	Na2O	K2O
Plagioclase (andesine, An30–50)	0.23	55.48	0.04	26.86	0.65	0.00	0.07	9.84	5.69	0.29
Plagioclase (labradorite, An50–70)	0.23	51.77	0.03	28.90	0.52	0.01	0.03	12.09	4.43	0.20
Plagioclase (bytownite, An 70–90)	0.20	48.65	0.02	31.01	0.65	0.02	0.04	14.99	2.71	0.10
Clinopyroxene	0.185	51.12	0.40	2.21	8.27	0.47	15.04	20.90	0.30	0.00
Amphibole	0.05	41.84	2.45	11.83	12.27	0.23	13.84	11.36	2.54	0.57
Orthopyroxene	0.045	52.73	0.28	1.84	16.16	0.94	25.27	1.78	0.05	0.00
Magnetite	0.06	0.08	7.58	3.21	81.31	0.62	2.54	0.03	0.00	0.00
CMIN		48.32	0.68	20.30	8.15	0.19	4.80	12.56	3.05	0.16

The proportions are based on the average relative abundances of the mineral phases in Bagana lavas (Supplementary File S1). The clinopyroxene, amphibole, orthopyroxene, magnetite, and plagioclase bytownite and labradorite values were taken from crystal cores. The plagioclase andesine value was taken from a population 1 plagioclase rim.

A significant proportion of macrocryst cores in Bagana lavas are not in chemical equilibrium with the bulk rocks. Equilibrium between plagioclase and (normalized anhydrous) bulk-rock compositions was assessed using the Ab–An exchange equilibrium test of Putirka (2008), implemented using Thermobar (version 1.0.13, Wieser et al., 2022). Plagioclase crystal compositions from 2005 and 2012 lavas were compared with average bulk rock compositions from the respective eruptions (~56 wt% SiO₂ for 2005, ~57 wt% SiO₂ for 2012). These equilibrium tests were performed assuming a pressure of 400 MPa, based on the thermobarometry calculations and comparisons with experiments presented below. Equilibrium tests show that none of the population 1 plagioclase cores, or population 2 plagioclase cores >An₆₇, are feasibly in equilibrium with the bulk rocks (Fig. 10b).

Basalts have not been reported from Bagana and are rare among Pliocene-Quaternary lavas on Bougainville (Rogerson, 1989). The nearest Pliocene-Quaternary volcanic centre to Bagana to have erupted basaltic lavas is the Numa-Numa volcano to the north (location 3 on Fig. 1). A significant proportion of Bagana population 1 and 2 plagioclase crystal cores (An_68–80) are feasibly in equilibrium with basaltic (~52 wt% SiO₂) compositions like that found at Numa-Numa (sample NGV748 of Rogerson, 1989, Supplementary Fig. S6). Plagioclase rims are feasibly in equilibrium with the bulk rocks (Fig. 10b), but also more evolved (~60 wt% SiO₂) compositions represented by the 2016 ash (Supplementary Fig. S6). This result may reflect the limited sensitivity of the equilibrium test. Our calculated groundmass compositions only serve as an approximation of the erupted melt composition and do not account for any dissolved volatile phases that may have been present in the melt from which plagioclase rims crystallized. Hence, our calculated groundmass may not accurately represent the melt in equilibrium with plagioclase rims and we have not tested for equilibrium between plagioclase rims and calculated groundmass.

Assuming K_D(Fe-Mg)^cpx-liq = 0.28 ± 0.08 (Putirka, 2008), the majority of clinopyroxene are plausibly in equilibrium with 2005–12 lava bulk-rock compositions (Fig. 11a), excluding rare clinopyroxene cores with Mg# >80 found only in glomerocrysts. It should be noted that the wide range of permitted K_D(Fe-Mg)^cpx-liq equilibrium values may limit the sensitivity of this equilibrium test. Melts in equilibrium with amphibole (amphibole equilibrium melts or AEMs) were calculated from measured amphibole core compositions using the equations provided by Humphreys et al. (2019). The calculated AEMs (57–65 wt% SiO₂) show that amphibole crystallized from melts with generally more evolved compositions than the bulk lavas (Fig. 7), with similar SiO₂, MgO, CaO, and K₂O to our calculated groundmass.

Chemical disequilibrium between macrocryst cores and bulk-rock compositions suggests that these crystals represent an accumulated cargo. Therefore, the 2005–12 Bagana lava bulk compositions most likely reflect mixtures of chemically evolved carrier melts and a mafic crystal cargo, consistent with highly porphyritic lavas at arc volcanoes worldwide (e.g. Reubi & Blundy, 2008; Cassidy et al., 2015; Namur et al., 2020; Reubi & Müntener, 2022; Weber et al., 2023; Zemeny et al., 2023).

In the absence of melt inclusion or groundmass glass analyses, our calculated groundmass compositions and the AEMs represent the best available proxy for the range of melt compositions that could accumulate a crystal cargo to form Bagana lavas. The chemical variation recorded by the AEMs suggests that carrier melts feeding Bagana eruptions may have a wide range of compositions spanning andesite to dacite. The crystal content of the 2016 ash has not been determined, but its more evolved composition, which generally overlaps the AEMs (Fig. 7), and slightly higher incompatible trace element content (Fig. 8) compared with the lavas suggest that it may represent a true melt composition or perhaps a more melt-rich, crystal-poor eruption.

We suggest that the ash, groundmass and AEM compositions were produced by fractionation of the Bagana lava mineral assemblage (plagioclase + clinopyroxene + amphibole + magnetite + orthopyroxene), in varying proportions, from basaltic precursor magmas, perhaps reflected by the presence of high-An population 1 plagioclase cores. The occurrence of rare olivine crystals in older lavas (Bultitude et al., 1978), in a glomerocryst associated with high Mg# clinopyroxene crystals in our Bagana samples (Fig. 6f) and in a hornblende gabbro plutonic inclusion sample (Bultitude, 1979) suggests that olivine may also be a fractionating phase during the early evolution of Bagana magmas, consistent with typical arc basalt crystallization sequences (Nandedkar et al., 2014; Marxer et al., 2021).

Bagana lavas erupted before 1989 contain the same mineral assemblage and similar proportions of groundmass and crystals (33–52 vol% crystals, Bultitude et al., 1978) to 2005–12 lavas. We suggest that the bulk-rock compositional variation observed in Bagana lavas over time (Fig. 7) is likely primarily controlled by variations in carrier melt composition, the proportion of carrier melt relative to crystals, and by compositional variations within the crystal cargo.

Compared with crystallization experiments on arc basaltic andesite (< 53 wt % SiO₂) starting compositions that produced a plagioclase-clinopyroxene-amphibole assemblage (Pichavant et al., 2002; Melekhova et al., 2017), Bagana lavas contain lower proportions of amphibole (< 14 vol% vs > 20 vol% in the experiments) and higher proportions of plagioclase and clinopyroxene. Hence, the crystal cargo in Bagana lavas may not be accumulated in approximately cotectic proportions, as observed at other arc volcanoes such as Calbuco (Namur et al., 2020). Our Bagana lavas have similar or lower incompatible element contents (Rb, Ba, K₂O) than historical eruptions with lower SiO₂ (Figs. 7 and 8), which demonstrates that bulk compositional variations are not simply related by fractional crystallization, which would result in increasing incompatible element concentrations with increasing SiO₂. Variations in incompatible element concentrations between eruptions with near identical SiO₂ (e.g. Rb in 1959–75 and 2005, Fig. 8) most likely reflect differences in the incompatible element content of the carrier melt.

Magmatic processes recorded by the crystal cargo

Population 1 plagioclase crystals are characterized by sharp decreases in An content (~30 mol%) between cores and rims and often display slightly rounded sieve-textured cores. The An content of plagioclase is influenced by crystallization temperature, pressure, melt composition, and H₂O content (Namur et al., 2012, and references therein). The sharp decrease in An between cores and rims in population 1 plagioclase crystals is inconsistent with growth during fractional crystallization and indicates a shift in crystallization conditions (e.g. Namur et al., 2020). Plagioclase An correlates positively with melt H₂O, CaO/Na₂O, and temperature (e.g. Sisson & Grove, 1993; Panjasawatwong et al., 1995; Ustunisik et al., 2014). Therefore, the lower-An rims could have crystallized from a more evolved, lower CaO/Na₂O, lower temperature melt, or a melt with lower H₂O content.

Sharp normal zoning in plagioclase crystals at other arc volcanoes has been explained by decompression, which results in H₂O loss from a carrier magma on ascent and subsequent crystallization of lower-An plagioclase rims during pre-eruptive storage at shallower levels or in the conduit (e.g. Coote & Shane, 2016; Namur et al., 2020; Moshrefzadeh et al., 2023; Weber et al., 2023). Sieve textures in plagioclase crystals have been shown experimentally to be produced by decompression (Nelson & Montana, 1992). Sieve textures may also be formed due to mixing with a more calcic or hotter melt (Tsuchiyama, 1985); however, this would be expected to result in crystallization of higher-An rims (reverse zoning). The experiments of Tsuchiyama (1985) also show that entrainment of calcic plagioclase into a more evolved melt, in equilibrium with relatively more sodic plagioclase than the entrained crystal, does not produce sieve textures. We do not have direct constraints on the H₂O contents of melts in the Bagana plumbing system, but note that entrainment of plagioclase in a H₂O rich plagioclase-undersaturated melt could also contribute to plagioclase dissolution (e.g. Ferlito et al., 2014).

The coarse sieve texture of Bagana population 1 plagioclase crystal cores and the lack of reverse zoning closely resembles the textures formed in the decompression experiments of Nelson & Montana (1992). Thus, we suggest that the combination of coarse sieve textured cores and abrupt decrease in An between cores and rims in our samples is best explained by decompression, during transport of high-An plagioclase cores towards the surface. Population 1 plagioclase rims are feasibly in equilibrium with the bulk rocks and 2016 ash (Supplementary Fig. S6), whereas cores are feasibly in equilibrium with basaltic compositions (Supplementary Fig. S6), suggesting that the rims crystallized from more evolved melts. Therefore, the lower An content of population 1 plagioclase rims compared with cores could reflect crystallization from a more evolved melt, which entrained and transported the high-An cores, H₂O loss from this melt upon ascent towards the surface, or the combined effects of both processes.

Population 2 plagioclase crystals are characterized by cores of variable An content, surrounded by finely oscillatory zoned mantles that extend to crystal rims. Some population 2 plagioclase crystals display high-An sieved cores matching population 1 plagioclase (Fig. 4). We suggest that these plagioclase cores were also subject to decompression, as opposed to mixing with a hotter or more mafic magma, since zones immediately surrounding the sieved cores are less calcic (Figs. 4 and 10d). The overall variation in population 2 plagioclase core compositions could reflect crystallization from magmas with different melt compositions (between basalt and andesite) and/or H₂O contents, temperature, or possibly at different pressures within the plumbing system.

Oscillatory zoning in plagioclase is often considered to develop due to frequent changes in magmatic environment experienced by a growing crystal. For example, movement of a crystal within a melt-rich magma body via convection can induce changes in local temperature, pressure, melt composition, and H₂O (e.g. Ginibre et al., 2002; Ustunisik et al., 2014). Alternatively, Blundy et al. (2006) suggested that oscillatory zoning may develop due to the opposing effects of decompression of a hydrous magma, and latent heat release due to decompression induced crystallization. The thicker (compared with typical oscillatory zones), occasionally sieved, higher-An bands found in the mantles of some population 2 plagioclase (Supplementary Fig. S4) most likely record recharge with a slightly more mafic or more hydrous magma or prolonged periods of growth at higher temperature. The lower An content of these bands (~An₇₀) compared with high-An plagioclase cores (An_72–86) is consistent with recharge. The variation in zoning patterns observed in population 2 plagioclase mantles, such as the occasional presence of thicker, higher-An bands, suggests that these crystals grew in distinct magma bodies or were subject to different conditions during growth. Population 3 plagioclase are unzoned and compositionally match the rims of population 1 and 2 plagioclase crystals. We suggest that population 3 plagioclase crystallized from the same melts as population 1 and 2 rims during magma ascent.

Clinopyroxene in Bagana lavas is characterized by concentric zoning that is superimposed on sector zoning in some crystals. Based on experimental evidence, the presence of sector zoning in euhedral clinopyroxene is consistent with crystallization at low-moderate degrees of undercooling (ΔT = 18–32°C, Kouchi et al., 1983, Masotta et al., 2020). Compositional zoning is most clearly defined by fluctuations in Mg# and Al₂O₃ and consistent compositional changes from core to rim are not observed between crystals (Fig. 12). Both point analyses and profiles across concentric zones demonstrate that Al₂O₃, TiO₂ and Na₂O broadly increase with decreasing clinopyroxene Mg# (Figs. 11 and 12). Covariation in Al, Ti, and Na with Mg is consistent with charge-balanced substitutions controlling the composition of clinopyroxene, whereby incompatible elements are enriched with increasing degree of undercooling (e.g. Kouchi et al., 1983; Mollo et al., 2018; Ubide et al., 2019; Masotta et al., 2020; MacDonald et al., 2023). Hence, the fluctuations in clinopyroxene Mg#, Al₂O₃ and TiO₂ across concentric zones could reflect crystal growth at variable degrees of undercooling. Alternatively, the compositional shifts could reflect changes in local melt composition around growing crystals (e.g. Elardo & Shearer, 2014). Cr₂O₃ in Bagana clinopyroxene is uniformly low, almost always below 0.1 wt% across a range of Mg# (Supplementary File S1). Recharge with a more mafic magma, which produces clinopyroxene zones with elevated Cr₂O₃ (e.g. Streck et al., 2002; Streck, 2008; Ubide & Kamber, 2018), is, therefore, unlikely to have contributed to the formation of concentric zonation in clinopyroxene crystals. Variations in the degree of undercooling and local melt composition experienced by clinopyroxene crystals could occur due to convection in a thermally or chemically stratified magma reservoir (e.g. Elardo & Shearer, 2014; Cao et al., 2023). The planar growth faces and euhedral form of clinopyroxene crystals, indicative of crystallization in a melt-rich environment (e.g. Holness et al., 2019), is consistent with this suggestion.

The lack of consistent core-to-rim zoning patterns between our clinopyroxene crystals suggests that these crystals may have formed in multiple magma bodies and/or were subject to variable thermochemical conditions during crystallization. We interpret the narrow lower-Mg# (~72) rims captured in some clinopyroxene crystal profiles (Fig. 12c, Supplementary Fig. S4) as growth from carrier melts during final ascent.

Based on isothermal decompression experiments, amphibole breakdown rims form due to H₂O loss from the coexisting melt upon magma ascent to low pressures outside the amphibole stability field (e.g. Rutherford & Hill, 1993; Rutherford & Devine, 2003; Browne & Gardner, 2006). These experiments simulate the decompression of magma during ascent from a deeper (~6–8 km) storage region where amphibole is stable to the surface to feed dome-building and Plinian eruptions. The constant-rate decompression experiments of Rutherford & Hill (1993), using a Mount St Helens dacite starting material, produced ~50 μm-thick breakdown rims over a duration of 20 days. In contrast, the constant-rate decompression experiments of Browne & Gardner (2006), using a dacite from Redoubt volcano, produced only 15 μm-thick breakdown rims in ~20 days. The compositions of the starting materials for these experiments (~63–64 wt% SiO₂, 1.8–2 wt% MgO) are close to our calculated groundmass compositions for Bagana lavas and overlap with the AEMs, suggesting that the results of these experiments are likely broadly applicable to Bagana. Therefore, we infer that the 54–105 μm thickness of Bagana amphibole breakdown rims implies that magmas are likely to have spent tens of days outside the amphibole stability field prior to eruption.

The presence of crystals joined along near planar boundaries in some glomerocrysts (Fig. 6) is consistent with an origin via synneusis, the aggregation of isolated crystals in a melt-rich or convecting environment (Vance, 1969; Schwindinger & Anderson, 1989; Holness et al., 2017, 2019). Variation in zoning patterns displayed by adjacent clinopyroxene crystals within the same glomerocrysts (Fig. 6) supports the idea that isolated crystals with different growth histories were aggregated together. In contrast, some glomerocrysts contain interlocking textures and interstitial plagioclase filling pore space between clinopyroxene crystals, closely resembling textures observed in partially to fully solidified crystal mushes (e.g. Holness et al., 2022; Winslow et al., 2022). We suggest that the glomerocrysts in Bagana lavas represent a combination of crystal aggregates formed via synneusis and entrained crystal mush fragments. Bultitude (1979) reported a single holocrystalline hornblende gabbro plutonic inclusion containing interstitial amphibole in Bagana lavas, providing further evidence for the presence of crystal mush zones in the Bagana plumbing system.

Magma storage conditions

To our knowledge, there are no published estimates for the crustal thickness beneath Bougainville. Seismic studies of the nearby Solomon Island arc suggest crustal thicknesses of ~25–35 km (Ku et al., 2020). Based on the similar tectonic history of Bougainville and the Solomon Islands (e.g. Holm et al., 2016), the crust beneath Bagana may be a similar thickness. A global model based on seismic data predicts a crustal thickness of 37 ± 6 km beneath Bougainville (Szwillus et al., 2019).

By combining thermobarometric and hygrometric calculations with constraints from experimental petrology, we can estimate the magma storage conditions (pressure/depth, temperature, melt H₂O content) under which the crystal cargo in Bagana lavas formed. Clinopyroxene crystallization temperatures and pressures, calculated using the machine learning-based thermobarometer of Jorgenson et al. (2022), are shown in Fig. 13a,b. The majority of clinopyroxene crystallized at 970–1100°C, ~90–370 MPa, with a few returning higher pressures of ~430–570 MPa. Uncertainties vary for individual pressure–temperature estimates, typically ±30–50°C and ± 100–150 MPa, with the higher-pressure estimates showing greater uncertainties of ±250 MPa. There is no systematic difference in crystallization temperatures and pressures between cores and rims or between clinopyroxene from 2005 and 2012 within uncertainty (Fig. 13a,b, see Supplementary Fig. S7 for further discussion of uncertainties). Due to the absence of constraints on the crustal structure beneath Bougainville, we convert pressures to depths assuming a uniform crustal density of 2700 kg/m³ (Wieser et al., 2022). Based on this assumption, we infer that clinopyroxene crystallized over a depth range of ~3–21 km (Supplementary Fig. S8).

Fig. 13

Calculated crystallization temperatures and pressures for clinopyroxene from a) 2005 and b) 2012, using the thermobarometer of Jorgenson et al. (2022). c) Calculated crystallization temperatures and pressures for amphibole, using the thermometer of Ridolfi & Renzulli (2012) and barometer of Ridolfi (2021). All pressures were calculated using the equation P = (P1b + P1c)/2 from Ridolfi (2021). d) Amphibole crystallization pressure vs H₂O content of melt in equilibrium with amphibole, calculated using the hygrometer of Ridolfi (2021).

We calculated amphibole crystallization temperatures and pressures using the barometer of Ridolfi (2021) and thermometer of Ridolfi & Renzulli (2012) (Fig. 13c). Ridolfi (2021) proposes five equations (1a to 1e), calibrated for different pressure ranges. An algorithm is used to determine which equation result should be taken as the final pressure estimate. Applying the algorithm as recommended, we used the formula P = (P1b + P1c)/2 for the majority of our amphibole analyses, returning calculated pressures of ~300–500 MPa. The algorithm suggested that equations P = P1d and P = P1e were most suitable for the remaining analyses, which output higher calculated pressures of up to 1400 MPa (Supplementary Fig. S9). We consider these pressures implausible as they correspond to mantle depths, assuming ~3.8 km per 100 MPa based on a crustal density of 2700 kg/m³ and typical island arc crustal thicknesses of <35 km (e.g. Melekhova et al., 2019). Therefore, we chose to calculate pressures for all amphibole analyses using the same equation, P = (P1b + P1c)/2. The resulting pressure estimates range from 300–500 MPa (Fig. 13c) and are associated with model uncertainties of ±39–59 MPa (±12% of calculated values). These pressures correspond to depths of ~11–19 km (±2 km) (Supplementary Fig. S8). Calculated amphibole crystallization temperatures range from 910–970°C ± 23.5°C (Fig. 13c). The pressure–temperature estimates for amphibole from 2005 and 2012 lavas overlap within uncertainty. Melts in equilibrium with amphibole have estimated H₂O contents of ~4–9 wt% (Fig. 13d), towards the higher end of the typical range for arc magmas (Wallace, 2005; Plank et al., 2013). Assuming a crustal thickness of 25–43 km beneath Bagana (see above), we consider it most plausible that clinopyroxene and amphibole crystallization occurs in the mid-upper (3–21 km depth) crustal regions of the plumbing system.

Our thermobarometry and hygrometry results are consistent with experimental petrology studies of similar arc magmas. The lower calculated crystallization temperatures for amphibole compared with clinopyroxene are consistent with typical arc basalt crystallization sequences (Nandedkar et al., 2014; Marxer et al., 2021). Based on crystallization experiments on basaltic andesite (SiO₂ = 53 wt%, MgO = 4.2 wt%, melt H₂O 6.0–6.3 wt%) starting compositions from the Lesser Antilles arc, which are similar to the least evolved Bagana lavas (Fig. 7), a plagioclase (An_76–85)-clinopyroxene (Wo_40-43En_41-43Fs_16–18, Mg# 70–72)-amphibole ± magnetite assemblage is stable at ~990–940°C and 240–400 MPa (Pichavant et al., 2002; Melekhova et al., 2017), equivalent to depths of ~9–15 km. The compositions of the plagioclase produced in these experiments match Bagana population 1 (and some higher An population 2) plagioclase cores. The compositions of the experimental clinopyroxene and amphibole are also similar to clinopyroxene and amphibole found in our Bagana lavas. In these experimental studies, a clinopyroxene (Wo_44-45En₄₂Fs_13–14, Mg# 75–76) and high-An_78–86 plagioclase (± magnetite) assemblage, without amphibole, is stable at temperatures >1000 °C, in line with our predicted clinopyroxene crystallization temperatures.

Fractional crystallization experiments on a basaltic (SiO₂ = 49 wt%, MgO = 5.4 wt%, melt H₂O = 4.6 wt%) starting composition by Marxer et al. (2021) produced a plagioclase (An₈₇)-clinopyroxene-amphibole-magnetite assemblage, with mineral compositions closely matching the clinopyroxene, amphibole, and population 1 plagioclase cores in our Bagana lavas, at 990 °C and 400 MPa. Their experiments also produced a plagioclase (An₈₄, matching Bagana population 1 plagioclase cores) + amphibole + magnetite assemblage from a basaltic andesite (SiO₂ = 52.5 wt%, MgO = 2.9 wt%, melt H₂O = 6.6 wt%) starting composition at 930 °C, 400 MPa, which more closely matches our predicted amphibole crystallization temperatures (Fig. 13c).

Phase equilibria experiments on andesitic-dacitic (60–64 wt% SiO₂) starting materials showed that the assemblage plagioclase (An_45–62) + orthopyroxene + clinopyroxene + magnetite is stable at pressures <130 MPa, whereas amphibole is unstable and breaks down (Rutherford & Hill, 1993; Rutherford & Devine, 2003). Plagioclase compositions from these experiments closely match our population 1 and 2 plagioclase rims and our unzoned population 3 plagioclase crystals (Fig. 10a). This supports our suggestion that low-An plagioclase rims and population 3 plagioclase crystallized from evolved melts during later stages of magma ascent. Experimental constraints also suggest that orthopyroxene mesocrysts, which are in equilibrium with evolved melts represented by the AEMs (assuming K_D (Fe-Mg)^opx-liq = 0.28 ± 0.04, Beattie, 1993), are likely to have crystallized at pressures <130 MPa during magma ascent.

Our thermobarometry and hygrometry results, combined with experimental constraints, suggest that the assemblage high-An plagioclase (corresponding to cores of population 1 and some population 2 crystals), clinopyroxene, amphibole, and magnetite crystallized from parental basaltic or basaltic andesite (<53 wt% SiO₂) melts with high H₂O contents (>4 wt%). Our temperature estimates and the experimental constraints are consistent with initial crystallization of clinopyroxene and high-An plagioclase from the parental melts, with amphibole joining the assemblage upon cooling below 1000 °C. Our pressure estimates and comparison with experiments suggests that crystallization occurs predominantly at pressures of 130–570 MPa, corresponding to ~5–21 km depth (Fig. 13, Supplementary Fig. S8, with the 5 km upper limit based on the crystallization depth of the majority of clinopyroxene). The greater range of crystallization pressures, extending to shallower depths, recorded by clinopyroxene compared with amphibole most likely reflects the instability and breakdown of amphibole at lower pressures/depths. Hence, some clinopyroxene crystals (or clinopyroxene cores) could have crystallized at similar depths to amphibole, while others may have crystallized entirely (or continued to grow) at shallower levels of the plumbing system where amphibole is not stable. Low-An plagioclase rims, unzoned population 3 plagioclase crystals and orthopyroxene crystallized from evolved andesitic-dacitic melts (>60 wt% SiO₂) at low pressures <130 MPa, corresponding to depths <5 km, accompanying the formation of amphibole breakdown rims.

Mantle-derived magma supply to Bagana

We infer that the crystal cargo in Bagana lavas records magmatic processes in a mid-upper crustal storage region (~5–21 km depth), fed by parental magmas derived from either a lower crustal region (>21 km depth) of the plumbing system or directly from the mantle source. Our thermobarometric modelling finds no evidence for high-pressure crystallization and we lack samples of primitive magmas at Bagana, e.g. basaltic lavas, mafic enclaves. This leaves the early stages of magmatic evolution and the deepest parts of the plumbing system poorly constrained. We use our new radiogenic isotope and trace element data for Bagana lavas to infer the characteristics of parental magmas feeding the Bagana plumbing system, which has not been addressed by previous geochemical studies (Bultitude et al., 1978; Rogerson, 1989).

Bagana lavas have high Nb/Yb relative to N-MORB (Fig. 8e), indicative of an enriched mantle source composition (Pearce & Peate, 1995). High Th/Yb in arc lavas is interpreted to reflect the influence of sediment melts from the subducting slab, while high Ba/Th is indicative of a slab-derived aqueous fluid component (e.g. Woodhead et al., 2001; Pearce et al., 2005). Bagana lavas have both elevated Th/Yb and Ba/Th compared with MORB (Fig. 8d, e), suggesting that the mantle wedge beneath Bagana has been modified by both subducting sediment melt and slab fluid components.

Bagana Pb isotopic compositions plot in the Pacific mantle domain (after Kempton et al., 2002; Fig. 9). This Pb isotopic signature is consistent with the addition of slab fluid and sediment components to an Indian affinity mantle wedge with an isotopic composition similar to MANUS basin MORB, as suggested by Woodhead et al. (1998) for New Britain lavas. Alternatively, this Pb isotope signature may have been inherited from pre-20 Ma subduction of the Pacific plate along the Melanesian trench northeast of Bougainville. Bagana lavas have higher ⁸⁷Sr/⁸⁶Sr and lower ¹⁴³Nd/¹⁴⁴Nd than MANUS and Woodlark basin MORB and Solomon Sea plate basalts, consistent with the addition of high ⁸⁷Sr/⁸⁶Sr, low ¹⁴³Nd/¹⁴⁴Nd subducting Solomon Sea plate sediments (Fig. 9). The high bulk-rock δ¹⁸O of Bagana lavas (6.7–8.3‰) and ash (7.6‰) compared with MORB (5.7 ± 0.2‰, Bindeman, 2008) may reflect the addition of high-δ¹⁸O (>9‰) sedimentary material, either from the subducting slab or through crustal assimilation (e.g. Bindeman, 2008; Jeffery et al., 2013; Deegan et al., 2021). Given the lack of visible alteration on our lava samples, we do not consider these positively enriched δ¹⁸O values to reflect hydrothermal alteration (c.f. Heap et al., 2022).

Overall, our isotopic and trace element data suggest that parental magmas supplying Bagana are derived from an enriched mantle source modified by both slab fluids and subducted sediments. Our Bagana lava samples are relatively evolved, hence bulk-rock trace element ratios may not accurately reflect the mantle source composition due to modification by crystal fractionation or accumulation processes. In addition, we cannot rule out crustal assimilation without a more extensive isotopic dataset and information on the composition of the crust beneath Bougainville.

A new model for the Bagana magma plumbing system

Combining our observations from Bagana bulk rock and mineral chemical compositions, mineral textures, thermobarometry, and evidence from experimental petrology in other arc systems, we propose the following general model for the Bagana magma plumbing system (Fig. 14). Partial melts of an enriched mantle source, modified by slab fluid and subducted sediment melt components, enter the crust beneath Bagana. Basaltic to basaltic andesite (<53 wt% SiO₂) parental magmas with high H₂O contents (>4 wt%) feed a magma storage region at ~5–21 km depth, and crystallize an assemblage of clinopyroxene + high-An plagioclase ± amphibole ± magnetite over a temperature interval of ~1100–900 °C. The parental magmas may be derived from a lower crustal (>21 km depth) region of the plumbing system, where mantle melts are stored and undergo chemical evolution, though little evidence of this domain is preserved in erupted lavas with the possible exception of rare olivine and high Mg# (>80) clinopyroxene crystals.

Fig. 14

Schematic (not drawn to scale) illustrating a general model for the Bagana magma plumbing system, during persistent effusive activity as sampled by 2005–12 lavas.

The chemical composition of erupted Bagana lavas varies according to differences in carrier melt composition and the proportion and composition of entrained crystal cargoes. The andesitic to dacitic carrier melts, as recorded by the AEM compositions, are the residues of crystallization of the high-An plagioclase + clinopyroxene + amphibole + magnetite assemblage.

Individual crystals within the cargo record different growth, storage, and ascent histories in the plumbing system. Population 1 plagioclase crystals have high-An cores with sieve textures and narrow low-An rims. The sieve textures and low-An rims form during decompression-induced degassing of carrier melts during ascent from crustal storage towards the surface. Some population 2 plagioclase crystals also have high-An cores with sieve textures, but thick oscillatory-zoned mantles, indicative of stalling and further growth before final ascent. Other population 2 plagioclase crystals, which show oscillatory zoning but lack high-An cores, are inferred to crystallize entirely within this intermediate environment. The same processes (upward transport of cores and further growth during stalling, or crystallization entirely within an intermediate environment) could also account for the concentrically zoned clinopyroxene crystals. Population 2 plagioclase crystals with oscillatory zoned mantles and concentrically zoned clinopyroxene crystals with variable zoning patterns are suggestive of crystallization under changing thermochemical conditions in multiple different magma bodies. The planar euhedral morphology of these population 2 plagioclase and clinopyroxene crystals indicates that these magma bodies were melt-rich environments.

Crystals could feasibly aggregate to form glomerocrysts at any level of the plumbing system, during initial crystallization and storage, upward transport in a common carrier melt, stalling, or final ascent to the surface. Therefore, adjacent crystals within glomerocrysts may have grown from the same melt in a given region of the plumbing system, or formed from different melts in different parts of the plumbing system and become aggregated during transport towards the surface. The presence of some glomerocrysts with textures typical of partially to fully solidified crystal mushes among the crystal cargo in Bagana lavas suggests that melt-rich portions of the plumbing system are lined by or separated by mushy regions, which can be disrupted by ascending magmas.

Andesitic-dacitic carrier melts ascend from the ~5–21 km depth storage region towards the surface, entraining varied crystal cargoes by interacting with melt-rich and mushy magma bodies en-route. During ascent, degassing-induced crystallization results in the growth of lower-An plagioclase crystals (population 3) from carrier melts alongside population 1 and 2 crystals’ lower-An rims. Experimental constraints suggest that crystallization of low-An plagioclase (+ orthopyroxene + additional clinopyroxene growth) and breakdown of amphibole are most likely to occur between ~5 km depth and the surface. Thick amphibole breakdown rims suggest that magmas spend tens of days outside the amphibole stability field prior to eruption. It is possible that crystals may be deposited on conduit walls during ascent, then scraped off and transported to the surface by subsequent batches of magma, adding to the diversity of crystal cargoes found in individual lava flows.

The characteristic long-term eruptive behaviour of Bagana, as observed in the 2000–2014 interval sampled by this study, comprises prolonged effusion of lava, and can plausibly be sustained by a steady state regime where magmas near-continually ascend from the crustal plumbing system. The scarcity of high frequency volcano-tectonic earthquakes detected at Bagana by RVO seismometers lends support to this idea. The continual degassing of rising magmas, and/or those residing deeper in the plumbing system, can explain Bagana’s persistent gas emissions (McCormick Kilbride et al., 2019).

The processes causing rare transitions to explosive behaviour at Bagana are poorly constrained. One possibility is plugging of the conduit, which impedes degassing and leads to build up of overpressure, resulting in explosions (e.g, Stix et al., 1993; Preece et al., 2016). Gas escape may also be suppressed by mineralisation sealing fractures in the summit dome (e.g. Heap et al., 2021). Alternatively, explosive eruptions may be fuelled by changes in the volatile supply from deeper parts of the plumbing system to the uppermost levels (e.g. Liu et al., 2020; Edmonds et al., 2022). Our 2016 ash samples have more evolved bulk-rock compositions than Bagana lavas, allowing us to tentatively suggest that explosive eruptions may be supplied by higher bulk SiO₂, or possibly more melt-rich/crystal poor magmas. Petrographic and geochemical analyses of explosive eruption products could provide valuable insights into the state of the plumbing system during these events and the processes controlling transitions in activity style at Bagana.

What type of volcano is Bagana?

Wadge et al. (2018) described Bagana as a ‘rare example of a very youthful, polygenetic, andesite volcano’. The young age of Bagana is based on stratigraphic relationships and estimates of mean extrusion rate derived from a series of digital elevation models constructed over different time intervals (Wadge et al., 2018). We cannot assess the potential youth of Bagana, though the textural and chemical data we present herein could form the basis for future studies of reservoir residence times, crystal growth rates, and magma ascent timescales. However, our work does support the idea of Bagana as a polygenetic volcano. Complex crystal cargoes, a range, albeit limited, of whole-rock and potential melt compositions, and mineral crystallization throughout a vertically extensive storage region are all consistent with a more complex plumbing system than is associated with low-volume monogenetic volcanism (Smith & Németh, 2017). Satellite observations of inter-eruptive pauses of comparable duration to Bagana’s eruptions (McCormick Kilbride et al., 2019) and dramatic shifts to explosive eruptions (Global Volcanism Program, 2023b, 2023c) are also consistent with a polygenetic nature.

If Bagana has grown within centuries (<500 years, Wadge et al., 2018), its plumbing system architecture (Fig. 14) and patterns of eruptive activity may be characteristic of the early stages of the lifecycle of composite arc volcanoes. However, it is unclear if Bagana constitutes an independent magmatic and volcanic system or is genetically related to the neighbouring volcano Billy Mitchell (see Figs. 1 and 2), an andesitic to dacitic pyroclastic shield volcano (Global Volcanism Program, 2023d). The last eruption of Billy Mitchell was the Plinian BM2 event, dated to 370 ± 19 years BP, which generated ~3.5 km³ of airfall tephra and a 10 km³ ignimbrite (Mckee et al., 1990). Bagana deposits consistently lie on top of BM2 deposits and these stratigraphic relationships raise the possibility that Bagana is a parasitic cone of the Billy Mitchell system. Caldera collapse after the BM2 event may have suppressed magma ascent directly beneath Billy Mitchell, with the younger volcanism that built the Bagana edifice being fed by an inclined conduit, which reached the surface 6 km to the southwest (Wadge et al., 2018).

There is no direct seismic or geodetic evidence to support the existence of such a conduit, but the notion is plausible. Many arc volcano edifices are laterally offset relative to the centroid location of their subvolcanic plumbing systems (Ebmeier et al., 2018; Lerner et al., 2020). Intriguingly, other young polygenetic mafic volcanoes worldwide appear to be parasitic cones related to larger or more well-established magmatic systems. Ngauruhoe is a 900 m-high basaltic andesite cone built of explosive eruption products and highly porphyritic lavas over ~2500 years, the youngest of nine cones in the Tongariro complex, New Zealand (Hobden et al., 2002). Cerro Negro is a small basaltic cone on the northwest flank of the El Hoyo-Las Pilas complex in Nicaragua, built from ash, scoria and lava flows in ~150 years (McKnight & Williams, 1997). Volcán de Izalco is a 650 m-high basaltic-andesitic stratovolcano composed of tephra deposits and lava flows, 4 km to the south of the Santa Ana complex in El Salvador and has grown in 225 years (Carr & Pontier, 1981). Izalco lies trenchward of the parent system Santa Ana, a spatial relationship repeated in three other closely spaced (<5 km) pairs of volcanoes in the Central American volcanic arc (Cerro Quemado–Santa Maria, Toliman–Atitlan, Acatenango–Fuego, Halsor & Rose, 1988). These volcanoes share some, but not necessarily all, of the following features: narrow, largely mafic compositional ranges; complex crystal cargoes indicating disequilibrium with the carrier melt; rapid edifice building in hundreds to a few thousand years; a substantial volume fraction of their deposits being blocky lava flows; and an association with a nearby more established magmatic system. Future work investigating a link between Bagana and Billy Mitchell, whether based on petrological and geochemical studies of their products or geophysical detection of a shared plumbing system, could be valuably extended by a systematic global analysis of young, mafic, polygenetic volcanoes and their parent systems to evaluate whether these relatively rare locations can provide useful insights into the long-term evolution of magmatic plumbing systems, edifice construction and volcanic hazards in arc settings.

CONCLUSIONS

Our bulk-rock and mineral chemical data for 2005–16 Bagana lavas and ash provide important new insights into sub-volcanic processes and magma storage conditions at this persistently active but little studied volcano. The bulk rock compositions of Bagana lavas are controlled by variations in composition of the andesitic-dacitic carrier melts and proportion and composition of the plagioclase + clinopyroxene + amphibole + orthopyroxene + magnetite crystal cargo.

Incompatible trace element and radiogenic isotopic compositions of Bagana lavas suggest that parental magmas are derived from an enriched mantle source modified by both slab fluid and subducted sediment melt components. Constraints from thermobarometry, hygrometry and experimental petrology suggest that H₂O-rich (>4 wt%) basaltic and basaltic andesite (<53 wt% SiO₂) parental magmas supply a magma storage region at ~5–21 km depth beneath Bagana. Parental magmas crystallize an assemblage of high-An plagioclase + clinopyroxene ± amphibole ± magnetite over a temperature interval of ~1100–900 °C, producing andesitic-dacitic residual melts. The plumbing system consists of multiple melt-rich magma bodies, where crystals can grow under changing thermochemical conditions, in addition to mushy regions.

Andesitic-dacitic carrier melts segregate and ascend from the ~5–21 km depth storage region towards the surface, variably sampling melt-rich magma bodies and mushes en-route, entraining a diverse crystal cargo. Degassing during magma ascent reduces carrier melt H₂O contents, leading to crystallization of low-An plagioclase and the formation of breakdown rims on amphibole, most likely between ~5 km depth and the surface. Breakdown rims on amphibole of >50 μm thickness imply that magmas spend tens of days at pressures outside the amphibole stability field prior to eruption. Our findings suggest that typical activity at Bagana, which comprises prolonged effusion of lava over months to years, can be sustained by near-continuous ascent, accompanied by degassing, of magmas from the crustal plumbing system.

The impacts of the July 2023 explosions highlight the need for improved monitoring of this persistently active volcano to protect local communities. Through identifying magma storage depths and pre-eruptive processes, our study provides a valuable framework for interpretation of future monitoring data. Future petrographic and geochemical studies of explosive eruption products could be compared with our findings from Bagana lavas to identify potential differences in magma storage conditions or changes in pre-eruptive processes during less frequent explosive activity. Our characterization of the Bagana plumbing system provides a foundation for future work to assess any potential link to the nearby caldera-forming Billy Mitchell volcano, which has important implications for understanding future volcanic hazards on Bougainville and more generally the evolution and growth of crustal magmatic plumbing systems and arc volcano edifices.

ACKNOWLEDGEMENTS

We thank Jon Fellowes for his assistance with EPMA analyses at the University of Manchester. Chris Harris kindly provided oxygen isotope analyses. We thank Wally Johnson, Jon Woodhead and Janet Hergt for numerous useful insights into volcanism on Bougainville and across the wider region. Finally, we are deeply grateful to the communities of Gotana, Torokina, and Patsikopa for their warm hospitality and invaluable support during challenging fieldwork and sample collection on Bougainville. We would like to thank two anonymous reviewers and editor Adam Kent for their thorough reviews and insightful comments that significantly improved this manuscript.

FUNDING

This work was supported by the Alfred P. Sloan Foundation’s support of the Deep Carbon Observatory Deep Earth Carbon Degassing program (DECADE); the University of Manchester; the Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics (COMET); Clare College, University of Cambridge. DAN acknowledges support from a NERC IRF (NE/T011106/1).

DATA AVAILABILITY

The data underlying this article are available in the article and in its online supplementary material.