Magmatic Processes at La Soufrière de Guadeloupe: Insights From Crystal Studies and Diffusion Timescales for Eruption Onset

Introduction

Pumice and scoria contain crystals that record magmatic processes at depth and provide information on the magma plumbing system conditions prior to a volcanic eruption, including pressure, temperature, fO₂ and timescales of magmatic processes (e.g., Cashman et al., 2017; Cooper, 2019; Costa et al., 2020). Many processes can trigger an eruption including mafic intrusions, magma mixing and differentiation and volatile transport (e.g., Cooper 2019 and references therein). Eruption triggers are often recorded by crystals in the reservoir through re-equilibration reactions with the intruding magma, resulting in diffusion gradients visible as compositional zoning (Druitt et al., 2012; Kilgour et al., 2014; Albert et al., 2019; Costa et al., 2020). Focusing on the outermost crystal zone gives an indication of the time elapsed between the most recent magma input and the subsequent eruption, which quenches the magma and effectively stops diffusion processes. Diffusion studies using olivine, plagioclase and pyroxene can be used to derive a wide range of timescales (i.e., days to millennia), across a range of volcano types and compositions, relating to a range of magmatic processes (e.g., Morgan et al., 2004; Costa et al., 2008; Chamberlain et al., 2014; Kilgour et al., 2014; Dohmen et al., 2017).

Magma reservoirs are now often interpreted as long-lived, relatively cold storage zones with mush lenses comprised of melt and crystals formed as a result of repeated intrusions from depth (e.g., Mahood, 1990; Nakamura, 1995; Bachmann et al., 2002; Zellmer et al., 2003; Moretti et al., 2013; Brown et al., 2014; Cooper and Kent 2014; Annen et al., 2015; Bergantz et al., 2015; Edmonds et al., 2016; Bergantz et al., 2017; Cashman et al., 2017; Spera and Bohrson, 2018; Carrara et al., 2019; Moretti et al., 2019). Crystals form a continuous interlocking framework due to close crystal contacts through which melt is distributed in varying volumetric proportions. Melt and/or fluid migration, amalgamation and mixing makes the system inherently unstable allowing remobilization (Costa et al., 2008; Mader et al., 2013; Cooper and Kent, 2014; Rubin et al., 2017; Sparks and Cashman, 2017; Jackson et al., 2018; Cooper 2019).

An eruption is an event resulting from a sequence of magmatic processes which form eruptible magma, and occur on different timescales, such as: magma genesis, ascent and storage, fractional crystallization, volatile infiltration and degassing, magma mixing and assimilation (Moretti et al., 2013; Edmonds et al., 2016; Carrara et al., 2019). Though crystal-mush disaggregation may occur as a gradual process, numerical modeling suggests continued injection leads to mush collapse involving small areas of the magma reservoir, due to intrusion geometry controlled by density differences between the intruding magma and host mush (Bergantz et al., 2015; Schleicher et al., 2016; Schleicher and Bergantz 2017; Bergantz et al., 2017; Carrara et al., 2019). This means not all the crystals interact with the same amount of intruding magma, with some crystals experiencing multiple recharge events before the system is fully ‘unlocked’ and destabilized prior to an eruption (Cheng et al., 2020).

Magmatic crystals act as effective time capsules containing information key to understanding the mush system structure, recharge, mush remobilization, and destabilization timescales and dynamics (e.g., Costa and Dungan 2005; Davidson et al., 2007; Ubide and Kamber, 2018; Costa et al., 2020). These processes change the chemical environment within which crystals grow, and so disturb the system’s equilibria leading to chemical gradient relaxation in the crystals by element diffusion as the system re-equilibrates (Morgan et al., 2004; Druitt et al., 2012; Barker et al., 2016). Investigation of disequilibria textures and kinetics is used to calculate the rate at which a system departs from equilibrium and reaches a new equilibrium (Costa and Dungan 2005; Costa et al., 2008; Saunders et al., 2012; Kilgour et al., 2014; Costa et al., 2020). This technique, termed diffusion chronometry, allows us to link the history and timing of magmatic processes to our current understanding of the magmatic system mush system and future periods of volcanic unrest (Kilgour et al., 2014; Pankhurst et al., 2018; Costa et al., 2020). Relating diffusion timescales to different eruption styles may allow us to improving hazard assessment and eruption forecasting for different eruption styles by interpreting the diffusion timescales in terms of periods of volcanic unrest.

The wealth of high-resolution diffusion timescale data acquired over the past decade for a variety of volcanic systems shows different timescale ranges are recorded by different minerals, with olivine, often used for basaltic systems, recording the shortest timescales and pyroxene and plagioclase recording longer timescales used for more silicic systems (Barker et al., 2016; Hartley et al., 2016; Flaherty et al., 2018; Costa et al., 2020). Orthopyroxene in particular is used over a broad range of temperature conditions (750–1,100°C; Dohmen et al., 2016), has a well constrained diffusion coefficient (D_Fe-Mg) and allows calculation of timescales on both a short and long scale, ranging from <1 day to 100,000’s of days (Dohmen et al., 2017; Costa et al., 2020).

In this paper we present a study on orthopyroxene diffusion chronometry to investigate magmatic process history and timing in a basaltic-andesitic, subduction zone setting which is characterized by a variety of eruption styles. But do the different eruption styles relate to the same mush system disaggregation rate. To address this, we study four La Soufrière de Guadeloupe explosive basaltic-andesite to andesitic eruptions: 1657 ± 20 Cal. CE, 1010 ± 10 Cal. CE, ∼341 Cal. CE and 5680 ± 164 Cal. BCE. These eruptions range from the La Soufrière de Guadeloupe’s first (5680 Cal. BCE) to most recent (1657 Cal. CE) explosive magmatic eruptions, covering a variety of eruption styles (Plinian, Vulcanian and Strombolian; Komorowski et al., 2005; Legendre, 2012; Komorowski et al., 2012; Komorowski et al., 2013).

La Soufrière de Guadeloupe is an ideal volcano to explore the correlation between mush disaggregation timescales and eruption size and style. This stratovolcano has a stable bulk composition (basalt to basaltic-andesite) of magma through the disparate events studied, making comparisons between eruptions and the magmatic processes preceding them possible (Komorowski et al., 2005; Legendre, 2012). Finally, the key hypothesis and conceptual model that La Soufrière de Guadeloupe is underlain by a mush system based on large variability in crystal populations, lack of extensive and systematic magma mixing textures and no clear time evolution in eruption product composition, allows us to link timescales to magmatic processes (Moretti et al., 2020).

Our study is timely due to recent, increasing unrest over the last 25 years with a significant risk of an eruption at La Soufrière de Guadeloupe. Approximately 70,000 people live and work in Southern Basse-Terre in the path of previous pyroclastic density currents and debris flows (Komorowski, 2008; Leone et al., 2019). This highlights the need to fully understand the magmatic system, including the timescales of magmatic processes related to unrest before eruption in order to improve early-warning systems, forecast models, crisis response management for different scenarios, and long-term risk reduction planning.

Geological Setting and Eruptive History

Guadeloupe is composed of several islands in the central Lesser Antilles arc. The arc formed due to the westward subduction of the Atlantic plate beneath the Caribbean plate at ∼2 cm/year (Feuillet et al., 2002, Feuillet et al., 2011; Figure 1A). The eastern islands (Grand-Terre and Marie-Galante) are associated with the external, older arc and the western island (Basse-Terre) is associated with the active, inner arc (Macdonald et al., 2000; Komorowski et al., 2005; Boudon et al., 2008; Kopp et al., 2011; Figure 1B). La Soufrière de Guadeloupe volcano is located in Southern Basse-Terre which is composed of seven volcanic complexes built over the last 3 Ma and shows a temporal evolution from North to South (Komorowski et al., 2005; Samper et al., 2007; Boudon et al., 2008).

