# Magma decompression rate calculations with EMBER: A user-friendly software to model diffusion of H 2 O, CO 2 and S in melt embayments

Guillaume Georgeais, Kenneth T. Koga, Yves Moussallam, Estelle F. Rose-Koga

► To cite this version:

Guillaume Georgeais, Kenneth T. Koga, Yves Moussallam, Estelle F. Rose-Koga. Magma decompression rate calculations with EMBER: A user-friendly software to model diffusion of H 2 O, CO 2 and S in melt embayments. Geochemistry, Geophysics, Geosystems, 2021, 10.1029/2020GC009542. hal-03265994

HAL Id: hal-03265994

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

Submitted on 21 Jun 2021

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

# Magma decompression rate calculations with EMBER: A user-friendly software to model diffusion of H₂O, CO₂ and S in melt embayments

Guillaume Georgeais¹, Kenneth T. Koga¹, Yves Moussallam²,³, Estelle F. Rose-Koga¹

¹Université Clermont Auvergne, CNRS, IRD, OPGC, Laboratoire Magmas et Volcans, F-63000 Clermont-Ferrand, France

²Lamont-Doherty Earth Observatory, Columbia University, New York, NY 10027, USA

³Department of Earth and Planetary Sciences, American Museum of Natural History, New York, NY 10024, USA

Corresponding author: Guillaume Georgeais; [guillaume.georgeais@uca.fr](mailto:guillaume.georgeais@uca.fr)

**Keywords:** volcanic eruption, magma degassing, conduit processes, volatile element exsolution, GUI

## Highlights:

* EMBER is a new, freely available, GUI software that models magma ascent rates for basalt to rhyolite with H₂O, CO₂ and S diffusion profiles

* EMBER was validated by reproducing previous published literature data

* For mafic eruptions, we found a notable correlation between maximum recalculated decompression rates and eruption magnitude or plume height

19

20

21

22

23

24
ABSTRACT

Magma decompression rate is one of the most important parameters in controlling eruption dynamics. One way to determine decompression rate is by fitting a volatile elements diffusion profile to a concentration gradient in crystal-hosted embayments. Previous studies have used a variety of diffusion models, limiting the possibility for inter-study comparison. Here, we introduce EMBER (EMBayment-Estimated Rates), a standalone versatile tool that models diffusion of volatile elements along melt embayments. Our model relies on the *pdepe* function of MATLAB to calculate diffusion profiles of H₂O, CO₂ and S through the finite difference method. EMBER uses a grid search seeking out the best fits for decompression rates, initial dissolved concentration of each studied volatile and initial exsolved gas content, while setting three constants: temperature along the ascent and pressure at the beginning and end of the ascent. Our model can compute the rate for basaltic, intermediate, and rhyolitic compositions. We applied EMBER to previous studies to evaluate and validate our model.

We then re-processed "homogeneously" the raw data from the literature for a comparison. In other words, the same protocol was used for each diffusion profiles removing the literature-specific strategies used to constrain unknown parameters. With this comparison, we found a statistically significant positive correlation between maximum magma decompression rates and explosivity of the related eruption. EMBER is expected to help increase the number of volatile diffusion in embayments studies aiming at constraining magma decompression and ascent rates and to facilitate inter-study comparisons.
46

# 1. INTRODUCTION

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

Magma ascent rate is a fundamental physical parameter in determining the behavior of a volcanic eruption. Magma ascent rate can be estimated in several ways, during ascent and eruption by geophysical methods, or after the eruption using geochemical and petrological methods. One example of a geophysical method is the analysis of volcanic earthquakes, in which the progressive migration of tremors from a depth to the surface is interpreted as the upward migration of ascending magma [e.g., Aki and Koyanagi., 1981; Scandone and Malone., 1985; Klein et al., 1987; Tryggvason., 1994]. Geochemical methods to determine magma ascent rate are numerous but most of them are only applicable to limited melt compositions. In andesitic magma for instance, the thickness of amphibole breakdown rims has been used to determine ascent rates based on the experimentally-determined rate of the break-down reaction [e.g., Carey and Sigurdsson., 1985; Rutherford and Hill., 1993; Geschwind and Rutherford., 1995; Browne and Gardner., 2006]. Another example is the use of diffusion profiles (e.g., Fe-Mg or H₂O) at the rims of minerals in mantle xenoliths as they re-equilibrate with their carrier melt on their way to the surface [e.g., Mackwell and Kohlstedt., 1990; Klügel et al., 1997; Klügel., 1998; Kohlstedt and Mackwell., 1998; Le Voyer et al., 2014]. Knowing the diffusion coefficient of the species of interest, one can derive an estimate of the magma ascent rate for these eruptions [e.g., Demouchy and Mackwell., 2006; Sparks et al., 2006; Rutherford., 2008]. Another method, based on nucleation theory, links the bubble number density in erupted products (measurable from 2D and 3D observations) to the magma decompression rate [Toramaru., 1989, 1995, 2006]. This theoretical relationship has been reproduced experimentally but is not applicable to natural samples, once bubbles start to coalesce [e.g., Martel and Iacono-Marziano., 2015].

69

70

One method, which holds the potential to be widely applicable for a wide range of magma compositions and decompression rates, is the analysis of the diffusion profiles of
volatile elements along melt embayments (also called re-entrants). Embayments are crystal-
hosted elongated melt pockets of various shapes and sizes, opened to the outside melt [e.g.,
Anderson., 1991]. Their formation mechanism is similar to that of melt inclusions [e.g., Faure
and Schiano., 2005], through either crystallization around a defect or dissolution of the host
crystal, with the exception that they remain connected to the surrounding bubbly melt. They
have been studied in quartz [e.g., Liu et al., 2007; Myers et al., 2018], in plagioclases [e.g.,
Humphreys et al., 2008] and olivine crystals [e.g., Lloyd et al., 2014; Ferguson et al., 2016;
Moussallam et al., 2019], in quenched material ranging from basalt to rhyolite in composition.
During magma ascent, the volatile content of the melt surrounding the embayment will
decrease, maintaining equilibrium with the exsolved gas phase. The limited volume of melt
embayments, however, often prevents bubble formation inside, resulting in embayment melt
being super-saturated in volatiles compared to the surrounding melt. This difference in
chemical potential leads to a concentration gradient and to diffusive transport of volatile
species from the interior to the mouth of the embayment resulting in a diffusion profile. If
such a profile is preserved in natural samples, it can then be inverted to derive a
decompression rate assuming that the elongated tubular shape of the embayment and the
relative impermeability of the host crystal led to unidirectional (1D) diffusion and prevented
any advective melt motion in the embayment. The dependency of volatile diffusion in
embayments on multiple parameters (melt composition, temperature, pressure, degassing
path, decompression rate) makes the interpretation of natural diffusion profiles a challenging
endeavour. A numerical model adapted to a wide range of magma compositions is therefore
needed to generate multiple diffusion profiles with known parameters and find the synthetic
profile which most closely reproduces the measurement.
Several such numerical models have been developed in the last decades (Table 1). The first published code was developed in FORTRAN 77 [Liu et al., 2007]. Assuming a certain temperature, pressure, initial concentration, exsolved gas, degassing path and decompression rate, the code generated time-dependent profiles of H₂O and CO₂ concentrations until the fragmentation pressure is reached. The process was repeated several times with different decompression rates to find the best fit. A second model was presented by Humphreys et al., [2008] using the COMSOL multiphysics software to model only H₂O profiles. Their model imposes the final concentration at the mouth of the embayment from the start of the calculation and allows the software to run the diffusion calculation. Lloyd et al., [2014] developed a model calculating simultaneously, for a range of decompression rates, three volatile element profiles at once: H₂O, CO₂ and S. In a following study, an improved development of a MATLAB code by Ferguson et al., [2016] took H₂O, CO₂ and S into account and considered not only a range of decompression rates but also the initial concentration of each volatile element as well as the exsolved gas content at the beginning of the ascent (M₀) for basaltic compositions. The addition of this new parameter: (M₀), the pre-existing (already exsolved volatile content in equilibrium with magma at the onset of a magma ascent), proved to have a significant impact on the modelled profile and made the grid search more complex and the result better constrained [Ferguson et al., 2016]. Another study subsequently built the FORTRAN 77 model from Liu et al., [2007] in MATLAB with updated diffusion coefficients for rhyolitic compositions and a best fit search algorithm [Myers et al., 2018]. This code was later updated taking the pressure at which degassing stops as a free parameter [Myers et al., 2021]. Another code, written in Rstudio and restricted to basaltic melts, took into account all the aforementioned parameters with fixed inputs and was made openly available [Moussallam et al., 2019]. One of the most recent model is tuned to intermediate magma compositions, contains a specific S solubility relation and a general H₂O
120

diffusion coefficient relation for intermediate magma, and looks for the decompression rate
and initial pressure [Newcombe et al., 2020]. The latest model to date is written in Matlab and
uses grid searches to find decompression rate, initial concentration of H₂O and S and exsolved
gas content [Moussallam et al., 2021].

121

122

123

124

At present, seven different data processing methods to model volatile diffusion profiles
in embayments exist; each published study uses its specific model written on four different
platforms (FORTRAN 77, COMSOL, MATLAB and RStudio). Other than the code of
Moussallam et al., [2019], none are directly downloadable without a specific request to the
authors. Each code considers different input parameters, different volatile species and is tuned
to a specific melt composition. This lack of consistency is an issue for inter-study comparison
and the lack of open software access can be an impediment to a large number of new studies
on natural products.

125

126

127

128

129

130

131

132

The aim of this article is to provide the community with a user-friendly and cross-
operating system MATLAB code that is able to constrain decompression rates from volatile
diffusion in melt embayments for rhyolitic to basaltic melt compositions, and for as wide a
range of starting conditions as possible. We then retroactively analyze all volatile diffusion
profiles from the literature using EMBER. Our results help identify potential discrepancies in
published decompression rates, notably for the Mt St Helens 1980's eruption, and provide an
easily comparable, self-consistent summary of decompression rates obtained from volatile
element diffusion in melt embayments published to date. Our software, EMBER, calculates
results likely comparable to those by the DIPRA software [Girona and Costa., 2013], which
can extract timescales from diffusion zoning in olivine crystals, and is also a widely
distributed MATLAB program.

133

134

135

136

137

138

139

140

141

142

143
## 2. Code architecture

### 2.1 Diffusion model

Volatile element diffusion in embayments can be regarded as a 1D process, because of their elongated, tube-like geometries and the incompatibility of the elements in their mineral host [Ferguson et al., 2016; Humphreys et al., 2008; Liu et al., 2007; Lloyd et al., 2014; Moussallam et al., 2019, 2021; Myers et al., 2018, 2021; Newcombe et al., 2020]. While this is a commonly accepted assumption, further studies are needed to assess the impacts of the three-dimensional shape of embayments to volatiles diffusion [deGraffenried and Shea., 2020]. If it proves to be relevant, EMBER will need to be updated accordingly to compute both 1D and 3D diffusion. The evolution of the concentration gradient is therefore described by Fick's Second Law (Eq.1):

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

$$ C = C_i \quad \text{at } x > 0, t = 0 \qquad (2) $$

$$ C = C_{sat}(t) \quad \text{at } x = 0, t > 0; \quad \frac{dc}{dx} = 0 \quad \text{at } x = X, t > 0 \qquad (3) $$

