Melt embayments record multi-stage magma decompression histories

Behnaz Hosseini *, Madison Myers
Department of Earth Sciences, Montana State University, Bozeman, MT 59717-3480, USA

A R T I C L E I N F O

Keywords:
Magma decompression
Melt embayments
Volatiles
Numerical modeling
Decompression experiments

A B S T R A C T

Over the last decade, the melt embayment has proven its merit as a robust petrological tool capable of recording
magma decompression rates for explosive eruptions. However, the models developed and applied to extract this
information from embayments have not accounted for the complexity and nonlinearity of magma flow in the
conduit. We present Embayment Decompression in Two Stages (EDiTS): a numerical model for extracting magma
decompression rates from measured volatile diffusion profiles preserved in crystal-hosted embayments,
approximating magma acceleration using two constant-rate decompression paths. This model solves for three
unknown parameters: initial (deeper) and final (shallower) decompression rates, as well as the pressure where a
transition occurs. We successfully benchmark EDiTS against existing numerical diffusion models, and use
controlled multi-stage decompression experiments on natural quartz-hosted embayments to test the ability of our
model to recover known decompression paths. We find that EDiTS is able to closely approximate the known two-
stage path in the mixed-volatile (H2O + CO2) experiment, while a constant-rate modeling approach is unable to
simultaneously fit H2O and CO2 gradients. However, in the H2O-saturated experiment, there is no unique so-
lution to the resulting gradient, with both constant-rate and two-stage models reproducing the measured profile,
and EDiTS notably overestimating the known total ascent time by several hours. Using decompression experi-
ments, we show that constant-rate models can provide misleadingly good fits to embayment H2O gradients
produced by more complex decompression histories, and thus the measurement and modeling of multiple
diffusing species, when available, can provide crucial constraints. We then apply EDiTS to re-evaluate mixed-
volatile embayment datasets from explosive silicic arc and caldera-forming eruptions from five volcanic centers
(Yellowstone, WY, USA; Valles, NM, USA; Long Valley, CA, USA; Taupo, NZ; Mount St. Helens, WA, USA). In
contrast to the minutes to hours of total ascent time extracted from embayment volatile profiles using constant-
rate models, our two-stage model resolves slower initial ascent times that span 3.5–11 h. Final ascent rates are
1–2 orders of magnitude faster than the initial extracted rates, in agreement with theoretical conduit flow model
predictions. Reassessment of embayments from the May 18th, 1980 eruption of Mount St. Helens results in an
initial stage of ascent consistent with the timing of magma arrival at the surface from the seismically-inferred
storage region (7–9 km) ~3.5 h after the initial blast, and a final stage of ascent (<1–5 min) in close agree-
ment with time-integrated bubble number densities. Our combined numerical and experimental results, and
reevaluation of natural datasets, suggest that, with the application of advanced models, the melt embayment can
provide a more nuanced picture of magma decompression timescales from the deep conduit to the surface.

1. Introduction

Magma ascent through the volcanic conduit is a dynamic process,
governed by complex feedbacks between intrinsic magma properties (e.
g., viscosity, density, crystallinity, vesicularity) and external controls (e.
g., conduit wall shearing, conduit geometry, etc.) (e.g., Gonnermann
and Manga, 2013). One critical factor driving acceleration and explosive
fragmentation of high-viscosity rhyolitic magmas is the interaction

between magma decompression, volatile exsolution and melt-coupling,
buoyancy, and resulting overpressure (e.g., Sparks, 2003). Recent nu-
merical conduit models have implemented time-dependent, non-steady-
state solutions to account for this dynamic feedback between magma
flow, conduit properties, and reservoir pressure (e.g., Melnik and
Sparks, 2006; Wong and Segall, 2019; Fig. 1a). Advances in numerical
techniques applied to modeling pyroclast textures (i.e., bubbles, crys-
tals) have followed suit, with incorporation of time-dependent

* Corresponding author.
E-mail address: behnazhosseini@montana.edu (B. Hosseini).

Contents lists available at ScienceDirect

Journal of Volcanology and Geothermal Research

journal homepage: www.journals.elsevier.com/journal-of-volcanology-and-geothermal-research

https://doi.org/10.1016/j.jvolgeores.2024.108211
Received 4 April 2024; Received in revised form 18 September 2024; Accepted 20 October 2024

Journal of Volcanology and Geothermal Research 456 (2024) 108211 

Available online 23 October 2024 
0377-0273/© 2024 The Author(s). Published by Elsevier B.V. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ). 

mailto:behnazhosseini@montana.edu
www.sciencedirect.com/science/journal/03770273
https://www.journals.elsevier.com/journal-of-volcanology-and-geothermal-research
https://doi.org/10.1016/j.jvolgeores.2024.108211
https://doi.org/10.1016/j.jvolgeores.2024.108211
https://doi.org/10.1016/j.jvolgeores.2024.108211
http://crossmark.crossref.org/dialog/?doi=10.1016/j.jvolgeores.2024.108211&domain=pdf
http://creativecommons.org/licenses/by/4.0/


decompression in models aimed at reproducing crystal size distributions
(CSD), microlite number densities (MND), and bubble number densities
(BND) observed in natural samples (Andrews and Befus, 2020; Haji-
mirza et al., 2021). These models have expanded upon earlier efforts
that applied a step-function (Hammer and Rutherford, 2002) or
constant-rate (Toramaru, 2006) boundary condition to extract decom-
pression rates. This time-integrated approach to modeling pyroclast
textures has allowed for seemingly decoupled, speedometer-specific
decompression rates spanning >4 orders of magnitude to be more
closely reconciled, as recently demonstrated for the 2008 eruption of
Chaitén in Chile and the 1991 eruption of Mount Pinatubo in the
Philippines (Hajimirza et al., 2021; Harris et al., 2024). The same
treatment, however, has not been extended to all petrological tools used
to infer magma decompression rates, limiting the comparability of re-
sults between different ascent rate meters and our ability to resolve
nuanced decompression histories. One such recorder frequently lever-
aged to estimate magma decompression rates (14 published studies
applied to 27 eruptions from 14 volcanic centers; Table 1) but still
confined to the realm of constant-rate modeling is the melt embayment
(open melt pocket or inclusion).

Over the last decade, embayments have gained a foothold among
petrological recorders used to estimate magma decompression rates in
explosive volcanic eruptions (Liu et al., 2007; Humphreys et al., 2008;
Lloyd et al., 2014; Ferguson et al., 2016; Myers et al., 2018, 2021;
Moussallam et al., 2019; Newcombe et al., 2020; Geshi et al., 2021;
Befus et al., 2022; Saalfeld et al., 2022; Zuccarello et al., 2022; Elms
et al., 2023; Harris et al., 2024). The premise behind this technique is
that as magma ascends, dissolved volatiles (H2O, CO2, S, Cl, F) diffuse
through melt embayment channels in response to decompression-
induced degassing of the surrounding melt and the consequent forma-
tion of a bubble at the outlet of the embayment (Fig. 1b). This generates

a chemical potential gradient between the embayment interior and
outlet, where the presence of the bubble implies a boundary layer of
melt with a volatile concentration fixed at the solubility-predicted,
equilibrium value. Upon eruption and explosive fragmentation,
embayment melt rapidly quenches to glass, locking in compositional
gradients that are measured using microanalytical techniques (e.g.,
nano secondary ion mass spectrometry: Lloyd et al., 2014, Newcombe
et al., 2020; Raman spectroscopy: Myers et al., 2021, Zuccarello et al.,
2022; Fourier transform infrared spectroscopy, Myers et al., 2018) or
other novel methods (e.g., calibrated back-scattered electron grayscale
intensities, Humphreys et al., 2008) and modeled numerically to extract
decompression rates. Although the embayment has proven its versatility
in recording decompression for a range of melt compositions (rhyolite:
Liu et al., 2007; rhyo-dacite: Humphreys et al., 2008; basaltic-andesite:
Lloyd et al., 2014; basalt: Ferguson et al., 2016), phenocryst hosts
(quartz: Saalfeld et al., 2022; plagioclase: Humphreys et al., 2008; py-
roxene: Myers et al., 2021; olivine: Moussallam et al., 2019), and rates
(0.0001–1 MPa s− 1), and has been subjected to rigorous numerical
(deGraffenried and Shea, 2021; Wei et al., 2024) and experimental
(Hosseini et al., 2023) evaluation, the execution of embayment
modeling using a constant-rate boundary condition has been met with
variable success. In particular, some studies have been unable to model
measured embayment profiles comprising multiple diffusing species
(H2O, CO2, S) directly from the inferred magma storage region,
requiring the application of two separate constant decompression rates
to best fit volatile profiles within individual embayments (e.g., Lloyd
et al., 2014; Myers et al., 2018; Geshi et al., 2021). Further, some studies
have applied modified (i.e., not storage) starting conditions to account
for negligible CO2 in embayments compared to co-erupted melt in-
clusions (e.g., Myers et al., 2018, 2021; Saalfeld et al., 2022; Elms et al.,
2023; Fig. 2; Table 1). These starting conditions assume either (i) some

Fig. 1. (a) Illustrative model showing the evolution of gas volume fraction (dark gray curve) and overpressure (light gray curve) for a rhyolitic magma approaching
the fragmentation front, modified from Gonnermann and Manga (2013). (b) Schematic of a two-stage decompression-diffusion model, where an initial slow stage of
decompression allows for partial re-equilibration of the embayment melt to lower volatile concentrations prior to a faster final decompression in the upper conduit.
Abbreviations are EMB – embayment, MI – melt inclusion, QTZ – quartz, and B – bubble. White arrows indicate the direction of diffusion from the interior to the
outlet of the embayment. (c) Pressure-time diagram illustrating how a two-stage model (EDiTS) better approximates a natural accelerating decompression path
compared with a standard constant-rate, continuous decompression model.

B. Hosseini and M. Myers Journal of Volcanology and Geothermal Research 456 (2024) 108211 

2 



degree of isobaric crystallization occurred prior to eruption, driving CO2
from the system via second boiling or (ii) exsolution of CO2 occurred
through a minor, initial decompression stage, potentially buoyancy-
driven (e.g., Myers et al., 2021; Saalfeld et al., 2022) (Fig. 2). Impor-
tantly, both scenarios can produce a moderate-high H2O and no CO2
starting condition distinctive from the melt inclusion population. A third
(iii) frequently applied ‘modified’ starting condition is the pressure
associated with the maximum volatile content preserved in the
embayment interior. The application of this initial condition assumes
that a period of slow decompression or stalling allowed for re-
equilibration of the embayment interior to a moderate H2O and no-
low CO2 condition (e.g., Myers et al., 2018; Saalfeld et al., 2022; Elms
et al., 2023) (Fig. 2). This scenario was also called upon by Lloyd et al.
(2014) to explain the consistently low embayment interior CO2 contents
of Fuego embayments (<200 ppm) compared to the much wider range
of co-erupted melt inclusion concentrations (100–750 ppm) (Table 1).
This same observation holds across all examined eruptions, where
embayment interiors retain <60 % of their initial (melt inclusion) CO2
concentrations (the exception being the Mesa Falls Tuff; Table 1). The
uncertainty in the starting volatile concentrations and corresponding
pressures required to reproduce measured profiles has non-negligible
effects on calculated decompression rates. Recent modeling of