FIGURE 1. (A) A schematic map showing the Lesser Antilles in relation to the subduction zone position where the Atlantic plate moves under the Caribbean plate. Guadeloupe is located in the in the central Lesser Antilles. (Adapted from Allard et al., 2014; Pichavant et al., 2018). (B) Schematic map of Guadeloupe which is made up of two islands: Grand-Terre associated with the older, inactive, outer arc and Basse-Terre associated with the active inner arc. La Soufrière de Guadeloupe, the current active vent, is located in south Basse-Terre. (Adapted from Allard et al. (2014) and Pichavant et al. (2018). (C) Annotated photo cross-section (Photo: Metcalfe, 2020) showing the key features of La Soufrière and the surrounding vents. This highlights the composite nature of La Soufrière de Guadeloupe, with some of the different sequences of lava flows and collapses labeled (Blanc, 1983; Carlut et al., 2000; Komorowski et al., 2005; Samper et al., 2009). The current dome is associated with the 1530 Cal. CE eruption crater.

La Soufrière de Guadeloupe is the most recent edifice of the Grande Découverte-Soufrière (GDS) volcanic complex, which was built during the last 0.455 million years over the older Sans Toucher composite volcano’s southern flanks (∼0.8–0.5 Ma) (Boudon et al., 1988; Komorowski et al., 2005; Samper et al., 2009; Legendre, 2012; Figures 1B,C). The GDS complex is composed of three composite volcanoes, Grande Découverte (active from 0.455 Ma to 44330 Cal. BCE), Carmichaël (active from 44330 Cal. BCE to 10870 Cal. BCE) and La Soufrière de Guadeloupe (active from ∼7140 Cal. BCE to present).

La Soufrière de Guadeloupe is a basaltic-andesitic system and is characterized by multi-phase activity including lava flows, domes, scoria cones, explosive eruptions with sustained and collapsing plumes, and explosive hydrothermal and phreatic activity (Legendre, 2012; Komorowski et al., 2012; Komorowski et al., 2013; Figure 1C). Legendre (2012) reports numerous, previously unknown eruptions (labeled GDS 1 to GDS 15), including five Plinian explosive eruptions (GDS 15, GDS 14, GDS 9, GDS 8, GDS 3). Here we provide the eruptive history necessary to contextualize our results (Table 1). Detailed stratigraphy, radiocarbon dating, and eruptive deposit analysis are available in Legendre (2012).

TABLE 1. Synthesis of the eruption record for La Soufrière of Guadeloupe in the last 9,160 years adapted from Komorowski et al. (2005), Legendre (2012), Komorowski et al. (2012), Komorowski et al. (2013) where more details on carbon dating is available.

The proto-Soufrière edifice, built from lava flows and dome-building eruptions, partially collapsed at 7140 ± 18 years Cal. BCE (∼9,160 years ago) producing a basal debris avalanche deposit, the first volcanic deposit of La Soufrière de Guadeloupe (Komorowski et al., 2005; Boudon et al., 2007; Komorowski, 2008; Legendre, 2012). Following this collapse, effusive activity resumed with lava flow and dome extrusion dated via K-Ar record an age of ∼6 ± 2 kyr, though the flows and domes are not well sampled/preserved (Samper et al., 2009; Figure 1C). Activity culminated at 5680 ± 164 Cal. BCE with the first Plinian eruption preserved in the geological record (GDS 15).

The 5680 Cal. BCE eruption generated a sustained eruption column reaching an estimated height of 10–25 km (deposits are poorly preserved), followed by plume collapse resulting in pyroclastic density currents (Legendre, 2012; Komorowski et al., 2012; Komorowski et al., 2013). Basaltic scoria deposits are found 7 km from the vent (Legendre, 2012). This is the most violent magmatic eruption from this system, estimated at Volcanic Explosivity Index (VEI) 4 and with a minimum estimated erupted volume of ∼0.1 km³. This was followed by a similar sized Plinian (GDS 14) dated at 5387 ± 84 Cal. BCE (Legendre, 2012; Komorowski et al., 2012; Komorowski et al., 2013).

Following this explosive period, an effusive, dome-building period began with eruptions associated with partial edifice collapse, debris-avalanches and violent laterally directed dome or cryptodome explosions. Four main events are observed in this period: 4366 ± 134 Cal. BCE (GDS 13), 3300 ± 110 Cal. BCE (GDS 12), 1870 ± 71 Cal. BCE (GDS 11), and 1370 ± 100 Cal. BCE (GDS 10). During this dome building period, no explosive magmatic events are identified due to an incomplete stratigraphic record. This is a common problem in tropical settings such as Guadeloupe, where the climate generates a high mechanical weathering rate between 800–4,000 t/km²/year (Rad et al., 2013).

Explosive magmatic activity resumed in 1080 ± 10 Cal. BCE with two events interpreted as VEI 4 eruptions and minimum erupted volumes of ∼0.1 km³ (Table 1; GDS 9 and GDS 7). Following these events, between 720 ± 10 Cal. BCE and 450 ± 10 Cal. CE, explosive magmatic activity migrated 0.8–1.4 km SE from the main La Soufrière de Guadeloupe vent, forming the older Echelle (GDS 6) and younger Citerne (GDS 5) scoria cones (Figure 1C). Boudon et al. 1988 and Komorowski et al. 2005 report an age of ∼341 Cal. CE for laharic debris-flow deposits reworking the Echelle deposits shortly after their emplacement. These two Strombolian-style eruptions are interpreted as VEI 2 events, with minimum erupted volumes of 0.001 km³ and an estimated erupted column height of ∼1–15 km. The explosive and effusive erupted products from these two vents are basaltic compared to the andesitic products at the La Soufrière de Guadeloupe vent. Echelle lava flows show evidence of magma mingling textures involving a more silicic magma (Boudon et al., 1988).

By 450 ± 10 Cal. CE magmatic activity had migrated back to La Soufrière de Guadeloupe beginning with an effusive phase (GDS 4). The following event in 1010 ± 10 Cal. CE is the most recent Plinian eruption (GDS 3). This eruption is categorized as a VEI 4 explosive eruption, with a minimum estimated eruptive volume of ∼0.1 km³ and an estimated eruptive column height of 10–25 km, which dispersed eruptive products 7 km south-west from the vent. This explosive eruption was associated with a dome-building phase that produced a laterally-directed explosion (blast), high-energy and dome collapse pyroclastic density currents.

The most well-studied and archetypal eruption is the 1530 Cal. CE eruption (GDS 2) which began with partial edifice collapse, followed by a paroxysmal sub-Plinian phase and finished with emplacement of the current lava dome (Komorowski et al., 2005; Boudon et al., 2008, Komorowski, 2008; Legendre, 2012; Figure 1C). The 1530 Cal. CE eruption is interpreted as a VEI 3 sub-Plinian eruption with a short-lived column reaching, which transported juvenile material, and associated column-collapse and scoriaceous pyroclastic density-currents up to 6 km from the vent (Komorowski, 2008; Komorowski et al., 2012; Komorowski et al., 2013). The combined pyroclastic DRE erupted volume is ∼0.06 km³ (Boudon et al., 2008; Komorowski, 2008; Legendre, 2012; Esposti Ongaro, 2020). The 1530 Cal CE eruption is used as a calibrated scenario for a future magmatic eruption, as it is representative of La Soufrière de Guadeloupe’s complete eruptive past (excluding the directed blast phenomena).

The 1530 Cal. CE eruption was, until recently, known as the most recent magmatic eruption (Komorowski et al., 2005; Boudon et al., 2008; Komorowski, 2008; Pichavant et al., 2018). However, new field data indicate a more recent, small magmatic eruption (GDS 1) occurred in 1657 ± 20 Cal. CE (Legendre, 2012; Komorowski et al., 2012; Komorowski et al., 2013). This eruption was a small, explosive Vulcanian eruption (VEI 2-3), with a minimum estimated erupted volume of ∼0.001 km³.

Six phreatic/hydrothermal explosions have occurred since 1657 Cal. CE (1690, 1797–1798, 1809–1812; 1836–1837, 1956 and 1976–1977; Feuillard et al., 1983; Komorowski et al., 2005; Hincks et al., 2014). The 1976–1977 phreatic eruption resulted in the establishment of a monitoring system including a seismic network, global navigation satellite system receivers (GNSS), meteorological networks, periodic extensometric surveys of lava dome fractures, fumarole gas and thermal spring sampling, and thermal monitoring (Komorowski et al., 2005; Villemant et al., 2005; Villemant et al., 2014; Jessop et al., 2020; Moretti et al., 2020).

The monitoring data following the 1976–1977 phreatic eruption, show that the level of degassing and seismic activity decreased, and a quiescent period began by 1984 lasting until May 1992 with the start of seismic swarms and a new degassing phase at summit fumaroles (Zlotnicki et al., 1992; Komorowski et al., 2001; Komorowski et al., 2005; OVSG-IPGP, 1999–2020). Degassing has continued to evolve with pre-existing fumaroles reactivating and new fumarolic vents appearing. Notably, late 1997 and early 1998 marked the onset of highly acidic chlorine-rich fluids degassing, continuing till present day, along with the hydrothermal fluid C, He and Cl isotopic ratios. This indicates persistent magma reservoir degassing, supplying fluids and heat to the shallow hydrothermal system (Komorowski et al., 2005; Ruzié et al., 2012; Ruzié et al., 2013; Villemant et al., 2014; Moretti et al., 2020; OVSG-IPGP, 1999–2020). Since 1992 there have been twenty felt volcanic earthquakes, with ten of these occurring since 2013, culminating in 2018 with the strongest recorded event since 1976 (M 4.1; Komorowski et al., 2005; OVSG-IPGP, 1999–2020; Moretti et al., 2020). Recent years have also seen an acceleration in the opening rate of several summit fractures and the area of heated ground at the summit increasing in size (Jessop et al., 2020). Passive ground degassing (CO₂) has occurred in areas near the summit vent, as well as on the flanks of the dome (OVSG-IPGP, 1999-2020; Moretti et al., 2020; Jessop et al., 2020).

The Magmatic System at La Soufrière de Guadeloupe

In an earlier study, Touboul et al. (2007) proposes a reservoir structure for La Soufrière de Guadeloupe with a shallower crystallized intrusion and a deeper zoned magma chamber, resulting from multiple magmatic recharge events and fractional crystallization. This zoned magma reservoir would account for heterogeneities such as mingling textures observed in the 1530 Cal. CE eruption. However, these heterogeneities are not consistent through the eruption series studied. Further to this, Touboul et al. 2007 observed that, despite showing large scale mingling textures, the 1530 Cal. CE eruption shows mineral phases which are relatively compositionally homogeneous, which is hard to resolve in a single, zoned magma reservoir.

Fractional crystallization is interpreted as the dominant process in a magma reservoir, which should lead to continuous melt composition at the surface (Bonnefoi et al., 1995; Geist et al., 1995; Peccerillo et al., 2003; Macdonald et al., 2008; Dufek and Bachmann 2010; Sliwinski et al., 2017). Eruption products may also be expected to reflect the magma chamber’s temporal evolution through fractional crystallization. However, though the compositional trends observed in the La Soufrière de Guadeloupe suite do most likely represent a degree of fractional crystallization, the youngest eruption products are not always the most evolved.

Large magma bodies can cool rapidly if the intruding, recharge magma mixes with cooler stored melts, assimilates with the host rock or is connected to a developed hydrothermal system. Maintaining a single magma reservoir feeding the long-lived Grande Découverte-Soufrière complex would require high average rates of magma injections to balance heat loss to the crust (Bachmann and Huber 2016; Degruyter et al., 2016; Cashman et al., 2017). However, it is also important to consider the response of the magma chamber to repeated magma recharges and eruptions as the frequency and duration have a huge effect on how the magma chamber cools (Degruyter et al., 2016).

More recently, results from petrological modeling of the 1530 Cal. CE eruption products indicate mafic magma intrudes into a shallow-depth, andesitic magma storage zone. This zone extends from a minimum depth of 5–6 km–8.5 km; this configuration allows efficient magma mixing and hybrid magma genesis (Pichavant et al., 2018). However, the current conceptual model proposed by Moretti et al. (2020) can better resolve the eruption product features and timescales of magmatic processes observed, based on the following evidence: the lack of extensive and systematic magma mixing textures and deeper mafic enclaves in the majority of eruptions (excluding the ∼341 Cal. CE and 1530 Cal. CE eruptions), the crystal populations record a large variability in chemical perturbations shown in mineral zonation patterns, and the short diffusion timescales and eruption product composition show no clear time evolution (Figure 2). This system has most likely varied through time, with the variability in the proportion of eruptible melt explaining why the magmatic eruptions of the last 10,000 years are characterized by limited erupted volumes. Though the mush system extent is unclear, the crystal, glass and whole rock compositions for Echelle and La Soufrière de Guadeloupe lie on the same compositional trend. This indicates the La Soufrière de Guadeloupe mush system could extend beneath the secondary, adventive Echelle and Citerne monogenetic cones (Boudon et al., 1988; Komorowski et al., 2005; Legendre, 2012).

FIGURE 2. Schematic cross section depicting the mush system beneath La Soufrière de Guadeloupe (adapted from Cashman et al., 2017; Edmonds et al., 2016; Moretti et al., 2020, pressures and temperatures from Pichavant et al., 2018). Based on bulk rock compositions, the crystal mush zone is hypothesized to extend from La Soufrière de Guadeloupe to the Echelle and Citerne scoria cones. Basaltic magma is injected from depth into the system at temperatures ranging from 975 to 1,025°C (Pichavant et al., 2018), this magma interacts with the deep mush system, allowing heat and fluids to ascend to shallow depths. Recharge events mainly affect the deep parts of the mush system, where regular magma interaction is recorded by zoning and compositional changes in the crystals. The mush system is comprised of crystal-rich and melt-rich lenses. Melt can episodically move up in the system leading to mixing and amalgamation, resulting in lenses of varying compositions and crystal populations through the system (Sparks and Cashman, 2017). Crystal and melt evolution in the mush zone creates compositional heterogeneity. Schematic crystals are shown with their zoning patterns and amount of resorption depending on the frequency and volume of magma interaction. Important processes occurring in the La Soufrière de Guadeloupe mush system include mafic recharge, magma mixing, fractional crystallization and the ascent of heat, gas and fluids from depth into the hydrothermal system. The mush system and the extensively developed hydrothermal system allow buffering of processes that lead to eruptions, such as recharge and mixing, as shown by the 2018 period of unrest, until recharge reaches a critical threshold.

The most recent model, proposed by Pichavant et al. 2018, is a heterogenous crystal mush system extending from depths of 5–8.5 km, based on melt volatile contents and the absence of amphibole, and is considered the source of heat and deep magmatic fluids rising into the shallow hydrothermal system. Heterogeneously distributed magma bodies make-up a mush system, although largely crystalline. They are also comprised of varying proportions of melt, interlocking crystals and exsolved volatiles that form different horizontal sill-like lenses, connected vertically by dykes (Bachmann and Bergantz, 2004; Hildreth, 2004; Huber et al., 2009; Bachmann and Huber, 2016; Sparks and Cashman, 2017). In any of the mush lenses, only a small volume percent of the melt is eruptible. Volcanic behavior is controlled by the whole transcrustal mush system which is laterally or vertically distributed (Cashman et al., 2017; Sparks and Cashman, 2017).

The mush system is inherently unstable due to compaction and shear leading to melt-rich layer segregation (Sparks and Cashman, 2017). This promotes melt migration between mush lenses and up through the system, leading to mixing, amalgamation and compositional evolution as the melt reacts with progressively cooler mush (Sparks and Cashman, 2017; Jackson et al., 2018). The system that comprises zones characterized by different degrees of cooling, rheology and volumetric melt fraction, is likely to experience regular intrusions from a deep, less evolved hotter source of varying size, location, composition, and rheology. Variation in the way the deep intrusions interact with the mush system, and its different zones, results in a heterogeneous composition through the system and a variation in the timescales recorded by crystals residing in it (Figure 2). As described above, less evolved hot magma intruding from depth is assumed to result in mush remobilization. However, given that the presence of fluids in magmatic systems lowers the melting temperature of rock favoring melt formation, Moretti et al. 2013 suggested that large amounts of fluid in a system allow crystal mushes to remobilize without the need for a less evolved hot magma intruding from depth. This is an important point to consider for the water-rich magmas of La Soufrière of Guadeloupe.

Methods

Sample Collection and Preparation

The samples used for this study were collected between 2001 and 2019, from both distal and proximal outcrops around La Soufrière de Guadeloupe (Figure 3; Table 2), and are described fully, with stratigraphy and carbon dates in Komorowski et al. (2005, Komorowski, 2008, Boudon et al. 2008 and Legendre (2012).

FIGURE 3. Schematic map of the south of Basse-Terre, showing the main volcanic complexes present (Adapted from Komorowski et al., 2005). Sampling locations are marked and correspond to the stratigraphic logs, excluding the ∼341 Cal. CE eruption (C). A mix of proximal and distal outcrops were sampled, depending on the quality and accessibility of the outcrop. Stratigraphic logs show the range of erupted products recorded at each location, in each case the sample used is labeled as well as any dated palaeosol horizons. (A) – 1657 Cal. CE (SOU 0913D) sampled from a fallout unit found between ash layers in Ravine Goyavier. This sample was found above a palaeosol horizon dated at 245 ± 30 years BP. (B) – 1010 Cal. CE (SOU 0111A and SOU 0924A) sampled from pumice fallout units found in sequences of ash and pumice flow units in Ravine Bonne Espérance. (C) - ∼341 Cal. CE (SOU 1904A) was sampled from a slump off the Echelle monogenetic cone, so no stratigraphic position is available. (D) – 5680 Cal. BCE (SOU 1101viA) sampled from a scoria fall out unit in a sequence of pumice and scoria fallout units at Chantier du SDIDS. This sample was taken from directly below a palaeosol horizon dated at 6380 ± 40 years BP. [Adapted from Komorowski et al. (2005) and Legendre (2012)].

TABLE 2. Summary table for the eruptions showing the key features for each eruption and sample location.

The 1657 Cal. CE eruption (GDS 1) was sampled in 2009 and 2013 from a proximal 15 and 20 cm thick fallout unit less than 1 km from La Soufrière de Guadeloupe (Figure 3, sampling location A; Table 2). The fallout unit was observed in the middle of a sequence of ash layers directly above a palaeosol unit, carbon dated at 245 ± 30 years BP (1656 ± 26 years Cal CE; Legendre, 2012). The light gray, dense andesite pumice is rich in plagioclase and pyroxene phenocrysts and was sampled across the layer at constant depth into the unit to ensure the sample is representative of the whole unit.

The 1010 Cal. CE eruption (GDS 3) was sampled in 2001 and 2009 from distal pyroclastic sequences, including fall units and column collapse pumiceous pyroclastic density-currents deposits (Figure 3, sampling location B; Table 2). Both units were sampled across the entire layer at constant depth. These samples are light beige to light gray andesitic pumice and are rich in plagioclase and pyroxene.

The ∼341 Cal. CE eruption (GDS 5) was sampled in 2019 from a proximal outcrop near the Echelle cone. Individual black, crystal rich, scoria clasts, 10–40 cm in size were picked from a recent landslide on the road that crosses the Echelle scoria cone not far from its base, allowing collection of the freshest scoria clasts (Figure 3, sampling location C; Table 2). Though it is unclear where the samples precisely fit in the eruptive sequence, they likely represent juvenile material that was erupted in the early phases of the eruption.

The 5680 Cal. BCE (GDS 15) eruption sampled from a distal pyroclastic sequence in 2011. The unit is composed of moderately to strongly vesicular, dark gray to black scoria and was sampled across the layer at constant depth (Figure 3, sampling location D; Table 2). The sampled scoria is rich in plagioclase and pyroxene.

Scoria and pumice samples for the 1657 Cal. CE, 1010 Cal. CE and 5680 Cal. BCE eruptions were treated with the same methods. Following washing and drying, samples were hand sieved in fractions from >32 mm down to <63 μm. Clasts from larger fractions (>32 mm and 16–32 mm fractions) were taken for crushing. Handpicking clasts ensured any alteration was removed prior to crushing. The clasts, with alteration removed, were sectioned and crushed to powders (100 μm) for inductively coupled plasma mass spectrometry (ICP-MS). The ∼341 Cal. CE eruption sample is composed of individual large scoria clasts which were treated individually. Any alteration was removed with only fresh sections of clasts taken for ICP-MS and crystal studies.

This study focuses on orthopyroxene, which was present in all samples and has well constrained diffusion coefficients (Dohmen et al., 2016). Crystals were randomly hand-picked from the 500 μm sieved size fraction, using a binocular microscope, with only highly altered or broken crystals discarded. 30–40 whole, euhedral crystals per eruption were mounted in epoxy grain mounts with the long (c) axes parallel to the grain mount surface in order to examine diffusion along the b-axis. In some minerals, diffusivities vary along different crystal axes; however, in orthopyroxenes there is no variation in diffusion along the c and b crystal axes meaning crystal orientation can be done by eye (e.g., Kilgour et al., 2014; Solaro et al., 2020). Fe and Mg fractionate between the M1 and M2 non-equivalent octahedral cation sites in orthopyroxene (Ganguly and Tazzoli, 1994), where the distances between the M1 and M2 sites along the b and c-crystallographic axes are relatively close, in comparison to the a-axis (Ganguly and Tazzoli, 1994; Schwant et al., 1998; Cherniak and Dimanov, 2010). Dohmen et al., 2016 report diffusion along the a-axis is ∼3.5 times faster making crystal orientation important. In order to discount any effects of crystal sectioning, only crystals with zones of equal size were imaged.

Each grain mount was polished using silicon carbide polishing pads (P-800, P-1200, P-2400 and P-4000) and finished using 0.3 μm grade aluminum powder. Finally, grain mounts were carbon coated in preparation for scanning electron microscopy (SEM) and Electron Micro-Probe analysis (EMPA).

Data Acquisition and Analysis Methods

Powders (100 µm) from crushed clasts were analyzed by ICPMS for major elements at Institut de Physique du Globe de Paris (IPGP) using the ICP-QMS 7900 Agilent.

Grain mounts were examined with back-scattered electrons (BSE) to image local chemical variations using the SEM EVO MA10 Zeiss at IPGP, following the method of Kilgour et al., 2014. Conditions of 15 kV accelerating voltage, 2 nA probe current and a working distance between 10 and 13.5 mm were used. BSE gray scale images are used as pyroxene greyscale variation is strongly correlated to Fe-Mg content (e.g., Saunders et al., 2012; Allan et al., 2013; Figure 4). Crystal zoning was categorized and ∼30 pyroxene crystals for each eruption, with clear congruent zoning, were selected for diffusion modeling. These crystals were imaged using high contrast and high magnification with the final images showing 50–100 μm of the crystal (Figures 4A,B).

FIGURE 4. (A) Overview of a crystal from the ∼341 Cal. CE eruption, the crystal is euhedral with melt inclusions and large oxide inclusions. (B) Close up of the zoning with the area chosen for diffusion studies shown by the white box. This zone is congruent and on closer inspection this dark zone has a thin, pale zone within it. This may be a result of melt composition variation which results in small changes in the zoning patterns. This is also reflected in the ImageJ gray scale profile. (C) Profile produced by ImageJ, a steep decline in greyscale at the boundary of the sharp zone is seen, this then increases again in the separation between the dark zones before decreasing again. The sudden rise in gray value at the end of the profile is an effect of the SEM creating a bright, white rim at the edge of the crystal. (D) EMPA transect across the crystal rim, confirming the crystal composition and the correlation between Fe-Mg and greyscale value. A similar shape is observed in both the ImageJ and EMPA profiles.

We used 10 crystals per eruption for EMPA chemical profiles (Figure 4D). Transects were completed using a SX-100 at University Pierre and Marie Curie (UPMC) calibrated with natural mineral samples with beam conditions of 15 kV accelerating voltage, 10 nA sample current and a focused beam of 1 μm. As well as confirming the majority of crystals are orthopyroxene, transects were taken along the crystal rims. This method provides quantitative analysis along the profile, allowing calibration of greyscale intensity to zoning composition (Allan et al., 2013).

Fe-Mg diffusion zones have a strong negative linear relationship with the BSE greyscale value, meaning greyscale values are proportional to Fe-Mg composition. This means variations in Fe-Mg composition records diffusion processes. (Figures 4C,D; Allan et al., 2013). For each crystal, three to six profiles on the same crystal zone were taken, and the timescales averaged. Where possible, images were taken of crystal zones on both sides of the same crystal face and transects completed. The method of Morgan et al. (2004) was followed when converting BSE greyscale images and EMPA compositional profiles. The main assumption is that the initial profile between two zones follows a step function, which is modified by diffusion to form sigmoidal concentration gradients. An initially steep profile becomes broader with time, in relation to the element inter-diffusion rate, which is stopped by eruptive quenching. This means maximum timescales are calculated (Morgan et al., 2004; Barker et al., 2016).

To obtain greyscale intensity values and profile quantification, we used the free software package ImageJ (Figure 4C), allowing greyscale profile quantification. Images were rotated where needed so the crystal boundary was perpendicular to the diffusion direction. This also prevents profiles from becoming anomalously stretched (Allan et al., 2013; Kilgour et al., 2014). To reduce the noise, the greyscale profiles were averaged over the columns of a portion of the image. The averaging areas were relatively consistent between all samples, but were adapted depending on the image and zone observed.

Diffusion Model

When converting the pyroxene crystal zoning greyscale profile into a timescale, we assume crystal growth occurred rapidly. Therefore, the diffusion profile is only related to chemical relaxation, as confirmed by the sigmoidal profiles observed. If the curvature of the mineral boundary is small and the initial mineral compositions and its surroundings are homogeneous, the diffusion process is approximately one dimensional and is expressed as per Costa et al. (2020),

where C is the concentration of a particular chemical species, x is the perpendicular distance away from the boundary, t is the time since the diffusion process started (i.e., the time elapsed since the crystal came into contact with fresh magma) and D is the diffusion coefficient that is modeled by (e.g. Dohmen et al., 2016),

which depends strongly on temperature, T, activation energy, E_A, and the ideal gas constant, R. Diffusion of the different chemical species in the mineral proceeds on timescales that depend on each species diffusivity. An initially sharp “step-like” boundary may exist between mineral and melt. However, as diffusion progresses, molecules of given species, initially in the mineral, move into the melt which smooths the mineral-melt boundary.

Diffusion-like processes such as Eq. 1 have an intermediate asymptotic for which an analytical solution is well known (e.g., Barenblatt and Vázquez, 2003; Costa et al., 2020) so that, after sufficient time, the concentration profile takes the form,

where C_melt and C_mineral refer to the bulk concentration in the melt and mineral, respectively, and µ is the profile midpoint. In this work, we used greyscale intensity as a proxy for these concentrations, so that C_melt should be read as the greyscale intensity in the part of the image corresponding to the melt. When the magma containing crystals is erupted (after t = τ), the diffusion process halts, and we say that τ is the diffusion timescale. We denote the set p = {C_melt, C_mineral, μ, D, τ} as the model parameters.

Optimization of Model Parameters

Each greyscale profile was fitted to the diffusion profile (Eq. 3) in order to obtain model parameters. The model parameters are determined by either propagating the model forward in time until the solution (model) matches the profile (data), as achieved by Girona and Costa (2013) and their software in the analysis of olivine by DIPRA, or adjusting parameters until the model “fits” the data, typically done by eye (e.g., Kilgour et al., 2014; Fabbro et al., 2017). However, manual adjustment is prone to bias as there is no metric for judging goodness-of-fit, i.e., whether the current solution is globally optimal. Due to the strong non-linearity of Eq. 3 and the interdependence of model parameters, particularly D, small changes to the model may require substantial changes in these parameters. Such pitfalls are avoided by using well-tested, robust numerical optimization schemes such as non-linear least-squares regression (e.g., Press et al., 2007). The metric that determines goodness-of-fit is often termed the misfit function and is based on the residual between model solution, f (p) and data, d,

Eq. 4 is a measure of how closely the model fits the available data points, which is how close the model comes to the data points. Hence, our goal is to minimize its value by adjusting the model parameters using the optimization schemes mentioned above. The optimal solution is obtained when Eq. 4 is minimised.

The data were manually selected from the profiles extracted from SEM images and only these selected data were input into the fitting algorithm (Figure 4). The suitability and quality of this selection was checked as part of the fitting quality process (see below). We only selected profiles that had clear “plateau” to the left and right of the diffusion zone (Figure 5, Supplementary Data S1). Several (3–6) profiles were extracted per sample, which later allowed us to check the consistency of profiles from the same sample and check for quality.

FIGURE 5. SEM greyscale intensity as a function of traverse location extracted from the SEM image for a sample from the 1010 Cal. CE eruption (blue dots), the data selected for fitting by the model (red dots), and the model solution (black line). Also shown are the fitted parameters (timescale, diffusion coefficient, D, and effective temperature), the theoretical profile midpoint (μ–gray line), and the “goodness-of-fit” (R). The gray shaded area represents the mean ±1 SD of the selected data.

As an initial estimate for p, we chose C_melt and C_mineral as the minimum and maximum greyscale values from the selected data, µ as the midpoint of the data range, and τ = 2 × 10⁶ s which is considered as a typical timescale value, given a typical diffusion coefficient (Dohmen et al., 2016). For the diffusion coefficient, the most important factor to consider is temperature, which has a large effect on Fe-Mg diffusion and therefore timescales. Even a small temperature variation results in a large change in the diffusion timescale; in the case of La Soufrière de Guadeloupe, a temperature variation of ±20°C results in a broad range. For example, a timescale of 50 days for the 1657 Cal. CE eruption crystals 0913D-B-05 is derived at 975°C, while at 995°C yields 30 days and at 955°C yields 78 days (Cherniak and Dimanov, 2010; Allan et al., 2013; Fabbro et al., 2017). For each eruption, the temperature, T₀, was determined using an orthopyroxene-melt geothermometer from Putirka (2008). Orthopyroxene crystal rim compositions and adjacent groundmass glass compositions were analyzed using EMPA measurements and used as melt-crystal pairs in the Putirka (2008) geothermometer (Eq. 28a in Putirka, 2008). Using glass-rim pyroxene pairs which are in equilibrium (pyroxene K_D = 0.23–0.35) ensures the temperature calculated is representative of the area over which diffusion occurred (Putirka, 2008). The geothermometer used here yields a 25–30°C error, which is the smallest magnitude for model error expected. Multiple temperature calculations are averaged for each eruption to produce the temperatures used in the final diffusion calculations (Table 2).

Another factor to consider is the magma oxidation state, which has a moderate effect on diffusivities.

Negative ΔNi-Ni-O (NNO) (Oxygen fugacity, expressed in relation to the NNO buffer) will produce small diffusion coefficients and, therefore, longer calculated timescales (Allan et al., 2013; Fabbro et al., 2017). Oxygen fugacity (fO₂) was taken from Pichavant et al. (2018) who calculated NNO+0.8 for the La Soufrière de Guadeloupe magma when the resident magma and arrival recharge magma are in equilibrium. This value is based on geothermometry and experiments on basaltic and andesitic eruption products. The 1530 Cal. CE eruption groundmass glass and melt inclusion values have a wide range of compositions but are most similar to 1010 Cal. CE and 1657 Cal. CE eruptions (Figure 6C; Pichavant et al., 2018). The 5680 Cal. BCE and ∼341 Cal. CE eruptions are in the same whole rock and groundmass glass compositional trend and show some compositional overlap with the more recent eruptions. This indicates NNO + 0.8 is a valid fO₂ to use for all the eruptions, particularly as the fO₂ has a minimal effect on the timescale calculation in comparison to the effect of temperature.

FIGURE 6. (A) Total alkali-silica diagram (Streckeisen et al., 2002) for ICP-MS whole rock data and groundmass glass compositions from the range of eruptions studied. Whole rock data is basaltic to andesitic in composition and the groundmass glass is more evolved, ranging from andesitic to rhyolitic. The older eruptions, 5680 Cal. BCE. and ∼341 Cal. CE show less evolved compositions in comparison to the younger eruptions, 1010 Cal. CE. and 1657 Cal. CE. (B) Ternary diagram showing the proportions of pyroxene endmembers Wollastonite (Wo), Enstatite (En) and Forsterite (Fo). The orthopyroxenes studied are enstatite rich and show little compositional variation between eruptions. (C) Groundmass glass SiO2 variation diagrams for (i) Al₂O₃, (ii) CaO, (iii) Na₂O, (iv) TiO₂, (v) MgO, (vi) FeO, (vii) K₂O and (viii) MnO for the eruptions studied. The 1530 Cal. CE eruption is also included to provide a more complete picture of La Soufrière de Guadeloupe (Boudon et al., 2008; Pichavant et al., 2018). All these glass compositions were determined to be in equilibrium with the nearest crystal, with the K_D (Mg/Fe)_melt/(Mg/Fe)_mineral between 0.23–0.35, the accepted equilibrium values for melt-pyroxene crystals (Putirka, 2008).

We took the value predicted by Eq. 2 using the melt-matrix temperature for that sample, T₀, and fO₂. While we placed no constraints on the values of C_melt, C_mineral, μ, and τ, we imposed upper and lower bounds on D according to T₀ ± 30 K, as Putirka (2008). We found the calculated temperatures (see Table 2) were well within these bounds, which provides confidence in the imposed conditions.

We have developed a procedure for automatically and robustly fitting the model parameters using well-known and robust numerical routines contained within the Scipy module in Python. The non-linearity of Eq. 3 requires a good estimate of the solution to be entered as a starting condition for the iterative procedure, which proceeds until either convergence is reached, or the maximum number of iterations is achieved. In the latter case, no valid model solution would be output so we either tried to improve the data selection (e.g., by extending the range), or rejected this profile from further analysis. In most cases, we added further profiles with longer plateaus on either side of the diffusion zone, from the same sample and iterated this process until the timescales converged.

Data Quality and Model “Goodness-of-Fit” Assessment

To check whether the data were of sufficient quality and quantity to be fitted using the model, we devised a pair of quality assessment criteria. First, we check how closely the model fits the data using the coefficient of determination,

where SS_res is the sum of squared residuals and SS_tot is the total sum of squares. If the model closely reproduces the data points, then R² is close to 1. We arbitrarily chose a threshold of 0.85 and exclude profiles with R² less than this value.

Secondly, to assess whether the selected data are representative of the general model shape, that is their range and shape were sufficient to capture the “sigmoidal” curve of the model (Figures 4, 5), we ensured the calculated midpoint (μ) lay within one standard deviation of the mean of the data points. When µ is too extreme compared to this test, the model solution poorly captures the shoulders and plateau of the profile, which leads to over–or under-estimated timescales. If a selected data profile failed either of these criteria, it was flagged to either be replaced by changing the range of data points from that profile, or to be replaced by a new profile.

Finally, we looked at the agreement between the set of timescales calculated for each sample. We calculated the mean and standard deviation of the timescales from the 3–6 profiles per sample. Profiles with extreme timescales were eliminated and new profiles added as necessary, until the set of all timescales converged on a common value.

Population Timescales

Given the random nature of the crystal sampling process, the population of timescales for each eruption is representative of that eruption. This is an important step to allow for comparison between eruption timescale populations. We may look at how these populations are distributed to determine the expected timescale for that eruption. The time between magma upwelling or recharge events is a random variable, whose cumulative frequency distribution is described by a gamma distribution (Cowan, 1998),

where y are the data points, alpha and beta are the shape and rate parameters, respectively, and γ(z) is the gamma function. The expected, i.e. typical, value of the distribution is given by Cowan (1998). We now show how such a distribution is fitted to timescale data to reveal, in addition to the expectation, many statistical features of the data including standard error estimates of the parameters.

Maximum likelihood estimators (Cowan, 1998) are constructed for the parameters in Eq. 6 by writing a likelihood function as,

Maximizing the likelihood function (or alternatively minimizing the negative log-likelihood function) leads to the following pair of equations,

where Ψ⁽ⁿ⁾ is the nth-order polygamma function. These equations were solved to find α and ß. We note the latter equation predicts that E [y] = $\overline{y}$ = 1/n ∑ y_i, the arithmetic mean of the data. For the parameter variance, we note the covariance matrix of the model parameters is given by,

where θ = [α, β] is the vector of model parameters and E [H (θ)] is the expected value of the Hessian (matrix of 2nd derivatives) of θ. This expression is termed Ficher’s information matrix (Cowan, 1998). The square-roots of the diagonal terms of Eq. 9 are the standard errors on [α, β]. For a sufficiently large sample size (typically several tens of samples), the standard errors become normally distributed under the central-limit theorem and we may write, for the population timescale standard error,

A final step is to check how well the distribution fits the population of timescales, for which we used the Kolmogorov-Smirnov goodness-of-fit test (Dodge, 2008), as formulated in Python’s Scipy module. This test returned a probability (p-value) allowing us to evaluate whether or not to accept the null-hypothesis (that the Gamma distribution does not represent the data) at a given confidence level, which we took to be 0.05 (5%, cf. Gibbons and Pratt, 1975). That is, p-values greater than 0.05 passed this goodness-of-fit test. We found the distributions to all eruptions had p-values well above 0.05.

Results

Geochemistry

Volcanic rocks from Basse-Terre show a large compositional range, from basalt to dacite, with the least evolved compositions observed in the Mont Caraïbes suite in Southern Basse-Terre (Bissainte 1995; Boudon et al., 2008; Samper et al., 2009). La Soufrière de Guadeloupe shows a narrower magma compositional range, from basalt to andesite (Boudon et al., 2008; Samper et al., 2009).

Though juvenile clasts from the 1530 Cal. CE eruption show characteristic light and dark bands (Boudon et al., 2008; Pichavant et al., 2018), the four eruptions studied here show no evidence for compositional banding. The clasts from each eruption used for ICP-MS analysis appear homogeneous in composition. The younger eruptions, 1657 Cal. CE and 1010 Cal. CE eruptions are composed of light gray clasts; in comparison, the older eruptions have much darker clasts. The ∼341 Cal. CE eruption clasts used are black juvenile scoria and the 5680 Cal. BCE eruption clasts are dark gray.

Whole rock analysis shows geochemical variation between basalt and andesite for the eruption suite studied (Figure 6A). Through time, the whole rock silica content becomes more silicic, and does not correlate with eruption style. The ∼341 Cal. CE Strombolian eruption is the most mafic composition with a basalt and basaltic-andesite composition observed (Figure 6A). In comparison, younger eruptions are more evolved with 1530 Cal. CE and 1657 Cal CE eruptions composed of basaltic-andesite and andesite (Figure 6A; Boudon et al., 2008). The 1010 Cal. CE Plinian eruption has the most evolved composition and is only composed of andesite. However, generally there is little variation in the more recent eruption compositions.

Across all the eruptions the main mineral phases are pyroxene and plagioclase with minor oxide phases. No amphibole or olivine are observed. Orthopyroxene is more abundant than clinopyroxene and are present in a ∼4:1 ratio. Both pyroxene and plagioclase contain melt inclusions and oxide inclusions; orthopyroxenes are described in detail below. Plagioclase often has large fluid inclusions; particularly in older samples. Oscillatory zoning is observed on the plagioclase rims and patchy zoning is observed in the plagioclase cores (Vance, 1965; Fabbro et al., 2017).

The orthopyroxene enstatite compositions across all eruptions range from En₄₈ to En₆₈. Overall, this a relatively narrow range in compositions, with the 5680 Cal. BCE, 1010 Cal. CE and 1657 Cal. CE eruptions showing the same orthopyroxene compositions (Figure 6B). The 5680 Cal. BCE eruption shows the largest orthopyroxene compositional range with enstatite compositions covering the full range i.e., En₄₈ to En₆₈. The ∼341 Cal. CE eruption has the highest enstatite compositions, overlapping only with the 5680 Cal. BCE eruption (Figure 6B). The 1010 Cal. CE eruption has a small enstatite compositional range and shows a similar composition to the 1657 Cal. CE eruption. These two eruptions show a large compositional overlap with the 5680 Cal. BCE eruption, but have a lower enstatite content than the ∼341 Cal. CE eruption (Figure 6B). The enstatite contents reported here are lower than those used by Dohmen et al., 2016 for calibration of the Fe-Mg interdiffusion coefficient which was calibrated for En₉₁. There is some speculation that there is weak dependence of the diffusion on Fe content (Dohmen et al., 2016), which could offset any timescale results calculated for enstatite contents lower than En₉₁; however, based on our current knowledge this diffusion coefficient is acceptable to use for lower En contents (e.g., Costa et al., 2013; Solaro et al., 2020).

The groundmass glass compositions cover a wider range from 57–77 wt% SiO₂ (andesitic to rhyolitic compositions; Figures 6B,C). The groundmass glass compositions of individual eruptions cannot be easily distinguished from major element binary plots (Figure 6C). The 5680 Cal. BCE shows the least evolved compositions but also the largest range in composition from basalt to rhyolite (57–74 wt% SiO₂). This eruption also has the highest Al₂O₃ content (15–24 wt%) and highest CaO content (3 – 9 wt%). The ∼341 Cal. CE eruption also has a large range in composition similar to the 5680 Cal BCE eruption, ranging from andesite to rhyolite in composition (60–74 wt% SiO₂).

Similar to the whole rock analysis, the younger eruptions show more evolved and narrower glass composition range and have large compositional overlap. The 1010 Cal. CE eruption has the narrowest compositional range and is dacitic in composition (67–70 wt% SiO₂). All major elements for this eruption have a narrow compositional range, for example Al₂O₃ varies over only 1 wt% (Figure 6B). The 1530 Cal. CE eruption is characterized by the most evolved groundmass glass, with compositions ranging from dacitic to rhyolitic (65–77 wt% SiO₂; Boudon et al., 2008; Pichavant et al., 2018). Despite the mingling textures observed in 1530 Cal. CE eruption products, the groundmass glass appears relatively homogenous covering similar range of compositions to the 1657 Cal. CE eruption. The 1657 Cal. CE eruption groundmass glass is also dacitic to rhyolitic in composition (68–76 wt% SiO₂) and shows a narrow range of major element. compositions.

The glass compositions were used with crystal rim compositions to calculate system temperatures for each eruption, with an error of ±30°C on each temperature (Putirka, 2008; Table 2). The highest temperature calculated was 1,025°C for the 1010 Cal. CE eruption, followed by 1010°C for the 5680 Cal. BCE and ∼341 Cal. CE eruptions, the lowest temperature of 975°C was calculated for the youngest eruption of 1657 Cal. CE eruption. For the 1530 Cal. CE eruption, temperatures are reported between 880 and 970°C; however, an equilibrium temperature of 900°C was used for diffusion calculations reported in Pichavant et al. (2018) and Bourgeoisat (2018).

Crystal Zoning

All crystals used for this study were taken from the 500 μm sieved size fraction, and large compositional variations are not observed. This means the crystals must be grouped based on differences in textures, including zoning and other disequilibrium textures. We identify six different crystal populations, which based on the observed textures, have experienced different histories (Figure 7).

FIGURE 7. Schematic representation of the five different crystal types with zoning, along with high magnification schematic of the crystal rims showing representative zoning for the crystal type. Enstatite compositions are shown for crystal rims and cores. (A) Type-1 crystals are simple and have a single sharp zone. (B) Type-2 crystals are also simple with single zone but the zones are diffuse and a sharp boundary is not observed. (C) Type-3 crystals have multiple zones, there are generally two or three zones which are sharp and often vary in composition. (D) Type-4 crystals are more complex and have been partially resorbed, the crystals still have a distinct core with an undulating resorption surface and a thick rim. The crystal core shows patchy zoning in some samples, with the core size depending on the amount of resorption the crystal has under gone. The resorption rim has multiple zones, which are significantly finer than multiple zones in the simply zoned crystals. Zoning is more complex and records a longer history than the simple crystals. (E) Type-5 crystals are fully resorbed crystals with no remaining core. They are characterised by patchy zoning in the center of the crystal. The crystal rims have similar styles of zoning as Type-4, with multiple zones of variable thickness. (F) Type-6 crystals are unzoned.

There are three types of crystals with similar features and histories, crystals Types-1-3 are simple crystals with few features and often host melt inclusions and oxides. These three crystal types are distinguished based on differences in the patterns of chemical zoning as they show little compositional variation (Figure 7). Type-1 crystals are simple crystals characterized by one simple, sharp, 2–15 μm wide zone (Figure 7A). In the majority of crystals, this zoning is reverse (i.e., Mg-rich rim and a relatively Fe-rich core). Type-2 crystals show one diffuse zone that is 10–80 μm wide, wider than the Type-1 crystals sharp single zones (Figure 7B). The majority of crystals exhibit reverse zoning. Type-3 crystals have multiple, sharp 2–15 μm wide zones, with the outer zone in the majority of crystals being a reverse zone (Figure 7C).

The other two types identified are more complex crystals, with various resorption textures and similar features and histories. The resorbed areas of these complex crystals have higher enstatite contents than the simple crystals, though across the whole dataset overlap in enstatite contents are observed (Figure 7). Type-4 crystals are partially resorbed crystals with thick, magnesium-rich rims and iron-rich cores with undulating core-rim boundaries and patchy zoning present in the rim (Figure 7D). The remaining cores show similar enstatite contents to the simple crystals, while the rims have higher enstatite contents (Figure 7). Rims range from ∼100 to ∼215 μm in thickness, with some crystals showing only a small core remaining. The crystal edges show complex zoning, with fine multiple zones on the scale of 1–6 μm, often with up to six zones present; reverse zoning is often present on the outer edge of the crystal. Type-5 crystals are fully resorbed with no core remaining and have high enstatite contents, comparable to the rims of Type-4 (Figure 7E). These crystals have patchy zoning close to the crystal core and multiple zones at the crystal rim, similar to the zoning observed in the partially resorbed crystals; reverse zoning is often present on the outer edge of the crystal. Type-4 and Type-5 crystals can be considered endmembers, experiencing various degrees of resorption.

As well as a large proportion of zoned crystals, significant unzoned crystal populations are observed in each eruption, with 30–40% of the crystals mounted for analysis unzoned. These crystals are Type-6 crystals, which show no visible zoning, few features and no compositional variation, with enstatite contents in the same range as the simple crystals and Type-4 cores (Figure 7F).

For the 5680 Cal. BCE eruption, all six crystal types are observed (Figure 8). The zones used for calculation are interpreted as reversely zoned with a light gray core and dark gray rim in backscattered electron (BSE) images, correlating to a relatively Fe-rich core and Mg-rich rim. This indicates a change from lower temperatures during core growth to higher temperatures during rim growth. Looking at the breakdown of crystal proportions in this eruption shows Type-1 crystals are the most common, simply zoned crystals with 11% of crystals showing this zoning type. Type-2 and Type-3 crystals are both uncommon with only 4% of crystals observed with these zoning types. This eruption shows a large proportion of crystals with a variable degree of resorption, with Type-4 and Type-5 crystals making up 48% of the crystals. There is also a substantial population of unzoned crystals, 38% of crystals are of Type-6.

FIGURE 8. Histogram showing the number of each crystal type per eruption. Each eruption is made up of a different proportion of each crystal type, with not all crystal types present in the different eruptions. The 1657 Cal. CE eruption is composed of entirely simple crystals showing no resorption, with unzoned crystals making up the largest crystal population. The 1010 Cal. CE eruption is also entirely composed of simple crystals and is predominantly composed of crystals with diffuse zoning. This crystal type is the most abundant in the 1010 Cal. CE eruption, with only one or two crystals observed in other eruptions. The ∼341 Cal. CE eruption is also mostly composed of simple crystals with small populations of fully resorbed and partially resorbed crystals. The 5680 Cal. BCE eruption has the largest populations of resorbed crystals and only small populations of simple crystals. This highlights that the complex, resorbed crystals are only found in the older eruptions.

The ∼341 Cal. CE eruption has crystal populations most similar to the 5680 Cal. BCE eruption (Figure 8). Type-1 crystals make up 39% of the crystals, 75% are reversely zoned with light gray cores and dark gray rims. Type-2 crystals are not present in this eruption and Type-3 crystals represent only 5%. In comparison to the 5680 Cal. BCE eruption, for the ∼341 CE eruption Type-4 represent only 4% and Type-5 represents 17% of the crystals. Again, there is a significant population of unzoned Type-6 crystals, making up 36% of the crystals.

The 1010 Cal. CE eruption is composed only of simple crystals, no partially or fully resorbed crystals are present (Figure 8). Type-1 crystals make up only 8% of the crystals with Type-2 crystals making up 46% of the crystals, this crystal type is the most abundant in the 1010 Cal. CE eruption in comparison to the other eruptions. The remainder of crystals are Type-3 crystals (14%) and of Type-6 crystals (32%) with no zoning. All zoned crystals show reverse zoning with a light gray, Fe-rich core and a dark gray, Mg-rich rim.

The 1657 Cal. CE eruption has the largest population of Type-6 crystals, with 65% showing no visible zoning. There are 32% of Type-1 crystals and 3% of Type-3 crystals; all the zoned crystals show reverse zoning (Figure 8). This eruption shows no complex crystals with resorption features.

Timescale Modeling Results

Diffusion timescales across the core-rim boundaries for the eruptions studied give a timescale range from <1 to 3,052 days (Table 2). Our method has investigated how the timescale populations are distributed, highlighting clear differences between eruptions, including variations in the maximum values, range of values and expected timescales.

The 5680 Cal. BCE Plinian eruption records short timescales with a range from <1 to 169 days and the expected (mean) timescale calculated as 22.8 ± 0.43 days (Figure 9). The ∼341 Cal. CE Strombolian eruption records a large range of timescales from <1 to 3,052 days. The modelled data yield an expected timescale value of 59.4 ± 0.40 days. This expected value is not comparable to any other eruption studied, though the range of values is comparable.

FIGURE 9. Distribution of timescales for the 1657 Cal. CE, 1010 Cal. CE, ∼341 Cal CE and 5680 Cal. BCE eruptions. Blue bars indicate the empirical cumulative density frequencies (CDF), red curves indicate the theoretical distributions fitted to these and the green curves represent the theoretical probability density frequencies. On each plot we indicate the mean and median timescales predicted by the theoretical distributions and the p-values (goodness-of-fit) between model and data.

The 1010 Cal. CE Plinian eruption records the longest timescales from <1 to 2,462 days. These data have an expected timescale value of 361.0 ± 0.40 days, considerably higher than calculated for any other eruption (Figure 9). This eruption is considerably different to the 5680 Cal. BCE Plinian eruption, despite the expectation that these eruptions are comparable. This is a direct consequence of the characteristic crystal populations in the 1010 Cal. CE eruption which shows a large proportion of Type-2 crystals with diffuse zoning, not present in the 5680 Cal. BCE eruption (Figure 8).

The 1657 Cal. CE Vulcanian eruption records the shortest range of timescales from <1 to of 84 days. The expected timescale was calculated as 18.8 ± 0.37 days, the shortest expected timescale value calculated for the eruptions studied (Figure 9). This expected value is comparable to the 5680 Cal. BCE expected value of 22.8 days, with most crystals for both these eruptions recording timescales in the range of 3 weeks.

The two Plinian eruption timescales are not comparable, with the 5680 Cal. BCE eruption, showing short timescales and the 1010 Cal CE eruption, showing relatively long timescales. This also highlights the different eruption styles are not distinguishable by their timescales, with the 1657 Cal. CE vulcanian and 5680 Cal. BCE Plinian eruptions showing the most comparable timescales, despite being different eruption styles.

Discussion

How Do We Interpret the Diffusion Timescales?

The zoning observed in the orthopyroxene crystals is a gradient in composition and results from magma mixing, mingling and crystallization causing melt differentiation (Davidson et al., 2007; Ginibre et al., 2007; Streck, 2008; Ubide and Kamber, 2018; Costa et al., 2020). These processes create chemical gradients in the crystals and/or melts, resulting in element diffusion as the system re-equilibrates during high temperature storage. Crystal zoning is generally the result of both crystal growth and diffusion, but the natural system is often simplified and zoning is attributed solely to diffusion, particularly when sigmoidal profiles are observed (Costa et al., 2020; Solaro et al., 2020). Experimental data on the 1530 Cal. CE deposits suggest magmas reside within a temperature range of 825–875°C (Pichavant et al., 2018; Figure 2). The temperatures calculated from the groundmass glasses in this study (975–1,024°C) indicate the intruding magma prior to each eruption was hotter than the existing magma reservoir. Hotter magma intruding into a cooler system will produce the reverse zoning observed in the orthopyroxene populations.

When interpreting the diffusion timescale data, the entire range of timescales calculated provides information on how an intrusion has moved through and interacted with the system. Numerical simulations of a crystal free basaltic liquid intruding through the base of a crystal rich pile shows magma interactions occur in a ‘mixing bowl’ (Bergantz et al., 2015; Schleicher et al., 2016; Bergantz et al., 2017; Schleicher and Bergantz 2017; Carrara et al., 2019). The mixing bowl forms through failure along frictional contacts of settled crystals, moving up through the system to form a wedge-shaped area above the intrusion location (Bergantz et al., 2015). The crystals experience a wide range of chemical environments due to the shape of the mixing bowl which results in varying thermo-chemical effects on the crystals. The mixing bowl results in crystals encountering new melt at different times, with crystals only recording an event when their environment is perturbed (Bergantz et al., 2015; Schleicher and Bergantz, 2017; Cheng et al., 2020). The time between melt injection into the system and crystals interacting with melt in the mixing bowl is defined as the delay time, and can explain the spread of data observed in timescale results (Cheng et al., 2020). The delay time assumes an intrusion is occurring from the system base and depends on many unknown factors such as crystal shape and size, viscosity, intrusion location, duration and both existing system and intrusion geometry (Cheng et al., 2020).

Crystals at the base of the wedge-shaped mixing bowl, and therefore close to a new intrusion, experience a perturbation of their environment before crystals that are located above, so the former record the shortest delay time. As these crystals spend the longest time interacting with the new magma, they record the longest timescales derived for a given eruption. There is also likely a delay between the onset of the intrusion and time when crystals start recording a change. However, the maximum timescale calculated for each eruption is still considered a minimum approximation for the onset of the intrusion. Crystals recording short timescales are located farthest from the initial intrusion location. They interact with the intrusion for the shortest period and so have a long delay time. By the time the expected timescale is reached, most crystals in the mixing bowl are interpreted to have interacted with the intrusion. At this point, as most crystals have interacted with the intrusion, this could imply that the system has a high enough proportion of eruptible melt to allow destabilization and eruption. Crystals in the mixing bowl which do not interact with the magma are quenched and unzoned on eruption, with the population size depending on the intrusion’s residence time (Cheng et al., 2020). Alternatively, syn-eruptive crystallization may also explain the populations of unzoned crystals (Type-6).

The mixing bowl model by Bergantz et al. 2015 is developed for the injection of basalt into an olivine-rich mush. For the more evolved system of La Soufrière de Guadeloupe, characterized by greater viscosity contrasts, this may result in different mixing dynamics. However, no mafic enclaves or mingling textures are observed for the explosive eruptions discussed here, suggesting, despite potential viscosity contrasts, that mixing has been relatively complete, and that the mixing bowl model is plausible.

These processes are reflected in the gamma distribution, which describes all events that share the same properties (i.e., are generated by the same process) and models them as random processes. In this case, the gamma distribution is derived from the individual likelihood of all intrusion events across the eruption's crystal populations. The timescale distribution reflects an intrusion producing a mixing bowl where the crystals come into contact with the intrusion at different times, with the expected timescale being the typical value of the distribution. The expected timescale/typical distribution value shows the average time it took a crystal for a specific eruption to come into contact with the intrusion. Therefore, the gamma distribution allows the probability of a crystal experiencing an event at a given time to be estimated. In the case of La Soufriére de Guadeloupe, crystals have a smaller probability of recording a short delay time/longer timescale (on the order of a year; e.g., 1010 Cal. CE), and a larger probability of recording a longer delay time/shorter timescale (on the order of days; e.g., 1657 Cal. CE). This could relate to several parameters including: the intrusions interaction with the system, the composition of the system (including volatile content) and system geometry (crystals recording a longer delay time/shorter diffusion timescale farthest from the intrusion are more likely to be erupted than those closer to the system base which interact with the intrusion first).

Though there is some debate as to whether Fe-Mg diffusion in orthopyroxene records magmatic processes on short timescales (on the order of days), and if such timescales are related to diffusion, many studies do report short timescales related to diffusion processes (Druitt et al., 2012; Singer et al., 2016; Dohmen et al., 2017; Costa et al., 2020). Fast diffusing elements, such as Fe-Mg, record shorter timescales and are likely to only record the final thermal or chemical perturbation as these elements often re-equilibrate between events (Costa et al., 2020). In this study the same mineral type and elements are used to calculate timescales which reduces uncertainties arising when comparing slow diffusing verses fast diffusing elements. Further investigating the importance of growth vs. diffusion by modeling anisotropic element diffusion in two dimensions, modeling multiple elements with varying diffusivities and investigating crystal-melt partitioning relationships, would allow us to truly understand these short timescales (Shea et al., 2015; de Maisonneuve et al., 2016; Costa et al., 2020).

Timescales and Magmatic Processes

Which Mush System Processes Are the Timescales Related To?

The repose time between eruptions provides important information about the state of the system prior to an eruption. Repose time is the time elapsed between eruptions, for example, the repose time for the 1657 Cal. CE eruption is 127 years, which is the time elapsed between the 1530 Cal. CE eruption and subsequent 1657 Cal. CE eruption. For most eruptions this is straightforward to calculate from Table 1; however, as there is some uncertainty on when activity began at La Soufrière de Guadeloupe, the repose time for the 5680 Cal. BCE eruption is difficult to calculate. The first volcanic deposit attributed to La Soufrière de Guadeloupe is the basal debris avalanche deposited in 7140 ± 18 Cal. BCE (∼9,160 years ago) resulting from partial collapse of the proto-Soufrière edifice (Komorowski et al., 2005; Boudon et al., 2007; Komorowski, 2008; Legendre, 2012). A lava flow originating at the La Soufrière de Guadeloupe edifice base, flowed SW from the multi-phase collapse structure, and is dated at 6,000 ± 2,000 years by K-Ar (Samper et al., 2009). However, the date eruptive activity resumed after this partial collapse, how many eruptions occurred, and the size of the volcano prior to the 5680 Cal. BCE eruption, are unknown. Looking at known events prior to the 5680 Cal. BCE eruption gives a minimum repose time estimate of 311 years and a mean repose time of ∼611 ± 414 years.

This shows for the 1657 Cal. CE, ∼341 Cal. CE and 5680 Cal. BCE eruptions, the diffusion timescales increase, as repose time increases (Figure 10A). The 1530 Cal. CE eruption is also resolved into this trend, as a minimum timescale estimate is reported as 20–30 days and a maximum of 90 days (Bourgeoisat, 2018; Pichavant et al., 2018). Short expected timescales of magmatic processes, particularly for the 5680 Cal. BCE and 1657 Cal. CE eruptions correlate to short repose times. This suggests these crystals interacted with a magma intrusion that was hotter than the existing mush, which moved rapidly through the mush system, perhaps due to a high melt proportion already in the system (Figures 11A, 1B). Other processes which would allow an intrusion to move rapidly through the system include partial degassing at depth or if the intrusion was large enough, and melt induced fracturing processes (Burgisser and Bergantz 2011; Huber et al., 2011; Costa et al., 2013; Moore et al., 2014; Bergantz et al., 2015; Cassidy et al., 2015; Davydova et al., 2018). An intrusion moving rapidly through a system could account for the large un-zoned crystal populations observed (Type-6), which have had limited interaction with an intrusion and were entrained shortly before eruption and quenching (Bergantz et al., 2015; Bergantz et al., 2017; Cheng et al., 2020). Alternatively, Type-6 crystals may result from syn-eruptive crystallisation.

FIGURE 10. Expected timescale variation diagrams for (A) Repose Time, (B) Volcanic Explosive Index (VEI). Repose time for the 5680 Cal. BCE eruption is based a collapse event at 7160 Cal. BCE, the closest datable event. *1530 Cal. CE eruption data from Boudon et al. (2008); Bourgeoisat (2018); Pichavant et al. (2018).

FIGURE 11. Schematic sketches of the mush lenses beneath La Soufriére de Guadeloupe and Echelle with representative crystals shown at approximate locations within the mush zone. (A) The La Soufriére de Guadeloupe mush systemprior to eruption, this study suggests intrusions begin a between a minimum of 84 and 3,052 days prior to eruption. (B) In the case of 1657 Cal. CE and 5680 Cal. BCE eruptions, the expected value shows that themajority of the system could be affected by an intrusion 19 and 23 days prior to eruption, respectively. The intrusion moves through the system rapidly, in the case of the 5680 Cal. BCE eruption crystals are erupted from deep in the mush system which re-equilibrate during smaller intrusion events which do not destabilise the mush system. The 1657 Cal. CE eruption only samples simple crystals and unzoned crystals. (C) The ∼341 CE eruption expected value shows that the majority of the systemcould be affected by an intrusion 59 days before eruption. This intrusion has also moved relatively rapidly through the system and variation in intrusion location, results in destabilisation and eruption from the Echelle vent. Both complex and simple crystal populations are erupted. (D) Processes resulting in the 1010 Cal. CE eruption occurred slower than the other eruptions with the expected value showing the majority of the system could be affected by an intrusion 361 days prior to eruption. The very long delay time suggests that magma may have stalled in the mid-upper mush system, crystals with the diffuse zoning have spent the longest time in contact with the intrusion. (E) System destabilisation results in eruption of populations of simply zoned and unzoned crystals, perhaps from the upper lenses.

The 1010 Cal. CE eruption is not resolved into this trend, with the last major, magmatic eruption recorded as a dome eruption in 450 Cal. CE giving a similar repose time to the 5680 Cal. BCE eruption (Figure 10A). The long expected timescale may be assumed to correlate with a longer repose time, with system reactivation and production of eruptible melt taking longer following a recharge event. The lack of major explosive activity prior to the 1010 Cal. CE eruption suggests that the region of the magmatic mush system which produced the eruption was cold and inactive, and therefore hard to remobilize due to the presence of crystalline mush lenses with little eruptible melt.

1010 Cal. CE is the only eruption with a large crystal population with diffuse zoning (Type-2) which results in the long-expected timescales and large range of diffusion timescales (Figure 12). This indicates the system remobilized slowly and the crystals have spent a relatively long period interacting and equilibrating with new magma, recorded in Type-2 crystals. Magma intruding into a colder system with a higher proportion of crystalline material and a lower proportion of melt may promote slower remobilization (Figures 11A, Figures 1D,E). In the case of the 1010 Cal. CE eruption, the repose time shows only one dome eruption occurred prior to the 1010 Cal. CE eruption and the crystal populations show no evidence for multiple intrusion events. The 1010 Cal. CE eruption geochemistry also shows the smallest compositional variation, suggesting prior to this eruption the longer timescales allowed more time for homogenization (Figure 6; Figure 11A, Figures 1D,E). A process allowing an intrusion to rise slowly is magma degassing and stalling in shallow areas of the mush system (Ebmeier et al., 2013; Lu and Dzurisin, 2014; Cassidy et al., 2015; Rasmussen et al., 2018). Another control on this diffusion timescale may include the time taken for the temperature of a deep cold-storage region to reach the solidus and following this, the time taken for melt to move and remobilize the shallower reservoir, as shown by reactive transport modeling (Jackson et al., 2018).

FIGURE 12. Maximum, minimum and average timescale values for the different crystal populations recorded in each eruption. For single points, only one crystal was identified. The largest range of values are observed in the simple crystals Type-1 and Type-2; however, this may be a reflection of population size. Type-1 is found in all eruptions and records the shortest timescales. Type-2, is mainly found in the 10101 Cal. CE eruption and records the longest timescales and a large range of timescales. Type-3 has the fewest number of crystals and shows the smallest range of values.While the 1657 Cal. CE, ∼341 Cal. CE and 5680 Cal. BCE eruptions show similar averages, the 1010 Cal. CE eruption shows a higher average and minimal overlaps with the other eruptions.

The ∼341 Cal. CE and 1010 Cal. BCE eruptions have similar maximum diffusion timescales but different expected timescales (Table 2). This indicates while intrusions may begin at similar times prior to an eruption, the diffusion timescale distribution and resulting expected timescales, suggest the intrusions interact differently with the system. A transition in conduit location to Echelle or the eruption of a monogenetic cone may result in different expected timescales of magmatic processes preceding eruptions. In a shallow storage region, crystals interact with melt shortly before eruption and could explain the crystals recording short timescales observed, particularly in Echelle and in other monogenetic systems (Ruprecht et al., 2007; Johnson et al., 2008; Denis et al., 2013; Brenna et al., 2018; Figure 11A,C).

The 1530 Cal. CE eruptions distinct, characteristic compositional banding observed in hand specimens, indicates late-stage magma-mixing processes have occurred in the La Soufrière de Guadeloupe system (Touboul et al., 2007; Boudon et al., 2008; Pichavant et al., 2018). Hand-specimens for the eruptions used in the present study show no evidence for magma mixing, though field evidence for complex mingling textures are observed in the ∼341 Cal. CE lava flow, with mixing of a small-volume of evolved magma with a dominant basaltic magma (Boudon et al., 1988). This suggests mixing and hybridization occurred in the 1530 Cal. CE and ∼341 Cal. CE eruptions, and these processes were not involved in all eruptions. In particular, this brings into question the link between the temporally close 1530 Cal. CE and 1657 Cal. CE eruptions, suggesting these eruptions are most likely unrelated magmatic events. Alternatively, following the 1530 Cal. CE eruption, magma may have stored and fully homogenized at a shallow depth and erupted in the 1657 Cal. CE eruption, explaining the differences in textures observed.

Are the Timescales Related To Eruptive Styles?

The lack of precise eruption volume data available for La Soufrière de Guadeloupe eruptions, as a result of poor preservation of the deposits in the field, means the magnitude or intensity as defined by Pyle (2000) cannot be calculated (except for the better exposed and recent 1530 Cal. CE eruption which has an estimated minimum magnitude of 3.2, Komorowski, 2008; Komorowski et al., 2012; Komorowski et al., 2013). Instead, a Volcanic Explosivity Index (VEI) estimation, calculated from an approximate eruption volume is used (Newhall and Self, 1982; Legendre, 2012; Figure 10A). VEI gives a semi-quantitative scale but provides the best approximation of the eruption size for these events. Excluding the 1010 Cal. CE outlier eruption, VEI shows a weak correlation of decreasing expected timescale with increasing VEI (Figure 10B).

An expected working hypothesis, assuming the magmatic system was in a similar state, is that the two Plinian eruptions (5680 Cal. BCE and 1010 Cal. CE) result from similar processes, and therefore should record similar timescales of magmatic processes. However, the 1010 Cal. CE eruption expected timescale is considerably longer than the other eruptions, and along with the tenuous link between timescales of magmatic process and eruption style, suggests the diffusion timescales do not relate to eruption style. Instead, they provide an indication of the state of the magma system prior to eruption, which is recorded by the crystal populations. Indeed, variation in crystal abundance and compositional zoning is a direct result of mush system heterogeneity that reflects variations in the intrusion size, location, composition, and rheology, as well as processes such as the mixing and migrating of melt-rich areas.

Across the eruptions, the same crystal populations show similar timescales ranges (Figure 12). Simple, Type-1 crystals show the largest diffusion timescale range and Type-2 crystals with diffuse zoning recording the longest diffusion timescales observed across all the eruptions. The older more complex crystals show similar diffusion timescale ranges as the simple crystals and are not easily distinguished. The complex crystals (Type-4 and Type-5), which show resorption features, multiple zones and higher enstatite compositions than simple crystals (e.g., Type-1), suggest the crystals experienced many environmental changes prior to eruption. These environmental changes may include: changes in temperature from intrusions, changes in pressure, gas fluxing or transfer to a different magmatic environment (Martel, 2012; Moretti et al., 2013; Solaro et al., 2020). In contrast, simple crystals with one zone (Type-1, Type-2), indicate the system did not experience multiple environmental changes prior to eruption. This could perhaps explain the contrasting diffusion timescales between Plinian eruptions, with intrusions interacting differently in a cold mush system vs. one that is hot, driven by frequent melt injections.

The complex crystals found in the older eruptions have a long history with wide re-equilibration rims, higher enstatite compositions and multiple zones (Type-4 and 5). This indicates that these complex crystals have resided in a different part of the mush in contrast to the simple crystals. The presence of these complex crystals suggests the 5680 Cal. CE and ∼341 Cal. CE eruptions have sampled regions of the mush system which have not been sampled by the more recent 1010 Cal. CE and 1657 Cal. CE eruptions. The observations that these complex crystals have the same core compositions as the simple crystals, but also have higher enstatite contents in the rims and reverse zoning, suggest these crystals have resided in a region experiencing frequent increases in temperature relating to the intrusion of hot magma. Volcanic systems are interpreted to grow incrementally through intrusion emplacement from the system’s base (Annen, 2014). Crystals in regions close to the locus of intrusion emplacement would interact more frequently with intrusions, promoting regular equilibrations and resorption feature formation. Crystals farther from an intrusion are less frequently disturbed and so will undergo fewer re-equilibrations. However, this is dependent on the system geometry and intrusion location (Annen et al., 2006; Annen et al., 2015). On this basis, the La Soufrière de Guadeloupe system will have different magmatic environments with crystals recording different processes. The variability in crystal populations from different eruptions shows that they have sampled different magmatic environments (Figure 2; Figure 11). Complex crystals have been sampled by large eruptions (e.g., 5680 Cal. BCE) or in eruptions which evacuated large parts of the mush system, had significant eruption volumes, and/or fast dynamic magma ascent with possibly a less viscous melt rheology (∼341 Cal. CE).

How do the La Soufrière de Guadeloupe Timescales Fit Into a Global Context?

La Soufrière de Guadeloupe is a basaltic-andesitic system which shows, in comparison to similar systems, a wide range of timescales for magmatic processes. An estimate of the average diffusion timescale for systems in the same compositional range (50–59 wt% SiO₂) as La Soufrière de Guadeloupe, is 330 days. Though the expected timescales calculated in this study cannot be directly compared to the averages from other studies, it is still important to put La Soufrière de Guadeloupe in a global context. The 1657 Cal. CE, 1530 Cal. CE, ∼341 Cal. CE and 5680 Cal. BCE eruptions show expected timescales on the order of tens of days, which is short in a global context (Figure 13). The 1010 Cal. CE Plinian eruption shows timescales on the order of hundreds of days which is at the upper end of global diffusion timescales in the same compositional range. Despite this, the majority of diffusion timescales calculated in the present study are short compared to similar systems. This implies basaltic-andesitic to andesitic systems can rapidly produce large-scale eruptions on short timescales, from deep magmatic processes, which depending on the depth and hydrothermal system, may not necessarily be detectable by monitoring network and surface activity.

FIGURE 13. Global database of diffusion timescales (adapted from Costa et al., 2020), for relatively evolved composition eruptions. The majority of studies are done on pyroxene with two examples of olivine studies. This shows the wide range of timescales that can be calculated through diffusion methods, from hours to centuries. Systems incomparable to La Soufriére are marked with an asterisk (*), these systems generally record the longest timescales and are more typically observed in evolved caldera forming systems. The1530 C.E., 1657 Cal. CE and 5680 Cal. BCE (whole rock estimated) eruptions record some of the shortest timescales for systems in a similar compositional range, with timescales most similar to Bezymianny and Novarupta, Alaska, which show similarly short timescales. (Davydova et al., 2018; Boudon et al., 2008; Pichavant et al., 2018; Singer et al., 2016; Ruth, et al., 2018; Costa et al., 2013; Weber et al., 2019; Chamberlain et al., 2014; Costa and Chakraborty, 2004; Fabbro et al., 2013; Barker et al., 2016; Flaherty et al., 2018; Kilgour et al., 2014; Gamble et al., 1999; Saunders et al., 2014; Fabbro et al., 2013; Bourgeoisat, 2018; Solaro et al., 2020; Sundermeyer et al., 2020).

Timescale data from the Dominican magmatic systems, the island directly to the south of Guadeloupe, are available for two large Plinian eruptions, the Layou (51 ka) and Roseau (33 ka) eruptions (Solaro et al., 2020). Although magmatic systems in Dominica share similarities with those of La Grande Découverte-Soufrière de Guadeloupe, the two eruptions reported by Solaro et al., 2020 are voluminous silicic eruption with the Roseau Tuff having a DRE of 58 km³. The Layou and Roseau record eruption records timescales of 1.65–33.7 years and 1.06–61.39 years respectively, for the last event responsible for magma remobilization (Figure 13; Solaro et al., 2020). These diffusion timescales on the orders of decades are not comparable to those calculated for La Soufrière de Guadeloupe, with larger volumes of magma involved in the eruptions studied on Dominica. The main reason is that over the 10,000 years of its activity the La Soufrière magmatic system did not produce such large magma volumes. However, previous work on older phases of the Grande Découverte volcanic system (Boudon et al., 1988; Komorowski et al., 2005; Boudon et al., 2008) has shown the existence of several large Plinian eruptions associated with much larger magmatic systems. The short timescales calculated for La Soufrière de Guadeloupe are not observed at Dominica. However, future studies looking at smaller eruptions on Dominica may find similar timescales to La Soufrière de Guadeloupe.

Some systems show magmatic processes occurring on timescales comparable to La Soufriére de Guadeloupe, such as the 1912 eruption of Novarupta (Alaska), which shows similarly short diffusion timescales. A smaller range of timescales are calculated for this eruption, but the average value, 45 days, is comparable to the expected values calculated for the 1657 Cal. CE, ∼341 Cal. CE and 5680 Cal. CE eruptions (Singer et al., 2016). Despite the similarities in the timescales, the Novarupta system is much larger and is more evolved than La Soufriére de Guadeloupe, limiting the comparability.

Another system with magmatic processes occurring on timescales comparable to La Soufrière de Guadeloupe is Ruapehu volcano (New Zealand). The average timescales for eruptions ranging from 1969 to 1996 are on the order of hundreds of days, slightly longer than the majority of expected timescales calculated for La Soufrière de Guadeloupe. The minimum diffusion timescales calculated across the range of Ruapehu eruption are not as short as the minimum timescales in the range calculated for La Soufrière de Guadeloupe. Though Ruapehu is more comparable to La Soufrière de Guadeloupe than Novarupta, the eruption styles are still not entirely comparable, with timescales calculated for small, sub-Plinian, phreatomagmatic eruptions at Ruapehu.

A more comparable system is the Bezymianny volcanic system (Kamchatka, Russia), which is an andesitic complex (Figure 13). Despite showing a higher rate of magma production, Bezymianny records analogue timescales to La Soufrière de Guadeloupe, spanning from days to hundreds of days (Davydova et al., 2018; Figure 13). The range in timescales is related to two characteristic zoning types observed in the Bezymianny eruptive products, formed by two separate magmatic processes. In the first scenario, recorded by crystals with outer high-magnesium rims, deep injections from depth would likely trigger large scale magma interactions with no long periods of magma stalling in the mush system, thereby generating short timescales of days to months prior to eruption (Costa et al., 2013; Moore et al., 2014; Cassidy et al., 2015; Davydova et al., 2018). The intruding magma may only partially degas at depth, allowing the magma to rise rapidly through the mush system with a short residence time (Costa et al., 2013; Moore et al., 2014; Cassidy et al., 2015; Davydova et al., 2018). As well as being a mechanism for large Plinian eruptions, this is also interpreted as an important process in Vulcanian explosions (Costa et al., 2013; Cassidy et al., 2015) such as the 1657 Cal. CE eruption. In the second scenario, recorded in crystals with high-magnesium rhythmic zoning, recurring small volume melt injections into the upper mush system result in a longer time to destabilize the mush system and would generate timescales of weeks to years (Davydova et al., 2018). Similar processes are reported for other systems with long timescales related to multiple, smaller volumes injected over a longer period (Zellmer et al., 1999; Zellmer et al., 2003; Davydova et al., 2018; Fabbro et al., 2017; Flaherty et al., 2018). This scenario would generate more extensive degassing, and reduce the ascent rate to the shallower mush system, and require a longer time to reach a critical combined recharge threshold to trigger an eruption (Ebmeier et al., 2013; Lu and Dzurisin, 2014; Cassidy et al., 2015; Rasmussen et al., 2018).

Applying this to the 1010 Cal. CE eruption crystal populations, suggests this mechanism may occur only in the upper mush lenses and produce simple crystal populations before mush system destabilization. One possible interpretation is that a large intrusion moving rapidly through the system with little time to degas could be a triggering mechanism for an explosive Plinian eruption (Costa et al., 2013). However, the longer timescales of magmatic processes recorded by the 1010 Cal. CE Plinian eruptions indicate the magma intrusion moved slowly through the system with a longer time to degas, but still generated a large Plinian eruption. This scenario raises questions surrounding our understanding of factors controlling and modulating the explosivity of eruptible melt in the mush system.

Can We Link Timescales to Unrest at La Soufrière de Guadeloupe?

A long-standing unrest at La Soufrière de Guadeloupe started in 1992 (Komorowski et al., 2005; Villemant et al., 2014; Moretti et al., 2020), but inferring the likelihood of eruption timescales and mush dynamics from unrest signals is difficult as no magmatic eruption has been observed at this system so there are no monitoring data from a magmatic eruption available. At other well-monitored analogue volcanoes, the causative links between the mush system processes, recorded unrest signals, and diffusion timescales from corresponding erupted magmatic crystals are more easily explored. In several studies, including recent eruptions at Mt. Ruapehu (1969, 1971, 1977, 1995), Mount St. Helens (1980–1986) and Eyjafjallajökull (2010), rapidly increasing volcanic seismicity is often shown to correlate to shorter timescale data and shorter residence times (Saunders et al., 2012; Kilgour et al., 2014; Hurst et al., 2018; Pankhurst et al., 2018).

At first approximation we may relate the diffusion timescales to the amount and rate of heat needed to unlock the mush system, leading to two alternate scenarios linking the expected timescales to unrest duration and evolution prior to eruption. These two scenarios relate to different interactions between multiple magma reservoirs of heat and fluids at varying depths in the mush system and so should result in different geophysical and geochemical signals recorded at the surface.

In the first scenario, short timescales on the order of days to weeks (i.e., 1657 Cal. CE, 1530 Cal. CE and 5680 Cal. BCE) could relate to similar short periods of unrest. Here, a deep source degasses supplying heat and fluids upward before magma is injected into the shallower system without a significant long-term temperature change recorded. An abrupt increase in unrest could mark the point at which most crystals have interacted with the intrusion and may relate to the expected timescale, provided a time lag occurs between deep processes and crystals recording a perturbation in their environment. For example, these associations have been inferred at Bezymianny volcano (Kamtchaka), where high S/HCl and CO₂/HCl ratios were observed prior to the December 2009 eruption (López et al., 2013; Davydova et al., 2018).

In the second scenario, long timescales on the order of a year (i.e., 1010 Cal. CE) suggest similar long unrest periods prior to the eruption could be observed. This unrest, recorded over months or years, may show a gradual increase in seismic activity and thermal anomalies, along with geochemical markers of the rising magma batch, peaking in a shallow degassing signature following its emplacement and the thermal (and chemical) exchange with the residing body. Behaviour corresponding to this scenario was observed prior to the May 2008 Bezymianny eruption, characterised by a decrease of CO₂/H₂O, S/HCl, CO₂/S and CO₂/HCl ratios (López et al., 2013; Davydova et al., 2018).

La Soufrière de Guadeloupe is characterized by a more extensive and prolonged hydrothermal system with more effective scrubbing of magmatic gases than Bezymianny volcano, testified by lower fumarolic outlet temperatures (∼100°C vs ∼300°C) and the presence of methane (Moretti and Stefánsson, 2020). This limits any comparison of these two volcanoes, with the La Soufrière de Guadeloupe hydrothermal system that potentially screens signals related to the processes of mush unlocking that could be precursory to an eruption.

Nevertheless, looking at major episodes of phreatic activity and accelerated unrest, interpreted as aborted eruptive events related to magma emplacements at La Soufrière de Guadeloupe, gives important insight into the evolution of unrest. Major unrest episodes include the 1976–77 phreatic eruption (Feuillard et al., 1983) and the well-monitored 2018 accelerated unrest episode that culminated in the April M4.1 seismic event (Moretti et al., 2020). This episode was interpreted as a failed phreatic event related to a small-volume magmatic intrusion at depth causing a magmatic fluid pulse which overheated and over-pressurized the hydrothermal system, enhancing the supply of hot, oxidized fluids from depth resulting in peaks of CO₂/H₂O, CO₂/CH₄ and CO₂/He ratios (Moretti et al., 2020). The geochemical features of this unrest event and the rapid onset with quick ramping-up of volcano-tectonic seismicity between 1 February and 28 April, is compatible with the short-expected timescales of mush system destabilization in scenario 1.

The longest timescales recorded (e.g., 1010 Cal. CE) suggest instead ongoing long-standing unrest (scenario 2) could last years before an eruption and may commence at a level which could potentially go unnoticed. For example, the gradually increasing distal volcanic-tectonic seismic events and deep long-period seismicity, that were observed at Mt Spurr 1992, and Kilgour et al., 2014; White and McCausland, 2019), could correspond to long-standing mush destabilization, Understanding how episodes of accelerating unrest (scenario 1), as observed in 2018, superimpose over long-standing unrest (scenario 2) and how they develop with respect to the timescales for mush unlocking and magmatic processes reported in this study, is important when preparing for a future eruption. This is highlighted by the lack of correlation between timescale and eruption style, which suggests sudden periods of accelerated unrest can result in large eruptions. In this respect, future sudden increase in seismicity could indicate escalating mush destabilization and build-up to a large eruption.

Conclusion

By investigating a range of eruption styles, we provide an indepth analysis of the diffusion timescales of magmatic processes occurring at La Soufriére de Guadeloupe and provide new insights in the processes occurring in the mush system. The method constitutes a significant advance for the calculation of orthopyroxene diffusion timescales, eradicating biases induced by fitting the profiles by eye, and optimising data quality by using well-tested and robust numerical schemes and “goodness-of-fit” analyses. The six different crystal population distributions, found across the four various types of eruptions of La Soufriére de Guadeloupe, indicate different eruptions are fed by different mush system lenses, which have distinct histories. We found distinct timescales for similar eruption styles, suggesting the diffusion timescales do not allow us to discriminate between eruptive styles. In detail, we determined expected values of 22.8 ± 0.43 days for the 5680 Cal. BCE Plinian eruption, 59.4 ± 0.40 days for the ∼341 Cal. CE Strombolian eruption, 361.0 ± 0.40 days for the 1010 Cal. CE Plinian eruption and 18.8 ± 0.37 days for the 1657 Cal. CE Vulcanian eruption. The 5680 Cal. BCE and 1657 Cal. CE eruption short timescales correlate to short repose suggesting a magma intrusion hotter than the existing mush moved rapidly through the mush system due to the presence of a magma system with a high proportion of melt. The 1010 Cal. CE eruption long-expected timescale and large range of timescales indicates the system remobilised comparatively slowly with new magma interacting with the system slowly. The majority of timescales calculated in this study are short when compared to global data sets calculated for similar systems. This implies basaltic-andesitic to andesitic volcanoes can rapidly produce large-scale eruptions. Paramount for hazard assessment and crisis response, the lack of a correlation between eruption explosive intensity (VEI) and timescales that also applies to short timescales, indicates that a future eruption of La Soufriére de Guadeloupe could broadly span a range from low to high explosivity. These results underscore the necessity to further: improve the reliability of detecting and interpretingmultiparameter monitoring data as eruption precursors, expand eruption forecast modelling, develop probabilistic expert judgement for crisis response, as well as enhance risk reduction and societal resilience.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author. All of our numerical routines are made publicly available via github: https://github.com/djessop/mineral_diffusion_timescales/tree/main/SEM_traverse_data.

Author Contributions

AM prepared the samples, acquired the geochemical and petrological data, processed and analyzed all the data, calculated the timescales information, wrote the draft of the manuscript and drafted the figures. J-CK, YL, SM, and AM undertook fieldwork, and sampling. J-CK and YL carried out the reconstruction of the eruptive past and physical volcanology. GK co-processed the backscattered electron contrast data and assisted with timescale calculations and experience from analogue volcanoes. DJ formulated and performed the statistical analyses of the data which were coded in Python. RM Provided current monitoring data on La Soufrière de Guadeloupe in the context of the recent unrest. All authors discussed the data, wrote and revised the manuscript and the figures.

Funding

The PhD scholarship of AM was co-funded by the “Make our Planet Great Again” initiative of the French Government (Campus France) and the IPGP Ecole Doctorale. We thank IPGP for providing a 5-weeks visiting research position for GK in January 2020 as well as for general funding to the Observatoires Volcanologiques et Sismologiques (OVS), the INSU-CNRS for funding provided by Service National d'Observation en Volcanologie (SNOV), a Tellus-Aleas project to SM (2020, Insights into the dynamics of the active hydrothermal system of La Soufrière de Guadeloupe from past eruptions), and the Ministère pour la Transition Écologique et Solidaire (MTES) for financial support. This work has been supported by the Clervolc (UCA-LMV), the AO-IPGP 2018 project “Depth to surface propagation of fluid-related anomalies at La Soufrière de Guadeloupe volcano (FWI): timing and implications for volcanic unrest” (coord: R. Moretti), the project “Vers la Plateforme Régionale de Surveillance Tellurique du futur” (PREST) co-funded by INTERREG Caraïbes V for the European Regional Development Fund, and the European Union’s Horizon 2020 research and innovation program, under grant agreement No 731070 (EUROVOLC project). This work was supported by the Institut National des Sciences de l’Univers (INSU-CNRS) for radiocarbon dating (AO artemis, LMC14– CEA Saclay), and the CASAVA Project (PI: J-CK) funded by Agence National de la Recherche, 2009–2015 (ANR-09-RISK-02).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

The handling editor declared a past co-authorship with one of the authors RM.

Acknowledgments

This work is part of AM’s doctoral dissertation at IPGP. The authors thank the OVSG-IPGP team for logistical support and help with data collection. We are grateful for help and information from Laurence and Eric Barret (Vert Intense). We thank G. Boudon for assistance in the field and discussions over the last 3 decades on the geology of La Soufrière de Guadeloupe and the reconstruction of its eruptive past. We are grateful to the Parc National de Guadeloupe for allowing us to undertake research and obtain geological samples. We are grateful for assistance in preparing samples and for analysis from Samia Hidalgo, Pierre Burckel, Stephan Borenborensztajn, Nicolas Rividi, Michel Fialin, and Nicole Metrich. Radiocarbon dates were obtained from Beta Analytical. This study contributes to the IdEx Université de Paris ANR-18-IDEX-0001. It is IPGP contribution number 4189 and is Laboratory of excellence ClerVolc contribution number 440.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2021.617294/full#supplementary-material.

References