Impact Factor 2.689 | CiteScore 3.3
More on impact ›

Original Research ARTICLE

Front. Earth Sci., 19 November 2020 |

A System Dynamics Approach to Understanding the deep Magma Plumbing System Beneath Dominica (Lesser Antilles)

  • 1Université de Paris, Institut de Physique du Globe de Paris, CNRS, Paris, France
  • 2Institut des Sciences de la Terre de Paris (ISTeP), UMR 7193, CNRS-Sorbonne Université, Paris, France
  • 3Institute of Geophysics and Tectonics, School of Earth & Environment, University of Leeds, Leeds, United Kingdom
  • 4Institut des Sciences de la Terre d’Orléans (ISTO), UMR 7327, Université d'Orléans-CNRS/INSU-BRGM, Orléans, France

To understand the dynamics of magmatic systems, one must first seek to characterize the time-dependent behavior of magma storage and ascent. Herein, we do this through a combination of the Crystal System Approach and careful study of Fe-Mg interdiffusion in orthopyroxene. This allows us to trace the pre-eruptive dynamics of magma plumbing systems, both in space and time. We apply this novel approach on two large silicic eruptions (about 3–5 km3 DRE/eruption) that occurred in the central part of Dominica Island (Lesser Antilles Arc): the eruptions of Layou (∼51 ka) from Morne Diablotins, and Roseau (∼33 ka) from Morne Trois Pitons-Micotrin. For the Roseau eruption, two magmatic environments (MEs) are identified on the basis of orthopyroxene composition, with a dominant reverse-zoning pattern from 50 to 54 to 54–59 mol% enstatite (En), indicating interaction with hotter magma. For the Layou eruption, three MEs are observed as represented by three populations of pyroxenes: En47-51, En51-53 and En53-58. The normal-zoning pathway from En51-53 to En47-51 is significantly registered by crystals, interpreted as convective mixing in a zoned reservoir. The reverse-zoning pathway from En47-51 to En51-53 and also En53-58 is also significantly present, supporting the mixing within the zoned reservoir but also suggesting mixing with a hotter magma, possibly stored in another part of a sub-volcanic mush. The crystal and glass compositions (melt inclusion and matrix glass) from both studied eruptions suggest heating and mixing between different magma pockets located within the mush that were the dominant process for mobilizing eruptible magma. In parallel, we constrain the associated pre-eruptive timescales by modeling the diffusive relaxation of Fe-Mg chemical gradients that originated within the zonation of the same orthopyroxene crystals. Diffusion modeling was considered along the b-axis of 66 zoned orthopyroxene crystals for these two eruptions, at a magmatic temperature of 850 ± 25°C. In light of these results, we propose that the Layou and the Roseau magma reservoirs were rejuvenated and heated by ∼25–50°C about 10 years prior to eruption by the injection of an underplating, hotter magma, creating the observed dominant reverse-zoning patterns of the erupted orthopyroxenes. We thus have evidence that silicic mush can be re-mobilized over timescales of decades prior to eruption, as previously suggested for Santorini and Taupo volcanoes.


How magma is stored at depth and how magma bodies assemble prior to eruption are two questions that are still widely debated. The architecture and associated dynamics of magma plumbing systems for silicic volcanoes have been in the last decade envisaged as transcrustal magma systems (Cashman et al., 2017) following the pioneering work of Bachmann and Bergantz (2003).

In the process of sequential accretion of sills by which large crustal magma reservoirs form (Annen, 2009; Bachmann and Huber, 2016; Blundy and Annen, 2016), the resident magma can cool and partially crystallize in between successive magma pulses to become a crystal mush (Marsh, 2006). Crystallinity stays in the 50–60% range and the solidus is not reached quickly because of the release of latent heat during crystallization (Marsh, 1989). This type of highly-crystallized crustal reservoir raises questions regarding magma mobility and the timescale over which magma can be remobilized through the rejuvenation of the mush to become eruptible melt. Several mechanisms have been proposed to enable mush rejuvenation and possible melt extraction: self-mixing (Couch et al., 2001), gas sparging (Bachmann and Bergantz, 2003; Bachmann and Huber, 2016; Huber et al., 2010), or convective stirring (Huber et al., 2009). In any case, underplating mafic magmas bringing heat and/or volatiles are considered as major players in the reactivation of mush. Heat and volatile transfers change the rheological properties of the mush (buoyancy, crystallinity), favoring its mobility (formation of melt connections, volatile foams). The efficiency of the enthalpy (±mass) exchange with newly injected recharge magma of higher enthalpy necessary to unlock a crystal mush highly depends not only on the size but also on the composition and temperature of the intruding magma body. Estimated timescales suggest that the reactivation timescale is much longer than the following homogenization timescale (Huber et al., 2012).

Detailed petrological investigations of the compositions of crystals, their melt inclusions and the residual melt provide strong constraints on the dynamics of the magmatic plumbing system in terms of architecture and timescales. Although such studies are major undertakings, they are highly informative (Ruth et al., 2018; Albert et al., 2019). Using the method developed by Palm (2005), the pre-eruptive dynamics of recent eruptions of Mt. Etna has been explored using the sequential variations in olivine crystal compositions, represented schematically in a box model and called Crystal System Analysis (CSA) (Kahl et al., 2011; Kahl et al., 2015). Coupling CSA with diffusion timescales can be used to show variable evolution of magmatic interactions prior to eruption and showed variable magmatic pathway systematics between eruptions. In recent decades, intracrystalline diffusion modeling has been increasingly used to study timescales of pre-eruptive processes, built on the pioneering works of Costa et al. (2003), Costa and Chakraborty (2004) and Chakraborty (2008). Because of the various diffusion rates of elements in crystals (Costa et al., 2020), this technique allows investigation of timescales spanning minutes (Charlier et al., 2012) to thousands of years (Allan et al., 2017) and, being applicable to single crystals, it can yield timescales of specific processes. Several studies, covering various volcanic settings and involving different compositions of magmas and eruptive dynamics have already been conducted with different element/crystal pairs, including olivine (Costa and Chakraborty, 2004; Costa and Morgan, 2010; Kahl et al., 2011; Kahl et al., 2015), plagioclase (Zellmer et al., 2003; Zellmer et al., 2003; Ginibre et al., 2004; Ruprecht and Wörner, 2007; Saunders et al., 2010; Shcherbakov et al., 2011; Druitt et al., 2012), pyroxenes (Morgan et al., 2004; Allan et al., 2013; Barker et al., 2016), magnetite (Boudon et al., 2015) and quartz (Matthews et al., 2012; Chamberlain et al., 2014).

In the Lesser Antilles Arc, the volume of magma involved in most of the eruptions is relatively low along the volcanic arc (generally <0.5 km3) as evidenced by many historical and prehistorical eruptions (Lindsay et al., 2005). Dominica is an exception: it is the sole island of the arc with several active volcanoes. Together with the neighboring islands of Guadeloupe and Martinique, it concentrates half of the active volcanoes of the arc within its central part, only accounting for 1/5 of its total length (Figure 1A). The pioneer works on the volcanic activity of Dominica recognized a voluminous silicic eruption, dated between 50 and 30 ka and called “the Roseau Tuff” (Sigurdsson, 1972; Carey and Sigurdsson, 1980). Approximatively 58 km3 Dense Rock Equivalent (DRE) of magma has been first estimated for this eruption, making it as the most voluminous eruption of the arc. Though recent works now recognized several eruptions, in particular in the central part of Dominica (Lindsay et al., 2005; Smith et al., 2013; Boudon et al., 2017), with necessarily lower emitted volumes, these volumes still contrast with those recognized for the other active volcanoes. For instance, the ongoing eruption of Soufrière Hill’s of Montserrat only emitted 1 km3 of magma (Wadge et al., 2014) and the famous devastating eruption of Montagne Pelée (Martinique) in 1902–1905 was responsible of the emission of only 0.2–0.3 km3 (Boudon, 1983). Dominica volcanic activity thus offers the unique opportunity to investigate the magma dynamics in the plumbing system prior to repeated large silicic eruptions, only recognized on this island.