where D is the diffusion coefficient of the studied volatile species; C, the concentration of studied species, x, the distance from the mouth of the embayment, X, the distance at the interior of the embayment, Csat, the saturation concentration calculated at the mouth of the embayment by using a solubility model, Cᵢ, the initial concentration along the embayment and t, the time. These diffusion equations show that diffusion coefficients depend heavily on H₂O concentration(CH₂O), temperature and pressure. These three parameters, and therefore the diffusion coefficient, all change during ascent. Diffusion coefficients are hence calculated along every point of the embayment from the start to the end of the calculation.
166

The boundary condition at the interior is defined by an absence of mass flux (Neumann condition) (Eq.2) and the volatile concentration at the mouth of the embayment is fixed but varies with respect to time t (Dirichlet boundary condition). The initial volatile concentration $C_{init}$ along the embayment is constant, and the value should be the concentrations of volatile elements at the initiation of magma ascent. Volatile concentration at the mouth of the embayment is set by the solubility of each volatile species along a pressure related path of the magma ascent, and these constraints are entered into EMBER as a text file. For example, EMBER's default setup reads output files from SolEx [Witham et al., 2012], and VolatileCalc [Newman and Lowenstern., 2002] depending on magma types.

167

168

169

170

171

172

173

174

175

Solubility of gas species depends on T, P, magma composition and exsolved gas content ($M_0$). EMBER does not calculate the solubility of volatiles and uses an “external” solubility model like VolatileCalc [Newman and Lowenstern., 2002], Solex [Witham et al., 2012] or any other model that the user chooses to calculate the degassing paths. Since the value of $M_0$ is initially unknown, a typical calculation is done by choosing a “target exsolved gas content” thus by setting a corresponding degassing path. EMBER works with seven solubility files accounting for $M_0$ values of 0, 0.1, 0.2, 0.4, 0.8, 1.6, 3.2 wt. %, which are supplied by the user. EMBER interpolates those seven degassing paths to find the degassing path associated with the targeted exsolved gas content. It should be noted that, when VolatileCalc-generated degassing paths are used, EMBER imposes the reference point concentrations of 0.01wt. % H$_2$O, 1 ppm CO$_2$ and 1 ppm S at P=1 bar because VolatileCalc does not model the degassing path down to one bar.

176

177

178

179

180

181

182

183

184

185

186

187

EMBER calculates all model diffusion profiles of the volatile elements following the grid search which range is defined by the user, before iteratively comparing measured and calculated profiles. First, profiles are calculated by varying three parameters for each calculation loop: constant decompression rate, initial volatile content and exsolved gas

188

189

190
content, for each element. Parameters such as temperature, melt composition, initial pressure prior to ascent and pressure at which decompression stops are specified together with the volatile species of interest (H₂O, CO₂, S) and they are fixed for all calculations. The program calculates diffusion profiles using finite difference formulation solved with the *pdepe* function, an ordinary differential equation (ODE) solver from MATLAB.

Diffusion coefficients used in the model calculations are based on the following equations and calculated for each step of the ascent and at each point of a given profile. Water diffusion coefficient ($D_{H_2O}$) in m²·s⁻¹ for basaltic melt is given by

$$D_{H_2O} = \exp[-11.924 - 1.003 \ln(C_{H_2O})] * \exp\left[\frac{-\exp(11.836-0.139 \ln(C_{H_2O}))}{RT}\right] \quad (4)$$

where $C_{H_2O}$ is water concentration in wt. %; R, the gas constant in J·mol⁻¹·K⁻¹ and T, the temperature in Kelvin [Freda et al., 2003]. For melt with intermediate composition (55 < SiO₂ < 70 wt. %), the water diffusion coefficient ($D_{H_2O_t}$ in m²·s⁻¹) is calculated using a combination of equations (5a-i), from Ni and Zhang., [2018], as previously used by Newcombe et al., [2020]:

$$D_{H_2O_t} = D_{H_2O_m} \left(1 - \frac{dX_{OH}}{2dX}\right) + D_{OH} \frac{dX_{OH}}{2dX}, \quad (5a)$$

$$D_{H_2O_m} = D_0 \exp(aX), \quad (5b)$$

$$D_{OH} = \text{const}, \quad (5c)$$

$$\frac{dX_{OH}}{2dX} = \frac{1-2X}{\sqrt{4X(X-1)(1-\frac{4}{K})+1}}, \quad (5d)$$

$$\ln K = 2.6Y_{SI} - \frac{4339Y_{SI}}{T}, \quad (5e)$$

$$a = -94.07 + 74.112Y_{SI} + \frac{198508-166674Y_{SI}}{T}, \quad (5f)$$
211
$$ \ln D_{OH} = -16.78 - 37.428Y_{SI} - \frac{39250-27576Y_{SI}}{T}, \quad (5g) $$

212
$$ \ln\left(\frac{D_{OH}}{D_0}\right) = -56.09 - 115.93Y_{SI} + 160.54\sqrt{Y_{SI}} - \frac{3970\sqrt{Y_{SI}}}{T}, \quad (5h) $$

213
$$ \ln D_0 = 8.02 - 31Y_{SI} + 2.348PY_{SI} + \frac{121824Y_{SI}-118323\sqrt{Y_{SI}}-(10016Y_{SI}-3648)P}{T}, \quad (5i) $$

214 with $Y_{SI}$, the mole fraction of Si among all cations, $D_0$, a diffusion parameter in $m^2 \cdot s^{-1}$, $K$, 215 the equilibrium constant, $X$, the mole fraction of each species, $P$ in GPa, $T$ in Kelvin, $a$, a 216 dimensionless parameter, $D_{OH}$ the OH diffusivity in $m^2 \cdot s^{-1}$, $D_{H_2O_m}$, the molecular H₂O 217 diffusivity in $m^2 \cdot s^{-1}$ and $D_{H_2O_t}$, the total H₂O diffusivity in $m^2 \cdot s^{-1}$. 218

For rhyolitic melt, the water diffusion coefficient (in $\mu m^2 \cdot s^{-1}$) is given by,

219
$$ D_{H_2O} = C_{H_2O} \exp\left(10.49 - \frac{10.661}{T} - \frac{1.772P_t}{T}\right) \quad (6a) $$

220 for low ($\le 2$ wt. %) water melt, where $P_t$ is the pressure in MPa at time $t$ [Zhang and 221 Behrens., 2000]. However, if the water content is high (from 2 to 8 wt. %), diffusion 222 coefficient is calculated by,

223
$$ D_{H_2O} = X \exp(m) \left\{ 1 + \exp \left[ 56 + m + X \left( -34.1 + \frac{44620}{T} + \frac{57.3P_t}{T} \right) - \sqrt{X} \left( 0.091 + \frac{4.77*10^6}{T^2} \right) \right] \right\} \quad (6b) $$

225 where $m = -20.79 - 5030/T - 1.4 P_t/T$, and $X = (C_{H_2O}/18.015)/[C_{H_2O}/18.015 + (100 - C_{H_2O})/32.49]$ [Zhang and Behrens., 2000]. The diffusion coefficient of CO₂ ($D_{CO_2}$), 226 in $m^2 \cdot s^{-1}$, is given by,

228
$$ D_{CO_2} = \exp\left[-14.34 - \frac{17360-0.6527P_t}{T} + \left(-0.7171 + \frac{1436.8}{T}\right) C_{H_2O}\right] \quad (7) $$
for basalts to rhyolites [Nowak et al., 2004; Zhang et al., 2007]. The program also calculates diffusion coefficient of S ($D_s$), in $m^2 \cdot s^{-1}$:

$$D_s(i) = \exp\left(-8.21 - \frac{27692 - 651.6 C_{H_2O}}{T}\right) \quad (8)$$

[Zhang et al., 2007]. It is to be noted that the sulfur diffusivity determined by Zhang et al., [2007] is applicable only to basaltic melts under reduced conditions [Zhang et al., 2007].

## 2.2 Best fit determination

The second part of the code compares a measured concentration profile to a series of calculated ones (with N set of parameters) to find the best fit and its associated parameters (i.e. grid search). EMBER favours the grid search over other more efficient optimization methods (e.g. gradient search method), because with the mixture of diffusion of dissimilar species, it is not clear a priori if there will always exist a unique solution. While computationally intensive (and inefficient), we consider it is more suitable to take an approach of “calculate-all”. The best fit is determined by comparing the dimensionless normalized residual square error ($NErr_j$) for each j of N sets of parameters,

$$NErr_j = \sum_i \frac{1}{u_j^2} (C_{ij} - C_{s_{ij}})^2 \quad (9)$$

where $u_j$ is the analytical measurement uncertainty applied at the highest volatile concentration measured in the embayment; usually at the interior of the embayment (i.e: an uncertainty of 5% of $C_{s_{ij}_max}$), i denotes discrete points along the distance to the embayment mouth, $C_{ij}$ is the calculated concentration and $C_{s_{ij}}$ is the measured concentration. Among the generated profiles, the one that gives the minimum $NErr_j$ is the best fit. One may calculate $NErr_j$ by considering the “weighting and scaling of the error” of each volatile (further details can be found in part 4.5). By doing so, the program first calculates $NErr_j$ values without
dimensions thanks to $u_j$. The term $u_j$ effectively constrains the “weight” of three diffusion profiles (H₂O, CO₂, and S) as measurements with better uncertainty (in relative term) have a stronger influence to the sum total of three $NErr$. EMBER then scales the $NErr$ values between each volatiles before adding them. If one chooses not to weight and scale each volatile species’ error, $u_j = 1$.

EMBER assesses statistical variations of $NErr$ by a Monte Carlo simulation accounting for the uncertainty of the measurement $C \pm \sigma$, using a Gaussian distribution random function to generate m iterations of possible profiles within measurement uncertainties. Therefore, for a parameter set (i.e. case j), there is m number of $NErr$ computed by comparing a calculated profile against m numbers of the randomly generated profiles $C_{s_i} + \varepsilon$ (Figure 1, Step 1). A similar approach is taken for the uncertainty of $x$. The result of this simulation is represented by a mean ($Nerr_m$) and its related 1-$\sigma$ confidence interval for a parameter set j (Figure 1, Step 2).

## 2.3 Propagation of errors

The uncertainties on the fit parameters $C_i$, $dPdt$ and $M_0$ (to a lesser extent) are assessed by the statistical distributions of the parameters that are extracted from the calculated profiles $y$ that satisfy the following condition (Figure 1, Step 3):

$$ \Pr(a - \sigma_1 \le y \le a + \sigma_2) = 68\% \qquad (10) $$

Here, $a$ is the value for which $Nerr_m$ is the smallest, and $\sigma_1$ and $\sigma_2$ respectively the lower and upper boundaries of the confidence interval, calculated using the *prctile* MATLAB function. For every j calculated profiles (with a set of parameters), there is a corresponding $NErr$. Equation (10) states that there are cases of calculated profiles for which $NErr$ is within the expected uncertainty. EMBER extracts the corresponding parameters (dPdt and $C_i$) and
report the distribution of values as mean and 1-σ confidence interval, which are the best fit
parameter and an associated uncertainty (Figure 1. Step 4 and 5). In consequence,
measurement uncertainties must be entered by the user in EMBER because of the Monte
Carlo error propagation. An uncertainty on the distance between two measurement spots is set
by default at ±2 µm and can be changed to whichever value the user requires, down to 0 µm if
needed.

280

281 **3. Examples of model outputs**

282 **3.1 Organization of the graphical user interface**

