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

Abigail Metcalfe, Séverine Moune, Jean-Christophe Komorowski, Geoff Kilgour, David E Jessop, Roberto Moretti, Yoann Legendre

<figure><img src="image_1.png" /></figure>

## To cite this version:

Abigail Metcalfe, Séverine Moune, Jean-Christophe Komorowski, Geoff Kilgour, David E Jessop, et al. Magmatic Processes at La Soufrière de Guadeloupe: Insights From Crystal Studies and Diffusion Timescales for Eruption Onset. Frontiers in Earth Science, 2021, 9, 10.3389/feart.2021.617294. hal-03382479

HAL Id: hal-03382479

https://uca.hal.science/hal-03382479

Submitted on 18 Oct 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Check for
updates

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

Abigail Metcalfe¹*, Séverine Moune¹,²,³, Jean-Christophe Komorowski¹, Geoff Kilgour⁴, David E. Jessop¹,²,³, Roberto Moretti¹,² and Yoann Legendre⁵

¹Institut de Physique du Globe de Paris, CNRS, Université de Paris, Paris, France, ²Observatoire Volcanologique et Sismologique de Guadeloupe, Institut de Physique du Globe de Paris, Gourbeyre, France, ³CNRS, IRD, OPGC Laboratoire Magmas et Volcans, Université Clermont Auvergne, Clermont-Ferrand, France, ⁴Wairakei Research Centre, GNS Science, Taupō, New Zealand, ⁵Bureau de Recherches Géologiques et Minières, Guadeloupe, France

Signals of volcanic unrest do not usually provide insights into the timing, size and style of future eruptions, but detailed analysis of past eruptions may uncover patterns that can be used to understand future eruptive behavior. Here, we examine basaltic-andesitic to andesitic eruption deposits from La Soufrière de Guadeloupe, covering a range of eruption styles, ages and magnitudes. Our work is timely given unrest at La Soufrière de Guadeloupe has increased over the last 25 years. We constrain the timescales of magmatic processes preceding four eruptions: 1657 Cal. CE (Vulcanian), 1010 Cal. CE (Plinian), ~341 Cal. CE (Strombolian) and 5680 Cal. BCE (La Soufrière de Guadeloupe's first known Plinian eruption). Using crystal-specific analyses of diffusion in orthopyroxenes, we calculate the timescale occurring between the last recharge/mixing event in the magma reservoir and the eruption. We use backscattered electron images, coupled with EMPA of the outermost crystal rim, to derive magmatic timescales. We model the timescale populations as random processes whose probability distributions provide expected ("mean") timescales and the associated standard errors for each eruption. This provides a new statistical method for comparing magmatic timescales between disparate eruptions. From this, we obtain timescales of magma storage at La Soufrière de Guadeloupe ranging from 34.8 ± 0.4 days to 847 ± 0.4 days, with no clear distinction between eruption style/size and timescales observed. Based on these data, magmatic interaction timescales are a poor predictor of eruption style/size. This study shows that magmatic processes prior to eruption can occur on relatively short timescales at La Soufrière de Guadeloupe. Further to this basaltic-andesitic to andesitic volcanoes can rapidly produce large-scale eruptions on short timescales. These relatively short timescales calculated for volcanic processes at this system constitute a critical new data set and warrant an urgency in enhancing modeling and interpretation capabilities for near-real time monitoring data. These integrated efforts will improve early warning, eruption forecasting and crisis response management for different scenarios, as well as planning for long-term risk reduction.

**OPEN ACCESS**

**Edited by:**

Fabio Arzilli,
The University of Manchester,
United Kingdom

**Reviewed by:**

Euan James Forsyth Mutch,
University of Maryland, United States
Eugenio Nicotra,
University of Calabria, Italy

*Correspondence:

Abigail Metcalfe
metcalfe@ipgp.fr

**Specialty section:**

This article was submitted to
Volcanology,
a section of the journal
Frontiers in Earth Science

**Received:** 14 October 2020

**Accepted:** 28 January 2021

**Published:** 15 April 2021

**Citation:**

Metcalfe A, Moune S, Komorowski J-C, Kilgour G, Jessop DE, Moretti R and Legendre Y (2021) Magmatic Processes at La Soufrière de Guadeloupe: Insights From Crystal Studies and Diffusion Timescales for Eruption Onset. *Front. Earth Sci.* 9:617294
doi: 10.3389/feart.2021.617294
Keywords: diffusion timescales, orthopyroxene, magma "mush", unrest, monitoring, eruption style, forecasting, La Soufrière de Guadeloupe

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<sub>Fe-Mg</sub>) 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).

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

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).
<figure><img src="image_4.png" /></figure>

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

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
<table><caption>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.</caption><thead><tr><th>Statistical age</th><th></th><th>Vent</th><th>Dominant eruption style</th><th>VEI</th><th>Erupted volume (km³)</th><th>Erupted column height (km)</th><th>Estimated dome volume (km³)</th><th>Flank destabilization</th><th>Dome destabilisation</th><th>Directed dome blast</th></tr></thead><tbody><tr><td rowspan="2">240 years BP</td><td>1657 ± 20 CE</td><td>Soufrière</td><td>GDS 1</td><td>Vulcanian</td><td>2-3</td><td>0.001</td><td>1-15</td><td>?</td><td>?</td><td>?</td></tr><tr><td>320 years BP</td><td>1530 ± 10 CE</td><td>Soufrière</td><td>GDS 2</td><td>Sub-plinian</td><td>3</td><td>0.09</td><td>3-15</td><td>0.08 ± 0.02</td><td>x</td><td>x</td></tr><tr><td rowspan="3">1,050 years BP</td><td>1010 ± 10 CE</td><td>Soufrière</td><td>GDS 3</td><td>Plinian</td><td>4</td><td>0.2</td><td>10-25</td><td>0.17 ± 0.05</td><td>?</td><td>?</td><td>x</td></tr><tr><td>1,570 years BP</td><td>450 ± 10 CE</td><td>Soufrière</td><td>GDS 4</td><td>Dome</td><td>?</td><td>?</td><td>?</td><td>0.1</td><td>x</td><td>x</td></tr><tr><td>1,700 years BP</td><td>341 CE</td><td>Citerne</td><td>GDS 5</td><td>Strombolian</td><td>2</td><td>0.001</td><td>1-15</td><td>0.1-0.16</td><td>?</td><td>?</td></tr><tr><td rowspan="2">1,700 years BP</td><td>341 CE</td><td>Echelle</td><td>GDS 6</td><td>Strombolian</td><td>2</td><td>0.001</td><td>1-15</td><td>0.1</td><td>?</td><td>?</td><td>?</td></tr><tr><td>2,460 years BP</td><td>720 ± 10 BCE</td><td>Soufrière</td><td>GDS 7</td><td>Dome</td><td>?</td><td>0.2</td><td>?</td><td>?</td><td>x</td><td>x</td></tr><tr><td rowspan="2">?</td><td>?</td><td>Soufrière</td><td>GDS 8</td><td>Plinian</td><td>4</td><td>0.1</td><td>10-25</td><td>?</td><td>?</td><td>?</td></tr><tr><td>2,900 years BP</td><td>1080 ± 10 BCE</td><td>Soufrière</td><td>GDS 9</td><td>Plinian</td><td>4</td><td>0.1</td><td>10-25</td><td>0.29 ± 0.09</td><td>x</td><td>x</td></tr><tr><td>3,060 years BP</td><td>1370 ± 100 BCE</td><td>Soufrière</td><td>GDS 10</td><td>Dome</td><td>?</td><td>?</td><td>?</td><td>0.26 ± 0.09</td><td>x</td><td>?</td></tr><tr><td>3,495 years BP</td><td>1870 ± 170 BCE</td><td>Soufrière</td><td>GDS 11</td><td>Dome</td><td>?</td><td>?</td><td>?</td><td>0.19 ± 0.07</td><td>x</td><td>?</td></tr><tr><td>4,485 years BP</td><td>3300 ± 110 BCE</td><td>Soufrière</td><td>GDS 12</td><td>Dome</td><td>?</td><td>?</td><td>?</td><td>0.1</td><td>x</td><td>?</td></tr><tr><td>5,550 years BP</td><td>4366 ± 134 BCE</td><td>Soufrière</td><td>GDS 13</td><td>Dome</td><td>?</td><td>?</td><td>?</td><td>0.1</td><td>x</td><td>?</td></tr><tr><td>6,300 years BP</td><td>5387 ± 84 BCE</td><td>Soufrière</td><td>GDS14</td><td>Plinian</td><td>4</td><td>0.1</td><td>10-25</td><td>?</td><td>?</td><td>?</td></tr><tr><td>6,770 years BP</td><td>5680 ± 164 BCE</td><td>Soufrière</td><td>GDS15</td><td>Plinian</td><td>4</td><td>0.1</td><td>10-25</td><td>?</td><td>?</td><td>?</td></tr></tbody></table>