FIGURE 1. Dominica and studied ignimbritic eruptions. (A) Lesser Antilles arc, with Dominica in the central part of the arc. (B) Dominica and location of the sampled outcrops (stars) for Layou and Roseau ignimbrites, after Boudon et al. (2017). Red: deposits from Morne Trois Pitons–Micotrin eruptive center, with dark red deposits corresponding to Roseau eruption. Blue: deposits from Morne Diablotins eruptive center, with dark blue deposits corresponding to Layou eruption. (C) Layou Plinian fallout deposit. (D) Roseau fallout deposit. (E) Layou pyroclastic density current deposit. (F) Roseau pyroclastic density current deposit.

Here, the aim is to identify and characterize the pre-eruptive disrupting events leading to two successive large silicic eruptions (∼3–5 km3 DRE; Boudon et al., 2017) in the central part of Dominica, Layou (∼51 ka) and Roseau (∼33 ka), and their associated timescales. For that purpose, we present an original study using orthopyroxenes that combine a variant of the Crystal System Analysis method and intracrystalline diffusion modeling. More precisely, we performed a detailed petrological investigation of pumices ejected during the first eruptive Plinian phase of these two major eruptions. We focus our analysis on orthopyroxenes, which crystallize over a wide range of temperature, pressure and oxygen fugacity. Thus, orthopyroxene crystals may be good tracers of evolving crustal storage conditions. This work builds upon a revised chronostratigraphy of the eruptions (Boudon et al., 2017), experimental investigations of the pre-eruptive magma storage conditions (Solaro et al., 2019), and a detailed study of major, trace, and volatile element compositions of melt inclusions (Balcone-Boissard et al., 2018). The results are discussed for their implications in terms of the more general framework of the magma plumbing system architecture and dynamics beneath Dominica. Furthermore, we assess the rates of silicic magma rejuvenation prior to voluminous explosive volcanic eruptions.

Geological Context and Previous Work

Geological Context

The Lesser Antilles Arc results from the westward subduction of the North-South American plates beneath the Caribbean plate, at an average rate of 2 cm/yr (Wadge, 1984). The volcanic arc is 850 km long, and extends from the island of Saba in the north to Grenada in the south, with an eastward convexity. Dominica is located in the central part of the arc (Figure 1).

North of Dominica, the arc splits into two branches: the Outer or Older arc (Eocene to mid-Oligocene) to the north-east and an Inner or Recent arc (Miocene to Recent) trending north-west. This jump in the northern sector has been ascribed to the ancient subduction of an aseismic ridge in the center of the arc causing a rupture in the subducting plate and a displacement of volcanic activity to the west in the early-Miocene (∼20 Ma; Bouysse and Westercamp, 1990). South of Dominica, the two branches merge, creating a single volcanic arc active since the Miocene. The active parts of the arc (Saba to Dominica in the North, and Dominica to Grenada in the South) have twelve, currently-active volcanic centers, the best-known being Soufrière (Guadeloupe), Montagne Pelée (Martinique), Soufrière Hills (Montserrat) and Soufrière (St.-Vincent). Conversely to the other islands, Dominica is characterized by four, potentially-active volcanic centers, from North to South: Morne aux Diables, Morne Diablotins, Morne Trois Pitons-Micotrin and Morne Plat Pays Volcanic Complex (Lindsay et al., 2005; Figure 1B).

Five ignimbritic eruptions are currently recognized on Dominica: Grand Bay, Roseau, Grand Fond, Layou and Grande Savane ignimbrites (Lindsay et al., 2005; Smith et al., 2013; Howe et al., 2014; Boudon et al., 2017; Figure 1B). These deposits have been recently re-investigated by Boudon et al. (2017) who derive new 14C ages for the Layou, Roseau and Grand Fond ignimbrites respectively at ∼51, ∼33 and ∼24 ka; each eruption involved magma volumes of the order of 3–5 km3 DRE. On the basis of geological arguments and different geochemical signatures observed in Roseau magmas compared to Layou ones, the Roseau ignimbrite originates from Morne Trois Pitons–Micotrin. The Layou ignimbrite likely originated from the Morne Diablotins eruptive center. The deposit of each eruption is subdivided into a basal Plinian fallout (<10% of total erupted volumes) (Figures 1C,D), overlain by thick (up to 100 m) pumiceous pyroclastic density current deposits that vary in grain size, welding and thickness distribution (Figures 1E,F).

Previous Works: Petrology of Pumice Clasts and Pre-eruptive Magma Storage Conditions


Whole rock, melt inclusions, and residual glass compositions (major and trace elements) have been studied together with the crystal contents of the natural samples of the Roseau and Layou eruptions (Boudon et al., 2017; Balcone-Boissard et al., 2018). Magma of the Roseau eruption is intermediate between andesite and dacite (62–63 wt% SiO2, ∼4.5–5.0 wt% Na2O + K2O) whereas magma from the Layou eruption is dacitic (>64 wt% SiO2, ∼5 wt% Na2O + K2O). For both eruption deposits, the residual glass is rhyolitic (76.5–78.0 wt% SiO2, 5.5–7.0 wt% Na2O + K2O). The samples are porphyritic with ∼30 ± 5 vol% phenocrysts (Boudon et al., 2017) consisting of plagioclase, orthopyroxene, amphibole, clinopyroxene, Ti-magnetite and ilmenite embedded in a microlite-free residual glass (Solaro et al., 2019). Modal composition of the magma of each eruption shows that plagioclase is the most abundant phenocryst (70%), followed by orthopyroxene (15–17%), amphibole (up to 17% for Layou, though <1% for Roseau) and, in lesser proportions, clinopyroxene and Fe-Ti oxides. Rare quartz phenocrysts are also recognized in the Layou deposit. Plagioclase crystals show oscillatory zoning patterns and resorption textures, with anorthite content (An, expressed in mol%) spanning a broad range between An42 to An90 (Solaro et al., 2019), with some high An spikes. Amphiboles in Layou are ferri-magnesio-hornblende with Mg# [Mg/(Mg + Fetot)] ranging from 0.52 to 0.60. Clinopyroxenes are augites with Mg# ranging from 0.48 to 0.71. Fe-Ti oxides are titanomagnetites (Mag71–75; Mag = magnetite mol%) and ilmenites (Ilm87; Ilm = ilmenite, mol%). The orthopyroxene compositions, the focus of this study, are presented hereafter.

Melt Inclusions

Only glassy Melt Inclusions (MIs) with no post-entrapment modifications have been investigated, as detailed in Balcone-Boissard et al. (2018). MIs are ubiquitous in orthopyroxenes (sizes up to ∼50 µm). In plagioclase, scarce MIs are smaller than 40 µm. The MIs are all rhyolitic regardless of their host crystal, with some variations in Na2O + K2O and SiO2 contents. MIs are among the most water-rich yet recorded (up to 8 wt% H2O). Carbon dioxide contents are generally low (<650 ppm), although MIs trapped at high pressures contain higher CO2 contents (up to 3,500–4,000 ppm). Some low-pressure (<100 MPa) inclusions have elevated CO2 contents (up to 1,100–1,150 ppm), suggestive of fluxing of shallow magmas with CO2-rich fluids. CO2-trace element systematics indicate that melts were volatile-saturated at the time of entrapment (Balcone-Boissard et al., 2018).