embayments subjected to experimentally-controlled decompression
revealed that applying shallower-than-known starting conditions yiel-
ded decompression rates up to 5× faster than the experienced rate, with
equally good or better model fits (Hosseini et al., 2023). Similar results
were shown numerically by deGraffenried and Shea (2021), where using
the maximum H2O content in an embayment as the starting condition,
when the known value is greater, compounded errors from geometry by
an additional factor of 2–10×. This finding highlights the challenges and
limitations of a constant-rate modeling approach combined with un-
certain starting conditions. In these cases, it is likely that magmas
experienced more dynamic, multi-stage decompression histories and
thus require a new modeling framework that can account for these
complexities.

Here we present, validate, and apply a new MATLAB-based numer-
ical model (Embayment Decompression in Two Stages, or EDiTS) that
allows for decoupling and resolution of lower and upper conduit
decompression rates recorded by volatile (H2O ± CO2) profiles in em-
bayments while maintaining the usability and extensibility of available
constant-rate decompression models (Fig. 1b-c). This approach ap-
proximates the presumed accelerating path of magma ascent, driven by
a positive feedback between decompression, volatile exsolution, buoy-
ancy, and overpressure, with two separate stages of decompression

Table 1
Compilation of embayment studies, including number of embayments analyzed, maximum melt inclusion and embayment CO2 contents, maximum percentage of CO2
retained in the embayment interior (compared to melt inclusion concentrations), and number of good model fits in the original studies (including percentage of total).

Eruption No. EMB Analyzed Max MI CO2 (ppm) Max EMB CO2 (ppm) Max CO2 Retained (%) No. Good Model Fits

RHYOLITIC TO DACITIC Bandelier[1]*
(1.61 and 1.26 Ma)

22 216 0 0 9 (41 %)

Oruanui[2]*
(26.5 ka)

9 200 100 50 1 (11 %)

Huckleberry[2]*
(2.08 Ma)

9 800 450 56 4 (44 %)

Bishop[2]*
(760 ka)

13 225 275 2 (15 %)

Aira[3]

(30 ka)
>27 250 0 0 –

Okataina[4]

(21.9 ka–1314 CE)
21 126 0 0 21 (100 %)

Santorini[5]

(3.6 ka)
11 <200 0 0 10 (91 %)

Pinatubo[6]

(1991)
11 <20 0 0 6 (55 %)

Mesa Falls[7]

(1.3 Ma)
40 669 558 83 –

Mount St. Helens[8]*
(1980)

3 – – – 3 (100 %)

BASALTIC TO ANDESITIC Fuego[9]

(1974)
4 726 173 24 1 (25 %)

Etna[10]

(2013, 2015, 2018)
15 – – – –

Ambae[11]

(2018)
10 3899 1672 43 –

Kilauea[12]

(1500 CE, 1600 CE, 1959)
4 1200 300 25 4 (100 %)

Seguam[13]

(1977)
5 850 150 18 5 (100 %)

Note: all model fits based on a constant-rate modeling approach. (− ) indicates no data collected or otherwise available. (*) indicates eruptions revisited in this study.
1 Saalfeld et al. (2022).
2 Myers et al. (2018).
3 Geshi et al. (2021).
4 Elms et al. (2023).
5 Myers et al. (2021).
6 Harris et al. (2024).
7 Befus et al. (2022).
8 Humphreys et al. (2008).
9 Lloyd et al. (2014).
10 Zuccarello et al. (2022).
11 Moussallam et al. (2019).
12 Ferguson et al. (2016).
13 Newcombe et al. (2020).

B. Hosseini and M. Myers Journal of Volcanology and Geothermal Research 456 (2024) 108211 

3 



(Fig. 1b-c). Currently the code is formulated for silicic systems, from
which a greater number of embayments have been measured (>80 %;
Table 1) but can be modified for andesitic to basaltic systems (and the
inclusion of S) using published diffusivities and solubilities (Georgeais

et al., 2021 and references therein). We first benchmark EDiTS against
existing numerical embayment decompression-diffusion models, before
validating it using experimentally-controlled decompression of natural
quartz-hosted embayments. Finally, we apply EDiTS to existing

Fig. 2. Scenarios that can produce embayment H2O + CO2 contents modified (t2) from melt inclusion (storage) concentrations (t1), shown (a) schematically and (b-c)
graphically. (i) Isobaric crystallization drives second boiling of CO2, resulting in low to no CO2 in embayments despite some amount in melt inclusions. (ii) Buoyancy-
driven decompression similarly results in CO2 exsolution in a saturated system and loss from the embayment, while retaining moderate-high H2O. (iii) A period of
slow decompression or stalling during magma ascent can result in both lower H2O in, and complete loss of CO2 from, the embayment compared to melt inclusions. In
these scenarios, applying a constant-rate model with melt inclusion starting conditions may lead to poor fits to H2O ± CO2 gradients.

B. Hosseini and M. Myers Journal of Volcanology and Geothermal Research 456 (2024) 108211 

4 



embayment datasets (H2O ± CO2) from highly explosive silicic erup-
tions from five volcanic centers (Yellowstone, WY, USA; Valles, NM,
USA; Long Valley, CA, USA; Taupo, NZ; Mount St. Helens, WA, USA) to
reassess timescales of magma decompression and ascent. Through this
multi-faceted approach, we aim to investigate whether (1) embayments
can faithfully record, and models can be developed to extract, more
complex and realistic magma decompression paths, and (2) a two-stage
embayment modeling framework provides magma decompression rates
in better agreement with independent estimates for historic eruptions
fromwell-monitored volcanoes (e.g., May 18th, 1980Mount St. Helens).

2. Model set-up and validation

2.1. Numerical model development and benchmarking

The available 1-, 2-, and 3-D models for recreating volatile diffusion
profiles within embayments apply numerical (finite-difference) solu-
tions to the diffusion equation (e.g., unnamed model (1D): Myers et al.,
2018; Hosseini et al., 2023; EMBER (1D): Georgeais et al., 2021; un-
named model (1D, 3D): deGraffenried and Shea, 2021; Embayment-
Diffuser (1D, 2D): Befus et al., 2022). Here, we will focus on the 1D
model set-up (majority of embayment models), since a 1D solution forms
the basis for EDiTS. In these models, H2O and CO2 (e.g., Myers et al.,
2018; Befus et al., 2022; Hosseini et al., 2023) and, in some cases, S (e.g.,
Lloyd et al., 2014; Ferguson et al., 2016; Georgeais et al., 2021; Zuc-
carello et al., 2022) diffuse through a melt-filled embayment in response
to the formation of a bubble at the channel outlet, which establishes a
concentration step function between the bubble and surrounding melt.
The presence of this bubble is a result of decreasing volatile solubilities
and exsolution induced by magma decompression. Volatile element
solubilities are well-established for a range of melt compositions and are
pressure-, temperature-, and fluid composition-dependent (e.g., New-
man and Lowenstern, 2002; Witham et al., 2012). Diffusivities have
been experimentally determined and are pressure, temperature, H2O-
concentration and, in the case of S, oxygen fugacity dependent (Nowak
and Behrens, 1997; Zhang and Behrens, 2000; Freda et al., 2003, 2005;
Nowak et al., 2004; Ni and Zhang, 2008; Zhang and Ni, 2010; Zhang
et al., 2010). At the crystal-melt interface, a no-flux boundary condition
(Neumann boundary) is imposed, justified by the incompatibility of the
diffusing volatile elements in the host crystal, such that diffusive flux is
uni-directional along the length of the embayment. The boundary con-
dition at the bubble-melt interface is updated at each decompression
time-step (Dirichlet boundary), with the volatile concentration assumed
to be in local equilibrium with the exsolved fluid present in the melt
surrounding the crystal. We note that this assumption of an equilibrium
degassing pathway is a current limitation of embayment modeling.
Previous decompression experiments have shown that bubble nucle-
ation delay and disequilibrium degassing in crystal-poor rhyolitic melts
can generate significant overpressures during magma ascent (e.g.,
Mangan and Sisson, 2000; Cichy et al., 2011), and recent numerical
studies have demonstrated the non-negligible effect of the assumed
degassing mechanism on calculated embayment decompression rates
(up to one order of magnitude; deGraffenried and Shea, 2021). Future
embayment models would thus benefit from a greater understanding of
bubble nucleation mechanisms in rhyolitic melt, as well as the devel-
opment of general models describing pressure-H2O concentration re-
lationships during disequilibrium degassing (e.g., deGraffenried and
Shea, 2021).

Model input parameters include (1) initial dissolved volatile con-
centration in the melt (CH2O, wt%; CCO2, ppm; CS, ppm) which is typi-
cally determined from co-erupted melt inclusions or otherwise a
‘modified’ starting condition as discussed above (Fig. 2); (2) initial
pressure (Pi; MPa), which is related to volatile concentrations via solu-
bility models (e.g., VolatileCalc, Newman and Lowenstern, 2002; Mag-
maSat, Ghiorso and Gualda, 2015); (3) temperature (T; ◦C), which is
often held constant, an assumption that is likely more appropriate in

silicic rather than basaltic systems (Mastin and Ghiorso, 2001); (4) the
mass of pre-exsolved fluid; and (5) embayment length from the interior
crystal-glass interface to the bubble-glass interface (x; μm). To date, all
published codes simulate a constant-rate, continuous decompression
path using some or all of the above model assumptions and input
parameters.

Our MATLAB-based two-stage decompression-diffusion model,
EDiTS, is applicable to H2O ± CO2 diffusing through rhyolitic melt
embayments, and incorporates 1D time- and concentration-dependent
diffusion defined by Fick’s second law:

∂Ci

∂t =
∂

∂x

(

Di
∂Ci

∂x

)

(1)

where C is the concentration of the diffusing species i in the melt (wt% or
ppm), D is the diffusivity of the volatile species in silicate melt (μm2 s− 1)
as a function of the total dissolved H2O concentration, x is distance (μm),
and t is time (s). Solving Eq. (1) numerically requires defining appro-
priate initial (t = 0, 0≤x≤L) and boundary (t > 0, x = 0, L) conditions at
the embayment outlet (x= 0) and interior (x= L), given in ourmodel by:

Ci0≤x≤L, t=0 = uniform (2a)

Cix=0, t>0 = Cisolubility (2b)

∂Cix=L

∂x = 0 (2c)