Highlighted in bold are eruptions for which new data is presented in this paper. Many eruptions, particularly the large explosive Plinian events are often multiphase with phases of explosive activity, effusive dome eruptions and flank collapses. CE: common Era, BP: before present = before 1950 CE).

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, advective Echelle and Citerne monogenetic cones (Boudon et al., 1988; Komorowski et al., 2005; Legendre, 2012).

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

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,
<figure><img src="image_5.png" /></figure>

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

<figure><img src="image_6.png" /></figure>

**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><caption>TABLE 2 | Summary table for the eruptions showing the key features for each eruption and sample location.</caption><thead><tr><th rowspan="2">Eruption</th><th rowspan="2">Vent</th><th colspan="2">Sample GPS position</th><th rowspan="2">Unit sampled from</th><th rowspan="2">Eruption style</th><th rowspan="2">Expected timescale (days)</th><th rowspan="2">Median timescale (days)</th><th rowspan="2">Min. timescale (days)</th><th rowspan="2">Max. timescale (days)</th><th rowspan="2">Timescale std. error (days)</th><th rowspan="2">Crystals analyzed</th><th rowspan="2">Temp (±30°C)</th><th rowspan="2">Av. Bulk SiO<sub>2</sub> (wt %)</th><th rowspan="2">VEI</th></tr><tr><th>UTM N (m)</th><th>UTM O (m)</th></tr></thead><tbody><tr><td>1657 CE</td><td>Soufrière</td><td>641,466</td><td>1,773,774</td><td>Fallout unit</td><td>Vulcanian</td><td>34.8</td><td>23.4</td><td>0.4</td><td>153.8</td><td>0.37</td><td>41</td><td>975</td><td>58.1</td><td>2–3</td></tr><tr><td>1530 CE</td><td>Soufrière</td><td>643,318</td><td>1,774,407</td><td>Fallout unit</td><td>Sub-plinian</td><td>20.5*</td><td>10</td><td>10</td><td>90</td><td></td><td>24</td><td>900</td><td>57.5</td><td>3</td></tr><tr><td>1010 CE</td><td>Soufrière</td><td>638,266</td><td>1,769,743</td><td>Fallout unit</td><td>Plinian</td><td>847.9</td><td>418.9</td><td>1.7</td><td>4,503.1</td><td>0.40</td><td>34</td><td>1,025</td><td>59.7</td><td>4</td></tr><tr><td>341 CE</td><td>Echelle</td><td>663,019</td><td>1,769,738</td><td>Scoria clasts</td><td>Strombolian</td><td>102.3</td><td>31.42</td><td>0.7</td><td>644.1</td><td>0.43</td><td>33</td><td>1,010</td><td>51.2</td><td>2</td></tr><tr><td>5680 CE</td><td>Soufrière</td><td>637,417</td><td>1,770,276</td><td>Fallout unit</td><td>Plinian</td><td>40.5</td><td>20.1</td><td>1.6</td><td>462.2</td><td>0.42</td><td>31</td><td>1,010</td><td>4</td><td>4</td></tr><tr><td>BCE</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

The crystal number refers to the number of crystals diffusion timescales were calculated on *Mean timescale not expected timescale, calculated from available data Bourgeoisat (2018); Boudon et al. (2008); Pichavant et al. (2018).

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

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

$$ \frac{\partial C}{\partial t} - \frac{\partial}{\partial x} \left( D \frac{\partial C}{\partial x} \right) = 0, \qquad (1) $$

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

$$ D(T) = D_0 \exp \left( \frac{-E_A}{RT} \right), \qquad (2) $$

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,

$$ C(x, t) = C_{\text{melt}} + \frac{C_{\text{mineral}} - C_{\text{melt}}}{2} \text{erfc} \left( \frac{x - \mu}{2\sqrt{Dt}} \right), \qquad (3) $$

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_{\text{melt}}, C_{\text{mineral}}, \mu, D, \tau\}$ 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
<figure><img src="image_7.png" /></figure>

**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.
termed the misfit function and is based on the residual between model solution, f (p) and data, d,

$$ \mathcal{F}(\mathbf{p}, \mathbf{d}) = \| (f(\mathbf{p}) - \mathbf{d}) \| \quad (4) $$

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.

As an initial estimate for p, we chose C<sub>melt</sub> and C<sub>mineral</sub> as the minimum and maximum greyscale values from the selected data, μ as the midpoint of the data range, and τ = 2 × 10<sup>6</sup> 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<sub>0</sub>, 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<sub>D</sub> = 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<sub>2</sub>) 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

<figure><img src="image_8.png" /><figcaption>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.</figcaption></figure>

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<sub>2</sub> to use for all the eruptions, particularly as the fO<sub>2</sub> has a minimal effect on the timescale calculation in comparison to the effect of temperature.

We took the value predicted by Eq. 2 using the melt-matrix temperature for that sample, T<sub>0</sub>, and fO<sub>2</sub>. While we placed no constraints on the values of C<sub>melt</sub>, C<sub>mineral</sub>, μ, and τ, we imposed upper and lower bounds on D according to T<sub>0</sub> ± 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
<figure><img src="image_9.png" /><figcaption>**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).</figcaption></figure>

assessment criteria. First, we check how closely the model fits the data using the coefficient of determination,

$$R^2 = 1 - \frac{SS_{res}}{SS_{tot}} \qquad (5)$$

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 underestimated 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),

$$ \Gamma(y; \alpha, \beta) = \frac{\beta^{\alpha} y^{\alpha-1} \exp(-\beta y)}{\Gamma(\alpha)}, \qquad (6) $$