Magma Storage Conditions

Phase equilibrium experiments starting from glasses of Roseau and Layou bulk compositions have been performed in an internally-heated pressure vessel, spanning a temperature range of ∼800–900°C, a pressure range of 200–400 MPa (under water saturation and under-saturation), and fixed oxygen fugacity of ≈ NNO +0.6 (NNO = Nickel- Nickel oxide buffer) (Solaro et al., 2019). By experimentally reproducing the natural melt and phenocryst rim compositions, phenocryst assemblages, phase proportions, and water contents in MIs, the magma storage conditions have been determined to be close to a pressure of ∼400 MPa, a temperature of 850 ± 25°C with 7–8 wt% H2O dissolved in melt. These parameters will be used for the diffusion modeling below.

Material and Methods

Sample Types and Preparation

Studied pumices and crystals come from the field campaigns performed previously for chronostratigraphic, melt inclusions and experimental petrological studies (Boudon et al., 2017; Balcone-Boissard et al., 2018; Solaro et al., 2019). The samples consist of white, homogeneous pumice clasts collected from the basal Plinian fallout deposits of Layou and Roseau (Figures 1C,D). Whole rock samples were gently crushed to liberate crystals and sieved. Orthopyroxene (751 for Layou, 701 for Roseau) crystals were separated by hand-picking under binocular microscope and mounted in epoxy resin. Orthopyroxenes were oriented during mounting to obtain b-c sections needed for diffusion modeling along the b-axis of crystals.

Analytical Techniques

Textural Investigations

All crystals were investigated using Back-Scattered Electron (BSE) imaging via Scanning Electron Microscopy (SEM; ZEISS-Supra 55, Sorbonne Université, Paris) with an accelerating voltage of 20 kV, a beam current of 8 nA and a counting time in order to reduce noise impact on each image and thus better identify the crystal zonation (scan speed: 7; number of passages per scan: 8; Figure 2). The SEM images were processed using the free software ImageJ ( Selected zones on orthopyroxene crystals have been imaged under high resolution to obtain a high spatial resolution profile (Figures 2B C).


FIGURE 2. From crystal to greyscale profile. (A) High resolution backscattering image [Back-Scattered Electron (BSE) image] of an orthopyroxene oriented with the c-axis in N-S direction. The core and the rim parts of the crystal are identified; (B) Detailed BSE image of the main core-rim boundary; (C) Core-to-rim greyscale profile: averaged greyscale values over the selected area of the boundary (mean width: 20–50 µm) in (B) using the free software ImageJ; (D) Microprobe Mg# (Mg# = Mg/[Mg + Fetot)] core-to-rim traverse over the selected area of the crystal. Spatial resolution is of 2 µm. Uncertainty: 0.2 mol% on Mg; (E) Final calibrated profile and modeled profile resulting from the intercalibration of Back-Scattered Electron greyscale values on microprobe analyses, assuming a step function as the initial condition (dashed lines). Data are presented in Supplementary Tables S2, S3 in supplementary material.

Electron Probe Microanalyses

Major elements (Si, Ti, Al, Mn Fe, Mg, Ca, Na, K, P) in orthopyroxene crystals were analyzed by electron microprobe (Cameca SX-Five and SX-100, Service Camparis, Paris) with an accelerating voltage of 15 kV, a beam current of 10 nA and a focused beam. Counting time on peak and background was set at 10 s for all elements, except for Fe and Mg with 60 s for orthopyroxenes. Zoned orthopyroxene crystals have been analyzed along 2 µm step rim-to-rim or rim-to-core traverses (in regions previously picked out via high-resolution BSE imaging; Figure 2D).

Diffusion Modeling

The method (Figures 2A–E) is based on the modeling as a diffusion profile of the core-to-rim traverse (Figure 2E) obtained by coupling a high analytical resolution profile by Electron Probe Microanalyses (Figure 2D) and a high spatial resolution profile by SEM (Figure 2C), following the methodology developed by Allan et al. (2013) and Couperthwaite et al. (2020).

From each BSE image, a greyscale profile was extracted (Figure 2C), averaging over a selected area. In effect, this averages together many individual pixel traverses across the image, necessary to reduce thermal noise from the backscatter detector (Figure 2B). Profiles have been measured across the c-axis, in the b-crystal direction and away from the terminal faces of crystals, which are more sensitive to recrystallization (Figure 2A). The modeling was accomplished using 1-D modeling, an acceptable procedure as the profiles are very short (∼20–30 µm) compared to the scale of the crystal (∼300 µm), and diffusion is effectively normal to the long face with no additional components. As the diffusion coefficient of Ganguly and Tazzoli (1994) is composition-dependent, modeled profiles develop asymmetry that is more pronounced as the compositional contrast becomes greater. The modeling of such profiles can be eased by exploiting the property that 1-D diffusion profiles proportionally share symmetry, so that if the compositional boundary conditions are the same, the difference between a profile at two different times is merely a stretch factor along the mean diffusion direction. This also holds true in the case of composition-dependent diffusion, where the curve shape is determined by the degree of compositional contrast, but the curve width is dependent on time (Couperthwaite et al., 2020). By using a library of profiles calculated for particular values of compositional contrast, the appropriate curve shape to apply to a natural crystal can be selected from the library based upon compositional data, and then scaled in width to solve for time, leading to very rapid processing. The scaling factor needed is a simple expression of crystal composition, the diffusivity used for the library profiles, the diffusivity in the crystal and the ratio of diffusion timescales in the library and the timescale of the crystal. As the scaling factor can be determined from data, the timescale for the crystal can be simply solved as the only unknown. The internal consistency of timescales obtained by this modeling method is tested by modeling several profiles on the same core-rim boundary of the same crystal and verifying for the acquisition of comparable timescales.

As we undertake modeling without correction for oxygen fugacity, the largest source of uncertainty in the diffusion timescales retrieved by this method is derived from the temperature. Varying the input values of temperature by ±25°C in the Fe-Mg interdiffusion modeling for the same boundary implies change up to one order of magnitude in the calculated timescales (Allan et al., 2013). Here a Monte Carlo simulation was used to estimate the integrated time uncertainty from the different sources of uncertainty: final uncertainties on our reported timescales arise from uncertainties in diffusivity (D) calculations, which derive in turn from 1) uncertainties in measurements of diffusion coefficients (D0, Ea); 2) uncertainty in temperature (the Arrhenius relation gives an inverse exponential dependence on diffusion); 3) considerations about the resolution of greyscale values of the BSE images used for the profile calibration (pixel size, 1σ uncertainty on pixel size, 1σ uncertainty on greyscale values defining a plateau). The resolution of BSE images has usually a small effect on final uncertainties (∼1%; Costa and Morgan, 2010). Uncertainties on temperature, the pre-exponential factor (D0) and the activation energy (Ea) all have a logarithmic effect on final calculated timescales such that the final uncertainties are not Gaussian in time but asymmetric, with larger error bar on the long-time side and smaller error bar on the short-time side, with uncertainties approaching Gaussian behavior in log-time.

Diffusion Parameter

Here we have used the diffusion coefficient resulting from the parameterization of Ganguly and Tazzoli (1994), without correction for oxygen fugacity.

LogD(FeMg)cb= 9.54+2.6  XFe12530T   (1)