283 EMBER runs inside a user-friendly Graphical User Interface (GUI) upon execution
284 of the application file. The GUI is separated in two parts: the input section, dedicated to
285 generating the grid search and decompression conditions (left hand side Part1, panels a, b, c;
286 Fig.2), and the results section dedicated to displaying the results once calculation is
287 terminated (Part 2, panels a, b; Fig. 2). Panel 2c display a process log during the calculations.
288 Figure 3 shows the result of a simulation using the volatile profiles measured in an
289 embayment from the study of Ferguson et al., [2016] as an example. Unticking the
290 “Weighting and scaling of error” checkbox will remove the scaling of NErr. The *u*<sub>j</sub> parameter
291 from equation (8) is set to 1 in this case and NErr is not scaled before being added for the
292 cumulative error calculation. In such case, as all concentration are treated as weight percent
293 internally, the H<sub>2</sub>O profile will usually weight more than the others as it usually records much
294 larger variations in absolute concentration.

295

296 **3.2 Outputs of the calculation**
297

Figures generated by EMBER have two purposes: (1) to display the results of the calculation and (2) to track the variation of the best-fit estimation with parameters from the grid search such as M₀ or Cᵢ. Figure 3a shows the influence of the exsolved gas content to the resulting best-fit diffusion profile. Other parameters can be tested. For example, the number of volatile species fitted at once (single or up to 3, solid line or dashed line, respectively; Fig. 3b) also influences the resulting best-fit profile. Figure 3c shows the sensitivity of the cumulative error over decompression rate and exsolved gas content (M₀). Figure 3d shows the influence of the initial concentration on each volatile best fit error calculation.

298

299

300

301

302

303

304

305

Each figure generated is saved in a unique file (.fig) directly openable with EMBER. Along with the figures, EMBER also produces four Excel (.xls) files (.csv for Mac). They contain the input parameters, the diffusion profiles of each best case, for every decompression rate with the best Cᵢ and M₀, and a copy of the results for each volatile with the respective exsolved gas content.

306

307

308

309

310

311

# 4. Calculation and performance test

312

## 4.1 *A priori* requirements

313

As with all computational software, meaningful results in EMBER will only be achieved with appropriate dataset. Investigations should be limited to entirely glassy embayments, sampled from rapidly quenched deposits (e.g., <2 cm sized tephra), exhibiting a geometry that is close to that of a cylinder with constant radius. More complex geometries would void the core assumption of unidirectional diffusion and necessitate 3D diffusion modeling. [e.g., deGraffenried and Shea., 2020], which is not currently supported by EMBER.

The longer the embayment, the more likely it is to display a concentration plateau in the diffusion profile indicative of Cᵢ and the starting pressure of ascent (Pₘₐt). In absence of such condition, we recommend the user to determine Pₘₐt, Cᵢ and T from melt inclusions and geo-

314

315

316

317

318

319

320

321
thermobarometry studies. Finally, volatile concentration measurements on surrounding glass
are required to assess P<sub>end</sub>, the pressure of quenching.

While we used VolatileCalc and SolEx as the main degassing path generators in the
examples below, it should be noted that EMBER can also read degassing paths from any other
software as long as the input files comply with the required format (see the tutorial in the
additional instructions). Hence, the choice of degassing path software is ultimately up to the
user's preference.

## 4.2 Comparison with previous studies

We re-analysed the natural volatile diffusion profiles in embayments from previous
studies, to assess the quality of EMBER's decompression rate calculation: from the 1980 Mt
St Helens [Humphreys et al., 2008], 1974 Fuego [Lloyd et al., 2014], 1500, 1650 and 1959
Kīlauea [Ferguson et al., 2016], 27 ka Taupo [Myers et al., 2018], 767 ka Long Valley [Myers
et al., 2018], 2 Ma Yellowstone [Myers et al., 2018], 2017-2018 Ambae/Aoba [Moussallam et
al., 2019], Late Bronze Age (LBA) Santorini [Myers et al., 2021] and December 2018
Ambrym [Moussallam et al., 2021] eruptions. We used the uncertainties reported in each
study when provided. Also when provided, we directly used the specific grid search (e.g.
Ferguson et al., [2016] provides a range of M₀). Otherwise, we estimated the possible range of
grid search from the reported uncertainties for decompression rates and initial concentrations.
Degassing paths were calculated respecting P, T, and volatile contents used in the original
studies (typically, we use SolEx or VolatileCalc when authors specified it). Similarly, the re-
analysed results were calculated using the same set of values for M₀ as the original literature
data to which they are compared. In detail this was M₀ = 0.1, 1.6, 1.6, 0.8 wt. % respectively
for IkiE1, ReticE1, ReticE2 and KeaE1 [Ferguson et al., 2016], M₀ = 0.4, 0, 0, 3.2 wt. %
respectively for AF2, AD5, AE38 and PG11 [Moussallam et al., 2021] and M₀ = 0 wt. % for
all embayments from Ambae [Moussallam et al., 2019]. For all other studies of which the
exsolved gas content was not specified we used M₀ = 0 wt. %. In some cases [Myers et al.,
2018], we limited our tests to three samples from each investigated eruption, making sure to
test the samples with the highest and lowest estimation. All results except those from
Ambrym were calculated with the weighting and scaling of the error. The results are shown in
Figure 4.

There is a good agreement between EMBER results and most literature values. Data from the
1980's Mt St Helens eruption [Humphreys et al., 2008] show consistent disparity between the
published decompression rates and EMBER's. This difference is due to significantly different
boundary conditions of this study compared with the original. The model from the original
study [Humphreys et al., 2008] fixed the final concentration at the mouth of the embayment,
prior to the ascension. Since the diffusion phenomenon is gradient-dependent, imposing a
high concentration gradient at the beginning of the simulation should lead to faster volatile
diffusion and a higher decompression rate estimation. This is a critical simplification leading
to significant differences in the calculated decompression rate. Other minor differences
observed between our re-analysis and earlier decompression rates estimates come from (1) the
use of a different H₂O diffusion coefficient as well as a single episode of ascent in all our
calculation (as opposed to two-step ascent; [Lloyd et al., 2014]), (2) the calculation of
diffusion coefficient (e.g., taking into account the change of diffusivity due to H₂O diffusion
or choosing different equations of diffusivity) and (3) the weighting and scaling of error
associated with each volatile diffusion profile in all our calculations (Fig.4).

## 4.3 Monte Carlo simulation

We ran tests on three samples with variable number of Monte Carlo iterations (m=3, 11, 101,
501 and 1001) to determine the lowest number of the iterations to achieve stable assessment
of uncertainties: BT251, a rhyolite from Long Valley, ReE2, a reticulite from 1500 CE
Kīlauea lava fountain and AO02, a basalt from phase 2 of the 2017-2018 eruption of Ambae.
For BT251, the representation gradually acquires the shape of a pseudo-gaussian distribution
with increasing iteration counts. The median and mean values of all filtered decompression
rate do not change significantly with increasing iterations (Fig 5) and a stability of the data
distribution is quickly achieved (for m=11). For ReE1, the precision of the nanoSIMS
measurements and the low decompression rate leads to a stable calculation and a needle
shaped histogram even for m=3, suggesting that the results are well constrained. For AO02
the shape of the histogram is pointy, unimodal and skewing right with a long tail. With higher
dP/dt values (i.e., sharper diffusion profiles), modelled diffusion profiles become more and
more similar, making it harder for the software to find the best solution. This results in a
skewed solution histogram with a long tail and a more noticeable difference between mean
and median value. For such cases, stability is reached for m = 101. We recommend 101
iterations (EMBER default iteration value) as it was enough to reach a stable solution for our
most uncertain case.

387

388 4.4 How well does EMBER constrain M₀?

Calculation of the exsolved gas content has been introduced within the embayment
method with the calculations made by Ferguson et al., [2016]. They noted it had a significant
influence on their calculation of decompression rate, which was confirmed by later studies
and by EMBER (Fig 3a and Fig 6a). The effect is mostly a result of the exsolved gas content
having a large influence on the degassing trend. From our investigation of literature data, we
see that a change of M₀ from 0 to 3.2 wt. % can cause the decompression rate estimation to
increase by 10 fold (Fig 6 (a)). This variation is most pronounced on H₂O profiles in
embayments of basaltic composition (Fig 6 (a), (b), (c)). Except for these few cases, the
impact of $M_0$ is hardly significant enough for it to be accurately determined using only
diffusion profiles (the associated error frequently covering the entire range studied).

The change of decompression rate in response to variable $M_0$ illustrates the importance
of $P_{end}$ (Fig. 6 (d)) (which must be determined from measured data) and the critical influence
of the modelled degassing path.

The best fit $M_0$ values are directly affected by the choice of $P_{end}$ (Fig. 6 (d)) which must
be determined from measured data. Taking as an example, the calculation of one profile from
Ferguson et al., [2016] which ends at a specified quenching pressure $P_{end} = 2.75$ MPa. Using
the decompression path (from SolEx) for $M_0<1.6$ wt. %, the resulting $H_2O$ concentration at
$P=2.75$ MPa would be higher than $C_s(x = 0)$, the measured concentration at the mouth of the
embayment approximated with the glass measurement, by up to 0.20 wt. %. Synthesized
diffusion profiles in this case cannot reproduce the concentration at the mouth of the
embayment for $M_0<1.6$ wt. %. Because of this limitation, as shown in figure 3a, only a
degassing paths associated with the best fit $M_0$ value could match the whole diffusion
profile by fitting the concentration at the mouth. It should be noted that if a different $P_{end}$ is
chosen, it is possible that EMBER finds another best-fit $M_0$ (Fig 6 (d)). Therefore, we advise
the user to be cautious when choosing $P_{end}$ for decompression rate modelling, especially in
case constrained by one volatile species (Fig 6 (e), (f)).

## 4.5 Weighting and scaling of error

EMBER leaves the user the choice to, or not to, equally weigh and scale the relative
contribution of each volatile profile equally, when calculating the cumulative error. For
example, the study from Ferguson et al., [2016] gives equal weight ($u_j = 1$) to all three
volatile species, except on a few cases. By doing so, the best fit $NErr_j$ calculated from $H_2O$
profiles is several orders of magnitude higher than the best fit $NErr_j$ calculated for CO₂ or S
profiles. That difference is due to the range of concentration variation, as H₂O usually varies
within a few wt. % and CO₂ and S usually varies within thousands of ppm at best. Cumulative
error calculation, and subsequently the best fit determination, becomes heavily dependent on
the H₂O $NErr_j$ value, making the constraints brought by CO₂ or S profiles almost negligible.
Another approach is to weight the $NErr_j$ values either with analytical uncertainty [Myers et
al., 2018] or maximum measured concentration in the embayment [Newcombe et al., 2020] of
volatile species. Weighting the concentration with the error value nondimensionalizes the
profile and gives equal weight to the quality of fit of each volatile species. In EMBER, the
$NErr_j$ values are weighted by the analytical uncertainty on each volatile concentration. It
results in an even consideration of the constraints brought by each volatile species. There are
pros and cons for this choice, if one thinks the S degassing model to not be as accurate as the
H₂O degassing model for instance, one may choose not to weigh and scale the errors. The
user ultimately must make the choice and it should be reported.

## 5. Recalculation of decompression rates

In the previous section, we presented EMBER calculation results using initial conditions ($C_i$, $M_0$ if studied, range of $dPdt$, $P_{end}$, $P_{start}$) directly from the original literature studies to demonstrate the quality of EMBER and its ability to reproduce the former literature results. Now, we use EMBER to reprocess the raw data from the literatures, but this time using a uniform protocol that takes into account the following parameters: the presence of exsolved gas, the same set of formulas for diffusion coefficient calculation, a model with a single step ascent and an initial volatile element concentration, determined from the concentration plateau (Table 2). Detailed modelled profiles of each embayments can be found in the supplementary materials. We restricted the calculations to 3 profiles per eruption for Myers et
446