where y are the data points, alpha and beta are the shape and rate parameters, respectively, and y(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,

$$ \mathcal{L} = \prod_{i=1}^{n} \frac{\beta^{\alpha} y_{i}^{\alpha-1} \exp(-\beta y_{i})}{\Gamma(\alpha)} \qquad (7) $$

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

$$ \begin{gathered} f(\alpha) = n(\log\alpha - \log\bar{y}) + \sum logy_i - n\Psi^{(0)}(\alpha) = 0, \quad (8) \\ \alpha/\beta = \bar{y}, \end{gathered} $$

where Ψ(n) is the nth-order polygamma function. These equations were solved to find α and β. We note the latter equation predicts that E [y] = ¯y = 1/n Σ yᵢ, the arithmetic mean of the data. For the parameter variance, we note the covariance matrix of the model parameters is given by,

$$ \mathrm{var}(\theta) = (-E[H(\theta)])^{-1} \qquad (9) $$

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,

$$ \sigma_{\tau} = \left( (\sigma_{\alpha}/\alpha)^2 + (\sigma_{\beta}/\beta)^2 \right)^{\frac{1}{2}} \qquad (10) $$

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_{91}$. 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_{91}$; 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_2$ (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_2$). This eruption also has the highest $Al_2O_3$ 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_2$).

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_2$). All major elements for this eruption have a narrow compositional range, for example $Al_2O_3$ 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_2$; 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_2$) 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**).

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
<figure><img src="image_10.png" /></figure>

**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 (Continued)

making up 48% of the crystals. There is also a substantial population of unzoned crystals, 38% of crystals are of Type-6.

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 day to 4503 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 2 to 462 days and the expected (mean) timescale calculated as 40.5 ± 0.42 days (Figure 9). The ~341 Cal. CE Strombolian eruption records a large range of timescales from <1 to 644 days. The modeled data yield an expected timescale value of 102.3 ± 0.4 days. This expected value is not comparable to any other eruption studied, though the range of values is comparable.

The 1010 Cal. CE Plinian eruption records the longest timescales from 2 to 4503 days. These data have an expected timescale value of 847.9 ± 0.4 days, considerably higher than

**FIGURE 7** | 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.
<figure><img src="image_11.png" /></figure>

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

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 134 days. The expected timescale was calculated as 34.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 41 days, with most crystals for both these eruptions recording timescales in the range of 4–6 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
<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>

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 years; 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.

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

The 5680 Cal. BCE and ~341 Cal. CE 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
<figure><img src="image_13.png" /></figure>

**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><img src="image_14.png" /></figure>

**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 system prior to eruption, this study suggests intrusions begin a between a minimum of 134 days and 4503 days prior to eruption. (B) In the case of 1657 Cal. CE and 5680 Cal. BCE eruptions, the expected value shows that the majority of the system could be affected by an intrusion 34 and 41 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 destabilize 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 system could be affected by an intrusion 102 days before eruption. This intrusion has also moved relatively rapidly through the system and variation in intrusion location, results in destabilization 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 848 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 destabilization results in eruption of populations of simply zoned and unzoned crystals, perhaps from the upper lenses.

and resulting expected timescales, suggest the intrusions interact differently with system. A transition in conduit location to Echelle or the eruption of a monogenetic cone may result in the longer expected timescale recorded for the ~341 Cal. CE eruption. The ~341 Cal. CE eruption, and to a lesser extent some 1657 Cal. CE eruption crystals, also record the shortest timescales observed across the eruption suite (≤1

day), which are not represented by the overall expected timescale (Table 2). 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 and Wörner, 2007; Johnson et al., 2008; Denis et al., 2013; Brenna et al., 2018; **Figures 11A, 1C**).
<figure><img src="image_15.png" /></figure>

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

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 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,
<figure><img src="image_16.png" /></figure>