where Eq. (2a) describes an initially uniform volatile concentration C of
diffusing species i along the entire length of the embayment, Eq. (2b)
defines a pressure- and time-dependent boundary condition at the
bubble-melt interface (outlet) determined by the equilibrium solubility
of the diffusing species i, and Eq. (2c) defines a no-flux boundary con-
dition at the melt-crystal interface (interior). EDiTS incorporates the
empirical expressions of Liu et al. (2005) to model H2O and CO2 solu-
bility (see Eqs. 2a and 2b in Liu et al., 2005). For H2O and CO2 diffu-
sivities, we use equations given by Zhang and Behrens (2000),
applicable to a wide range of conditions (400–1200 ◦C, 0.1–810 MPa,
and 0.1–7.7 wt% total H2O), and Zhang et al. (2007), respectively.

The code and all relevant documentation are available on GitHub
and archived on Zenodo (https://zenodo.org/record/10042159). EDiTS
expands on the models of Myers et al. (2018) and Hosseini et al. (2023)
by implementing a new free parameter iteration structure. The original
model cycles through a range of unknown inputs including initial
exsolved fluid content (r; wt%), decompression rate (dP/dt; MPa s− 1),
and fragmentation pressure (Pf ; MPa), where diffusion ceases due to
rapid quenching of embayment melt. A χ2 misfit is calculated based on
the difference between the measured and modeled profiles with this
value determined iteratively for every combination of parameters:

χ2 =
∑L

x=1

(
Ci,xmeasured − Ci,xmodeled

)2

σ2i,x
(3)

Where Ci,xmeasured is the measured concentration of volatile species i at
position x in the compositional profile, Ci,xmodeled is the modeled volatile
concentration, and σi,x is the uncertainty of the measured concentration.
The best-fit profile is determined by the combination of parameters that
minimizes χ2 (Eq. 3) for fitting H2O and CO2 simultaneously. Here, we
apply the samemisfit calculation but instead cycle through the following
three user-defined, incremented parameters: initial (deeper; dP/dti) and
final (shallower; dP/dtf ) decompression rates, as well as the pressure
where this transition occurs (Pt). In this framework, the code now con-
sists of two constant-rate, continuous decompression stages where, at a
given transition pressure, the H2O ± CO2 diffusion profiles produced
during the first (initial) stage of decompression are used as inputs
(starting conditions) for the second (final) stage of decompression. The
tradeoff of this model set-up is that initial exsolved gas content and

B. Hosseini and M. Myers Journal of Volcanology and Geothermal Research 456 (2024) 108211 

5 

https://zenodo.org/record/10042159


quench pressure are no longer free model parameters and thus need to
be estimated independently. For example, the quench pressure can be
determined from the H2O ± CO2 concentration near the embayment
outlet or adhering matrix glass (e.g., Lloyd et al., 2014), and the
exsolved gas content can be estimated using other embayment models
(e.g., Myers et al., 2018). Fig. A1 provides a graphical flow-chart sum-
mary of the code architecture.

While there are no independent two-stage models to benchmark
EDiTS against, we are able to validate our code by using constant-rate,
continuous models to produce synthetic H2O and CO2 concentration
profiles characterized by known deeper and shallower decompression
rates (Fig. A2). Using the code of Hosseini et al. (2023), we first modeled
an initial decompression stage (dP/dti) from 150 MPa (3 wt% H2O and
~680 ppm CO2) to an intermediate pressure (Pt; 100, 90, or 80MPa) at a
specific rate (0.001, 0.005, or 0.009 MPa s− 1). The resulting H2O and
CO2 concentration profiles from this initial stage were then used as the
starting conditions for the final decompression stage

(
dP/dtf

)
from the

intermediate pressure to the quench pressure (Pf ; 30, 15, or 10 MPa) at a
specific rate (0.01, 0.05, or 0.09 MPa s− 1). The resulting ‘synthetic’
profiles were input as ‘measured’ H2O + CO2 profiles for recreation by
EDiTS. For the three distinct synthetic profiles that were produced using
different initial and final decompression rates and transition pressures,
we were able to successfully extract the known decompression paths
with χ2 misfits of 0 using EDiTS (Fig. A2; Table 2). This is an expected
result since both models use the same parameterizations, but never-
theless demonstrates that the two-stage code is performing in the
intended manner. We are confident that tests using synthetic volatile
profiles generated by other 1D models should lead to similarly low χ2
misfits in EDiTS, with any divergence resulting from the volatile solu-
bility and diffusivity equations defined within each model.

We also test how EDiTS performs in recreating synthetic H2O and
CO2 profiles produced by a single stage, constant-rate decompression.
Here we find that the model can recover initial and final rates that
bracket the experienced constant rate (i.e., a false positive result)
(Table 2). For instance, for synthetic H2O + CO2 profiles produced by a
constant decompression rate of 0.008 MPa s− 1 from 150 to 30 MPa,
EDiTS recovers a best-fit initial decompression rate of 0.004 MPa s− 1

from 150 to 80 MPa, and 0.01 MPa s− 1 from 80 to 30 MPa. We discuss
the implications of this result later.

2.2. High pressure-temperature decompression experiments

Our second validation approach is to apply EDiTS to model embay-
ment profiles produced during two-stage, mixed-volatile (H2O + CO2)
and pure H2O-saturated decompression experiments. Hosseini et al.
(2023) present a method to successfully re-equilibrate natural quartz-
hosted embayments to new starting pressures and fluid compositions
using a rapid-quench cold-seal pressure vessel at the University of
Oregon. Each experimental charge consists of Bishop ash (<60 μm)

doped with 10 wt% oxalic acid dihydrate (OAD), natural quartz crystals
with embayments, and nickel‑nickel oxide (NNO) buffer. The embay-
ment melt is re-equilibrated at high pressures (150 MPa) and tempera-
tures (780 ◦C), where the OAD breaks down to form an H2O-CO2 fluid
mixture. The charge is then subjected to controlled, constant-rate,
continuous decompression. Embayments record volatile diffusion pro-
files characterized by known decompression rates that were able to be
successfully recovered via numerical diffusion modeling. Following the
methodology outlined in Hosseini et al. (2023), we further develop the
technique by conducting multi-stage, continuous decompression ex-
periments, with the goal of testing the ability of EDiTS to recover known
two-stage decompression paths and transition pressures from re-
equilibrated natural embayments (Fig. A3; Table A1). One H2O-satu-
rated experiment (E4–2) was conducted in which we decompressed the
charge from 150 MPa at an initial rate of 0.006 MPa s− 1 before tran-
sitioning at 100 MPa to 0.02 MPa s− 1 and air quenching at 30 MPa (total
decompression time 3.28 h; Fig. 3; Fig. A3; Table A1). In the mixed-
volatile (H2O + CO2) saturated experiment (E3–2), we decompressed
the charge from 150 MPa at 0.009 MPa s− 1 before ramping up, at 80
MPa, to 0.06 MPa s− 1 and air quenching at 30 MPa (total decompression
time 2.39 h; Fig. 4; Fig. A3; Table A1). We additionally conducted one
stalling experiment in which the charge was decompressed from 150
MPa to 80MPa at a constant rate of 0.02MPa s− 1, stalled for three hours,
and then decompressed at the same constant rate until quench at 30 MPa
(total decompression time 1.66 h, total experiment time 4.66 h; Fig. A3;
Table A1). Three glassy embayments recovered from these experiments
(one from each) were prepared for analyses of H2O ± CO2 transects
using a Nicolet iN10 infrared microscope and integrated spectrometer at
Montana State University. The embayment from the H2O-saturated
experiment preserves a gradient from ~2.5 wt% H2O at the outlet that
plateaus at ~3.5 wt% H2O at the interior (Fig. 3; Table A2). Lower H2O
(2–3 wt%) characterizes the gradient in the embayment from the mixed-
volatile experiment, with CO2 ranging from 175 to 250 ppm (Fig. 4;
Table A2). The stalled embayment preserves a H2O gradient from 1 to
1.9 wt% with CO2 ranging from 35 to 175 ppm (Table A2). Full exper-
imental set-up and procedures, along with FTIR spectroscopy data
collection and processing, are outlined in Hosseini et al. (2023) and
provided in the Supplementary Material (Text A1-A2, Tables A1-A2, Fig.
A3).

We applied EDiTS to assess how accurately the known experimental
decompression paths could be extracted from preserved embayment
H2O ± CO2 concentration gradients, allowing dP/dti, dP/dtf , and Pt to
vary freely. Modeling the three recovered embayments from two-stage
experiments, one from a continuously decompressed H2O-saturated
experiment (E4–2), one from a continuously decompressed mixed-
volatile (H2O + CO2) experiment (E3–2), and the final from a stalled
mixed-volatile (H2O + CO2) experiment (E1–2), we note some signifi-
cant observations. For the H2O-saturated experiment (E4–2), we recover
equally good fits to the H2O gradient using the constant-rate (0.011 MPa

Table 2
Summary of results from applying constant-rate (Hosseini et al., 2023) and two-stage (EDiTS) models to experimental and synthetic H2O± CO2 profiles. ‘True positive’
results occur when a given model recovers known synthetic or experimental decompression conditions. ‘False positive’ results occur when a given model recovers non-
experienced synthetic or experimental decompression conditions.

Model

Constant-rate
(Hosseini et al., 2023)

Two-stage
(EDiTS)

Experimental Constant
rate

Validated by Hosseini et al. (2023) for both H2O-only and H2O + CO2

profiles.
This study. Can produce a false positive result for both H2O-only and
H2O + CO2 profiles.

Two-stage This study. Cannot simultaneously fit H2O + CO2 profiles. Can produce a
false positive fit to a H2O-only profile.

This study. Can successfully fit H2O + CO2 profiles with known
dP/dti, Pt , and dP/dtf . No unique fit to H2O-only profiles.

Synthetic Constant
rate

Perfect fits, N/A. This study. Can produce a false positive result, where dP/dti, Pt , and
dP/dtf are recovered with a misfit of 0.

Two-stage This study. Can produce a false positive fit to both H2O-only and H2O+ CO2

profiles that roughly averages the decompression path.
This study. Can successfully recover dP/dti, Pt , and dP/dtf with a
misfit of 0.

B. Hosseini and M. Myers Journal of Volcanology and Geothermal Research 456 (2024) 108211 

6 



s− 1; χ2 = 0.09) and two-stage (0.002 ➔ 0.012 MPa s− 1; χ2 = 0.08)
models (Fig. 3; Table 2). This highlights the previously-recognized
challenge of discerning unique embayment decompression paths using
only a single diffusing volatile specie (e.g., Lloyd et al., 2014; Hosseini
et al., 2023). We also observe that the constant-rate model recovers a
best-fit decompression rate that roughly averages the two-stage
decompression path, confirming that a constant-rate modeling
approach can blur the effects of acceleration or multi-stage ascent his-
tories. Translating these recovered decompression rates to ascent time
via:

tconstant =
Pi − Pf
dP/dt

(4a)

ttwo− stage =
Pi − Pt
dP/dti

+
Pt − Pf
dP

/
dtf

(4b)

where Pi, Pt, and Pf are the initial, transition, and final pressures and
dP/dt is the decompression rate, we find that the constant-rate model
(3.03 h) closely approximates the known total ascent time (i.e., exper-
iment duration, 3.28 h), while the two-stage model (8.56 h) poorly
approximates the known ascent time. However, the two-stage model
reasonably approximates the final, rapid stage of decompression from
100 to 30 MPa (97.2 mins recovered vs. 58.3 mins known), while this
time cannot be resolved using the constant-rate model. While a wide
range of initial decompression rates can fit the measured profile
(0.001–0.009 MPa s− 1 translating to 1.5–13 h of ascent time), the final
stage of decompression is much more constrained (good fits between
0.01 and 0.02 MPa s− 1 equating to 1.94 h–58 mins of ascent time).
Overall, these results highlight that most of the error in recovering the
experienced path from H2O profiles results from the poor constraints on
the initial, slower stage of decompression, where in this case EDiTS re-
covers an initial ascent time (6.9 h) a factor of three greater than that

experienced (2.3 h) (Fig. 3).
Conversely, for the continuously-decompressed two-stage mixed-

volatile experiment (E3–2), we are only able to fit the H2O and CO2
embayment profiles simultaneously using EDiTS (Fig. 4). For this
embayment, the constant-rate model recovers a best fit decompression
rate (characterized by an unacceptable fit to the CO2 profile) of 0.025
MPa s− 1, with a recovered ascent time (1.33 h) that underestimates the
experienced time (2.39 h). This code again recovers a decompression
rate that averages the known initial (0.009 MPa s− 1) and final (0.06
MPa s− 1) stages. On the other hand, the two-stage model recovers an
initial rate of 0.008 MPa s− 1 (0.009 MPa s− 1 expected) that transitions at
80 MPa (the imposed transition pressure) to 0.1 MPa s− 1 (0.06 MPa s− 1

expected; Fig. 4). This translates to a total ascent time of 2.57 h recov-
ered by the two-stage model, in close agreement with the known
decompression duration (2.39 h), as well as a very well-constrained final
stage of decompression (8.3 mins recovered vs. 13.8 mins expected).
These results suggest that modeling multiple diffusing species simulta-
neously provides crucial constraints for resolving the true decompres-
sion path, and more significantly, that a mixed-volatile system that has
experienced multi-stage or accelerating magma decompression may not
be successfully modeled using a constant-rate approach with ’melt in-
clusion’ starting conditions (Fig. 4). Finally, in applying the constant-
rate model of Hosseini et al. (2023), as well as EDiTS, to the stalled
embayment (E1–2), we are unable to simultaneously fit H2O and CO2
concentration gradients. This optimistically suggests that neither
constant-rate nor two-stage models may reproduce embayment volatile
profiles generated by more complex scenarios, such as hours of crustal
stalling.

A final important observation is that the saturation pressure pre-
served in the interior of each of the two continuously-decompressed (i.
e., not stalled) experimental embayments (3.55 wt% H2O = 80 MPa for
E4–2 and 3.16 wt%H2O+ 264 ppm CO2= 110MPa for E3–2; calculated
using the solubility model of Liu et al. (2005) in VESIcal: Iacovino et al.,

Fig. 3. Modeling results for embayment E4–2, shown in transmitted light in the upper left and annotated with the measured transect, with inset showing a cartoon
rendering of the host crystal for context. E4–2 was recovered from a H2O-saturated experiment that decompressed from 150 to 100 MPa at 0.006 MPa s− 1, and from
100 to 30 MPa at 0.02 MPa s− 1. Horizontal gray bar denotes the starting condition used to model this embayment (5.1 wt% H2O), determined by an H2O-saturated
equilibrium (dwell) experiment conducted by Hosseini et al. (2023) at 780 ◦C and 150 MPa. Measured transect shown as teal circles. Black curves indicate the
expected profile and light teal curves the recovered best-fit profile, both for the initial stage (dashed) and the final stage (solid) of decompression. Also shown is the
best-fit profile using the constant-rate model of Hosseini et al. (2023) (dark teal curve). Shaded gray curved region signifies the range of good model fits to the
measured profile at a fixed transition pressure (100 MPa). Analytical error shown as light gray vertical bars on measured points. Lower right χ2 misfit contour plot
illustrates the range of initial decompression rates (x-axis) at a fixed transition pressure (100 MPa) that produce good model fits, compared to much more constrained
final rates (y-axis). The misfits for the expected (orange star) and recovered (yellow star) model fits are shown. (For interpretation of the references to colour in this
figure legend, the reader is referred to the web version of this article.)

B. Hosseini and M. Myers Journal of Volcanology and Geothermal Research 456 (2024) 108211 

7 



2021) is within ~20–30 MPa of the transition pressure imposed in the
experiment (Figs. 3 and 4). Although this is not necessarily a surprising
result, as the increased final decompression rate will provide less time
for interior relaxation to occur, it does provide some useful constraints
on model input parameters. We therefore suggest that when applying
EDiTS to natural embayments the user should use the saturation pres-
sure preserved in the embayment interior plateau to inform the range of
transition pressures cycled through during a model run. This will also
help reduce the model run time and convergence on a best-fit solution,
since there is a significant tradeoff between initial decompression rate
and transition pressure which results in many equally good fits to the
measured profiles.

To again test the result of applying a multi-stage decompression code
to an embayment that experienced a known constant-rate path, we
applied EDiTS to a subset of embayments from the pure H2O- andmixed-
volatile (H2O + CO2) saturated, constant-rate experiments presented in
Hosseini et al. (2023). The H2O-saturated experiment (E8–1) was
decompressed from 150 to 30 MPa at 0.015 MPa s− 1 and the H2O + CO2
experiment (E3–1) decompressed from 150 to 30 MPa at 0.008 MPa s− 1.
For embayment E8–1, EDiTS recovers a best-fit profile characterized by
an initial decompression rate of 0.004 MPa s− 1 that transitions at 60
MPa to 0.015 MPa s− 1. Although the recovered final stage of decom-
pression equals the imposed rate (0.015 MPa s− 1), the initial slow stage
adds several hours to the total ascent time (e.g., ~6.5 h recovered vs. 2.2
h imposed) while maintaining a similarly low misfit (χ2 = 0.02 for two-

stage vs. 0.03 for constant-rate). We find a similar result for embayment
E3–1 from the mixed-volatile experiment; while EDiTS recovers a best-
fit profile defined by an initial decompression rate of 0.007 MPa s− 1 that
transitions at 50 MPa to 0.01 MPa s− 1 (~4.5 h of ascent time), the
experienced constant decompression rate was 0.008 MPa s− 1 from 150
to 30 MPa (~4.2 h of ascent time). Thus, similar to the synthetic test
results, we find that in both H2O-only and mixed-volatile systems, EDiTS
can produce a false positive result, meaning a best-fit diffusion profile
characterized by two stages of decompression can be recovered even
when the embayment experienced a constant rate of decompression
(Table 2). We argue, however, that this result is a red herring and does
not have any bearing on the application of EDiTS to natural systems,
since most natural magmas are unlikely to experience a constant, linear
decompression history from storage to fragmentation.

3. Model application

3.1. Revisiting existing embayment studies

To date, there have been 14 embayment studies published on rhyo-
litic to basaltic systems, representing 27 eruptions (14 pre-historic and
13 historic) and including over 200 embayments (Table 1). Below, we
present results from applying EDiTS to a subset of previously modeled
embayments (n = 5), reevaluating four pre-historic eruptions (all
caldera-forming, rhyolitic; Table 3), as well as reassessing the

Fig. 4. Modeling results for embayment E3–2, shown in transmitted light in the upper right and annotated with the measured transect, with inset showing a cartoon
rendering of the host crystal for context. E3–2 was recovered from a mixed-volatile (H2O + CO2)-saturated experiment that decompressed from 150 to 80 MPa at
0.009 MPa s− 1, and from 80 to 30 MPa at 0.06 MPa s− 1. Horizontal gray bars on upper (H2O) and lower (CO2) plots denote the starting conditions used to model this
embayment (3.4 wt% H2O and 425 ppm CO2), as determined by the mixed-volatile saturated dwell experiments of Hosseini et al. (2023), conducted at 780 ◦C and
150 MPa. Measured transects shown as green circles. Black curves indicate the expected profile, and light green curves the recovered best-fit profile, both for the
initial (dashed) and final (solid) stages of decompression. Also shown is the best-fit profile using the constant-rate model of Hosseini et al. (2023) (dark green curve).
Shaded gray curved region signifies the range of good model fits to the measured profile at a fixed transition pressure (80 MPa). Analytical error shown as light gray
vertical bars on measured points. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

B. Hosseini and M. Myers Journal of Volcanology and Geothermal Research 456 (2024) 108211 

8 



embayments (n = 3) from the closely monitored and well-documented
May 18th, 1980 eruption of Mount St. Helens (rhyo-dacitic) (Table 3;
Fig. 5). To allow for comparability of EDiTS results with those presented
in the original constant-rate modeling studies, we use the same starting
H2O ± CO2 and corresponding pressures, quench pressures, and
exsolved gas contents as the original studies.

3.1.1. Applications to pre-historic eruptions
The majority of embayments from pre-historic caldera-forming

eruptions have lacked measurable CO2, although co-erupted melt in-
clusions often contain a moderate amount (Table 1). The sparsity of CO2
in embayments from many eruptions (e.g., Bandelier, NM, USA; Bishop,
CA, USA; Oruanui, NZ) has been interpreted to result from one or both of
the following scenarios: either CO2 was driven out of the system prior to
eruption due to isobaric crystallization (e.g., second boiling), or the
system experienced an initial slow stage of ascent prior to eruption
permitting CO2 re-equilibration (Myers et al., 2018, 2021; Saalfeld et al.,
2022). We applied EDiTS to these systems to test previously presented
ideas, allowing for melt inclusion H2O ± CO2 concentrations to be used
as starting conditions. We remodeled measured volatile gradients from
five select embayments from four caldera systems (Yellowstone, WY;
Valles, NM; Long Valley, CA; Taupo, NZ), targeting samples with both
H2O and CO2 gradients present, when possible, for better modeling
constraints (Table 3). We find that for all systems, volatile gradients are
best fit with an initial deeper decompression stage equivalent to 3.5–11
h, a final shallower stage equivalent to hours to 10s of minutes, and
transition pressures (70–155 MPa) that agree closely with embayment
interior pressures (Table 3). This lends credence to the idea that for
eruptions with CO2-bearing melt inclusions, embayments with low-no
CO2 could be a record of longer (>10 h) initial decompression, allow-
ing partial to full re-equilibration to occur, as was suggested for 1974
eruption of Fuego (Lloyd et al., 2014), the Late Bronze Age eruption of
Santorini (Myers et al., 2021), and the 1.6 and 1.2 Ma eruptions of the
Bandelier Tuff (Saalfeld et al., 2022). Further, our two-stage model may
better constrain the final, explosive stage of each eruption, providing
times on the order of 10s of minutes, where the constant-rate model
produces total ascent times on the order of hours. However, because pre-
historic eruptions lack independent, real-time constraints on magma
ascent rates, we are limited in terms of howwe interpret and validate the
EDiTS results, and so, as a case study, we turn our assessment to a
modern, instrumented eruption that includes both independent petro-
logic (e.g., bubble number densities) and real-time observational
decompression rate estimates: the May 18th, 1980 Plinian eruption of
Mount St. Helens.

3.1.2. Applications to historic instrumented eruptions
Integrating petrologic records with real-time observations of magma

ascent and eruption remains a great ambition of volcano science (e.g.,
Rasmussen et al., 2018; Ruth et al., 2018; Caricchi et al., 2021). While
the fidelity of embayments to record magma decompression has been
well-established (e.g., Hosseini et al., 2023), the recovered rates have
largely remained decoupled from independent records (e.g., Fig. 6a). We
argue that this discordance between different decompression rate esti-
mates is in part an artifact of the model framework. We next demon-
strate the ability of EDiTS to reconcile embayment- and monitoring-
based magma ascent rates, as applied to the 1980 Plinian eruption of
Mount St. Helens (Washington, USA). OnMay 18th, 1980 at 0832 PDT, a
magnitude 5 earthquake triggered the failure of the north flank leading
to, in sequence, a debris avalanche, lateral blast, and eruptive column
(Voight et al., 1981; Kieffer, 1981; Eichelberger and Hayes, 1982).
Sustained Plinian activity continued for nine hours, with the eruption
column reaching a mean height of 16 km (Carey and Sigurdsson, 1985;
Carey et al., 1990). The record of strong ground shaking for the May
18th eruption reveals two distinct phases of activity, where the first
coincides with the M5 earthquake, landslide, and lateral blast, lasting
~15 min (Scandone and Malone, 1985). The onset of the second phaseTa

bl
e
3

Re
su
lts

fr
om

re
m
od
el
in
g
em

ba
ym

en
ts
fr
om

fo
ur

si
lic
ic
ca
ld
er
a-
fo
rm

in
g
(H
uc
kl
eb
er
ry
,B

an
de
lie
r,
Bi
sh
op
,O

ru
an
ui
)
an
d
on
e
ar
c
(M

ay
18
th
,1
98
0
M
ou
nt

St
.H

el
en
s)
er
up
tio

n
us
in
g
ED

iT
S,
an
d
co
m
pa
ri
so
n
w
ith

or
ig
in
al

co
ns
ta
nt
-r
at
e
m
od
el
in
g
re
su
lts

fr
om

:[
†]
M
ye
rs
et
al
.(
20
18
),
[‡
]
Sa
al
fe
ld
et
al
.(
20
22
),
an
d
[ο
]
H
um

ph
re
ys

et
al
.(
20
08
).
(-
)
in
di
ca
te
s
th
at
CO

2
w
as

no
tm

ea
su
re
d
in
th
e
or
ig
in
al
st
ud
y,
an
d

χ2
m
is
fit
s
w
er
e
no
tr
ep
or
te
d.

Co
ns
ta
nt
-r
at
e
m
od
el

Tw
o-
st
ag
e
m
od
el
(E
D
iT
S)

Er
up
tio

n
Sa
m
pl
e

St
ar
tin

g
pr
es
su
re

(M
Pa
)

H
2O

(w
t%

)
CO

2

(p
pm

)
D
ec
om

pr
es
si
on

ra
te

(M
Pa

s−
1 )

Q
ue
nc
h

pr
es
su
re

(M
Pa
)

A
sc
en
t

ra
te

(m
s-1
)

A
sc
en
t

tim
e

(m
in
s)

χ2
D
ec
om

pr
es
si
on

ra
te

(i
ni
tia
l)

(M
Pa

s−
1 )

D
ec
om

pr
es
si
on

ra
te

(fi
na
l)

(M
Pa

s−
1 )

Tr
an
si
tio

n
pr
es
su
re

(M
Pa
)

A
sc
en
tr
at
e

(i
ni
tia
l)

(m
s-1
)

A
sc
en
tr
at
e

(fi
na
l)

(m
s-1
)

A
sc
en
tt
im
e

(i
ni
tia
l)

(m
in
)

A
sc
en
tt
im
e

(fi
na
l)

(m
in
)

χ2

H
uc
kl
eb
er
ry

†
M
M
7
RE

no
.

10
15
0

4.
1

32
5

0.
02
7

17
1.
06

82
.2

0.
20

0.
00
2

0.
06
3

70
0.
08

2.
47

68
0

13
.8

0.
75

Ba
nd
el
ie
r‡

T1
–3

16
2

5.
5

57
0.
00
5

15
0.
18

50
0

0.
33

0.
00
4

0.
1

20
0.
16

1.
95

59
2

1.
8

0.
08

Ba
nd
el
ie
r‡

T2
–3

12
1

4.
5

11
7

0.
00
7

30
0.
28

21
5

4.
22

0.
00
1

0.
03
7

80
0.
04

1.
45

67
7

19
.2

0.
55

Bi
sh
op

†
BT

F9
_1
38

RE
no
.2

22
0

5.
4

47
5

0.
04
2

30
1.
53

83
.4

0.
9

0.
00
5

0.
06

15
5

0.
20

2.
35

21
7

43
0.
68

O
ru
an
ui

†
P1
97
0-
A
RE

no
.6

19
0

5.
0

32
5

0.
01
7

45
0.
67

14
2.
2

0.
34

0.
00
5

0.
08

75
0.
20

3.
14

38
3

6
0.
48

M
SH

ο
KV

C5
18
b

23
0

6.
4

-
1.
4

33
54
.8
9

2.
4

-
0.
00
9

0.
76

60
0.
35

29
.8
0

39
2

0.
6

0.
9

M
SH

ο
M
SH

1-
3

23
0

6.
4

-
1.
6

27
62
.7
3

2.
4

-
0.
01

0.
23

10
0

0.
39

9.
02

28
7

5.
4

4.
1

M
SH

ο
M
SH

1-
6

23
0

6.
4

-
1.
3

9
50
.9
7

3
-

0.
00
8

0.
65

90
0.
31

25
.4
8

37
9

1.
8

4.
8

B. Hosseini and M. Myers Journal of Volcanology and Geothermal Research 456 (2024) 108211 

9 



manifested as strong ground motion accompanied by an eruptive col-
umn ~3.5 h after the lateral blast; this was used to calculate magma
ascent on the order of 0.6–0.7 m s− 1 from 7 to 9 km depth (Fig. 6;
Scandone and Malone, 1985; Eichelberger, 1995). This sequence has
been interpreted as the eruption of a shallow magmatic body (lateral
blast) that triggered the ascent of deep, volatile-rich magma (Plinian
phase). Currently there is an orders-of-magnitude discrepancy between
the rapid embayment-based decompression rates (1.3–1.6 MPa s− 1

assuming CH2Oi= 6.4 wt% and 0.9–1.2 MPa s− 1 assuming CH2Oi= 4.6 wt
%) presented by Humphreys et al. (2008) using a constant-rate model,
and those estimated by independent real-time observations and seis-
micity (0.01–0.02 MPa s− 1; Endo et al., 1981, Scandone and Malone,
1985, Eichelberger, 1995) (Fig. 6a). We apply EDiTS to remodel
grayscale-calibrated H2O concentration gradients in these plagioclase-
hosted embayments (n = 3) from the Plinian (May 18th) phase of the
eruption (Humphreys et al., 2008), allowing for the incorporation of an
initial deeper stage of decompression (Table 3; Fig. 5).

Assuming the same input parameters presented by Humphreys et al.
(2008) and derived from Blundy and Cashman (2005) and Blundy et al.
(2006) (i.e., CH2Oi = 6.4 wt% H2O, Pi = 230 MPa, T = 880 ◦C), EDiTS
recovers best fit profiles characterized by an initial slow decompression
(0.008–0.01 MPa s− 1) followed by a final stage up to two orders of
magnitude faster (0.2–0.8 MPa s− 1). Combined with recovered transi-
tion pressures (60–100 MPa, 2.4–4 km), initial ascent rates and times

range from 0.3 to 0.4 m s− 1 and 4.8–6.5 h, respectively, with final ascent
rates of 9–30 m s− 1 translating to <1–5 min in the shallow conduit
(Fig. 5). These final rates are significantly faster than those determined
for the larger caldera forming eruptions above, but are in agreement
with numerical conduit model predictions (e.g., Mastin, 2005). More
importantly, our results bring petrologic estimates of magma ascent into
closer agreement with the monitoring record by supporting the presence
of an initial slower stage of magma decompression on the order of hours,
where previously the total embayment-derived ascent time was 2.4–3
min (Fig. 6). We additionally find reconciliation between the final
decompression rates recorded by the embayment population (0.2–0.8
MPa s− 1) and updated bubble number density estimates based on time-
integrated, heterogeneous bubble nucleation models (0.8 MPa s− 1;
Hajimirza et al., 2021; Fig. 6b). Although based on a limited dataset,
these results give confidence that embayments, when combined with
appropriate diffusion models, can provide valuable and complementary
insights into the rates of magma ascent for historic, well-monitored and
instrumented eruptions.

4. Broad implications

Our combined numerical and experimental results suggest that the
melt embayment can record, and models can be developed to recover,
multi-stage decompression histories that provide crucial insights into

Fig. 5. EDiTS applied to embayments (a) MSH1–6, (b) MSH1–3, and (c) KVC518b of Humphreys et al. (2008). Using the same high melt inclusion starting condition
(6.4 wt%, gray horizontal bar), the two-stage model recovers an initial stage of decompression spanning 4.78–6.54 h (light arrow and horizontal bar), and a final
stage <1–5 min (dark arrow and curve), whereas the original estimate was 2.4–3 min of total ascent time. The cartoon to the right of each panel (a-c) shows the
plagioclase crystal containing the analyzed embayment, along with the annotated grayscale transect (white lines). All grayscale transect data plotted in (a-c) are
sourced from Humphreys et al. (2008). (d) Pressure-time plot comparing results from the single-stage model of Humphreys et al. (2008)† (1.3–1.6 MPa s− 1) and
EDiTS applied to the same suite of Mount St. Helens embayments. EDiTS recovers consistent initial ascent rates ranging from 0.31 to 0.39 m s− 1, with final ascent
rates up to two orders of magnitude faster (9–30 m s− 1). Horizontal gray bars represent the range of Pinitial, Pintermediate, and Pfinal used for modeling MSH em-
bayments. Pintermediate is a free parameter in the model, and should be controlled dynamically by where H2O begins to rapidly exsolve. Model input parameters are
shown in the lower right of the plot, including a starting pressure of 230 MPa, final pressures between 9 and 33 MPa, an isothermal temperature of 880 ◦C, initial H2O
of 6.4 wt%, and melt density of 2500 kg m− 3.

B. Hosseini and M. Myers Journal of Volcanology and Geothermal Research 456 (2024) 108211 

10 



how quickly magma ascends to the surface during explosive volcanic
eruptions. Applying EDiTS to H2O + CO2 profiles in embayments from
both pre-historic and modern eruptions provides consistent results that
magma ascent from the storage region and through the deeper conduit
can be quite slow, with the shallower conduit characterized by faster
ascent times on the order of an hour to minutes. For the May 18th, 1980

eruption of Mount St. Helens, we find that EDiTS more effectively re-
solves the time of initial decompression in the deep conduit, which pairs
well with independent estimates of magma ascent (e.g., direct obser-
vations and seismicity). Final decompression in the shallow conduit
occurs on the order of minutes and is in close agreement with time-
integrated bubble number density decompression estimates, as was
similarly shown for the 1991 eruption of Mount Pinatubo (Hajimirza
et al., 2021; Harris et al., 2024). We argue that the embayment is the
only petrologic tool from which we are currently able to recover mul-
tiple decompression rates representative of a single magma’s decom-
pression path along different regions of the conduit. Previous studies
have presented evidence that initial magma ascent can be aseismic (i.e.,
too slow to produce detectable seismicity), as was the case for the 2009
eruption of Redoubt Volcano (Roman and Cashman, 2018). Our study
highlights that melt embayment geospeedometry, while largely a
retrospective approach to understanding magma ascent rates, can
potentially identify systems that experience slow magmatic ascent
leading to staging in the crust which may not be detectable in all cases
via seismic monitoring. Where in some systems this will allow for the
merging of monitoring and petrologic approaches for refining models of
magma ascent, in systems that lack monitoring datasets (i.e., pre-
historic eruptions), the embayment could provide important insights
into the behavior of magma ascent to inform on future monitoring
efforts.

Lastly, we highlight some outstanding challenges for the embayment
modeling community. Most notably, the non-uniqueness of best-fit so-
lutions in H2O-saturated systems, especially for approximating the
initial decompression rate, begs further investigation into whether
modeling a single diffusive species is a reliable approach to recovering
magma decompression rates, especially in the absence of any indepen-
dent decompression rate estimates. Secondly, the ability of EDiTS to
successfully model both synthetic and experimental embayment volatile
profiles produced by constant-rate decompression (i.e., a false positive
result) further highlights the challenges of embayment modeling and
interpretation. Again, however, we argue that because magma ascent is
dynamically a nonlinear process, this result should not impact the
application of EDiTS to natural systems. Finally, embayment modeling
assumptions (e.g., degassing mechanism; deGraffenried and Shea, 2021)
shown to have significant influence on calculated decompression rates
require further investigation and incorporation into future iterations of
embayment diffusion models. Numerical conduit models that consider
the influence of disequilibrium degassing of H2O and CO2 on predicted
decompression rates have been recently developed for basaltic magmas
(e.g., La Spina et al., 2016, 2021, 2022; Bamber et al., 2022, 2024). The
embayment modeling community would similarly benefit from the
incorporation of these complexities in future iterations of embayment
models for all melt compositions.

5. Conclusions

In this study, we have developed, validated, and applied a new two-
stage decompression-diffusion model for melt embayments: Embayment
Decompression in Two Stages (EDiTS). After successfully numerically
benchmarking EDiTS, we used high pressure-temperature decompres-
sion experiments on natural embayments to determine how faithfully
they record, and how reliably EDiTS recovers, the imposed decom-
pression paths. While EDiTS performs well in mixed-volatile (H2O +

CO2) systems, we note challenges in applying EDiTS to (1) H2O-satu-
rated systems and (2) embayment volatile profiles generated both syn-
thetically and experimentally by a constant-rate decompression (i.e.,
false positive results). To address (1), we suggest using the interior
embayment saturation pressures to constrain the range of model pa-
rameters cycled through. For (2), we argue that, because magma
dynamically does not experience constant, linear decompression, these
false positive results do not necessarily have any bearing on the appli-
cation of EDiTS to natural systems.

Fig. 6. Compilation of magma decompression rates for the May 18th, 1980
Plinian eruption of Mount St. Helens, based on petrological, observational, and
numerical approaches. (a) In the “Old View,” magma decompression rate es-
timates span five orders of magnitude between different petrologic recorders,
with no overlap between any approach. Further, we find no agreement between
petrologic estimates and other independent, direct observations. (b) Notably
EDiTS brings embayment-based ascent time estimates (4.8–6.5 h) into closer
agreement with observational data, which indicate the onset of the second
Plinian explosive pulse occurred 3.5 h following the unroofing of the system
due to an earthquake-triggered landslide and blast. Additionally, application of
new models incorporating time-integrated, heterogeneous bubble nucleation
into models of vesicle number densities lower original estimates by two orders
of magnitude, in good agreement with the final rapid stage of decompression
recovered from embayments using our two-stage model (Hajimirza et al.,
2021). Data sources are as follows: [1] Toramaru (2006), [2] Carey and
Sigurdsson (1985), [3] Humphreys et al. (2008), [4] Endo et al. (1981), [5]
Scandone and Malone (1985), [6] Eichelberger (1995), [7] Rutherford and Hill
(1993), [8] Hajimirza et al. (2021) and [9] Andrews and Befus (2020).

B. Hosseini and M. Myers Journal of Volcanology and Geothermal Research 456 (2024) 108211 

11 



Through application of EDiTS to existing embayment datasets, we
find that initial decompression in the deep conduit can occur over hours,
while the shallow conduit is characterized by much faster decompres-
sion on the order of minutes. For the May 18th, 1980 Plinian eruption of
Mount St. Helens in particular, EDiTS is able to reconcile embayment-
based estimates of magma decompression with both real-time observa-
tions and independent petrologic recorders. Altogether, our work em-
phasizes that embayments can preserve nuanced records of magma
ascent and should be leveraged in actively monitored settings for an
integrated understanding of magma movement from the deep conduit to
the surface. This is especially true for systems containing multiple
diffusing species (e.g., H2O, CO2, S). We thus recommend the develop-
ment of similar multi-stage decompression models applicable to and
parameterized for more intermediate to basaltic melt compositions. As
these are also more frequently active systems, the volcano community
would be positioned to pair embayment decompression rates with
robust, independent observations of magma ascent and eruption.

CRediT authorship contribution statement

Behnaz Hosseini:Writing – original draft, Visualization, Validation,
Software, Project administration, Methodology, Investigation, Funding
acquisition, Formal analysis, Data curation, Conceptualization. Madi-
son Myers: Writing – review & editing, Supervision, Software, Re-
sources, Project administration, Methodology, Investigation, Funding
acquisition, Conceptualization.

Declaration of competing interest

The authors declare that they have no known competing financial
interests or personal relationships that could have appeared to influence
the work reported in this paper.

Data availability

We have included GitHub and Zenodo links for the EDiTS model
presented in the paper: https://zenodo.org/record/10042159. We also
have included experimental and analytical data in the supplement.

Acknowledgements

We acknowledge and thank Jim Watkins for assistance with con-
ducting experiments at the University of Oregon and for fruitful dis-
cussions in the early stages of model development. The authors also
thank Fabio Arzilli and Rebecca deGraffenried for thorough and
constructive reviews that greatly improved this publication, and Ed
Llewellin for editorial handling. This work was supported by National
Science Foundation EAR Award 1922513 to M.M., and Mineralogical
Society of America‘s Grant for Student Research in Mineralogy and
Petrology and the Geological Society of America‘s Graduate Student
Research Grant to B.H.

Appendix A. Supplementary data

All data associated with this study are presented in this paper or
otherwise available in the supplemental document. Appendix A. Sup-
plementary Material contains detailed experimental and analytical
procedures and results, as well as schematics of the EDiTS model ar-
chitecture, numerical benchmarking results, and the experimental
capsule set-up. Supplementary data to this article can be found online at
[https://doi.org/10.1016/j.jvolgeores.2024.108211].

References

Andrews, B.J., Befus, K.S., 2020. Supersaturation Nucleation and Growth of Plagioclase:
A numerical model of decompression-induced crystallization. Contrib. Mineral.
Petrol. 175 (23), 1–20. https://doi.org/10.1007/s00410-020-1660-9.

Bamber, E.C., La Spina, G., Arzilli, F., de ’Michieli Vitturi, M., Polacci, M., Hartley, M.E.,
Petrelli, M., Fellowes, J., Burton, M., 2022. Basaltic Plinian eruptions at Las Sierras-
Masaya volcano driven by cool storage of crystal-rich magmas. Commun. Earth
Environ. 3 (1), 1–17. https://doi.org/10.1038/s43247-022-00585-5.

Bamber, E.C., La Spina, G., Arzilli, F., Polacci, M., Mancini, L., de ’Michieli Vitturi, M.,
Andronico, D., Corsaro, R.A., Burton, M.R., 2024. Outgassing behaviour during
highly explosive basaltic eruptions. Commun. Earth Environ. 5 (1), 1–16. https://
doi.org/10.1038/s43247-023-01182-w.

Befus, K.S., Ruefer, A.C., Allison, C.M., Thompson, J.O., 2022. Quartz-hosted inclusions
and embayments reveal storage, fluxing, and ascent of the Mesa Falls Tuff,
Yellowstone. Earth Planet. Sci. Lett. 601, 1–13. https://doi.org/10.1016/j.
epsl.2022.117909.

Blundy, J., Cashman, K., 2005. Rapid decompression-driven crystallization recorded by
melt inclusions from Mount St. Helens volcano. Geology 33 (10), 793–796. https://
doi.org/10.1130/G21668.1.

Blundy, J., Cashman, K., Humphreys, M., 2006. Magma heating by decompression-driven
crystallization beneath andesite volcanoes. Nature 443 (7107), 76–80. https://doi.
org/10.1038/nature05100.

Carey, S., Sigurdsson, H., 1985. The May 18, 1980 eruption of Mount St. Helens: 2.
Modeling of dynamics of the Plinian phase. J. Geophys. Res.: Solid Earth 90 (B4),
2948–2958. https://doi.org/10.1029/JB090iB04p02948.

Carey, S., Sigurdsson, H., Gardner, J.E., Criswell, W., 1990. Variations in column height
and magma discharge during the May 18, 1980 eruption of Mount St. Helens.
J. Volcanol. Geothermal Res. 43 (1–4), 99–112. https://doi.org/10.1016/0377-0273
(90)90047-J.

Caricchi, L., Townsend, M., Rivalta, E., Namiki, A., 2021. The build-up and triggers of
volcanic eruptions. Nat. Rev. Earth Environ. 2 (7), 458–476. https://doi.org/
10.1038/s43017-021-00174-8.

Cichy, S.B., Botcharnikov, R.E., Holtz, F., Behrens, H., 2011. Vesiculation and microlite
crystallization induced by decompression: A case study of the 1991-1995 Mt Unzen
eruption (Japan). J. Petrol. 52 (7–8), 1469–1492. https://doi.org/10.1093/
petrology/egq072.

deGraffenried, R., Shea, T., 2021. Using volatile element concentration profiles in
crystal-hosted melt embayments to estimate magma decompression rate:
Assumptions and inherited errors. Geochem. Geophys. Geosyst. 22, 1–18. https://
doi.org/10.1029/2021GC009672.

Eichelberger, J.C., 1995. Silicic volcanism: Ascent of viscous magmas from crustal
reservoirs. Ann. Rev. Earth Planet. Sci. 23, 41–63. https://doi.org/10.1146/annurev.
ea.23.050195.000353.

Eichelberger, J.C., Hayes, D.B., 1982. Magmatic model for the Mount St. Helens blast of
May 18, 1980. J. Geophys. Res. Solid Earth 87 (B9), 7727–7738. https://doi.org/
10.1029/JB087iB09p07727.

Elms, H.C., Myers, M.L., Nichols, A.R.L., Wallace, P.J., Wilson, C.J.N., Barker, S.J.,
Charlier, B.L.A., 2023. Pre-eruptive rhyolite magma ascent rate is rapid and
independent of eruption size: A case study from Ōkataina Volcanic Centre, Aotearoa
New Zealand. Bullet. Volcanol. 85 (4), 1–20. https://doi.org/10.1007/s00445-023-
01630-7.

Endo, E.T., Malone, S.D., Noson, L.L., 1981. The 1980 eruptions of Mount St. Helens,
Washington. US Geological Survey Professional Paper 1250, 93–107.

Ferguson, D.J., Gonnermann, H.M., Ruprecht, P., Plank, T., Hauri, E.H., Houghton, B.F.,
Swanson, D.A., 2016. Magma decompression rates during explosive eruptions of
Kı̄lauea volcano, Hawaii, recorded by melt embayments. Bull. Volcanol. 78 (10),
1–12. https://doi.org/10.1007/s00445-016-1064-x.

Freda, C., Baker, D.R., Romano, C., Scarlato, P., 2003. Water diffusion in natural potassic
melts. Geol. Soc. Lond. Spec. Publ. 213 (1), 53–62. https://doi.org/10.1144/GSL.
SP.2003.213.01.04.

Freda, C., Baker, D.R., Scarlato, P., 2005. Sulfur diffusion in basaltic melts. Geochim.
Cosmochim. Acta 69 (21), 5061–5069. https://doi.org/10.1016/j.gca.2005.02.002.

Georgeais, G., Koga, K.T., Moussallam, Y., Rose-Koga, E.F., 2021. Magma decompression
rate calculations with EMBER: A user-friendly software to model diffusion of H2O,
CO2, and S in melt embayments. Geochem. Geophys. Geosyst. 22 (7), 1–20. https://
doi.org/10.1029/2020GC009542.

Geshi, N., Yamasaki, T., Miyagi, I., Conway, C.E., 2021. Magma chamber decompression
during explosive caldera-forming eruption of Aira caldera. Commun. Earth Environ.
2 (1), 1–10. https://doi.org/10.1038/s43247-021-00272-x.

Ghiorso, M.S., Gualda, G.A., 2015. An H2O–CO2 mixed fluid saturation model compatible
with rhyolite-MELTS. Contrib. Mineral. Petrol. 169, 1–30. https://doi.org/10.1007/
s00410-015-1141-8.

Gonnermann, H.M., Manga, M., 2013. Dynamics of Magma Ascent in the Volcanic
Conduit. Modeling volcanic processes: The Physics and Mathematics of Volcanism,
p. 55.

Hajimirza, S., Gonnermann, H.M., Gardner, J.E., 2021. Reconciling bubble nucleation in
explosive eruptions with geospeedometers. Nat. Commun. 12 (1), 1–8. https://doi.
org/10.1038/s41467-020-20541-1.

Hammer, J.E., Rutherford, M.J., 2002. An experimental study of the kinetics of
decompression-induced crystallization in silicic melt. J. Geophys. Res. Solid Earth
107 (B1), ECV–8. https://doi.org/10.1029/2001JB000281.

Harris, M., Hosseini, B., Myers, M., Bouley, L., 2024. Reconciling petrologic magma
ascent speedometers for the June 12th, 1991 eruption of Mount Pinatubo. Volcanica.
Philippines 7 (1), 117–133. https://doi.org/10.30909/vol.07.01.117133.

B. Hosseini and M. Myers Journal of Volcanology and Geothermal Research 456 (2024) 108211 

12 

https://zenodo.org/record/10042159
https://doi.org/10.1016/j.jvolgeores.2024.108211
https://doi.org/10.1007/s00410-020-1660-9
https://doi.org/10.1038/s43247-022-00585-5
https://doi.org/10.1038/s43247-023-01182-w
https://doi.org/10.1038/s43247-023-01182-w
https://doi.org/10.1016/j.epsl.2022.117909
https://doi.org/10.1016/j.epsl.2022.117909
https://doi.org/10.1130/G21668.1
https://doi.org/10.1130/G21668.1
https://doi.org/10.1038/nature05100
https://doi.org/10.1038/nature05100
https://doi.org/10.1029/JB090iB04p02948
https://doi.org/10.1016/0377-0273(90)90047-J
https://doi.org/10.1016/0377-0273(90)90047-J
https://doi.org/10.1038/s43017-021-00174-8
https://doi.org/10.1038/s43017-021-00174-8
https://doi.org/10.1093/petrology/egq072
https://doi.org/10.1093/petrology/egq072
https://doi.org/10.1029/2021GC009672
https://doi.org/10.1029/2021GC009672
https://doi.org/10.1146/annurev.ea.23.050195.000353
https://doi.org/10.1146/annurev.ea.23.050195.000353
https://doi.org/10.1029/JB087iB09p07727
https://doi.org/10.1029/JB087iB09p07727
https://doi.org/10.1007/s00445-023-01630-7
https://doi.org/10.1007/s00445-023-01630-7
http://refhub.elsevier.com/S0377-0273(24)00203-8/rf0075
http://refhub.elsevier.com/S0377-0273(24)00203-8/rf0075
https://doi.org/10.1007/s00445-016-1064-x
https://doi.org/10.1144/GSL.SP.2003.213.01.04
https://doi.org/10.1144/GSL.SP.2003.213.01.04
https://doi.org/10.1016/j.gca.2005.02.002
https://doi.org/10.1029/2020GC009542
https://doi.org/10.1029/2020GC009542
https://doi.org/10.1038/s43247-021-00272-x
https://doi.org/10.1007/s00410-015-1141-8
https://doi.org/10.1007/s00410-015-1141-8
http://refhub.elsevier.com/S0377-0273(24)00203-8/rf0110
http://refhub.elsevier.com/S0377-0273(24)00203-8/rf0110
http://refhub.elsevier.com/S0377-0273(24)00203-8/rf0110
https://doi.org/10.1038/s41467-020-20541-1
https://doi.org/10.1038/s41467-020-20541-1
https://doi.org/10.1029/2001JB000281
https://doi.org/10.30909/vol.07.01.117133


Hosseini, B., Myers, M.L., Watkins, J.M., Harris, M.A., 2023. Are we recording? Putting
embayment speedometry to the test using high pressure-temperature decompression
experiments. Geochem. Geophys. Geosyst. 24 (6), 1–18. https://doi.org/10.1029/
2022GC010770.

Humphreys, M.C., Menand, T., Blundy, J.D., Klimm, K., 2008. Magma ascent rates in
explosive eruptions: Constraints from H2O diffusion in melt inclusions. Earth Planet.
Sci. Lett. 270 (1–2), 25–40. https://doi.org/10.1016/j.epsl.2008.02.041.

Iacovino, K., Matthews, S., Wieser, P.E., Moore, G.M., Bégué, F., 2021. VESIcal Part I: An
open-source thermodynamic model engine for mixed volatile (H2O-CO2) solubility in
silicate melts. Earth Space Sci 8 (11), 1–55. https://doi.org/10.1029/
2020EA001584.

Kieffer, S.W., 1981. Blast dynamics at Mount St Helens on 18 May 1980. Nature 291
(5816), 568–570. https://doi.org/10.1038/291568a0.

La Spina, G., Burton, M., de ’Michieli Vitturi, M., Arzilli, F., 2016. Role of syneruptive
plagioclase disequilibrium crystallization in basaltic magma ascent dynamics. Nat.
Commun. 7 (1), 1− 10. https://doi.org/10.1038/ncomms13402.

La Spina, G., Arzilli, F., Llewellin, E.W., Burton, M.R., Clarke, A.B., Vitturi, M.D.M.,
Polacci, M., Hartley, M.E., Di Genova, D., Mader, H.M., 2021. Explosivity of basaltic
lava fountains is controlled by magma rheology, ascent rate and outgassing. Earth
Planet. Sci. Lett. 553 (1), 1–11. https://doi.org/10.1016/j.epsl.2020.116658.

La Spina, G., Arzilli, F., Burton, M.R., Polacci, M., Clarke, A.B., 2022. Role of volatiles in
highly explosive basaltic eruptions. Commun. Earth Environ. 3 (1), 1–13. https://
doi.org/10.1038/s43247-022-00479-6.

Liu, Y., Zhang, Y., Behrens, H., 2005. Solubility of H2O in rhyolitic melts at low pressures
and a new empirical model for mixed H2O–CO2 solubility in rhyolitic melts.
J. Volcanol. Geotherm. Res. 143 (1–3), 219–235. https://doi.org/10.1016/j.
jvolgeores.2004.09.019.

Liu, Y., Anderson, A.T., Wilson, C.J., 2007. Melt pockets in phenocrysts and
decompression rates of silicic magmas before fragmentation. J. Geophys. Res.: Solid
Earth 112,1− 12. doi:10.1029/2006JB004500.

Lloyd, A.S., Ruprecht, P., Hauri, E.H., Rose, W., Gonnermann, H.M., Plank, T., 2014.
NanoSIMS results from olivine-hosted melt embayments: Magma ascent rate during
explosive basaltic eruptions. J. Volcanol. Geotherm. Res. 283, 1–18. https://doi.org/
10.1016/j.jvolgeores.2014.06.002.

Mangan, M., Sisson, T., 2000. Delayed, disequilibrium degassing in rhyolite magma:
Decompression experiments and implications for explosive volcanism. Earth Planet.
Sci. Lett. 183 (3–4), 441–455. https://doi.org/10.1016/S0012-821X(00)00299-5.

Mastin, L.G., 2005. The controlling effect of viscous dissipation on magma flow in silicic
conduits. J. Volcanol. Geotherm. Res. 143, 17–28. https://doi.org/10.1016/j.
jvolgeores.2004.09.008.

Mastin, L.G., Ghiorso, M.S., 2001. Adiabatic temperature changes of magma–gas
mixtures during ascent and eruption. Contrib. Mineral. Petrol. 141, 307–321.
https://doi.org/10.1007/s004100000210.

Melnik, O., Sparks, S., 2006. Transient Models of Conduit Flows During Volcanic
Eruptions. Geological Society of London. https://doi.org/10.1144/IAVCEI001.16.

Moussallam, Y., Rose-Koga, E.F., Koga, K.T., Médard, E., Bani, P., Devidal, J.-L., Tari, D.,
2019. Fast ascent rate during the 2017–2018 Plinian eruption of Ambae (Aoba)
volcano: A petrological investigation. Contrib. Mineral. Petrol. 174 (11), 1–24.
https://doi.org/10.1007/s00410-019-1625-z.

Myers, M.L., Wallace, P.J., Wilson, C.J., Watkins, J.M., Liu, Y., 2018. Ascent rates of
rhyolitic magma at the onset of three caldera-forming eruptions. Am. Mineral. 103
(6), 952–965. https://doi.org/10.2138/am-2018-6225.

Myers, M.L., Druitt, T.H., Schiavi, F., Gurioli, L., Flaherty, T., 2021. Evolution of magma
decompression and discharge during a Plinian event (late Bronze-Age eruption,
Santorini) from multiple eruption-intensity proxies. Bull. Volcanol. 83 (3), 1–17.
https://doi.org/10.1007/s00445-021-01438-3.

Newcombe, M.E., Plank, T., Barth, A., Asimow, P.D., Hauri, E., 2020. Water-in-olivine
magma ascent chronometry: Every crystal is a clock. J. Volcanol. Geotherm. Res.
398, 1–17. https://doi.org/10.1016/j.jvolgeores.2020.106872.

Newman, S., Lowenstern, J.B., 2002. VolatileCalc: A silicate melt–H2O–CO2 solution
model written in Visual basic for excel. Comput. Geosci. 28 (5), 597–604. https://
doi.org/10.1016/S0098-3004(01)00081-4.

Ni, H., Zhang, Y., 2008. H2O diffusion models in rhyolitic melt with new high pressure
data. Chem. Geol. 250 (1–4), 68–78. https://doi.org/10.1016/j.
chemgeo.2008.02.011.

Nowak, M., Behrens, H., 1997. An experimental investigation on diffusion of water in
haplogranitic melts. Contrib. Mineral. Petrol. 126 (4), 365–376. https://doi.org/
10.1007/s004100050256.

Nowak, M., Schreen, D., Spickenbom, K., 2004. Argon and CO2 on the race track in
silicate melts: A tool for the development of a CO2 speciation and diffusion model.
Geochim. Cosmochim. Acta 68 (24), 5127–5138. https://doi.org/10.1016/j.
gca.2004.06.002.

Rasmussen, D.J., Plank, T.A., Roman, D.C., Power, J.A., Bodnar, R.J., 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. https://doi.org/10.1016/j.
epsl.2018.01.001.

Roman, D.C., Cashman, K.V., 2018. Top–down precursory volcanic seismicity:
Implications for ‘stealth’ magma ascent and long-term eruption forecasting. Front.
Earth Sci. 6, 124. https://doi.org/10.3389/feart.2018.00124.

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

Rutherford, M.J., Hill, P.M., 1993. Magma ascent rates from amphibole breakdown: An
experimental study applied to the 1980− 1986 Mount St. Helens eruptions.
J. Geophys. Res. Solid Earth 98 (B11), 19667–19685. https://doi.org/10.1029/
93JB01613.

Saalfeld, M.A., Myers, M., deGraffenried, R., Shea, T., Waelkens, C., 2022. On the rise:
Using reentrants to extract magma ascent rates in the Bandelier Tuff caldera
complex, New Mexico, USA. Bull. Volcanol. 84, 1–21. https://doi.org/10.1007/
s00445-021-01518-4.

Scandone, R., Malone, S.D., 1985. Magma supply, magma discharge and readjustment of
the feeding system of Mount St. Helens during 1980. J. Volcanol. Geotherm. Res. 23
(3–4), 239–262. https://doi.org/10.1016/0377-0273(85)90036-8.

Sparks, R.S.J., 2003. Dynamics of Magma Degassing doi.org/10.1144/GSL.
SP.2003.213.01.02.

Toramaru, A., 2006. BND (bubble number density) decompression rate meter for
explosive volcanic eruptions. J. Volcanol. Geotherm. Res. 154 (3–4), 303–316.
https://doi.org/10.1016/j.jvolgeores.2006.03.027.

Voight, T., Blacken, H., Jana, R.J., Douglass, P.M., 1981. The 1980 eruptions of Mount St.
Helens, Washington. US Geological Survey Professional Paper 1250, 347–377.

Wei, Z., Ruefer, A.C., Pamukcu, A.S., Suckale, J., 2024. Deciphering clues regarding
magma composition encoded in quartz-hosted embayments and melt inclusions
through direct numerical simulations. J. Geophys. Res. Solid Earth 129 (4), 1–25.
https://doi.org/10.1029/2023JB028080.

Witham, F., Blundy, J., Kohn, S.C., Lesne, P., Dixon, J., Churakov, S.V., Botcharnikov, R.,
2012. SolEx: A model for mixed COHSCl-volatile solubilities and exsolved gas
compositions in basalt. Comput. Geosci. 45, 87–97. https://doi.org/10.1016/j.
cageo.2011.09.021.

Wong, Y.Q., Segall, P., 2019. Numerical analysis of time-dependent conduit magma flow
in dome-forming eruptions with application to Mount St. Helens 2004–2008.
J. Geophys. Res. Solid Earth 124 (11), 11251–11273. https://doi.org/10.1029/
2019JB017585.

Zhang, Y., Behrens, H., 2000. H2O diffusion in rhyolitic melts and glasses. Chem. Geol.
169 (1–2), 243–262. https://doi.org/10.1016/S0009-2541(99)00231-4.

Zhang, Y., Ni, H., 2010. Diffusion of H, C, and O components in silicate melts. Rev.
Mineral. Geochem. 72 (1), 171–225. https://doi.org/10.2138/rmg.2010.72.5.

Zhang, Y., Xu, Z., Zhu, M., Wang, H., 2007. Silicate melt properties and volcanic
eruptions. Rev. Geophys. 45 (4), 1–27. https://doi.org/10.1029/2006RG000216.

Zhang, Y., Ni, H., Chen, Y., 2010. Diffusion data in silicate melts. Rev. Mineral. Geochem.
72 (1), 311–408. https://doi.org/10.2138/rmg.2010.72.8.

Zuccarello, F., Schiavi, F., Viccaro, M., 2022. The eruption run-up at Mt. Etna volcano:
Constraining magma decompression rates and their relationships with the final
eruptive energy. Earth Planet. Sci. Lett. 597, 1–13. https://doi.org/10.1016/j.
epsl.2022.117821.

B. Hosseini and M. Myers Journal of Volcanology and Geothermal Research 456 (2024) 108211 

13 

https://doi.org/10.1029/2022GC010770
https://doi.org/10.1029/2022GC010770
https://doi.org/10.1016/j.epsl.2008.02.041
https://doi.org/10.1029/2020EA001584
https://doi.org/10.1029/2020EA001584
https://doi.org/10.1038/291568a0
https://doi.org/10.1038/ncomms13402
https://doi.org/10.1016/j.epsl.2020.116658
https://doi.org/10.1038/s43247-022-00479-6
https://doi.org/10.1038/s43247-022-00479-6
https://doi.org/10.1016/j.jvolgeores.2004.09.019
https://doi.org/10.1016/j.jvolgeores.2004.09.019
https://doi.org/10.1016/j.jvolgeores.2014.06.002
https://doi.org/10.1016/j.jvolgeores.2014.06.002
https://doi.org/10.1016/S0012-821X(00)00299-5
https://doi.org/10.1016/j.jvolgeores.2004.09.008
https://doi.org/10.1016/j.jvolgeores.2004.09.008
https://doi.org/10.1007/s004100000210
https://doi.org/10.1144/IAVCEI001.16
https://doi.org/10.1007/s00410-019-1625-z
https://doi.org/10.2138/am-2018-6225
https://doi.org/10.1007/s00445-021-01438-3
https://doi.org/10.1016/j.jvolgeores.2020.106872
https://doi.org/10.1016/S0098-3004(01)00081-4
https://doi.org/10.1016/S0098-3004(01)00081-4
https://doi.org/10.1016/j.chemgeo.2008.02.011
https://doi.org/10.1016/j.chemgeo.2008.02.011
https://doi.org/10.1007/s004100050256
https://doi.org/10.1007/s004100050256
https://doi.org/10.1016/j.gca.2004.06.002
https://doi.org/10.1016/j.gca.2004.06.002
https://doi.org/10.1016/j.epsl.2018.01.001
https://doi.org/10.1016/j.epsl.2018.01.001
https://doi.org/10.3389/feart.2018.00124
https://doi.org/10.1038/s41467-018-05086-8
https://doi.org/10.1038/s41467-018-05086-8
https://doi.org/10.1029/93JB01613
https://doi.org/10.1029/93JB01613
https://doi.org/10.1007/s00445-021-01518-4
https://doi.org/10.1007/s00445-021-01518-4
https://doi.org/10.1016/0377-0273(85)90036-8
https://doi.org/10.1016/j.jvolgeores.2006.03.027
http://refhub.elsevier.com/S0377-0273(24)00203-8/rf0275
http://refhub.elsevier.com/S0377-0273(24)00203-8/rf0275
https://doi.org/10.1029/2023JB028080
https://doi.org/10.1016/j.cageo.2011.09.021
https://doi.org/10.1016/j.cageo.2011.09.021
https://doi.org/10.1029/2019JB017585
https://doi.org/10.1029/2019JB017585
https://doi.org/10.1016/S0009-2541(99)00231-4
https://doi.org/10.2138/rmg.2010.72.5
https://doi.org/10.1029/2006RG000216
https://doi.org/10.2138/rmg.2010.72.8
https://doi.org/10.1016/j.epsl.2022.117821
https://doi.org/10.1016/j.epsl.2022.117821

	Melt embayments record multi-stage magma decompression histories
	1 Introduction
	2 Model set-up and validation
	2.1 Numerical model development and benchmarking
	2.2 High pressure-temperature decompression experiments

	3 Model application
	3.1 Revisiting existing embayment studies
	3.1.1 Applications to pre-historic eruptions
	3.1.2 Applications to historic instrumented eruptions


	4 Broad implications
	5 Conclusions
	CRediT authorship contribution statement
	Declaration of competing interest
	Data availability
	Acknowledgements
	Appendix A Supplementary data
	References