al., [2018], including the ones associated with the highest and lowest decompression rate recorded, and to 5 profiles for Myers et al., [2021].

447

448

## 5.1 The impact of calculation strategies

449

Our hope with developing EMBER is to minimize differences in modelling parameters (e.g., diffusion coefficients, error treatment, minimization strategy...) in future studies. It is important to realize however that differences in strategies will persist. For instance, whilst EMBER allow the user to set the exsolved gas content ($M_0$) as a free parameter to be determined by a grid-search, some users might prefer to impose an exsolved gas content based on independent constraint. This is the case for instance of the Moussallam et al., [2019] study of the 2017-2018 Ambae eruption where, in a second set of calculations, the authors used the difference between the volatile content of the melt inclusions and the starting volatile content of the embayments to estimate the likely amount of exsolved volatile in the system at the initiation of ascent.

450

451

452

453

454

455

456

457

458

459

Another strategy is to use the volatile content of melt inclusions as starting conditions for the diffusion model instead of the measured plateau values in the embayment interior. This was the strategy adopted by Myers et al., [2021] for the Minoan eruption of Santorini, where a fixed initial concentration of 5.2 wt. % H₂O (except for one embayment at 5.6 wt. % H₂O) was used for all model calculations. Such strategy was also used by Humphreys et al., [2008] for Mt St Helens calculations.

460

461

462

463

464

465

In this study, we recalculated all embayment profiles using the same strategy throughout, but we do not pass judgement on the validity of one strategy over another. Our aim is to present a dataset which is as comparable as possible. The data presented in Fig. 7 and Table 2 are hence all calculated leaving the exsolved volatile content ($M_0$) as unconstrained (i.e., as a part of the

466

467

468
469

grid search) and using the plateau values in volatile content recorded in embayment interiors
as the model starting conditions (Cᵢ).

470

471

This difference in starting assumptions can lead to significant differences in the resulting
decompression rate. Recalculating the embayment data from Myers et al., [2021] with our
protocol leads to decompression rates 3 to 20 times higher. It is therefore of paramount
importance that users of EMBER explicitly report their assumptions and starting conditions.
We recommend that future compilations continue to reprocess original data in a consistent
manner (as done here) in order to render inter-study comparison as coherent as possible.

472

473

474

475

476

477

479

## 5.2 Decompression rates vs eruption parameters

480

481

482

483

484

485

486

487

488

489

490

A recent literature compilation showed a clear relationship between magma decompression
rate and whether an eruption is explosive or effusive in character [Cassidy et al., 2018]. This
relationship deserves more scrutiny to establish if finer relationships between magma
decompression rates and explosivity exists in nature. The VEI (Volcanic Explosivity Index),
as an approximation of eruption explosivity, is related to eruption magnitude and/or plume
height [Newhall and Self., 1982]. To show the utility of EMBER for the studies of volcanic
explosivity, we tested the presence (or absence) of a correlation between decompression rate
and both eruption magnitude and plume height, with our reprocessed data. Previous studies
have found that the eruption magnitude is positively correlated with decompression rate
[Ferguson et al., 2016; Moussallam et al., 2019]. However, our compilation shows no clear
correlation between decompression rate and eruption magnitude (Fig. 7a).

491

492

493

We notice that an increase in decompression rate with eruption magnitude is noticeable for
basaltic magmas (Kīlauea, Ambae and Ambrym) but the correlation is weakly significant
(Pearson coefficient of 0.24 with a p-value of 0.35 and R²=0.47, which corresponds to a weak
494

positive correlation). In a single eruption, there can be varying flow regime creating a range of decompression rate within the conduit; a batch of magma potentially ascend faster than the others (e.g., [Martel et al., 1998; Gonnermann and Manga., 2007; Cassidy et al., 2015]). If we only consider the highest decompression rate of any given eruption however, the aforementioned weak correlation between decompression rate and magnitude becomes more significant for mafic eruptions (Figure 7a) (Pearson coefficient of 0.93 with a p-value of 0.01 and R²=0.86). The calculated decompression rates for rhyolitic eruptions on the other hand shows no correlation with eruption magnitude.

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

516

517

The eruption intensity is also assessed with the plume height, which is directly observed or calculated through several empirically-determined relations involving isopachs for explosive eruptions (e.g., [Woods and Wohletz., 1991; Mastin et al., 2009; Pyle., 2015]). The eruption intensity, estimated from the mass eruption rate, has been shown to be positively correlated with decompression rate [Barth et al., 2019; Ferguson et al., 2016; Newcombe et al., 2020]. Thus, the plume height is expected to correlate with the decompression rate. We used EMBER-recalculated decompression rates to test this hypothesis. For the two eruptions of Kīlauea for which the plume height is not constrained, we used the maximal height of the lava fountain instead. Ambrym's plume height data for the 2018 eruption is undetermined and therefore not added to Figure 7b. Results show no global correlation between these two parameters (Figure 7b). However, considering again the highest decompression rate recorded for each eruption, a strong positive correlation emerges for basaltic eruptions (with a Person coefficient of 0.76, a p-value of 0.005 and R²=0.88). Again, the calculated decompression rate for rhyolitic eruptions shows no correlation with eruption plume height. The maximum decompression rate for basaltic eruption therefore shows a statistically significant positive correlation with both eruption magnitude and plume height. A first order positive correlation
therefore exists between the maximal magma decompression rate and the explosivity of an eruption for basaltic eruptions.

520

521

## 6. Conclusion

* We present EMBER, a user-friendly GUI program that calculates decompression rates from H₂O, CO₂ and S concentration profiles along embayments of basaltic to rhyolitic compositions.

* We found that decompression rate calculations are particularly sensitive to variations of M₀ especially for the H₂O diffusion profile. Variations of Pₐnd are also accompanied by a variation of the best fit exsolved gas content, but not necessarily by a variation of associated decompression rate.

* We recalculated decompression rates from previous studies twice: first, to validate and test how well EMBER reproduced existing results using the parameters from the original studies, and secondly, to homogenize determined decompression rates applying the same protocol to the existing raw data from previous studies, in order to improve inter-study comparison.

* In the first case, recalculated decompression rates are in the same order of magnitude as original calculations but notable differences do occur such as for the 1980 Mt St Helens eruption which recalculated decompression rate are at 0.15-0.41 MPa/s, half of the previously reported values [Humphreys et al., 2008].

* In the second case, recalculated dataset shows no significant correlation between magma decompression rate and eruption magnitude when considering the entire dataset and shows a weak correlation when considering the subset of decompression rates of basaltic magma (Pearson coefficient of 0.24 with a p-value of 0.35 and
543

544

545

546

547

548

549

R²=0.47.) The correlation is significant when considering only the maximum decompression rates of each basaltic eruption (Pearson coefficient of 0.93 with a p-value of 0.01 and R²=0.86). Additionally, there is no significant correlation between decompression rate and plume height when considering the entire dataset. However, once again, a statistically significant trend appears when considering only the maximum decompression rate of the basaltic eruptions (with a Person coefficient of 0.84, a p-value of 0.007 and R²=0.88).

550

551

552

* Our results suggest for the first time, a significant positive correlation, between embayment-calculated maximum decompression rate and eruption explosivity parameters such as magnitude and plume height, for basaltic eruptions.

553
554 **7. Computer code availability**

555 • Title of the software: EMBER

556

557 • Developers: Guillaume Georgeais, Kenneth T. Koga, Yves Moussallam, Estelle F. Rose-Koga

558

559 • Hardware used: EMBER was run on a computer with 4 core (4.5GHz) and 16GB of RAM

560

561 • Link to the software, tutorial and test using data from [Ferguson et al., 2016] DOI: 10.6084/m9 figshare.13663811

562

563 • Software required: MATLAB Runtime Environment 2019b and up
564 DOI: https://www.mathworks.com/products/compiler/matlab-runtime.html

565

569 **Acknowledgement**

570 We would like to thank M. Humphreys for providing the raw data of the Mt St Helens

571 embayments. GG was supported by a PhD fellowship from the French Government

572 "Ministère de l'Enseignement Supérieur, de la Recherche et de l'Innovation". ERK

573 acknowledges partial funding from Laboratory of Excellence initiative n°ANR-10-LABX-

574 0006, the Région Auvergne and the European Regional Development Funds. Y.M

575 acknowledges funding from INSU and the Région Auvergne Rhone Alpes. This is Laboratory

576 of Excellence ClerVolc contribution number 463. We thank Madison Myers for sharing her

577 code during the review process of our manuscript. We would like to thank Madison Myers
578 and an Anonymous reviewer for their comments on the original manuscript and Marie

579 Edmonds for editorial handling.

580 **Bibliography**

581 Aki, K., & Koyanagi, R. (1981). Deep volcanic tremor and magma ascent mechanism under Kilauea, Hawaii. *Journal of Geophysical Research, 86*(B8), 7095.

582

583 https://doi.org/10.1029/JB086iB08p07095
584 Anderson, A. T. (1991). Hourglass inclusions: Theory and application to the Bishop Rhyolitic Tuff.

585 *American Mineralogist*, **76**(3-4), 530-547.

586 Barth, A., Newcombe, M., Plank, T., Gonnermann, H., Hajimirza, S., Soto, G. J., et al. (2019). Magma

587 decompression rate correlates with explosivity at basaltic volcanoes — Constraints from

588 water diffusion in olivine. *Journal of Volcanology and Geothermal Research*, **387**, 106664.

589 https://doi.org/10.1016/j.jvolgeores.2019.106664

590 Behrens, H., & Zhang, Y. (2001). Ar diffusion in hydrous silicic melts: Implications for volatile diffusion

591 mechanisms and fractionation. *Earth and Planetary Science Letters*, **192**, 363-376.

592 https://doi.org/10.1016/S0012-821X(01)00458-7

593 Browne, B., & Gardner, J. (2006). The influence of magma ascent path on the texture, mineralogy,

594 and formation of hornblende reaction rims. *Earth and Planetary Science Letters*, **246**(3-4),

595 161-176. https://doi.org/10.1016/j.epsl.2006.05.006

596 Carey, S., & Sigurdsson, H. (1985). The May 18, 1980 eruption of Mount St. Helens: 2. Modeling of

597 dynamics of the Plinian Phase. *Journal of Geophysical Research: Solid Earth*, **90**(B4), 2948-

598 2958. https://doi.org/10.1029/JB090iB04p02948

599 Cassidy, M., Cole, Paul. D., Hicks, K. E., Varley, N. R., Peters, N., & Lerner, A. H. (2015). Rapid and

600 slow: Varying magma ascent rates as a mechanism for Vulcanian explosions. *Earth and*

601 *Planetary Science Letters*, **420**, 73-84. https://doi.org/10.1016/j.epsl.2015.03.025

602 Cassidy, M., Manga, M., Cashman, K., & Bachmann, O. (2018). Controls on explosive-effusive volcanic

603 eruption styles. *Nature Communications*, **9**(1), 2839. https://doi.org/10.1038/s41467-018-05293-3

604 05293-3

605 deGraffenried, R., & Shea, T. (2020). Modeling diffusion in 1D within melt embayments: correcting

606 for 3D geometry. Presented at the AGU Fall Meeting 2020, AGU. Retrieved from

607 https://agu.confex.com/agu/fm20/meetingapp.cgi/Paper/676363
608 Demouchy, S., & Mackwell, S. (2006). Mechanisms of hydrogen incorporation and diffusion in iron-

609 bearing olivine. *Physics and Chemistry of Minerals*, 33(5), 347–355.