**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. The 1530 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).

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 CO2/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 years (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. Behavior corresponding to this scenario was observed prior to the May 2008 Bezymianny eruption, characterized by a decrease of CO2/H2O, S/HCl, CO2/S and CO2/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 CO2/H2O, CO2/CH4 and CO2/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 in-depth 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 optimizing 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 40.5 ± 0.42 days for the 5680 Cal. BCE Plinian eruption, 102.3 ± 0.43 days for the ~341 Cal. CE Strombolian eruption, 847.9 ± 0.4 days for the 1010 Cal. CE Plinian eruption and 34.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 remobilized 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 interpreting multiparameter monitoring data as eruption precursors, expand eruption forecast modeling, develop probabilistic expert judgment 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

## REFERENCES

Albert, H., Costa, F., Di Muro, A., Herrin, J., Métrich, N., and Deloule, E. (2019). Magma interactions, crystal mush formation, timescales, and unrest during caldera collapse and lateral eruption at ocean island basaltic volcanoes (Piton de la Fournaise, La Réunion). *Earth Planet. Sci. Lett.* 515, 187–199. doi:10.1016/j.epsl.2019.02.035

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

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

Allan, A. S. R., Morgan, D. J., Wilson, C. J. N., and Millet, M.-A. (2013). From mush to eruption in centuries: assembly of the super-sized Oruanui magma body. *Contrib. Mineral. Petrol.* 166 (1), 143–164. doi:10.1007/s00410-013-0869-2

Allard, P., Aiuppa, A., Beauducel, F., Gaudin, D., Di Napoli, R., Calabrese, S., et al. (2014). Steam and gas emission rate from La Soufrière volcano, Guadeloupe (Lesser Antilles): implications for the magmatic supply during degassing unrest. *Chem. Geology* 384, 76–93. doi:10.1016/j.chemgeo.2014.06.019
Annen, C. (2014). Modelling the growth of volcanoes, plutons, and magma reservoirs (Doctoral dissertation, Clermont Ferrand: Université Blaise Pascal).

Annen, C., Blundy, J. D., Leuthold, J., and Sparks, R. S. J. (2015). Construction and evolution of igneous bodies: towards an integrated perspective of crustal magmatism. *Lithos* 230, 206–221. doi:10.1016/j.lithos.2015.05.008

Annen, C., Blundy, J. D., and Sparks, R. S. J. (2006). The genesis of intermediate and silicic magmas in deep crustal hot zones. *J. Petrol.* 47 (3), 505–539. doi:10.1093/petrology/egi084

Bachmann, O., and Bergantz, G. W. (2004). On the origin of crystal-poor rhyolites: extracted from batholithic crystal mushes. *J. Petrol.* 45 (8), 1565–1582. doi:10.1093/petrology/egh019

Bachmann, O., Dungan, M. A., and Lipman, P. W. (2002). The Fish Canyon magma body, San Juan volcanic field, Colorado: rejuvenation and eruption of an upper-crustal batholith. *J. Petrol.* 43, 1469–1503. doi:10.1093/petrology/43.8.1469

Bachmann, O., and Huber, C. (2016). Silicic magma reservoirs in the Earth's crust. *Am. Mineral.* 101 (11), 2377–2404. doi:10.2138/am-2016-5675

Barenblatt, G. I., and Vázquez, J. L. (2003). Nonlinear diffusion and image contour enhancement. *arXiv preprint math/0302349* 6 (1). doi:10.1017/cbo9780511814921

Barker, S. J., Wilson, C. J. N., Morgan, D. J., and Rowland, J. V. (2016). Rapid priming, accumulation, and recharge of magma driving recent eruptions at a hyperactive caldera volcano. *Geology* 44 (4), 323–326. doi:10.1130/g37382.1

Bergantz, G. W., Schleicher, J. M., and Burgisser, A. (2017). On the kinematics and dynamics of crystal-rich systems. *J. Geophys. Res. Solid Earth* 122 (8), 6131–6159. doi:10.1002/2017jb014218

Bergantz, G. W., Schleicher, J. M., and Burgisser, A. (2015). Open-system dynamics and mixing in magma mushes. *Nat. Geosci.* 8 (10), 793–796. doi:10.1038/ngeo2534

Bissainte, M. (1995). Le complexe des Monts Caraïbes de Guadeloupe (Arc des Petites Antilles) Etude pétrographique, minéralogique et géochimique. Ph.D. thesis. Paris.

Blanc, F. (1983). Corrélations chronologiques et géochimiques des formations volcaniques du sud de la Basse-Terre de Guadeloupe (Petites Antilles). Début du cycle récent. 3e cycle Thesis, Grenoble: Univ. Sci. Medic. 171.

Bonnefoi, C. C., Provost, A., and Albarède, F. (1995). The 'Daly gap' as a magmatic catastrophe. *Nature* 378, 270–272. doi:10.1038/378270a0

Boudon, G., Dagain, J., Semet, M. P., and Westercamp, D. (1988). *Notice explicative de la carte géologique au 1/20 000ème du massif volcanique de La Soufrière (Département de la Guadeloupe, Petites Antilles)*. BRGM, CNRS, DRM, IPGP, 43.

Boudon, G., Komorowski, J.-C., Villemant, B., and Semet, M. P. (2008). A new scenario for the last magmatic eruption of La Soufrière of Guadeloupe (Lesser Antilles) in 1530 A.D. Evidence from stratigraphy radiocarbon dating and magmatic evolution of erupted products. *J. Volcanol. Geothermal Res.* 178 (3), 474–490. doi:10.1016/j.jvolgeores.2008.03.006

Boudon, G., Le Friant, A., Komorowski, J.-C., Deplus, C., and Semet, M. P. (2007). Volcano flank instability in the Lesser Antilles Arc: diversity of scale, processes, and temporal recurrence. *J. Geophys. Res.* 112 (B8), B08205. doi:10.1029/2006jb004674

Bourgeoisat, C. (2018). *Spatio-temporal dynamics of the magmatic reservoir at the origin of the 1530 AD eruption of La Soufrière of Guadeloupe*. Paris: Master's thesis.

Brenna, M., Cronin, S. J., Smith, I. E. M., Tollan, P. M. E., Scott, J. M., Prior, D. J., et al. (2018). Olivine xenocryst diffusion reveals rapid monogenetic basaltic magma ascent following complex storage at Pupuke Maar, Auckland Volcanic Field, New Zealand. *Earth Planet. Sci. Lett.* 499, 13–22. doi:10.1016/j.epsl.2018.07.015

Brown, R. J., Civetta, L., Arienzo, I., D'Antonio, M., Moretti, R., Orsi, G., et al. (2014). Geochemical and isotopic insights into the assembly, evolution and disruption of a magmatic plumbing system before and after a cataclysmic caldera-collapse eruption at Ischia volcano (Italy). *Contrib. Mineralogy Petrol.* 168, 1035. doi:10.1007/s00410-014-1035-1

Burgisser, A., and Bergantz, G. W. (2011). A rapid mechanism to remobilize and homogenize highly crystalline magma bodies. *Nature* 471 (7337), 212–215. doi:10.1038/nature09799

Carlut, J., Quidelleur, X., Courtillot, V., and Boudon, G. (2000). Paleomagnetic directions and K/Ar dating of 0 to 1 Ma lava flows from La Guadeloupe Island

(French West Indies): implications for time-averaged field models. *J. Geophys. Res.* 105 (B1), 835–849. doi:10.1029/1999jb900238

Carrara, A., Burgisser, A., and Bergantz, G. W. (2019). Lubrication effects on magmatic mush dynamics. *J. Volcanol. Geothermal Res.* 380, 19–30. doi:10.1016/j.jvolgeores.2019.05.008

Cashman, K. V., Sparks, R. S., and Blundy, J. D. (2017). Vertically extensive and unstable magmatic systems: a unified view of igneous processes. *Science* 355 (6331), 3055. doi:10.1126/science.aag3055

Cassidy, M., Cole, P. D., Hicks, K. E., Varley, N. R., Peters, N., and Lerner, A. H. (2015). Rapid and slow: varying magma ascent rates as a mechanism for Vulcanian explosions. *Earth Planet. Sci. Lett.* 420, 73–84. doi:10.1016/j.epsl.2015.03.025

Chamberlain, K., Morgan, D., and Wilson, C. (2014). Timescales of mixing and mobilisation in the Bishop Tuff magma body: perspectives from diffusion chronometry. *Contrib. Mineralogy Petrol.* 168 (1). doi:10.1007/s00410-014-1034-2

Cheng, L., Costa, F., and Bergantz, G. (2020). Linking fluid dynamics and olivine crystal scale zoning during simulated magma intrusion. *Contrib. Mineralogy Petrol.* 175 (6). doi:10.1007/s00410-020-01691-3

Cherniak, D. J., and Dimanov, A. (2010). Diffusion in pyroxene, mica and amphibole. *Rev. Mineralogy Geochem.* 72, 641–690. doi:10.2138/rmg.2010.72.14

Cooper, K. (2019). Time scales and temperatures of crystal storage in magma reservoirs: implications for magma reservoir dynamics. *Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci.* 377 (2139). doi:10.1098/rsta.2018.0009

Cooper, K. M., and Kent, A. J. (2014). Rapid remobilization of magmatic crystals kept in cold storage. *Nature* 506 (7489), 480–483. doi:10.1038/nature12991

Costa, F., and Chakraborty, S. (2004). Decadal time gaps between mafic intrusion and silicic eruption obtained from chemical zoning patterns in olivine. *Earth Planet. Sci. Lett.* 227 (3–4), 517–530. doi:10.1016/j.epsl.2004.08.011

Costa, F., and Dungan, M. (2005). Short time scales of magmatic assimilation from diffusion modelling of multiple elements in olivine. *Geology* 33 (10), 837. doi:10.1130/G21675.1

Costa, F., Andreastuti, S., Bouvet de Maisonneuve, C., and Pallister, J. S. (2013). Petrological insights into the storage conditions, and magmatic processes that yielded the centennial 2010 Merapi explosive eruption. *J. Volcanol. Geothermal Res.* 261, 209–235. doi:10.1016/j.jvolgeores.2012.12.025

Costa, F., Dohmen, R., and Chakraborty, S. (2008). Time scales of magmatic processes from modeling the zoning patterns of crystals. *Rev. Mineralogy Geochem.* 69 (1), 545–594. doi:10.2138/rmg.2008.69.14

Costa, F., Shea, T., and Ubide, T. (2020). Diffusion chronometry and the timescales of magmatic processes. *Nat. Rev. Earth Environ.* 1 (4), 201–214. doi:10.1038/s43017-020-0038-x

Cowan, G. (1998). *Statistical data analysis. Oxford science publications*. ISBN: Clarendon Press. 9780198501565.

Davidson, J. P., Morgan, D. J., Charlier, B. L. A., Harlou, R., and Hora, J. M. (2007). Microsampling and isotopic analysis of igneous rocks: implications for the study of magmatic systems. *Annu. Rev. Earth Planet. Sci.* 35 (1), 273–311. doi:10.1146/annurev.earth.35.031306.140211

Davydova, V. O., Shcherbakov, V. D., and Plechov, P. Y. (2018). The timescales of magma mixing in the plumbing system of Bezymianny volcano (kamchatka): insights from diffusion chronometry. *Mosc. Univ. Geol. Bull.* 73, 444–450. doi:10.3103/s0145875218050058

de Maisonneuve, C., Costa, F., Huber, C., Vonlanthen, P., Bachmann, O., and Dungan, M. (2016). How do olivines record magmatic events? Insights from major and trace element zoning. *Contrib. Mineralogy Petrol.* 171 (6). doi:10.1007/s00410-016-1264-6

Degruyter, W., Huber, C., Bachmann, O., Cooper, K. M., and Kent, A. J. R. (2016). Magma reservoir response to transient recharge events: the case of Santorini volcano (Greece). *Geology* 44 (1), 23–26. doi:10.1130/g37333.1

Denis, C. M. M., Demouchy, S., and Shaw, C. S. J. (2013). Evidence of dehydration in peridotites from Eifel Volcanic Field and estimates of the rate of magma ascent. *J. Volcanol. Geothermal Res.* 258, 85–99. doi:10.1016/j.jvolgeores.2013.04.010

Dodge, Y. (2008). “Kolmogorov-Smirnov test”. The concise encyclopedia of statistics. *Springer New York*, 283–287. doi:10.1007/978-0-387-32833-1_214

Dohmen, R., Faak, K., and Blundy, J. D. (2017). Chronometry and speedometry of magmatic processes using chemical diffusion in olivine, plagioclase and
pyroxenes. *Rev. Mineralogy Geochem.* 83 (1), 535–575. doi:10.2138/rmg.2017.83.16

Dohmen, R., Ter heege, J. H., Becker, H.-W., and Chakraborty, S. (2016). Fe-Mg interdiffusion in orthopyroxene. *Am. Mineral.* 101, 2210–2221. doi:10.2138/am-2016-5815

Druitt, T. H., Costa, F., Deloule, E., Dungan, M., and Scaillet, B. (2012). Decadal to monthly timescales of magma transfer and reservoir growth at a caldera volcano. *Nature* 482 (7383), 77–80. doi:10.1038/nature10706

Dufek, J., and Bachmann, O. (2010). Quantum magmatism: magmatic compositional gaps generated by melt-crystal dynamics. *Geology* 38 (8), 687–690. doi:10.1130/g30831.1

Ebmeier, S. K., Biggs, J., Mather, T. A., and Amelung, F. (2013). On the lack of InSAR observations of magmatic deformation at Central American volcanoes. *J. Geophys. Res. Solid Earth* 118 (5), 2571–2585. doi:10.1002/jgrb.50195

Edmonds, M., Kohn, S., Hauri, E., Humphreys, M., and Cassidy, M. (2016). Extensive, water-rich magma reservoir beneath southern Montserrat. *Lithos*, 252–253. doi:10.1016/j.lithos.2016.02.026

Esposti Ongaro, T. (2020). Komorowski, J.-C., Legendre, Y., Neri, A Modelling pyroclastic density currents from a subplinian eruption at La Soufrière de Guadeloupe (West-Indies, France). *Bull. Volcanol.* 82, 76. doi:10.1007/s00445-020-01411-6

Fabbro, G., Druitt, T., and Scaillet, S. (2013). Evolution of the crustal magma plumbing system during the build-up to the 22-ka caldera-forming eruption of Santorini (Greece). *Bull. Volcanol.* 75 (12). doi:10.1007/s00445-013-0767-5

Fabbro, G. N., Druitt, T. H., and Costa, F. (2017). Storage and eruption of silicic magma across the transition from dominantly effusive to caldera-forming states at an arc volcano (santorini, Greece). *J. Petrol.* 58 (12), 2429–2464. doi:10.1093/petrology/egy013

Feuillard, M., Allegre, C., Brandeis, G., Gaulon, R., Le Mouel, J., Mercier, J., et al. (1983). The 1975–1977 crisis of La Soufrière de Guadeloupe (F.W.I): a still-born magmatic eruption. *J. Volcanol. Geothermal Res.* 16 (3-4), 317–334. doi:10.1016/0377-0273(83)90036-7

Feuillet, N., Beauducel, F., and Tapponnier, P. (2011). Tectonic context of moderate to large historical earthquakes in the Lesser Antilles and mechanical coupling with volcanoes. *J. Geophys. Res.* 116 (B10). doi:10.1029/2011JB008443

Feuillet, N., Manighetti, I., Tapponnier, P., and Jacques, E. (2002). Arc parallel extension and localization of volcanic complexes in Guadeloupe, Lesser Antilles. *J. Geophys. Res.* 107 (B12), 3–1. doi:10.1029/2001JB000308

Flaherty, T., Druitt, T., Tuffen, H., Higgins, M., Costa, F., and Cadoux, A. (2018). Multiple timescale constraints for high-flux magma chamber assembly prior to the Late Bronze Age eruption of Santorini (Greece). *Contrib. Mineralogy Petrol.* 173 (9). doi:10.1007/s00410-018-1490-1

Gamble, J. A., Wood, C. P., Price, R. C., Smith, I. E. M., Stewart, R. B., and Waight, T. (1999). A fifty year perspective of magmatic evolution on Ruapehu Volcano, New Zealand: verification of open system behaviour in an arc volcano. *Earth Planet. Sci. Lett.* 170 (3), 301–314. doi:10.1016/s0012-821x(99)00106-5

Ganguly, J., and Tazzoli, V. (1994). Fe2+-Mg interdiffusion in orthopyroxene: retrieval from the data on intracrystalline exchange reaction. *Am. Mineral.* 79 (9–10), 930–937.

Geist, D., Howard, K. A., and Larson, P. (1995). The generation of oceanic rhyolites by crystal fractionation: the basalt-rhyolite association at volc n Alcedo, gal pagos archipelago. *J. Petrol.* 36, 965–982. doi:10.1093/petrology/36.4.965

Gibbons, J., and Pratt, J. (1975). P-values: interpretation and methodology. *The Am. Statistician* 29 (1), 20. doi:10.2307/2683674

Ginibre, C., Worner, G., and Kronz, A. (2007). Crystal zoning as an archive for magma evolution. *Elements* 3 (4), 261–266. doi:10.2113/gselements.3.4.261

Girona, T., and Costa, F. (2013). DIPRA: a user-friendly program to model multi-element diffusion in olivine with applications to timescales of magmatic processes. *Geochem. Geophys. Geosyst.* 14 (2), 422–431. doi:10.1029/2012gc004427

Hartley, M. E., Morgan, D. J., Maclennan, J., Edmonds, M., and Thordarson, T. (2016). Tracking timescales of short-term precursors to large basaltic fissure eruptions through Fe-Mg diffusion in olivine. *Earth Planet. Sci. Lett.* 439, 58–70. doi:10.1016/j.epsl.2016.01.018

Hildreth, W. (2004). Volcanological perspectives on long valley, mammoth mountain, and mono craters: several contiguous but discrete systems.

J. Volcanol. Geothermal Res. 136 (3-4), 169–198. doi:10.1016/j.jvolgeores.2004.05.019

Hincks, T., Komorowski, J. C., Sparks, S., and Aspinall, W. (2014). Retrospective analysis of uncertain eruption precursors at La Soufrière volcano, Guadeloupe, 1975–77: volcanic hazard assessment using a Bayesian Belief Network approach. *J. Appl. Volcanol.* 3 (1). doi:10.1186/2191-5040-3-3

Huber, C., Bachmann, O., and Dufek, J. (2011). Thermo-mechanical reactivation of locked crystal mushes: melting-induced internal fracturing and assimilation processes in magmas. *Earth Planet. Sci. Lett.* 304 (3-4), 443–454. doi:10.1016/j.epsl.2011.02.022

Huber, C., Bachmann, O., and Manga, M. (2009). Homogenization processes in silicic magma chambers by stirring and mushification (latent heat buffering). *Earth Planet. Sci. Lett.* 283 (1-4), 38–47. doi:10.1016/j.epsl.2009.03.029

Hurst, T., Kilgour, G., and Hamling, I. (2018). Magmatic triggering of earthquakes on distal faults as a potential medium-term warning signal from Ruapehu volcano. *Geophys. Res. Lett.* 45 (23), 12–776. doi:10.1029/2018GL080677

Jackson, M., Blundy, J., and Sparks, R. S. J. (2018). Chemical differentiation, cold storage and remobilization of magma in the Earth's crust. *Nature* 564 (7736), 405–409. doi:10.1038/s41586-018-0746-2

Jessop, D. E., Moune, S., Moretti, R., Gibert, D., Komorowski, J.-C., Robert, V., et al. (2020). A multi-decadal view of the heat and mass budget of a volcano in unrest: La Soufrière de Guadeloupe (French West Indies). *arXiv* [preprint]. doi:10.1007/s00445-021-01439-2

Johnson, E. R., Wallace, P. J., Cashman, K. V., Granados, H. D., and Kent, A. J. (2008). Magmatic volatile contents and degassing-induced crystallization at Volcan Jorullo, Mexico: implications for melt evolution and the plumbing systems of monogenetic volcanoes. *Earth Planet. Sci. Lett.* 269 (3-4), 478–487. doi:10.1016/j.epsl.2008.03.004

Kilgour, G. N., Saunders, K. E., Blundy, J. D., Cashman, K. V., Scott, B. J., and Miller, C. A. (2014). Timescales of magmatic processes at Ruapehu volcano from diffusion chronometry and their comparison to monitoring data. *J. Volcanol. Geothermal Res.* 288, 62–75. doi:10.1016/j.jvolgeores.2014.09.010

Komorowski, J.-C., Boudon, G., Semet, M., Beauducel, F., Anténor-Habazac, C., Bazin, S., et al. (2005). in *Volcanic hazard atlas of the lesser Antilles, seismic research uni*. Editors J. M. Lindsay, R. E. A. Robertson, J. B. Shepherd, and S. Ali (Trinidad and Tobago: The University of the West Indies), 65–102.

Komorowski, J.-C. (2008). *Du volcan au pyroclaste: une approche pluridisciplinaire de la compréhension des processus éruptifs et de l'évaluation des aléas. Habilitation à Diriger des Recherches (HDR), spécialité : sciences de la Terre, de l'environnement et des planète*. Paris: Université René Diderot Paris 7, 1–549. Mémoire d'Habilitation à Diriger des Recherches, soutenue le 30 juin.

Komorowski, J.-C., Boudon, G., Antenor-Habazac, C., Hammouya, G., Semet, M., David, J., et al. (2001). *L'activité éruptive et non-éruptive de la Soufrière de Guadeloupe : problèmes et implications de la phénoménologie et des signaux actuellement enregistrés (French)*. *Eruptive and non-eruptive activity from La Soufrière de Guadeloupe: problems and implications posed by the current phenomenology and monitoring signals. Workshop on volcanic hazards - lesser Antilles volcanoes: from processes to signals*. Paris: INSU Lesser Antilles Volcanic Hazard Workshop, 18–19.

Komorowski, J.-C., Legendre, Y., Barsotti, S., Esposti Ongaro, T., Jenkins, S., Baxter, P., et al. (2013). *Assessing long term hazards for La Soufrière of Guadeloupe volcano: insights from a new eruptive chronology, credible scenario definition, and integrated impact modelling*. Kagoshima, Japan: International Association of Volcanology and Chemistry of the Earth's Interior (IAVCEI), 19–24.

Komorowski, J.-C., Legendre, Y., Boudon, G., Barsotti, S., Esposti-Ongaro, T., Jenkins, S., et al. (2012). *A new Holocene eruptive chronology for the Soufrière de Guadeloupe volcano: implications for credible scenario definition as well as hazard and impact modelling*, Colima, Mexico: Cities on Volcanoes, IAVCEI, 7. 18–23.

Kopp, H., Weinzierl, W., Becel, A., Charvis, P., Evain, M., Flueh, E., et al. (2011). Deep structure of the central lesser antilles island arc: relevance for the formation of continental crust. *Earth Planet. Sci. Lett.* 304 (1–2), 121–134. doi:10.1016/j.epsl.2011.01.024

Legendre, Y. (2012). Reconstruction de l'histoire éruptive et scenarii éruptifs à La Soufrière de Guadeloupe : vers un modèle intégrale de fonctionnement du volcan. (French) [A high resolution reconstruction of the eruptive past and definition of eruptive scenario at La Soufrière of Guadeloupe]. Ph.D. thesis. Paris.
Leone, F., Komorowski, J-C., Gherardi-Leone, M., and Lalubie, G. (2019). Integrating spatial accessibility in the design of volcano evacuation plans in the French West Indies (Guadeloupe and Martinique). *J. Appl. Volcanol.* 8, 1–22. doi:10.1186/s13617-019-0089-1

López, T., Ushakov, S., Izbekov, P., Tassi, F., Cahill, C., Neill, O., et al. (2013). Constraints on magma processes, subsurface conditions, and total volatile flux at Bezymianny Volcano in 2007–2010 from direct and remote volcanic gas measurements. *J. Volcanol. Geotherm. Res.* 263, 92–107.

Lu, Z., and Dzurisin, D. (2014). *InSAR imaging of aleutian volcanoes: monitoring a volcanic arc from space*. Springer Berlin Heidelberg, 87–345.

Macdonald, R., Belkin, H. E., Fitton, J. G., Rogers, N. W., Nejbert, K., Tindle, A. G., et al. (2008). The roles of fractional crystallization, magma mixing, crystal mush remobilization and volatile-melt interactions in the genesis of a young basalt-peralkaline rhyolite suite, the greater olkaria volcanic complex, Kenya rift valley. *J. Petrol.* 49, 1515–1547. doi:10.1093/petrology/egn036

Macdonald, R., Hawkesworth, C. J., and Heath, E. (2000). The Lesser Antilles volcanic chain: a study in arc magmatism. *Earth-Science Rev.* 49, 1–76. doi:10.1016/s0012-8252(99)00069-0

Mader, H. M., Llewellin, E. W., and Mueller, S. P. (2013). The rheology of two-phase magmas: a review and analysis. *J. Volcanol. Geothermal Res.* 257, 135–158. doi:10.1016/j.jvolgeores.2013.02.014

Mahood, G. A. (1990). Second reply to comment of R.S.J. Sparks, H.E. Huppert and C.J.N. Wilson on "Evidence for long residence times of rhyolitic magma in the Long Valley magmatic system: the isotopic record in the precaldera lavas of Glass Mountain". *Earth Planet. Sci. Lett.* 99, 395–399. doi:10.1016/0012-821x(90)90145-n

Martel, C. (2012). Eruption dynamics inferred from microlite crystallization experiments: application to plinian and dome-forming eruptions of Mt. Pelee (Martinique, lesser antilles). *J. Petrol.* 53, 699–725. doi:10.1093/petrology/egr076

Moore, A., Coogan, L., Costa, F., and Perfit, M. (2014). Primitive melt replenishment and crystal-mush disaggregation in the weeks preceding the 2005–2006 eruption950′N, EPR. *Earth Planet. Sci. Lett.* 403, 15–26. doi:10.1016/j.epsl.2014.06.015

Moretti, R., Komorowski, J. C., Ucciani, G., Moune, S., Jessop, D. E., de Chabalier, J., et al. (2020). The 2018 unrest phase at La Soufrière de Guadeloupe (French West Indies) andesitic volcano: scrutiny of a failed but prodromal phreatic eruption. *J. Volcanol. Geothermal Res.* 393, 106769. doi:10.1016/j.jvolgeores.2020.106769

Moretti, R., Arienzo, I., Di Renzo, V., Orsi, G., Arzilli, F., Brun, F., et al. (2019). Volatile segregation and generation of highly vesiculated explosive magmas by volatile-melt fining processes: the case of the Campanian Ignimbrite eruption. *Chem. Geology* 503, 1–14. doi:10.1016/j.chemgeo.2018.10.001

Moretti, R., Arienzo, I., Orsi, G., Civetta, L., and D'Antonio, M. (2013). The deep plumbing system of ischia: a physico-chemical window on the fluid-saturated and CO2-sustained neapolitan volcanism (southern Italy). *J. Petrol.* 54, 951–984. doi:10.1093/petrology/egt002

Moretti, R., and Stefánsson, A. (2020). Volcanic and geothermal redox engines. *Elements* 16, 179–184. doi:10.2138/gselements.16.3.179

Morgan, D., Blake, S., Rogers, N., DeVivo, B., Rolandi, G., Macdonald, R., et al. (2004). Time scales of crystal residence and magma chamber volume from modelling of diffusion profiles in phenocrysts: vesuvius 1944. *Earth Planet. Sci. Lett.* 222 (3–4), 933–946. doi:10.1016/j.epsl.2004.03.030

Nakamura, M. (1995). Continuous mixing of crystal mush and replenished magma in the ongoing Unzen eruption. *Geol* 23, 807–810. doi:10.1130/0091-7613(1995)023<0807:cmocma>2.3.co;2

Newhall, C. G., and Self, S. (1982). The volcanic explosivity index (VEI) an estimate of explosive magnitude for historical volcanism. *J. Geophys. Res.* 87, 1231–1238. doi:10.1029/jc087ic02p01231

OVSG-IPGP (1999-2020). *Monthly reports on the activity of La Soufrière de Guadeloupe and on regional seismicity*. Gourbeyre: Guadeloupe http://www.igpg.fr/fr/ovsg/bulletins-mensuels-de-lovsg

Pankhurst, M. J., Morgan, D. J., Thordarson, T., and Loughlin, S. C. (2018). Magmatic crystal records in time, space, and process, causatively linked with volcanic unrest. *Earth Planet. Sci. Lett.* 493, 231–241. doi:10.1016/j.epsl.2018.04.025

Peccerillo, A., Barberio, M. R., Yirgu, G., Ayalew, D., Barbieri, M. W., and Wu, T. W. (2003). Relationships between mafic and peralkaline silicic magmatism in

continental rift settings: a petrological, geochemical and isotopic study of the Gedemsa volcano, central Ethiopian rift. *J. Petrol.* 44 (11), 2003–2032. doi:10.1093/petrology/egg068

Pichavant, M., Poussineau, S., Lesne, P., Solaro, C., and Bourdier, J.-L. (2018). Experimental parametrization of magma mixing: application to the ad 1530 eruption of La Soufrière, Guadeloupe (lesser antilles). *J. Petrol.* 59 (2), 257–282. doi:10.1093/petrology/egy030

Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P. (2007). *Numerical recipes: the art of scientific computing*. Third. ISBN: Cambridge University Press. 9780521880688.

Putirka, K. D. (2008). Thermometers and barometers for volcanic systems. *Rev. Mineralogy Geochem.* 69 (1), 61–120. doi:10.2138/rmg.2008.69.3

Pyle, D. M. (2000). "Sizes of volcanic eruptions," in *Encyclopedia of volcanoes*. Editors H. Sigurdsson, B. Houghton, H. Reimer, J. STix, and S. McNutt (San Diego: Academic Press), 263–269.

Rad, S., Rivé, K., Vittecoq, B., Cerdan, O., and Allègre, C. J. (2013). Chemical weathering and erosion rates in the Lesser Antilles: an overview in Guadeloupe, Martinique and Dominica. *J. South Am. Earth Sci.* 45, 331–344. doi:10.1016/j.jsames.2013.03.004

Rasmussen, D. J., Plank, T. A., Roman, D. C., Power, J. A., Bodnar, R. J., and Hauri, E. H. (2018). When does eruption run-up begin? Multidisciplinary insight from the 1999 eruption of Shishaldin volcano. *Earth Planet. Sci. Lett.* 486, 1–14. doi:10.1016/j.epsl.2018.01.001

Rubin, A. E., Cooper, K. M., Till, C. B., Kent, A. J. R., Costa, F., Bose, M., et al. (2017). Rapid cooling and cold storage in a silicic magma reservoir recorded in individual crystals. *Science* 356 (6343), 1154–1156. doi:10.1126/science.aam8720

Ruprecht, P., and Wörner, G. (2007). Variable regimes in magma systems documented in plagioclase zoning patterns: el Misti stratovolcano and Andahua monogenetic cones. *J. Volcanol. Geothermal Res.* 165 (3–4), 142–162. doi:10.1016/j.jvolgeores.2007.06.002

Ruth, D., Costa, F., Bouvet de Maisonneuve, C., Franco, L., Cortés, J., and Calder, E. (2018). Crystal and melt inclusion timescales reveal the evolution of magma migration before eruption. *Nat. Commun.* 9 (1). doi:10.1038/s41467-018-05086-8

Ruzié, L., Aubaud, C., Moreira, M., Agrinier, P., Dessert, C., Gréau, C., et al. (2013). Carbon and helium isotopes in thermal springs of La Soufrière volcano (Guadeloupe, Lesser Antilles): implications for volcanological monitoring. *Chem. Geology* 359, 70–80. doi:10.1016/j.chemgeo.2013.09.008

Ruzié, L., Moreira, M., and Crispi, O. (2012). Noble gas isotopes in hydrothermal volcanic fluids of La Soufrière volcano, Guadeloupe, Lesser Antilles arc. *Chem. Geology* 304-305, 158–165. doi:10.1016/j.chemgeo.2012.02.012

Samper, A., Quidelleur, X., Komorowski, J., Lahitte, P., and Boudon, G. (2009). Effusive history of the Grande Découverte volcanic complex, southern basse-terre (Guadeloupe, French west Indies) from new K-Ar cassignol-gillot ages. *J. Volcanol. Geothermal Res.* 187 (1–2), 117–130. doi:10.1016/j.jvolgeores.2009.08.016

Samper, A., Quidelleur, X., Lahitte, P., and Mollex, D. (2007). Timing of effusive volcanism and collapse events within an oceanic arc island: Basse-Terre, Guadeloupe archipelago (Lesser Antilles Arc). *Earth Planet. Sci. Lett.* 258 (1–2), 175–191. doi:10.1016/j.epsl.2007.03.030

Saunders, K., Blundy, J., Dohmen, R., and Cashman, K. (2012). Linking petrology and seismology at an active volcano. *Science* 336 (6084), 1023–1027. doi:10.1126/science.1220066

Saunders, K., Buse, B., Kilburn, M. R., Kearns, S., and Blundy, J. (2014). Nanoscale characterisation of crystal zoning. *Chem. Geology* 364, 20–32. doi:10.1016/j.chemgeo.2013.11.019

Schleicher, J. M., Bergantz, G. W., Breidenthal, R. E., and Burgisser, A. (2016). Time scales of crystal mixing in magma mushes. *Geophys. Res. Lett.* 43 (4), 1543–1550. doi:10.1002/2015gl067372

Schleicher, J. M., and Bergantz, G. W. (2017). The mechanics and temporal evolution of an open-system magmatic intrusion into a crystal-rich magma. *J. Petrol.* 58 (6), 1059–1072. doi:10.1093/petrology/egx045

Schwandt, C. S., Cygan, R. T., and Westrich, H. R. (1998). Magnesium self-diffusion in orthoenstatite. *Contrib. Mineral. Petrol.* 130 (3–4), 390–396.

Shea, T., Lynn, K. J., and Garcia, M. O. (2015). Cracking the olivine zoning code: distinguishing between crystal growth and diffusion. *Geology* 43 (10), 935–938. doi:10.1130/g37082.1
Singer, B. S., Costa, F., Herrin, J. S., Hildreth, W., and Fierstein, J. (2016). The timing of compositionally-zoned magma reservoirs and mafic 'priming' weeks before the 1912 Novarupta-Katmai rhyolite eruption. *Earth Planet. Sci. Lett.* 451, 125–137. doi:10.1016/j.epsl.2016.07.015

Sliwinski, J. T., Bachmann, O., Dungan, M. A., Huber, C., Deering, C. D., Lipman, P. W., et al. (2017). Rapid pre-eruptive thermal rejuvenation in a large silicic magma body: the case of the Masonic Park Tuff, Southern Rocky Mountain volcanic field, CO, USA. *Contrib. mineralogy Petrol.* 172 (5), 30. doi:10.1007/s00410-017-1351-3

Solaro, C., Balcone-Boissard, H., Morgan, D. J., Boudon, G., Martel, C., and Ostorero, L. (2020). A system dynamics approach to understanding the deep magma plumbing system beneath Dominica (Lesser Antilles). *Front. Earth Sci.* 8. doi:10.3389/feart.2020.574032

Sparks, R. S. J., and Cashman, K. V. (2017). Dynamic magma systems: implications for forecasting volcanic activity. *Elements* 13 (1), 35–40. doi:10.2113/gselements.13.1.35

Spera, F. J., and Bohrson, W. A. (2018). Rejuvenation of crustal magma mush: a tale of multiply nested processes and timescales. *Am. J. Sci.* 318, 90–140. doi:10.2475/01.2018.05

Stelten, M. E., Cooper, K. M., Vazquez, J. A., Calvert, A. T., and Glessner, J. J. G. (2015). Mechanisms and timescales of generating eruptible rhyolitic magmas at Yellowstone Caldera from zircon and sanidine geochronology and geochemistry. *J. Petrol.* 56, 1607–1642. doi:10.1093/petrology/egv047

Streck, M. J. (2008). Mineral textures and zoning as evidence for open system processes. *Rev. Mineralogy Geochem.* 69 (1), 595–622. doi:10.2138/rmg.2008.69.15

Streckeisen, A., Zanettin, B., Le Bas, M. J., Bonin, B., Bateman, P., Bellieni, G., et al. (2002). “Igneous rocks,” in *A classification and glossary of terms*. Editor R. W. Le Maitre (Cambridge University Press), 236.

Sundermeyer, C., Gätjen, J., Weimann, L., and Wörner, G. (2020). Timescales from magma mixing to eruption in alkaline volcanism in the Eifel volcanic fields, western Germany. *Contrib. Mineralogy Petrol.* 175 (8), 1–23. doi:10.1007/s00410-020-01715-y

Touboul, M., Bourdon, B., Villemant, B., Boudon, G., and Joron, J. (2007). 238U-230Th-226Ra disequilibria in andesitic lavas of the last magmatic eruption of Guadeloupe Soufriere, French Antilles: processes and timescales of magma differentiation. *Chem. Geology* 246 (3–4), 181–206. doi:10.1016/j.chemgeo.2007.09.009

Ubide, T., and Kamber, B. (2018). Volcanic crystals as time capsules of eruption history. *Nat. Commun.* 9 (1). doi:10.1038/s41467-017-02274-w

Vance, J. A. (1965). Zoning in igneous plagioclase: patchy zoning. *J. Geology* 73 (4), 636–651. doi:10.1086/627099

Villemant, B., Hammouya, G., Michel, A., Semet, M., Komorowski, J. C., Boudon, G., et al. (2005). The memory of volcanic waters: shallow magma degassing

revealed by halogen monitoring in thermal springs of La Soufrière volcano (Guadeloupe, Lesser Antilles). *Earth Planet. Sci. Lett.* 237 (3–4), 710–728. doi:10.1016/j.epsl.2005.05.013

Villemant, B., Komorowski, J. C., Dessert, C., Michel, A., Crispi, O., Hammouya, G., et al. (2014). Evidence for a new shallow magma intrusion at La Soufrière of Guadeloupe (Lesser Antilles). *J. Volcanol. Geothermal Res.* 285, 247–277. doi:10.1016/j.jvolgeores.2014.08.002

Weber, G., Arce, J. L., Ulianov, A., and Caricchi, L. (2019). A recurrent magmatic pattern on observable timescales prior to plinian eruptions from nevado de Toluca (Mexico). *J. Geophys. Res. Solid Earth* 124, 10999–11021. doi:10.1029/2019jb017640

White, R. A., and McCausland, W. A. (2019). A process-based model of pre-eruption seismicity patterns and its use for eruption forecasting at dormant stratovolcanoes. *J. Volcanol. Geothermal Res.* 382, 267–297. doi:10.1016/j.jvolgeores.2019.03.004

Zellmer, G. F., Blake, S., Vance, D., Hawkesworth, C., and Turner, S. (1999). Plagioclase residence times at two island arc volcanoes (Kameni Islands, Santorini, and Soufriere, St. Vincent) determined by Sr diffusion systematics. *Contrib. Mineralogy Petrol.* 136 (4), 345–357. doi:10.1007/s004100050543

Zellmer, G. F., Sparks, R. S. J., Hawkesworth, C. J., and Wiedenbeck, M. (2003). Magma emplacement and remobilization timescales beneath Montserrat: insights from Sr and Ba zonation in plagioclase phenocrysts. *J. Petrol.* 44, 1413–1431. doi:10.1093/petrology/44.8.1413

Zlotnicki, J., Boudon, G., and Le Mouël, J.-L. (1992). The volcanic activity of La Soufrière of Guadeloupe (lesser antilles): structural and tectonic implications. *J. Volcanol. Geothermal Res.* 49, 91–104. doi:10.1016/0377-0273(92)90006-y

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

Copyright © 2021 Metcalfe, Moune, Komorowski, Kilgour, Jessop, Moretti and Legendre. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.