where D is the Fe-Mg interdiffusivity, D0 = −9.54 (m2s−1) is the pre-exponential factor, XFe is the molar fraction of the ferrosilite component and T is the temperature in Kelvin. The Ganguly and Tazzoli (1994) equation is formulated for T = 500–800°C, oxygen fugacity of IW to IW + 0.8 (IW = Iron-Wüstite buffer) and XFe = 0.10–0.50. The question of considering or not the oxygen fugacity is subject to controversy (Supplementary Table S1). In their paper, Ganguly and Tazzoli (1994) conjectured that the oxygen fugacity dependence could follow the same form as that for olivine, with an exponent of 1/6, assuming the same mechanism of vacancy generation through iron oxidation, but had no further data upon which to constrain that hypothesis. This parameterization appears to be consistent with the experimental results of Dohmen et al. (2016), who argue for a minor dependence of the diffusion coefficient on oxygen fugacity (n = 1/20 instead of 1/6), although that study examined pyroxene very close to enstatite, and which is not a good match for these samples. A further study of Klugel (2001) in iron rich orthopyroxenes (in the temperature range 1,050–1,200°C and at QFM [quartz-fayalite-magnetite] oxygen buffer) reported a diffusivity of DFe-Mg = 3 × 10−19 m2s−1 at T = 1,130°C. If extrapolated to higher temperatures, the results of Ganguly and Tazzoli (1994) are in good agreement with those of Klugel (2001).

In summary, all modeled profiles and timescales have been here calculated for a fixed temperature of 850°C (±25°C) based upon the experimental work of Solaro et al. (2019) on these magma compositions. This yields diffusivities of the order of D ∼10–20 m2 s−1 (Ganguly and Tazzoli, 1994), without dependence of the diffusivity on oxygen fugacity. The initial condition is considered to be a step-wise profile (i.e., instantaneous growth at a fixed composition followed by diffusional relaxation of the interface), which yields maximum timescales if all broadening occurs via diffusion.

Compositional profiles for these eruptions have been interpreted as pure diffusion profiles with no or limited recorded growth, i.e., considering profiles as due purely to diffusion. This constrains maximum timescales, because if the profiles have a growth component, the required timescale is shortened. The fitted profiles have features typical of diffusive character—smooth compositional variation of sigmoid shape (Supplementary Figure 2). Oscillations can be seen toward the rim but are minor in amplitude (±0.2 mol% or less) and could plausibly be attributed to noise. Typically, growth-dominant profiles would display linear trends or curves on graphs of composition vs. distance, which would not be congruent nor consistent with a diffusion sigmoid calculated for the given initial compositional contrast.

Crystal System Analysis Approach

The CSA technique has been successfully applied to the study of Etna, in order to define the main active magmatic environments (ME) prior eruption from crystal compositional zoning, their connections (by heating, decompression, mixing, mingling, crystal settling, etc) and the pre-eruptive processes occurring at depth (Kahl et al., 2011; Kahl et al., 2015). Following this definition, a ME is not necessarily a single batch of magma, but may be different part of the same reservoir at depth. In the case of Etna, a variety of ME were described based on olivine-melt equilibria and their occurrence with time. We applied this method to orthopyroxenes (opx) that show variations in major-element composition (Enstatite content, XEn) that can be related to changes in the magmatic system. These changes can lead to either reverse zoning, typically associated with higher temperatures, higher pressures, and higher water fugacity, or normal zoning, typically associated with cooler and shallower environments.

The compositions of unzoned, single zoned and multiple zoned opx along rim-to-rim or core-to-rim traverses allow us to identify compositional plateaus (En mol% plateaus) (Supplementary Table S3). Compositional plateaus have been here defined as regions of the profiles with a minimum length of ∼15 µm (respect to 315–500 µm wide crystals) within which the En content of the crystal shows a constant value, with accepted variability of ±0.5 mol%. The composition of each plateau is then calculated as an average of all values defining the plateau. The occurrence of plateaus can be interpreted as a result of crystal growth in different, stable ME, with residence times in each environment long enough to be recorded by crystal growth, but also short enough that the compositional gradient has not been erased by diffusional equilibration. The repetitive occurrence of plateaus indicates that compositional zoning observed in a given crystal has not been produced by pure fractionation during growth, but by subsequent and fast changes during crystal growth in response to a change from one ME to another, i.e., by changing pressure, temperature, oxygen fugacity or water fugacity, or physically transferring a crystal to a different melt. In this way, assuming that crystal growth occurred sequentially from core to rim, cores represent the first ME recorded by crystals while outer rims represent the last ME. Analysis of compositional plateaus from core to rim allows us to track either movements of crystals and their carrying melts between different ME or modification of storage conditions by external processes (heating, decompression event, mixing, etc.) through the plumbing system until eruption. This procedure has to be repeated on all analyzed crystals in order to define crystal populations, proportions, and the frequently-travelled crystal transfer pathways, before drawing a CSA analysis diagram. This method allows a simple graphical representation that avoids much of the difficulties of combining and interpreting diverse and complicated zoning patterns between crystals of the same sample.

Here we applied this approach to opx to highlight magma transfers and the magnitudes of these transfers at depth between the different ME. This approach is based on the recognition of the diversity of opx composition interpreted as ME changes. Eventually, we can identify the disrupting events at depth responsible for magma remobilization up to eruption. That is also why a ME does not necessarily correspond to a physical lens of magma within the crust; this will be discussed hereafter.


The Diversity of Orthopyroxene Texture and Composition

A total of 1,452 opx were analyzed for the two eruptions (701 for Roseau, 751 for Layou), including their textural features and chemical compositions (Figure 3; Table 1; Supplementary Tables S2, S3).


FIGURE 3. Overview of opx textures and composition. Selected SEM images illustrating the opx zonations for Roseau (A: normal; B: reverse; C: multiple, reverse+normal) and Layou (D: normal; E: reverse; F: multiple, normal+reverse). Frequency histograms of the zoned and unzoned orthopyroxene compositions in Roseau (left side) and Layou (right side) expressed in terms of En (mol%). (G,J) unzoned crystals in grey; (H,K) cores of zoned crystals in blue for Layou and red for Roseau; (I, L) rims of zoned crystals, in dashed blue for Layou and dashed red for Roseau. Unzoned orthopyroxene compositions: single point analyses in the crystal core. Core and rim composition: averaged plateau compositions of the zoned crystals, with a minimum length of ∼15 μm. N: number of crystals.


TABLE 1. Statistic repartition of zonations in opx from Layou and Roseau eruption. Opx are separated in the 315–500 µm fraction.

Textural analysis by high-resolution BSE images shows that, for both eruptions, 85% of crystals of the 315–500 µm fraction are unzoned (Figures 3A–F; Table 1). Opx belonging to larger size classes (500 μm–1 mm) usually show only zonation linked to the melt evolution by fractional crystallization while smaller size classes (250–315 µm) present clear zonations, in the same proportions as the 315–500 µm fraction. Roseau opx show significantly more simple reverse zonation (higher En% in the rim than in the core; 69%—Figures 3B,C,H,I) compared to normal zoning (lower En% in the rim than in the core; 8%—Figure 3A). Some multiple zonations are present (23%) with a prevalence of the “reverse + normal” zonation (Figure 3I; Supplementary Figure S5). Layou opx show more scattered proportions of zonation types than observed in Roseau, despite the similar proportion of zoned crystals (15%): reverse zonations (44%—Figures 3E,F,K,L) are less abundant and normal zoning (28%—Figure 3D) more prevalent in Layou than for Roseau. Of the zoned Layou crystals, 28% show multiple zonation, but with a predominance of “normal + reverse” zonation (Figure 3L), contrary to Roseau.