610 https://doi.org/10.1007/s00269-006-0081-2

611 Faure, F., & Schiano, P. (2005). Experimental investigation of equilibration conditions during

612 forsterite growth and melt inclusion formation. *Earth and Planetary Science Letters*, 236(3),

613 882–898. https://doi.org/10.1016/j.epsl.2005.04.050

614 Ferguson, D. J., Gonnermann, H. M., Ruprecht, P., Plank, T., Hauri, E. H., Houghton, B. F., & Swanson,

615 D. A. (2016). Magma decompression rates during explosive eruptions of Kīlauea volcano,

616 Hawaii, recorded by melt embayments. *Bulletin of Volcanology*, 78(10), 71.

617 https://doi.org/10.1007/s00445-016-1064-x

618 Freda, C., Baker, D. R., Romano, C., & Scarlato, P. (2003). Water diffusion in natural potassic melts.

619 *Geological Society, London, Special Publications*, 213(1), 53–62.

620 https://doi.org/10.1144/GSL.SP.2003.213.01.04

621 Freda, Carmela, Baker, D. R., & Scarlato, P. (2005). Sulfur diffusion in basaltic melts. *Geochimica et*

622 *Cosmochimica Acta*, 69(21), 5061–5069. https://doi.org/10.1016/j.gca.2005.02.002

623 Gardner, J. E., Sigurdsson, H., & Carey, S. N. (1991). Eruption dynamics and magma withdrawal during

624 the Plinian Phase of the Bishop Tuff Eruption, Long Valley Caldera. *Journal of Geophysical*

625 *Research*, 96(B5), 8097. https://doi.org/10.1029/91JB00257

626 Geschwind, C.-H., & Rutherford, M. J. (1995). Crystallization of microlites during magma ascent: the

627 fluid mechanics of 1980–1986 eruptions at Mount St Helens. *Bulletin of Volcanology*, 57(5),

628 356–370. https://doi.org/10.1007/BF00301293

629 Girona, T., & Costa, F. (2013). DIPRA: A user-friendly program to model multi-element diffusion in

630 olivine with applications to timescales of magmatic processes. *Geochemistry, Geophysics,*

631 *Geosystems*, 14(2), 422–431. https://doi.org/10.1029/2012GC004427

632 Gonnermann, H. M., & Manga, M. (2007). The Fluid Mechanics Inside a Volcano. *Annual Review of*

633 *Fluid Mechanics*, 39(1), 321–356. https://doi.org/10.1146/annurev.fluid.39.050905.110207
634 Humphreys, M. C. S., Menand, T., Blundy, J. D., & Klimm, K. (2008a). Magma ascent rates in explosive eruptions: Constraints from H2O diffusion in melt inclusions. *Earth and Planetary Science Letters*, 270(1–2), 25–40. https://doi.org/10.1016/j.epsl.2008.02.041

637 Humphreys, M. C. S., Menand, T., Blundy, J. D., & Klimm, K. (2008b). Magma ascent rates in explosive eruptions: Constraints from H2O diffusion in melt inclusions. *Earth and Planetary Science Letters*, 270(1–2), 25–40. https://doi.org/10.1016/j.epsl.2008.02.041

640 Johnston, E., R.S.J., S., Phillips, J., & Carey, S. (2014). Revised Estimates for the Volume of the Late Bronze Age Minoan Eruption, Santorini, Greece. *Journal of the Geological Society, 171*. https://doi.org/10.1144/jgs2013-113

643 Klein, F. W., Koyanagi, R. Y., Nakata, J. S., & Tanigawa, W. R. (1987). The seismicity of Kilauea's magma system. *US Geol. Surv. Prof. Pap*, 1350(2), 1019–1185.

645 Klügel, A. (1998). Reactions between mantle xenoliths and host magma beneath La Palma (Canary Islands): constraints on magma ascent rates and crustal reservoirs. *Contributions to Mineralogy and Petrology*, 131(2), 237–257. https://doi.org/10.1007/s004100050391

648 Klügel, A., Hansteen, T. H., & Schmincke, H.-U. (1997). Rates of magma ascent and depths of magma reservoirs beneath La Palma (Canary Islands). *Terra Nova*, 9(3), 117–121. https://doi.org/10.1046/j.1365-3121.1997.d01-15.x

651 Kohlstedt, D. L., & Mackwell, S. J. (1998). Diffusion of hydrogen and intrinsic point defects in olivine. *Zeitschrift Fur Physikalische Chemie*, 207(1–2), 147–162.

653 Le Voyer, M., Asimow, P. D., Mosenfelder, J. L., Guan, Y., Wallace, P. J., Schiano, P., et al. (2014). Zonation of H2O and F Concentrations around Melt Inclusions in Olivines. *Journal of Petrology*, 55(4), 685–707. https://doi.org/10.1093/petrology/egu003

656 Liu, Y., Anderson, A. T., & Wilson, C. J. N. (2007). Melt pockets in phenocrysts and decompression rates of silicic magmas before fragmentation. *Journal of Geophysical Research*, 112(B6), B06204. https://doi.org/10.1029/2006JB004500
659 Lloyd, A. S., Ruprecht, P., Hauri, E. H., Rose, W., Gonnermann, H. M., & Plank, T. (2014). NanoSIMS
660 results from olivine-hosted melt embayments: Magma ascent rate during explosive basaltic
eruptions. *Journal of Volcanology and Geothermal Research*, 283, 1–18.
661 https://doi.org/10.1016/j.jvolgeores.2014.06.002

662 Mackwell, S. J., & Kohlstedt, D. L. (1990). Diffusion of hydrogen in olivine: Implications for water in
663 the mantle. *Journal of Geophysical Research*, 95(B4), 5079.
664 https://doi.org/10.1029/JB095iB04p05079

665 Martel, C., Pichavant, M., Bourdier, J.-L., Traineau, H., Holtz, F., & Scaillet, B. (1998). Magma storage
666 conditions and control of eruption regime in silicic volcanoes: experimental evidence from
667 Mt. Pelée. *Earth and Planetary Science Letters*, 156(1), 89–99.
668 https://doi.org/10.1016/S0012-821X(98)00003-X

669 Martel, Caroline, & Iacono-Marziano, G. (2015). Timescales of bubble coalescence, outgassing, and
670 foam collapse in decompressed rhyolitic melts. *Earth and Planetary Science Letters*, 412,
671 173–185. https://doi.org/10.1016/j.epsl.2014.12.010

672 Mason, B. G., Pyle, D. M., & Oppenheimer, C. (2004). The size and frequency of the largest explosive
eruptions on Earth. *Bulletin of Volcanology*, 66(8), 735–748. https://doi.org/10.1007/s00445-004-0355-9

673 Mastin, L., Van Eaton, A., & Lowenstern, J. (2014). Modeling ash fall distribution from a Yellowstone
supereruption. *Geochemistry, Geophysics, Geosystems*, 15.
674 https://doi.org/10.1002/2014GC005469

675 Mastin, L. G., Guffanti, M., Servranckx, R., Webley, P., Barsotti, S., Dean, K., et al. (2009). A
multidisciplinary effort to assign realistic source parameters to models of volcanic ash-cloud
transport and dispersion during eruptions. *Journal of Volcanology and Geothermal Research*,
*186*(1), 10–21. https://doi.org/10.1016/j.jvolgeores.2009.01.008

680 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
685 investigation. *Contributions to Mineralogy and Petrology*, 174(11), 90.
686 https://doi.org/10.1007/s00410-019-1625-z
687 Moussallam, Y., Médard, E., Georgeais, G., Rose-Koga, E. F., Koga, K. T., Pelletier, B., et al. (2021).
688 How to turn off a lava lake? A petrological investigation of the 2018 intra-caldera and
689 submarine eruptions of Ambrym volcano. *Bulletin of Volcanology*.
690 https://doi.org/10.1007/s00445-021-01455-2
691 Myers, M. L., Wallace, P. J., Wilson, C. J. N., Watkins, J. M., & Liu, Y. (2018). Ascent rates of rhyolitic
692 magma at the onset of three caldera-forming eruptions. *American Mineralogist*, 103(6), 952-
693 965. https://doi.org/10.2138/am-2018-6225
694 Myers, M. L., Druitt, T. H., Schiavi, F., Gurioli, L., & Flaherty, T. (2021). Evolution of magma
695 decompression and discharge during a Plinian event (Late Bronze-Age eruption, Santorini)
696 from multiple eruption-intensity proxies. *Bulletin of Volcanology*, 83(3), 18.
697 https://doi.org/10.1007/s00445-021-01438-3
698 Newcombe, M. E., Plank, T., Barth, A., Asimow, P. D., & Hauri, E. (2020). Water-in-olivine magma
699 ascent chronometry: Every crystal is a clock. *Journal of Volcanology and Geothermal
700 Research*, 398, 106872. https://doi.org/10.1016/j.jvolgeores.2020.106872
701 Newhall, C. G., & Self, S. (1982). The volcanic explosivity index (VEI) an estimate of explosive
702 magnitude for historical volcanism. *Journal of Geophysical Research*, 87(C2), 1231.
703 https://doi.org/10.1029/JC087iC02p01231
704 Newman, S., & Lowenstern, J. B. (2002). VolatileCalc: a silicate melt–H2O–CO2 solution model
705 written in Visual Basic for excel. *Computers & Geosciences*, 28(5), 597–604.
706 https://doi.org/10.1016/S0098-3004(01)00081-4
707 Ni, H., & Zhang, L. (2018). A general model of water diffusivity in calc-alkaline silicate melts and
708 glasses. *Chemical Geology*, 478, 60–68. https://doi.org/10.1016/j.chemgeo.2017.10.010
709 Nowak, M., & Behrens, H. (1997). An experimental investigation on diffusion of water in haplogranitic melts. *Contributions to Mineralogy and Petrology*, 126(4), 365-376.
710 https://doi.org/10.1007/s004100050256
711 Nowak, Marcus, 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. *Geochimica et Cosmochimica Acta*, 68, 5127-5138. https://doi.org/10.1016/j.gca.2004.06.002
712 Pyle, D. M. (2015). Sizes of Volcanic Eruptions. In *The Encyclopedia of Volcanoes* (pp. 257-264).
713 Elsevier. https://doi.org/10.1016/B978-0-12-385938-9.00013-4
714 Rutherford, M. J. (2008). Magma Ascent Rates. *Reviews in Mineralogy and Geochemistry*, 69(1), 241-
715 271. https://doi.org/10.2138/rmg.2008.69.7
716 Rutherford, Malcolm J., & Hill, P. M. (1993). Magma ascent rates from amphibole breakdown: An experimental study applied to the 1980-1986 Mount St. Helens eruptions. *Journal of Geophysical Research: Solid Earth*, 98(B11), 19667-19685.
717 https://doi.org/10.1029/93JB01613
718 Scandone, R., & Malone, S. D. (1985). Magma supply, magma discharge and readjustment of the feeding system of mount St. Helens during 1980. *Journal of Volcanology and Geothermal Research*, 23(3-4), 239-262. https://doi.org/10.1016/0377-0273(85)90036-8
719 Self, S. (2006). The effects and consequences of very large explosive volcanic eruptions. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, 364(1845), 2073-2097. https://doi.org/10.1098/rsta.2006.1814
720 Shreve, T., Grandin, R., Boichu, M., Garaebiti, E., Moussallam, Y., Ballu, V., et al. (2019). From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics. *Scientific Reports*, 9(1), 18868.
721 https://doi.org/10.1038/s41598-019-55141-7
722
733 Sparks, R. S. J., Baker, L., Brown, R., Field, M., Schumacher, J., Stripp, G., & L. Walters, A. (2006).