Core compositions, as well as the unzoned crystals, help us to constrain the equilibrium state of the magma before the disrupting event responsible for magma remobilization; core compositions of each opx are presented in Figure 3. The unzoned opx from Roseau have a range of composition between En50 and En54 with a clear mode at En51–53 (Figure 3A), whereas unzoned opx of Layou encompass a compositional range between En48 and En53, with a moderate mode at En51–53 (Figure 3D). For zoned opx, core compositions are centered at En51–53 for 80% of Roseau zoned crystals (Figure 3B) and for about 60% of Layou zoned crystals (Figure 3E), but with a second small peak at En49–50 (25%). Thus, unzoned opx and the core of zoned opx describe the same dominant composition domain of En51–53 for both eruptions (Figures 3A,B,D,E). Thus, the unzoned populations are a good proxy for the cores of the zoned populations for both eruptions. The cores have a common origin, and the rims represent dynamic events that happen subsequently, as detailed hereafter.

Orthopyroxene-Hosted Melt Inclusions

Measuring the composition and volatile content in MIs trapped in phenocrysts is one straightforward manner to infer magma storage depth when associated to matrix glass investigation. In addition to the MI study of Balcone-Boissard et al. (2018), a dedicated investigation of a possible correlation between the MI composition and the host opx has been undertaken on 15 MI-opx pairs (Figure 4). Correction of post-entrapment crystallization was considered on the basis of the KD (Fe-Mg)opx-melt (Putirka, 2008) and only represented between 1 and 2 vol% of the MI volume (Supplementary Table S5). For both eruptions, the measured MIs are trapped within opx that encompass the whole compositional range, En48–57 (Figure 4). Reporting the data together with existing data in the same samples (Balcone-Boissard et al., 2018) shows that there is no systematic variation of MI composition with host-opx core and rim composition (Figure 4).


FIGURE 4. Selected melt inclusion compositions as a function of host-orthopyroxene composition (En mol%). (A) SiO2; (B) CaO; (C) FeO; (D) Na2O.

Fe-Mg Diffusion Modeling in Orthopyroxene

Core-to-rim composition profiles were obtained and modeled for 44 selected zoned opx (21 for Roseau and 24 for Layou). Layou crystals included both normal- and reverse-zoned examples while only reverse-zoned crystals were studied for Roseau (given the dominance of reversely-zoned crystals in this eruption). Examples of data interpreted as diffusion profiles for Layou and Roseau are shown in Figure 5. As the purpose here is to constrain the last pre-eruptive processes, only the last change in composition - the outermost zone, excluding quench - has been used for timescale determination in each case.


FIGURE 5. Examples of two modeled diffusion profiles for Roseau (A) and Layou (B). Red/blue line on each Back-Scattered Electron (BSE) image (left) indicates the profile position [Electron Probe Microanalyses (EPMA) full profile] on the respective orthopyroxene (Roseau, Layou). Diffusion profiles (right) are reported as Mg# [Mg/(Mg + Fetot)], as a function of the distance from the center of the profile (in µm). Red rhombs: EPMA and BSE greyscale intercalibrated data; black line: diffusion model. Time (t) is given for each profile. The initial conditions are given by the dashed line (staircase profile).

Results obtained on simple zoned crystals (reverse and normal, for Layou and reverse only for Roseau) show a significant peak of Fe-Mg interdiffusion timescales at 3–10 years before eruption (Figures 6A,C; Table 2). For Roseau, the 17 crystals representing the histogram peak and showing timescales under 10 years (Figures 6A,B) allow us to calculate a timescale of ∼4–5 years prior to eruption, with an absolute uncertainty of −3/+ 11 years (Figure 6B). For Layou, the 14 crystals representing the histogram peak for timescales under 10 years (Figures 6C,D) lead to a calculated timescale of ∼3–4 years prior to eruption, with an absolute uncertainty of -2/+9 years (Figure 6D).


FIGURE 6. Timescale estimate from Fe-Mg interdiffusion modelling in orthopyroxenes. Frequency histogram of Fe-Mg interdiffusion timescales modelled at 850°C for the main core-rim boundary of (A) Roseau opx and (C) Layou opx; (B) and (D) individual ages and their absolute uncertainties based on the propagation of a temperature uncertainty of ±25°C for (B) Roseau and (D) Layou. The inset graphs inside (B) and (D) represent selected individual ages from the peak in interdiffusion ages.


TABLE 2. Timescales estimates for Layou and Roseau eruption. Sample labeling refers to the stratigraphy of Boudon et al. (2017) with DOM41a3 for Layou and DOM60d1a for Roseau. Timescales are given in year with the detail in uncertainty (σ+ and σ−). The type of zonation (simple or multiple) and the position within the multiple zoned opx is given. The corresponding plateaus in En content and the ME changes are detailed.


Pre-Eruptive Magma Dynamics and Identification of the Disrupting Event at Depth

Our starting conception of the architecture of the magmatic system beneath Dominica is a transcrustal magmatic system, comprising a series of magma bodies stored at different depths within the crust that may be interconnected by dykes and influenced by volatile fluxing (Bachmann and Bergantz, 2003; Cashman et al., 2017) or various regions of a large zoned magma body with different conditions (partially molten rocks, mush, cold storage zone; Edmonds et al., 2019). This creates a variety of ME characterized by differences in pressure, temperature, oxygen fugacity and volatile saturation. For both eruptions, the stable pre-eruptive conditions of the main magma body were constrained by Solaro et al. (2019) to be similar for both eruptions: 850 ± 25°C, ∼400 MPa (depth of ∼16 km, considering a local rock density of 2.45 g/cm3), a melt H2O content of ∼7–8 wt% (i.e., slightly H2O-undersaturated conditions) and an oxygen fugacity close to NNO +0.6. Such stable conditions of the main magma body are highlighted by a population of ∼85% of unzoned opx crystals in the 315–500 μm range. Here we discuss for both eruptions what the other ∼15% of zoned opx can tell us about the last magmatic events that triggered the eruption.

These ∼15% of zoned crystals display similar core compositions [in agreement with the stable pre-eruptive conditions determined by Solaro et al. (2019)] but a variety of rim compositions, and are thus evidence of transfers from a common initial environment to other ME, where compositional rims developed. For Roseau opx, the prevalence of Fe-Mg reverse zoning (69% of singly-zoned crystals) suggests few perturbations of the stored magma. For Layou, the complex zonation pattern (∼44% of reversely zoned opx, 28% of normal zoning, and 28% multiple zoning) comparatively testifies a heterogeneous perturbation of the main magma body during the pre-eruptive period. In the following section, these zoning patterns are interpreted in terms of changes in the physical parameters of the storage region, in order to trace the crystal history and the dynamics of the storage system. We then investigate the origin of the last event responsible for magma remobilization prior to eruption. The En content of opx is above all dependent on temperature: a temperature increase causes an increase in En (e.g., Martel et al., 1999; Solaro et al., 2019). Conversely, the En content of opx is less dependent on pressure. Thus, an increase or a decrease in En content may in some cases be interpreted as a heating or a cooling event, respectively, when magmas of different compositions or temperatures are mixing.

Magma Mixing and Reheating

For both eruptions, the relative homogeneity in the compositions of the melt inclusions trapped in opx showing large XEn variations from En47 to En57 (Figure 4) suggests a rather homogeneous and buffered silicate melt composition of low temperature (based on Zr content data as evidenced by Balcone-Boissard et al., 2018). In both eruptions, matrix glasses are rhyolitic but slightly less differentiated (silica- and alkali-poorer) than MIs (Balcone-Boissard et al., 2018). Magma mixing processes are able to create XEn increase, either by mixing a dacitic magma with a mafic recharge that produces a less differentiated matrix melt or by mixing different parts of the same mush displaying different degrees of melt differentiation. The hypothesis of a mixing with a more mafic magma has been discarded on the bases of i) the melt inclusion compositions, which reflect a rather homogeneous differentiated melt composition trapped within opx, whatever their En content, and also on the basis of ii) the rhyolitic matrix glasses that do not show textural evidences of mixing. The hypothesis of mixing different parts of the same mush is the most realistic, as i) the matrix glasses are slightly less differentiated than MIs, which is incompatible with a pure evolution by decompression-induced crystallization of the MI trapped in the opx and ii) MIs reflect an evolved melt formed at lower temperature than the matrix glass, as evidenced by the Zr saturation in MIs but not in the matrix glasses (Balcone-Boissard et al., 2018).