734 Dynamical constraints on kimberlite volcanism. *Journal of Volcanology and Geothermal*

735 *Research*, 155, 18–48. https://doi.org/10.1016/j.jvolgeores.2006.02.010

736 Toramaru, A. (2006). BND (bubble number density) decompression rate meter for explosive volcanic

737 eruptions. *Journal of Volcanology and Geothermal Research*, 154(3–4), 303–316.

738 https://doi.org/10.1016/j.jvolgeores.2006.03.027

739 Toramaru, Atsushi. (1989). Vesiculation process and bubble size distributions in ascending magmas

740 with constant velocities. *Journal of Geophysical Research*, 94(B12), 17523.

741 https://doi.org/10.1029/JB094iB12p17523

742 Toramaru, Atsushi. (1995). Numerical study of nucleation and growth of bubbles in viscous magmas.

743 *Journal of Geophysical Research: Solid Earth*, 100(B2), 1913–1931.

744 https://doi.org/10.1029/94JB02775

745 Tryggvason, E. (1994). Surface deformation at the Krafla volcano, North Iceland, 1982–1992. *Bulletin*

746 *of Volcanology*, 56(2), 98–107. https://doi.org/10.1007/BF00304105

747 Van Eaton, A., Herzog, M., Wilson, C., & McGregor, J. (2012). Ascent dynamics of large

748 phreatomagmatic eruption clouds: The role of microphysics. *Journal of Geophysical Research*

749 (*Solid Earth*), 117, 3203. https://doi.org/10.1029/2011JB008892

750 Witham, F., Blundy, J., Kohn, S. C., Lesne, P., Dixon, J., Churakov, S. V., & Botcharnikov, R. (2012).

751 SolEx: A model for mixed COHSCI-volatile solubilities and exsolved gas compositions in basalt.

752 *Computers & Geosciences*, 45, 87–97. https://doi.org/10.1016/j.cageo.2011.09.021

753 Woods, A. W., & Wohletz, K. (1991). Dimensions and dynamics of co-ignimbrite eruption columns.

754 *Nature*, 350(6315), 225–227. https://doi.org/10.1038/350225a0

755 Zhang, Y., & Behrens, H. (2000). H2O diffusion in rhyolitic melts and glasses. *Chemical Geology*,

756 169(1), 243–262. https://doi.org/10.1016/S0009-2541(99)00231-4

757 Zhang, Y., Xu, Z., Zhu, M., & Wang, H. (2007). Silicate melt properties and volcanic eruptions: SILICATE

758 MELT PROPERTIES. *Reviews of Geophysics*, 45(4). https://doi.org/10.1029/2006RG000216
759 Zhang, Y., Ni, H., & Chen, Y. (2010). Diffusion Data in Silicate Melts. *Reviews in Mineralogy and Geochemistry*, 72(1), 311–408. https://doi.org/10.2138/rmg.2010.72.8

760

761

762
763 **Figure caption**

764 **Figure 1**: Flow chart showing the calculation steps used in EMBER to determine the
765 parameters and associated errors resulting in the best fit of natural volatile diffusion profiles.
766 The outputs are a decompression rate and an initial (dissolved) volatile concentration as well
767 as their associated confidence interval. This process of calculation is repeated for each value
768 of exsolved volatile content ($M_0$). Estimation of the $M_0$ value and its related uncertainty is
769 done in a similar manner but by considering the range of $M_0$ variation as the 4th dimension of
770 the starting matrix.
771

771

772 **Figure 2**: A display example of EMBER's GUI after a cycle of calculations. Input parameters
773 are entered on the left panel (part 1) and results are displayed on the right (part 2). The input
774 section is divided in three main parts, the grid search definition (1a), the main parameters
775 and model selection (1b), and the parameters for the Monte-Carlo simulation (uncertainty
776 and number of iterations) (1c). The result section is in three parts: comparison between best
777 fit and natural diffusion profiles (2a), the result parameter of each best fit (2b) and the log
778 section to follow the evolution of the calculation (2c). Clicking on the "Display all figures"
779 once the program is done with the calculation will display best fits for each exsolved volatile
780 content value and for each studied volatile specie. It also displays 3 3D-plots, of the evolution
781 of the calculative error for each volatile(1) vs. the decompression rate and initial
782 concentration, (2) the evolution of the decompression rate vs. exsolved volatile content and
783 initial volatile concentration for each studied volatile specie and (3) the evolution of the
784 cumulative error vs. decompression rate and exsolved volatile content (Figure 3 ).
785

785

786 **Figure 3**: Subset of plots generated from a calculation on the 1500CE Kīlauea reticulite
787 ReaE1 studied previously by Ferguson et al., [2016]. (a) Influence of exsolved volatile
788 content $M_0$ on best fit diffusion profiles of H₂O. Initial parameters from the grid search are
789 listed in the caption. (b) Best fit diffusion profiles calculated by considering H₂O only (solid
790 line) or all studied species (dashed lines). (c) Evolution of the cumulated scaled log₁₀
791 (NErr_m) values over a range of decompression rate and exsolved volatile content. (d)
792 Evolution of the NErr_m value for each volatile over a range of decompression rate and
793 initial volatile content. The NErr_m values are in a log₁₀scale. The values of interest, linked to
794 the smallest error value, are highlighted by a blue, red or yellow dot for H₂O, CO₂ and S
795 respectively (d), and a black dot for the cumulative error (c).
796

796

797 **Figure 4**: Comparison between recalculated data with EMBER and previous estimation from
798 the literature. EMBER estimations are comparable with literature estimates with a few
799 systematic differences (see text for details). Dashed line indicates the 1:1 ratio and thin solid
800 lines indicates 2:1 and 0.5:1 ratios. Results from EMBER were acquired taking the input
801 conditions of the respective original studies. Calculations for samples from Lloyd et al.,
802 :[2014] were made using their original H₂O diffusion coefficient. All the results displayed are
803

804

805

calculated for M₀=0 excepted when original studies provided a range of M₀ [Ferguson et al., 2016]. EMBER calculations were made using the “weighting and scaling of error” option which is a source of difference with previously reported literature decompression rate.

806

807

808

809

810

811

812

**Figure 5:** Statistical repartition of all filtered decompression rate values with increasing iteration count of the Monte Carlo approach. The repartition gradually acquires a more pronounced shape. The red bars indicates the 1-σ confidence interval of those decompression rates while the blue, green and black ones respectively indicate the best-fit result, the median and the mean. Calculation is done with fixed initial concentration and exsolved volatile content. All displayed data are in MPa/s

813

814

815

816

817

818

819

820

821

**Figure 6:** Evolution of modeled relative decompression rates from H₂O (a), CO₂ (b) and S (c) as a function of the exsolved volatile content for a selection of embayments. Each case in a), b) or c) is calculated with a fixed Pₐ₇₉, reported in Table 2. For figures (d), (e) and (f), Pᵣₑf is the Pₐ₇₉ value displayed in Table 2 for ReE2 and BT251. (d) The initial exsolved volatile content needed to generate the best fit (Mₚₑf) increases for ReE2 when changing Pₐ₇₉ but not systematically for M413. (e) Variation of relative decompression rate modeled from the H₂O diffusion profile only with Pₐ₇₉. (f) Variation of relative decompression rate modeled from all studied species (H₂O, CO₂ and S for ReE2, and H₂O and CO₂ for M413) with Pₐ₇₉.

822

823

824

825

826

827

828

829

830

831

832

833

834

835

836

837

838

**Figure 7:** (a) Compilation of EMBER-recalculated decompression rate as a function of eruption magnitude. The bigger mark indicates the highest decompression rate estimate for each eruption. For basaltic magmas, there is a trend of increasing decompression rate related to increasing magnitude only when considering the maximum decompression rate values. The trend does not apply to rhyolitic magmas. (b) Compilation of EBER-recalculated decompression rate as a function of eruption plume height. The bigger mark indicates the highest decompression rate estimate for each eruption. With the exception of Ambrym, which plume maximal height is unknown and assessed by lava fountain height, there is a notable correlation between maximal decompression rate and plume height for basaltic magmas. The same caption applies for both figures. Magnitude estimates were calculated from either original article or respectively Shreve et al., [2019], Johnston et al., [2014], Self, [2006] and Mason et al., [2004] for Ambrym's 2018 eruption, Santorini's Minoan eruption, Taupo's Oruanui eruption and both Long Valley's Bishop Tuff and Yellowstone's Huckleberry Ridge Tuff. Plume height estimation comes from modelisations from Gardner et al., [1991], Mastin et al., [2014] and Van Eaton et al., [2012] respectively for Long Valley's Bishop Tuff, Yellowstone's Huckleberry Ridge Tuff and Taupo's Oruanui eruptions.

839

**Table caption**

840

841

842

843

844

845

**Table 1:** Comparison of previously published codes. EMBER is a fully available complete software that covers a large spectrum of magma composition: dPdt = decompression rate, Cᵢ = initial concentration, M₀ = exsolved volatile content, Pᵢ = pressure of quench, gs stands for “grid search”. Every study that model decompression rates from CO₂ and/or S uses the diffusion coefficients calculated from [Zhang et al., 2007]. The only two exceptions are the CO₂ diffusivity from [Liu et al.,
846 2007] which is calculated from [Behrens and Zhang, 2001] and the S diffusivity from [Ferguson et al.,

847 2016], which is calculated from [Freda et al., 2005].

848

849

850 **Table 2**: Re-calculated decompression rates with initial conditions and related eruption parameters.

851 Ascent rates were calculated with density values from each original studies using the equation

852 Asc_rate = (Decomp_rate/9.81 * Crustal_density) 10^6 with Decomp_rate in MPa/s, Asc_rate

853 in m/s and the crustal density in kg/m3

854

855

856
<table><caption>858 Table 1</caption><thead><tr><th>Article</th><th>Architecture</th><th>Volatiles studied</th><th>Magmatic composition</th><th>H<sub>2</sub>O diffusion coefficient</th><th>Input parameters</th><th>M<sub>0</sub> study</th><th>Availability</th></tr></thead><tbody><tr><td>Liu et al, 2007</td><td>FORTRAN 77</td><td>H<sub>2</sub>O</td><td>Rhyolitic</td><td>a</td><td>dPdt, C<sub>i</sub>, T, M<sub>0</sub>, P<sub>f</sub>, P<sub>0</sub> *</td><td>yes</td><td>On request</td></tr><tr><td>Humphreys et al, 2008</td><td>COMSOL</td><td>H<sub>2</sub>O</td><td>Rhyolitic</td><td>a ; b</td><td>dPdt, C<sub>i</sub>, P<sub>f</sub>, P<sub>0</sub>, T *</td><td></td><td></td></tr><tr><td>Lloyd et al, 2014</td><td>n.d</td><td>H<sub>2</sub>O,CO<sub>2</sub>,S</td><td>Intermediate</td><td>Self-determined empiric law</td><td>dPdt, P<sub>f</sub>, P<sub>0</sub>, C<sub>i</sub>, T *</td><td></td><td></td></tr><tr><td>Ferguson et al, 2016</td><td>MATLAB</td><td>H<sub>2</sub>O,CO<sub>2</sub>,S</td><td>Basaltic</td><td>c</td><td>gs(dPdt, C<sub>i</sub>, M<sub>0</sub>), P<sub>f</sub>, P<sub>0</sub>, T</td><td>yes</td><td></td></tr><tr><td>Myers et al, 2018</td><td>MATLAB</td><td>H<sub>2</sub>O,CO<sub>2</sub></td><td>Rhyolitic</td><td>a</td><td>gs(dPdt, P<sub>f</sub>, M<sub>0</sub>), C<sub>i</sub>, P<sub>0</sub>, T</td><td>yes</td><td>On request</td></tr><tr><td>Moussallam et al, 2019</td><td>RStudio</td><td>H<sub>2</sub>O</td><td>Basaltic</td><td>Constant (e)</td><td>dPdt, C<sub>i</sub>, P<sub>f</sub>, P<sub>0</sub>, T</td><td>yes</td><td>Publicly Available</td></tr><tr><td>Newcombe et al, 2020</td><td>MATLAB</td><td>H<sub>2</sub>O,CO<sub>2</sub>,S</td><td>Intermediate</td><td>d</td><td>dPdt, M<sub>0</sub>, P<sub>f</sub>, P<sub>0</sub>, T *</td><td>yes</td><td></td></tr><tr><td>EMBER (This study)</td><td>MATLAB</td><td>H<sub>2</sub>O,CO<sub>2</sub>,S</td><td>Basaltic, Intermediate and Rhyolitic</td><td>e ; d ; a</td><td>gs(dPdt, C<sub>i</sub>, M<sub>0</sub>), P<sub>f</sub>, P<sub>0</sub>, T</td><td>yes</td><td>Publicly Available</td></tr></tbody></table>