Investigating the nature of the zoning changes gives insights into the pre-eruptive magma dynamics, as displayed in Figure 7 and detailed hereafter. In agreement with the CSA method, the fluctuations in XEn within each single opx crystal have been divided into compositional groups in order to evaluate the magnitude of the heating or cooling event recorded by the crystals. The Roseau eruption shows 2 ME: MER1 (En50-54) and MER2 (En54-59), with the population dominantly displaying MER1 compositions. The dominant zoning is reverse and may be attributed to an increase of temperature. Though MER2 -MER1 normal zoning is present, it represents a statistically minor component of the population and is typically followed by reverse zoning back into MER2. This suggests a dominant reheating process or magma transfer from MER2 to MER1 with time. This reheating may be linked to the injection of a hotter and slightly less differentiated silicic magma into a rather homogeneous cooler magma reservoir (Figures 8A,B).


FIGURE 7. A variant Crystal Systems Analysis (CSA) diagram dedicated to opx, after the methodology of Kahl et al. (2015), demonstrating the direction and nature of sequential crystal zones of the opx from the Roseau (A) and Layou (B) eruptions. Each crystal is represented by a single arrow. Black arrows denote simple zoning and red arrows denote multiple zoning, from core to rim. The size of the magmatic environment boxes (ME1 to ME3, with subscript L for Layou and R for Roseau) are proportional to the number of unzoned crystals (grey) of the given composition (XEn) plus the number of zoned crystals (white) with cores of the given composition.


FIGURE 8. A dynamic magma model within the crust. Inset: Localization within the crust of the mush zone studied here. The color of the crystal is linked to magma temperature, from “cold” environment in blue to “hot” one in red, with intermediate temperature in orange; relative temperature estimates are on the basis of opx composition. Initial stage for (A) Roseau and (C) Layou. The dynamic model proposed on the basis of opx for (B) Roseau and (D) Layou.

For the Layou eruption, 3 ME are described (by order of abundance): MEL1 (En51-53), MEL2 (En47-51), MEL3 (En53-58). The volumetrically dominant body of MEL1 (En51-54) interacts by exchanging with a cooler, slightly smaller body of MEL2 (En47-51). This could imply alternating heating and cooling of a magma body, or exchange across an interface in a thermally stratified body, as proposed in Figure 8C. There is a significant interaction between the coldest MEL2 (En47-51) and the hottest MEL3 (En54-59), suggesting that the coolest part, likely the walls of the magma reservoir, was interacting with the hottest magma prior to eruption. As for Roseau, the hottest magma may be the injection of magma batches initially stored in an independent reservoir, approximatively at the same pressure/depth. Interestingly, there is only one direct transfer between MEL1 and MEL3, compared to significant transfers from MEL1 to MEL3 via MEL2. This suggests that MEL1 and MEL3 may have been more isolated in space and/or in time (Figure 8D). Overall, the CSA depicts a stratified or cooling system, with limited exchanges, prior to a significant disruption of the cooler endmember by a slightly hotter magma, likely occurring within a decade before eruption (Figures 6, 8; Table 2). In the CSA diagram, the normal zonation linked to the transfer from ME1L to ME2L, is related to the convection within the main reservoir of the crystals first stored in the hottest part and then in contact with the coolest part of the reservoir (Figures 8C,D). As both reverse and normal zonation have a timescale of the same order of magnitude (less than 10 years; Table 2), this would imply that the timescales calculated for both processes are too short for accurate discrimination by the method or that both processes are related. The injection of a new batch of hot magma may be responsible for the magma convection within the intruded reservoir.

The complexity of the pre-eruptive magma dynamics for the Layou eruption confirms the hypothesis of a system heterogeneous in temperature, with parts of the magma stored at lower temperature (allowing quartz crystallization; Solaro et al., 2019). In contrast, the CSA results suggest more direct paths and fewer magma transfers in the case of Roseau magmatic system.

Comparison With Previous Studies

Various studies have already used the Fe-Mg interdiffusion chronometer in opx to retrieve timescales of magma dynamics prior to eruption (Table 3). These include examples of eruptions across a large range of differentiated magma volumes.


TABLE 3. Comparison of diffusion timescales obtained for different large volume silicic eruptions.

Recent petrological studies have shown that the large magma chambers that feed Plinian eruptions at Santorini (Greece; Flaherty et al., 2018) are probably transient features on geological timescales. Magmas are now widely believed to be generated and incubated in bodies of crystal mush, from which they are extracted and accumulated in melt-dominated magma chambers prior to eruption (Bachmann and Huber, 2016). Conversely to what is generally expected for phenocrysts in magmas, at Santorini, silicic magmas were not the result of a slow cooling over geologically long time periods. By modeling, as here, the diffusive re-equilibration of chemical gradients in zoned opx (feldspars and cpx too), it has been shown that the phenocrysts grew within only a few hundred years, and possibly only a few decades, for the eruption of the Late Bronze Age eruption of Santorini (Flaherty et al., 2018).

The Fe-Mg interdiffusion in opx of the Oruanui eruption (Taupo volcano, New Zealand), involving a magma volume of 530 km3, helped to constrain the timing of the decompression event and extraction of melt occurring in the centuries before the eruption (Allan et al., 2013). For the same eruption the whole history of generation, accumulation and extraction of the Oruanui magma body has been discussed by combining diffusion modeling in opx, magnetites and olivines (Allan et al., 2017). The timescales of priming and accumulation of magma prior to the three youngest eruptions (<2.15 ka) of the Taupo volcanic center indicate that a rapid heating and priming of the crystal-mush, from which the final melt-dominant magma body is extracted, occurred in <120 years prior to eruption (Barker et al., 2016). For the Whakamaru eruption (∼340 ka, >1,000 km3 of magma), the diffusion modeling of Ti in quartz shows that timescales of the final recharge rejuvenating the mush reservoir are in the range of decades before eruption (Matthews et al., 2012).

By studying the geochemical signature of zircon and sanidines, the timescales of magma storage and differentiation of voluminous high-silica rhyolites at Yellowstone caldera (Wyoming) have been investigated (Stelten et al., 2015). Though this is not the same chronometry, the authors argue that eruptible rhyolites are generated by extracting melt and zircons from a long-lived mush of immobile crystal-rich magma. Residence times of the large-volume rhyolites (∼40–70 km3) are ≤1 kyr (conservatively <6 kyr), which suggests that large volumes of rhyolite can be generated rapidly by extracting melt from a crystal mush. As for Santorini, the Yellowstone magmatic system evidences that magma reservoirs spend most of their time in a largely crystalline uneruptible state and that eruptible magma bodies are ephemeral features.

These examples of Santorini, Taupo or Yellowstone, as Dominica, highlight that volcanic systems are characterized by thousands of years of volcanic repose: magma reservoirs spend most of their time in a largely crystalline uneruptible state. Consequently, eruptible magma bodies are ephemeral features, which do not necessarily require multi-decadal timescales for becoming eruptive by magma assembly and mush rejuvenation. Instead, large magma volumes can be mobilized to eruption over relatively short timescales, with no evident correlation between magma volume and timescales. Such highly dynamic volcanic systems bring important implications for ongoing monitoring and hazard mitigation, particularly in areas of high population density (Sparks and Cashman, 2017; Edmonds et al., 2019).


The coupled investigation of crystal and melt composition lead us to identify the disrupting event at depth responsible for the magma mobilization, leading to two major explosive eruptions in Dominica. Crystal zonation, interpreted using the revised CSA approach, coupled with diffusion modeling, provides a detailed view of the pre-eruptive architecture and dynamics of the crystal mush prior to eruption. This highlights the diverse complexity of the magmatic system with several magmatic environments and short remobilization timescales. In the last 10 kyrs, several Plinian eruptions in the central part of Dominica have been described and should be investigated to detail the plumbing system dynamics in time and space.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

CS acquired all data during her PhD, wrote the first draft. HBB and GB supervised the PhD, wrote the submitted version with new figures. DM helped during data acquisition and wrote the submitted version. CM wrote the submitted version. LO helped for data and figure management before submission.


This work is part of CS (financial support of the French MNRT doctoral Grant). We would like to thank the PPR fund of IPGP. This study was partly founded by the program Tellus-Aleas of INSU (CNRS) and ANR V-CARE-18-CE03-0010.

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.


This work is part of the PhD thesis of Solaro-Müller (2017). DM would like to thank IPGP for the award of a visiting lectureship which enabled this work to take place. Michel Fialin, Nicolas Rividi and Omar Boudouma are thanked for assistance during EPMA and SEM analyses. We would like to thank the editor V. Acocella, the asociated editor J. Lowenstern and the two reviewers, G. Kilgour and D. Ruth, for their constructive comments.

Supplementary Material

The Supplementary Material for this article can be found online at:


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

CrossRef Full Text | Google Scholar

Allan, A. S. R., Barker, S. J., Millet, M-A., Morgan, D. J., Rooyakkers, S. M., Schipper, C. I., et al. (2017). A cascade of magmatic events during the assembly and eruption of a super-sized magma body. Contrib. Mineral. Petrol. 172, 49. doi:10.1007/s00410-017-1367-8

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Annen, C. (2009). From plutons to magma chambers: thermal constraints on the accumulation of eruptible silicic magma in the upper crust. Earth Planet Sci. Lett. 284 (3–4), 409–416. doi:10.1016/j.epsl.2009.05.006

CrossRef Full Text | Google Scholar

Bachmann, O., and Bergantz, G. W. (2003). Rejuvenation of the fish canyon magma body: a window into the evolution of large-volume silicic magma systems. Geology 31 (9), 789–792. doi:10.1130/G19764.1

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Balcone-Boissard, H., Boudon, G., Blundy, J. D., Martel, C., Brooker, R. A., Deloule, E., et al. (2018). Deep pre-eruptive storage of silicic magmas feeding plinian and dome-forming eruptions of central and northern Dominica (Lesser Antilles) inferred from volatile contents of melt inclusions. Contrib. Mineral. Petrol. 173 (12), 101. doi:10.1007/s00410-018-1528-4

CrossRef Full Text | Google Scholar

Barker, S. J., Wilson, C. J., 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

CrossRef Full Text | Google Scholar

Blundy, J. D., and Annen, C. J. (2016). Crustal magmatic systems from the perspective of heat transfer. Elements 12, 115–120. doi:10.2113/gselements.12.2.115

CrossRef Full Text | Google Scholar

Boudon, G. (1983). La montagne Pelée, Martinique: évolution volcanologique. Mém. Soc. Géol. 163, 231–238.

Google Scholar

Boudon, G., Balcone-Boissard, H., Solaro, C., and Martel, C. (2017). Revised chronostratigraphy of recurrent ignimbritic eruptions in Dominica (Lesser Antilles Arc): implications on the behavior of the magma plumbing system. J. Volcanol. Geoth. Res. 343, 135–154. doi:10.1016/j.jvolgeores.2017.06.022

CrossRef Full Text | Google Scholar

Boudon, G., Balcone-Boissard, H., Villemant, B., and Morgan, D. J. (2015). What factors control superficial lava dome explosivity? Sci. Rep. 5, 14551. doi:10.1038/srep14551

CrossRef Full Text | Google Scholar

Bouysse, P., and Westercamp, D. (1990). Subduction of Atlantic aseismic ridges and late cenozoic evolution of the lesser Antilles island arc. Tectonophysics 175 (4), 349357–355380. doi:10.1016/0040-1951(90)90180-G

CrossRef Full Text | Google Scholar

Carey, S. N., and Sigurdsson, H. (1980). The Roseau ash: Deep-sea tephra deposits from a major eruption on Dominica, Lesser Antilles arc. J. Volcanology Geothermal Res. 7, 67–86. doi:10.1016/0377-0273(80)90020-7

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

CrossRef Full Text | Google Scholar

Chakraborty, S. (2008). Diffusion in solid silicates: a tool to track timescales of processes comes of age. Annu. Rev. Earth Planet Sci. 36, 153–190. doi:10.1146/

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Charlier, B. L. A., Morgan, D. J., Wilson, C. J. N., Wooden, J. L., Allan, A. S. R., and Baker, J. A. (2012). Lithium concentration gradients in feldspar and quartz record the final minutes of magma ascent in an explosive supereruption. Earth Planet Sci. Lett. 319, 218–227. doi:10.1016/j.epsl.2011.12.016

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Costa, F., Chakraborty, S., and Dohmen, R. (2003). Diffusion coupling between major and trace elements and a model for the calculation of magma chamber residence times using plagioclase. Geochem. Cosmochem. Acta 67 (12), 2189–2200. doi:10.1016/S0016-7037(02)01345-5

CrossRef Full Text | Google Scholar

Costa, F., and Morgan, D. (2010). “Time constraints from chemical equilibration in magmatic crystals,” in Timescales of magmatic processes: from core to atmosphere. Editors A. Dosseto, S. P. Turner, and J. A. Van Orman (Chichester, England: John Wiley & Sons), 125–159.

Google Scholar

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

CrossRef Full Text | Google Scholar

Couch, S., Sparks, R. S. J., and Carroll, M. R. (2001). Mineral disequilibrium in lavas explained by convective self-mixing in open magma chambers. Nature 411 (6841), 1037–1039. doi:10.1038/35082540

CrossRef Full Text | Google Scholar

Couperthwaite, F. K., Thordarson, T., Morgan, D. J., Harvey, J., and Wilson, M. (2020). Diffusion timescales of magmatic processes in the Moinui lava eruption at Mauna Loa, Hawai’I, as inferred from bimodal olivine populations. J. Petrol. 58. 1–19. doi:10.1093/petrology/egaa058

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Edmonds, M., Cashman, K. V., Holness, M., and Jackson, M. (2019). Architecture and dynamics of magma reservoirs. Philos. Trans. Royal Soc. A Math. Phys. Eng. Sci. 377 (2139), 20180298. doi:10.1098/rsta.2018.0298

CrossRef Full Text | Google Scholar

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. Mineral. Petrol. 173, 31. doi:10.1007/s00410-018-1490-1

CrossRef Full Text | Google Scholar

Ganguly, J., and Tazzoli, V. (1994). Fe 2+–Mg interdiffusion in orthopyroxene. Am. Mineral. 79 (9–10), 930–937. doi:10.2138/am-2016-5815

CrossRef Full Text | Google Scholar

Ginibre, C., Wörner, G., and Kronz, A. (2004). Structure and dynamics of the Laacher See magma chamber (Eifel, Germany) from major and trace element zoning in sanidine: a cathodoluminescence and electron microprobe study. J. Petrol. 45 (11), 2197–2223. doi:10.1093/petrology/egh053