a = [Zhang and Behrens, 2000]

c = [Zhang et al., 2010]

e = [Freda et al., 2003]

b = [Nowak and Behrens, 1997]

d = [Ni and Zhang, 2018]

*unconfirmed use of a grid search
<table><thead><tr><th rowspan="2"></th><th colspan="4">Decompression rate (Mpa/s)</th><th rowspan="2">Ascent rate (m/s)</th><th colspan="3">Initial concentration</th><th rowspan="2">M (wt%)</th><th rowspan="2">Time of<br>ascension (s)</th><th rowspan="2">Pstart (MPa)</th><th rowspan="2">Pend (MPa)</th><th rowspan="2">T(°C)</th><th rowspan="2">Magnitude</th></tr><tr><th>H2O based</th><th>CO2 based</th><th>S based</th><th>all species</th><th>H2O (wt%)</th><th>CO2 (ppm)</th><th>S (ppm)</th></tr></thead><tbody><tr><td colspan="16"><strong>Ferguson et al., 2016</strong></td></tr><tr><td>IkiE1</td><td>0.03 + 0.01</td><td>0.02 + 0.00</td><td>0.04 + 0.01</td><td><strong>0.04</strong> + <strong>0.00</strong></td><td>1.4</td><td>0.60 + 0.00</td><td>300 + 0</td><td>1390 + 0</td><td>0.1 + 0.1</td><td>2755</td><td>107</td><td>0.1</td><td>1192</td><td>1.5</td></tr><tr><td>ReticE1</td><td>0.08 + 0.03</td><td>0.09 + 0.00</td><td>0.08 + 0.01</td><td><strong>0.08</strong> + <strong>0.01</strong></td><td>2.9</td><td>0.55 + 0.00</td><td>184 + 0</td><td>1316 + 0</td><td>1.6 + 0.0</td><td>553</td><td>47</td><td>2.75</td><td>1163</td><td>2.5</td></tr><tr><td>ReticE2</td><td>0.07 + 0.03</td><td>0.03 + 0.01</td><td>0.40 + 0.00</td><td><strong>0.06</strong> + <strong>0.02</strong></td><td>2.1</td><td>0.57 + 0.00</td><td>170 + 0</td><td>1320 + 0</td><td>3.2 + 0.0</td><td>676</td><td>40</td><td>0.1</td><td>1163</td><td>2.5</td></tr><tr><td>KeaE1</td><td>0.28 + 0.10</td><td>0.11 + 0.01</td><td>0.35 + 0.16</td><td><strong>0.26</strong> + <strong>0.04</strong></td><td>9.3</td><td>0.52 + 0.00</td><td>265 + 0</td><td>1415 + 0</td><td>1.6 + 0.0</td><td>185</td><td>50</td><td>2</td><td>1160</td><td>2.6</td></tr><tr><td colspan="16"><strong>Moussallam et al., 2019</strong></td></tr><tr><td>AO2</td><td><strong>0.46</strong> + <strong>0.27</strong></td><td><strong>0.18</strong></td><td></td><td></td><td>16.0</td><td>1.32 + 0.00</td><td></td><td></td><td>0.1 + 0.1</td><td>110</td><td>50</td><td>0.1</td><td>1150</td><td>3.3</td></tr><tr><td>AO13</td><td><strong>0.12</strong> + <strong>0.03</strong></td><td></td><td></td><td></td><td>4.2</td><td>1.15 + 0.00</td><td></td><td></td><td>1.6 + 0.0</td><td>253</td><td>30</td><td>0.1</td><td>1150</td><td>3.3</td></tr><tr><td>AO38</td><td>0.17 + 0.03</td><td></td><td>0.39 + 0.01</td><td><strong>0.19</strong> + <strong>0.07</strong></td><td>6.7</td><td>1.15 + 0.00</td><td></td><td>403 + 0</td><td>0.8 + 0.5</td><td>237</td><td>45</td><td>0.1</td><td>1150</td><td>3.3</td></tr><tr><td colspan="16"><strong>Lloyd et al., 2014</strong></td></tr><tr><td>E1</td><td>0.24 + 0.01</td><td>0.03 + 0.00</td><td>0.01 + 0.00</td><td><strong>0.06</strong> + <strong>0.02</strong></td><td>2.2</td><td>4.00 + 0.00</td><td>100 + 0</td><td>1500 + 0</td><td>0.1 + 0.1</td><td>3333</td><td>202</td><td>2</td><td>1030</td><td>4.4</td></tr><tr><td>E2</td><td>1.07 + 0.29</td><td>0.10 + 0.00</td><td>0.13 + 0.03</td><td><strong>0.30</strong> + <strong>0.05</strong></td><td>11.1</td><td>2.20 + 0.00</td><td>150 + 0</td><td>1900 + 0</td><td>0 + 0.2</td><td>175</td><td>54.6</td><td>2</td><td>1030</td><td>4.4</td></tr><tr><td>E3</td><td>0.92 + 0.47</td><td>1.47 + 0.03</td><td>0.59 + 0.03</td><td><strong>1.10</strong> + <strong>0.22</strong></td><td>40.9</td><td>2.80 + 0.00</td><td>160 + 0</td><td>0 + 0.0</td><td>1.6 + 0.2</td><td>69</td><td>78</td><td>2</td><td>1030</td><td>4.4</td></tr><tr><td>E4</td><td>0.35 + 0.08</td><td>0.28 + 0.08</td><td>0.13 + 0.08</td><td><strong>0.44</strong> + <strong>0.10</strong></td><td>16.3</td><td>2.80 + 0.00</td><td>170 + 0</td><td>2070 + 0</td><td>0.1 + 0.0</td><td>498</td><td>221</td><td>2</td><td>1030</td><td>4.4</td></tr><tr><td colspan="16"><strong>Myers et al., 2018</strong></td></tr><tr><td>BTF7-1_RE_no_2</td><td><strong>0.021</strong> + <strong>0.402</strong></td><td><strong>0.801</strong></td><td></td><td></td><td>0.8</td><td>3.70 + 0.10</td><td></td><td></td><td>0.4 + 0.2</td><td>3238</td><td>78</td><td>10</td><td>740</td><td>8.5</td></tr><tr><td>BTF_F2-5_RE_no_1</td><td><strong>0.041</strong> + <strong>0.511</strong></td><td><strong>0.810</strong></td><td></td><td></td><td>1.6</td><td>3.80 + 0.05</td><td></td><td></td><td>1.6 + 0.2</td><td>1014</td><td>77</td><td>35</td><td>740</td><td>8.5</td></tr><tr><td>BTF8-2_RE_no_1</td><td><strong>0.273</strong> + <strong>0.861</strong></td><td><strong>0.561</strong></td><td></td><td></td><td>10.7</td><td>5.00 + 0.00</td><td></td><td></td><td>0.1 + 0.1</td><td>264</td><td>131</td><td>59</td><td>740</td><td>8.5</td></tr><tr><td>P1963-6_RE_no_1</td><td>0.007 + 0.502</td><td>0.0001 + 0.0004</td><td></td><td><strong>0.002</strong> + <strong>0.002</strong></td><td>0.1</td><td>1.80 + 0.00</td><td>90 + 0</td><td>0</td><td>0 + 0.0</td><td>3846</td><td>38</td><td>13</td><td>780</td><td>8</td></tr><tr><td>P2305-F_RE_no_1</td><td><strong>0.020</strong> + <strong>0.800</strong></td><td><strong>0.518</strong></td><td></td><td></td><td>0.8</td><td>3.40 + 0.10</td><td></td><td></td><td>3.2 + 0.3</td><td>900</td><td>68</td><td>50</td><td>780</td><td>8</td></tr><tr><td>P1970-A_RE_no_6</td><td>0.068 + 0.505</td><td>0.010 + 0.004</td><td></td><td><strong>0.033</strong> + <strong>0.072</strong></td><td>1.3</td><td>3.40 + 0.00</td><td>90 + 0</td><td>10</td><td>0.1 + 0.1</td><td>515</td><td>89</td><td>54</td><td>780</td><td>8</td></tr><tr><td>MM4_RE_no_12</td><td>0.009 + 0.504</td><td>0.001 + 0.000</td><td></td><td><strong>0.005</strong> + <strong>0.000</strong></td><td>0.2</td><td>2.00 + 0.00</td><td>390 + 0</td><td>0</td><td>0 + 0.0</td><td>4934</td><td>85</td><td>40</td><td>800</td><td>9</td></tr><tr><td>MM7_RE_no_10</td><td>0.099 + 0.503</td><td>0.007 + 0.000</td><td></td><td><strong>0.013</strong> + <strong>0.006</strong></td><td>0.5</td><td>2.50 + 0.00</td><td>170 + 0</td><td>0</td><td>3.2 + 0.2</td><td>467</td><td>71</td><td>25</td><td>800</td><td>9</td></tr><tr><td>MM4_RE_no_13</td><td>0.069 + 0.503</td><td>0.152 + 0.086</td><td></td><td><strong>0.105</strong> + <strong>0.081</strong></td><td>4.1</td><td>2.80 + 0.00</td><td>290 + 0</td><td>0</td><td>1.6 + 0.1</td><td>1272</td><td>98</td><td>10</td><td>800</td><td>9</td></tr><tr><td colspan="16"><strong>Humphreys et al., 2008</strong></td></tr><tr><td>KV518b-1</td><td><strong>0.16</strong> + <strong>0.00</strong></td><td><strong>0.80</strong></td><td></td><td></td><td>6.4</td><td>4.60 + 1.66</td><td></td><td></td><td>0</td><td>662</td><td>137</td><td>33</td><td>880</td><td>4.9</td></tr><tr><td>MSH1-3</td><td><strong>0.23</strong> + <strong>0.00</strong></td><td><strong>0.80</strong></td><td></td><td></td><td>9.3</td><td>6.50 + 0.00</td><td></td><td></td><td>0</td><td>481</td><td>137</td><td>27</td><td>880</td><td>4.9</td></tr><tr><td>MSH1-6</td><td><strong>0.82</strong> + <strong>0.00</strong></td><td><strong>0.84</strong></td><td></td><td></td><td>33.3</td><td>6.50 + 0.00</td><td></td><td></td><td>0</td><td>157</td><td>137</td><td>9</td><td>880</td><td>4.9</td></tr><tr><td colspan="16"><strong>Moussallam et al. (2021)</strong></td></tr><tr><td>AF2</td><td>2.64 + 1.55</td><td>5.87 + 0.14</td><td>2.96 + 0.00</td><td><strong>2.64</strong> + <strong>1.55</strong></td><td>99.5</td><td>1.48 + 0.03</td><td></td><td>220 + 10</td><td>3.2 + 0.0</td><td>26</td><td>68</td><td>0.1</td><td>1110</td><td>4</td></tr><tr><td>AD5</td><td>0.49 + 0.12</td><td>0.14 + 0.04</td><td>0.96 + 0.00</td><td><strong>0.53</strong> + <strong>0.14</strong></td><td>19.9</td><td>1.05 + 0.05</td><td></td><td>950 + 50</td><td>0.1 + 0.3</td><td>64</td><td>34</td><td>0.1</td><td>1110</td><td>4</td></tr><tr><td>AE38</td><td>0.11 + 0.03</td><td>0.01 + 0.00</td><td>0.00 + 0.00</td><td><strong>0.11</strong> + <strong>0.03</strong></td><td>4.2</td><td>1.05 + 0.00</td><td></td><td>950 + 233</td><td>0 + 0.2</td><td>335</td><td>37</td><td>0.1</td><td>1110</td><td>4</td></tr><tr><td>PG11</td><td>1.70 + 0.87</td><td>0.10 + 0.00</td><td>0.00 + 0.00</td><td><strong>1.70</strong> + <strong>0.86</strong></td><td>64.1</td><td>1.80 + 0.20</td><td></td><td>797 + 33</td><td>1.6 + 1.4</td><td>68</td><td>115</td><td>0.1</td><td>1110</td><td>4</td></tr><tr><td colspan="16"><strong>Myers et al., 2021</strong></td></tr><tr><td>18R1F</td><td><strong>1.70</strong> + <strong>0.86</strong></td><td><strong>0.75</strong></td><td></td><td></td><td>66.7</td><td>4.70 + 0.10</td><td></td><td></td><td>0</td><td>49</td><td>143</td><td>60</td><td>850</td><td>7.4</td></tr><tr><td>18R10</td><td><strong>1.70</strong> + <strong>0.56</strong></td><td><strong>0.40</strong></td><td></td><td></td><td>66.7</td><td>4.30 + 0.00</td><td></td><td></td><td>0</td><td>55</td><td>119</td><td>26</td><td>850</td><td>7.4</td></tr><tr><td>18R20</td><td><strong>0.43</strong> + <strong>0.85</strong></td><td><strong>0.85</strong></td><td></td><td></td><td>16.9</td><td>4.60 + 0.00</td><td></td><td></td><td>0</td><td>274</td><td>120</td><td>2</td><td>850</td><td>7.4</td></tr><tr><td>20R10</td><td><strong>0.47</strong> + <strong>0.81</strong></td><td><strong>0.81</strong></td><td></td><td></td><td>18.4</td><td>3.35 + 0.00</td><td></td><td></td><td>0</td><td>153</td><td>77</td><td>5</td><td>850</td><td>7.4</td></tr><tr><td>11R60</td><td><strong>2.93</strong> + <strong>1.31</strong></td><td></td><td></td><td></td><td>114.8</td><td>4.30 + 0.00</td><td></td><td></td><td>0</td><td>26</td><td>116</td><td>40</td><td>850</td><td>7.4</td></tr></tbody></table>
<caption>Figure 1.</caption>
k values of dPdt