CrossRef Full Text | Google Scholar

Howe, T. M., Lindsay, J. M., Shane, P., Schmitt, A. K., and Stockli, D. F. (2014). Re‐evaluation of the Roseau Tuff eruptive sequence and other ignimbrites in Dominica, Lesser Antilles. J. Quat. Sci. 29 (6), 531–546. doi:10.1002/jqs.2723

CrossRef Full Text | Google Scholar

Huber, C., Bachmann, O., and Dufek, J. (2012). Crystal-poor versus crystal-rich ignimbrites: a competition between stirring and reactivation. Geology 40 (2), 115–118. doi:10.1130/G32425.1

CrossRef Full Text | Google Scholar

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), 38–47. doi:10.1016/j.epsl.2009.03.029

CrossRef Full Text | Google Scholar

Huber, C., Bachmann, O., and Manga, M. (2010). Two competing effects of volatiles on heat transfer in crystal-rich magmas: thermal insulation vs defrosting. J. Petrol. 51 (4), 847–867. doi:10.1093/petrology/egq003

CrossRef Full Text | Google Scholar

Kahl, M., Chakraborty, S., Costa, F., and Pompilio, M. (2011). Dynamic plumbing system beneath volcanoes revealed by kinetic modeling, and the connection to monitoring data: an example from Mt. Etna. Earth Planet Sci. Lett. 308 (1–2), 11–22. doi:10.1016/j.epsl.2011.05.008

CrossRef Full Text | Google Scholar

Kahl, M., Chakraborty, S., Pompilio, M., and Costa, F. (2015). Constraints on the nature and evolution of the magma plumbing system of Mt. Etna Volcano (1991–2008) from a combined thermodynamic and kinetic modelling of the compositional record of minerals. J. Petrol. 56 (10), 2025–2068. doi:10.1093/petrology/egv063

CrossRef Full Text | Google Scholar

Klügel, A. (2001). Prolonged reactions between harzburgite xenoliths and silica-undersaturated melt: implications for dissolution and Fe-Mg interdiffusion rates of orthopyroxene. Contrib. Mineral. Petrol. 141 (1), 1–14. doi:10.1007/s004100000222

CrossRef Full Text | Google Scholar

Lindsay, J. M., Trumbull, R. B., and Siebel, W. (2005). Geochemistry and petrogenesis of late Pleistocene to recent volcanism in southern Dominica, lesser Antilles. J. Volcanol. Geotherm. Res. 148, 253–294. doi:10.1016/j.jvolgeores.2005.04.018

CrossRef Full Text | Google Scholar

Marsh, B. D. (1989). Magma chambers. Annu. Rev. Earth Planet Sci. 17, 439–474. doi:10.1146/annurev.ea.17.050189.002255

CrossRef Full Text | Google Scholar

Marsh, B. D. (2006). Dynamics of magmatic systems. Elements 2 (5), 287–292. doi:10.2113/gselements.2.5.287

CrossRef Full Text | Google Scholar

Martel, C., Pichavant, M., Holtz, F., Scaillet, B., Bourdier, J.-L., and Traineau, H. (1999). Effect of fO2 and H2O on andesite phase relations between 2 and 4 kbars. J. Geophys. Res. 104, 29453–29470. doi:10.1029/1999JB900191

CrossRef Full Text | Google Scholar

Matthews, N. E., Huber, C., Pyle, D. M., and Smith, V. C. (2012). Timescales of magma recharge and reactivation of large silicic systems from Ti diffusion in quartz. J. Petrol. 53 (7), 1385–1416. doi:10.1093/petrology/egs020

CrossRef Full Text | Google Scholar

Morgan, D. J., Blake, S., Rogers, N. W., 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), 933–946. doi:10.1016/j.epsl.2004.03.030

CrossRef Full Text | Google Scholar

Palm, W. J. (2005). System dynamics. New York, NY: McGraw-Hill Higher Education, 834.

Google Scholar

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

CrossRef Full Text | Google Scholar

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. Geotherm. Res. 165 (3), 142–162. doi:10.1016/j.jvolgeores.2007.06.002

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Saunders, K. E., Morgan, D. J., Baker, J. A., and Wysoczanski, R. J. (2010). The magmatic evolution of the Whakamaru supereruption, New Zealand, constrained by a microanalytical study of plagioclase and quartz. J. Petrol. 51, 1644. doi:10.1093/petrology/egq064

CrossRef Full Text | Google Scholar

Sigurdsson, (1972). Partly-welded pyroclast flow deposits in Dominica, Lesser Antilles. Bull. Volcanol. 36, 148–163.

Shcherbakov, V. D., Plechov, P. Y., Izbekov, P. E., and Shipman, J. S. (2011). Plagioclase zoning as an indicator of magma processes at Bezymianny Volcano, Kamchatka. Contrib. Mineral. Petrol. 162 (1), 83–99. doi:10.1007/s00410-010-0584-1

CrossRef Full Text | Google Scholar

Smith, A. L., Roobol, M. J., Mattioli, G. S., Fryxell, J. E., Daly, G. E., and Fernandez, L. A. (2013). The volcanic geology of the mid-arc Island of Dominica. Geol. Soc. Am. 496, 249. doi:10.1130/2013.2496

CrossRef Full Text | Google Scholar

Solaro, C., Martel, C., Champallier, R., Boudon, G., Balcone-Boissard, H., and Pichavant, M. (2019). Petrological and experimental constraints on magma storage for large pumiceous eruptions in Dominica island (Lesser Antilles). Bull. Volcanol. 81, 55. doi:10.1007/s00445-019-1313-x

CrossRef Full Text | Google Scholar

Solaro-Müller, C. (2017). Storage conditions and dynamics of magma reservoirs feeding the major pumiceous eruptions of Dominica (Lesser Antilles arc). PhD thesis. Paris (France): University Paris Diderot (Sorbonne Paris Cité), 330.

Google Scholar

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

CrossRef Full Text | Google Scholar

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 (8), 1607–1642. doi:10.1093/petrology/egv047.

Wadge, G. (1984). Comparison of volcanic production rates and subduction rates in the lesser Antilles and Central America. Geology 12 (9), 555–558. doi:10.1130/0091-7613(1984)12<555:COVPRA>2.0.CO;2

CrossRef Full Text | Google Scholar

Wadge, G., Voight, B., Sparks, R. S. J., Cole, P. D., Loughlin, S. C., and Robertson, R. E. A. (2014). An overview of the eruption of Soufrière Hills volcano, Montserrat, from 2000 to 2010. Geol. Soc. Lond. Mem. 39, 1–40. doi:10.1144/M39.1

CrossRef Full Text | Google Scholar

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 (8), 1413–1431. doi:10.1093/petrology/44.8.1413

CrossRef Full Text | Google Scholar

Keywords: orthopyroxenes, Fe-Mg interdiffusion, timescale, plumbing system, crystal system analysis, Dominica (Lesser Antilles), ignimbritic eruptions

Citation: Solaro C, Balcone-Boissard H, Morgan DJ, 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:574032. doi: 10.3389/feart.2020.574032

Received: 18 June 2020; Accepted: 14 October 2020;
Published: 19 November 2020.

Edited by:

Jacob B. Lowenstern, Volcano Disaster Assistance Program (USGS), United States

Reviewed by:

Dawn Ruth, United States Geological Survey, United States
Geoff Kilgour, GNS Science, New Zealand

Copyright © 2020 Solaro, Balcone-Boissard, Morgan, Boudon, Martel and Ostorero. 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.

*Correspondence: Hélène Balcone-Boissard,