<table><thead><tr><th></th><th>NErr<sub>1,1,m</sub></th><th>...</th><th>...</th><th>...</th><th>NErr<sub>1,k,m</sub></th></tr></thead><tbody><tr><th>NErr<sub>1,1,2</sub></th><td>...</td><td>...</td><td>...</td><td>...</td><th>NErr<sub>1,k,2</sub></th></tr><tr><th>NErr<sub>1,1,1</sub></th><td>...</td><td>...</td><td>...</td><td>...</td><th>NErr<sub>1,k,1</sub></th></tr><tr><th>...</th><td>...</td><td>...</td><td>...</td><td>...</td><th>...</th></tr><tr><th>...</th><td>...</td><td>...</td><td>...</td><td>...</td><th>...</th></tr><tr><th>...</th><td>...</td><td>...</td><td>...</td><td>...</td><th>...</th></tr><tr><th>NErr<sub>g,1,1</sub></th><td>...</td><td>...</td><td>...</td><td>...</td><th>NErr<sub>g,k,1</sub></th></tr></tbody></table>

m iterations

g values of C_i

STEP 1

- Mean for each NErr<sub>ij</sub> values over the m iterations

- Select the minimal NErr_m value

- Calculation of 1 sigma confidence interval over the NErr_m<sub>ab</sub> value

j=a

<table><tbody><tr><th>NErr_m<sub>1,1</sub></th><td>...</td><td>...</td><td>...</td><td>NErr_m<sub>1,k</sub></td></tr><tr><th>...</th><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><th>...</th><td>...</td><td>NErr_m<sub>a,b</sub></td><td>...</td><td>...</td></tr><tr><th>...</th><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><th>NErr_m<sub>g,1</sub></th><td>...</td><td>...</td><td>...</td><td>NErr_m<sub>g,k</sub></td></tr></tbody></table>

STEP 2

i=b

- Filtering of error values using equation 10

<table><thead><tr><th></th><th>j<sub>1</sub></th><th>j<sub>2</sub></th><th></th></tr></thead><tbody><tr><th>NErr_m<sub>1,1</sub></th><td>NErr_m<sub>1,2</sub></td><td>...</td><td>NErr_m<sub>1,k</sub></td></tr><tr><th>NErr_m<sub>2,1</sub></th><td>...</td><td>...</td><td>...</td></tr><tr><th>...</th><td>...</td><td>NErr_m<sub>i1,j1</sub></td><td>...</td></tr><tr><th>...</th><td>...</td><td>...</td><td>...</td></tr><tr><th>...</th><td>...</td><td>...</td><td>NErr_m<sub>i2,j2</sub></td></tr><tr><th>...</th><td>...</td><td>NErr_m<sub>i3,j1</sub></td><td>...</td></tr><tr><th>NErr_m<sub>g,1</sub></th><td>...</td><td>...</td><td>NErr_m<sub>g,k</sub></td></tr></tbody></table>

STEP 3

i<sub>1</sub>

i<sub>2</sub>

i<sub>3</sub>

- Extract the n values of dPdt and C<sub>i</sub> that were used to calculate all the filtered NErr_m<sub>ij</sub>

<table><thead><tr><th>from all j</th><th>from all i</th></tr></thead><tbody><tr><td>dPdt<sub>1</sub></td><td>C<sub>i1</sub></td></tr><tr><td>...</td><td>...</td></tr><tr><td>dPdt<sub>n</sub></td><td>C<sub>in</sub></td></tr></tbody></table>

STEP 4

- Result value of C<sub>i</sub> and dPdt

STEP 5

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

- 1 sigma confidence interval of C<sub>i</sub> and dPdt

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

Smallest value among the whole matrix with (a,b) as its coordinates

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

All NErr_m<sub>ij</sub> values that satisfy equation 10
<caption>Figure 2.</caption>
<figure><img src="image_5.png" /><figcaption>## Part 2: Results</figcaption></figure>
<caption>Figure 3.</caption>
<figure><img src="image_6.png" /><figcaption>a</figcaption></figure>

Comparison of H2O-estimated best fits for each M value

<figure><img src="image_7.png" /><figcaption>b</figcaption></figure>

Cumulative best fit at 0.08 MPa/s

<figure><img src="image_8.png" /><figcaption>c</figcaption></figure>

Cumulative error repartition over dPdt and M

<figure><img src="image_9.png" /><figcaption>d</figcaption></figure>

H2O Error repartition over dPdt at M=1.6 wt%

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

CO2 Error repartition over dPdt at M=3.2 wt%

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

S Error repartition over dPdt at M=0.2 wt%
<caption>Figure 4.</caption>
Decompression rate recalculated with EMBer (MPa/s)

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

Decompression rate from the litterature (MPa/s)
<caption>Figure 5.</caption>
<figure><img src="image_13.png" /></figure>

<table><thead><tr><th>number of iterations</th><th>3</th><th>11</th><th>101</th><th>501</th><th>1001</th></tr></thead><tbody><tr><th scope="row">BT251</th><td></td><td></td><td></td><td></td><td></td></tr><tr><th scope="row">Result</th><td>0.019</td><td>0.021</td><td>0.021</td><td>0.021</td><td>0.021</td></tr><tr><th scope="row">Median</th><td>0.039</td><td>0.023</td><td>0.025</td><td>0.027</td><td>0.027</td></tr><tr><th scope="row">Mean</th><td>0.042</td><td>0.022</td><td>0.026</td><td>0.027</td><td>0.028</td></tr><tr><th scope="row">Uncertainty</th><td>-0.006</td><td>-0.004</td><td>-0.006</td><td>-0.008</td><td>-0.01</td></tr></tbody><tfoot><tr><th scope="row"></th><td>+0.03</td><td>+0.004</td><td>+0.006</td><td>+0.008</td><td>+0.01</td></tr></tfoot></table>

<table><thead><tr><th>number of iterations</th><th>3</th><th>11</th><th>101</th><th>501</th><th>1001</th></tr></thead><tbody><tr><th scope="row">ReE1</th><td></td><td></td><td></td><td></td><td></td></tr><tr><th scope="row">Result</th><td>0.020</td><td>0.020</td><td>0.020</td><td>0.020</td><td>0.020</td></tr><tr><th scope="row">Median</th><td>0.020</td><td>0.020</td><td>0.030</td><td>0.020</td><td>0.030</td></tr><tr><th scope="row">Mean</th><td>0.022</td><td>0.023</td><td>0.027</td><td>0.023</td><td>0.028</td></tr><tr><th scope="row">Uncertainty</th><td>-0</td><td>-0</td><td>-0.01</td><td>-0</td><td>-0.01</td></tr></tbody><tfoot><tr><th scope="row"></th><td>+0.01</td><td>+0.02</td><td>+0.02</td><td>+0.03</td><td>+0.02</td></tr></tfoot></table>

<table><thead><tr><th>number of iterations</th><th>3</th><th>11</th><th>101</th><th>501</th><th>1001</th></tr></thead><tbody><tr><th scope="row">AO02</th><td></td><td></td><td></td><td></td><td></td></tr><tr><th scope="row">Result</th><td>0.218</td><td>0.218</td><td>0.159</td><td>0.159</td><td>0.159</td></tr><tr><th scope="row">Mean</th><td>0.248</td><td>0.218</td><td>0.278</td><td>0.218</td><td>0.218</td></tr><tr><th scope="row">Median</th><td>0.236</td><td>0.250</td><td>0.498</td><td>0.287</td><td>0.326</td></tr><tr><th scope="row">Uncertainty</th><td>-0.059</td><td>-0.059</td><td>-0.059</td><td>-0.059</td><td>-0.059</td></tr></tbody><tfoot><tr><th scope="row"></th><td>+0.089</td><td>+0.202</td><td>+0.296</td><td>+0.178</td><td>+0.237</td></tr></tfoot></table>
<caption>Figure 6.</caption>
M₀ study

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

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

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

P<sub>end</sub> study

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

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

<figure><img src="image_19.png" /></figure>
<caption>Figure 7.</caption>
<figure><img src="image_20.png" /></figure>

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