Guillaume Georgeais
email: guillaume.georgeais@uca.fr

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

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 H2O, CO2 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 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 2 O, CO 2 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.

INTRODUCTION

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., [START_REF] Aki | Deep volcanic tremor and magma ascent mechanism under Kilauea, Hawaii[END_REF][START_REF] Scandone | Magma supply, magma discharge and readjustment of the feeding system of mount St. Helens during 1980[END_REF][START_REF] Klein | The seismicity of Kilauea's magma system[END_REF][START_REF] Tryggvason | Surface deformation at the Krafla volcano, North Iceland, 1982-1992[END_REF]. 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., [START_REF] Carey | The May 18, 1980 eruption of Mount St. Helens: 2. Modeling of dynamics of the Plinian Phase[END_REF][START_REF] Rutherford | Magma ascent rates from amphibole breakdown: An experimental study applied to the 1980-1986 Mount St. Helens eruptions[END_REF][START_REF] Geschwind | Crystallization of microlites during magma ascent: the fluid mechanics of 1980-1986 eruptions at Mount St Helens[END_REF][START_REF] Browne | The influence of magma ascent path on the texture, mineralogy, and formation of hornblende reaction rims[END_REF]. Another example is the use of diffusion profiles (e.g., Fe-Mg or H 2 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., [START_REF] Mackwell | Diffusion of hydrogen in olivine: Implications for water in the mantle[END_REF][START_REF] Klügel | Rates of magma ascent and depths of magma reservoirs beneath La Palma (Canary Islands)[END_REF][START_REF] Klügel | Reactions between mantle xenoliths and host magma beneath La Palma (Canary Islands): constraints on magma ascent rates and crustal reservoirs[END_REF][START_REF] Kohlstedt | Diffusion of hydrogen and intrinsic point defects in olivine[END_REF][START_REF] Voyer | Zonation of H2O and F Concentrations around Melt Inclusions in Olivines[END_REF]. Knowing the diffusion coefficient of the species of interest, one can derive an estimate of the magma ascent rate for these eruptions [e.g., [START_REF] Demouchy | Mechanisms of hydrogen incorporation and diffusion in ironbearing olivine[END_REF][START_REF] Sparks | Dynamical constraints on kimberlite volcanism[END_REF][START_REF] Rutherford | Magma Ascent Rates[END_REF]. 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 [START_REF] Toramaru | Vesiculation process and bubble size distributions in ascending magmas with constant velocities[END_REF][START_REF] Toramaru | Numerical study of nucleation and growth of bubbles in viscous magmas[END_REF][START_REF] Toramaru | BND (bubble number density) decompression rate meter for explosive volcanic eruptions[END_REF]]. This theoretical relationship has been reproduced experimentally but is not applicable to natural samples, once bubbles start to coalesce [e.g., [START_REF] Martel | Timescales of bubble coalescence, outgassing, and foam collapse in decompressed rhyolitic melts[END_REF].

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 crystalhosted elongated melt pockets of various shapes and sizes, opened to the outside melt [e.g., [START_REF] Anderson | Hourglass inclusions: Theory and application to the Bishop Rhyolitic Tuff[END_REF]. Their formation mechanism is similar to that of melt inclusions [e.g., [START_REF] Faure | Experimental investigation of equilibration conditions during forsterite growth and melt inclusion formation[END_REF], 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., [START_REF] Liu | Melt pockets in phenocrysts and decompression rates of silicic magmas before fragmentation[END_REF][START_REF] Myers | Ascent rates of rhyolitic magma at the onset of three caldera-forming eruptions[END_REF], in plagioclases [e.g., [START_REF] Humphreys | Magma ascent rates in explosive eruptions: Constraints from H2O diffusion in melt inclusions[END_REF] and olivine crystals [e.g., [START_REF] Lloyd | NanoSIMS results from olivine-hosted melt embayments: Magma ascent rate during explosive basaltic eruptions[END_REF][START_REF] Ferguson | Magma decompression rates during explosive eruptions of Kīlauea volcano, Hawaii, recorded by melt embayments[END_REF][START_REF] Moussallam | Fast ascent rate during the 2017-2018 Plinian eruption of Ambae (Aoba) volcano: a petrological investigation[END_REF], 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.

Existing embayment volatile diffusion models

Several such numerical models have been developed in the last decades (Table 1). The first published code was developed in FORTRAN 77 [START_REF] Liu | Melt pockets in phenocrysts and decompression rates of silicic magmas before fragmentation[END_REF]. Assuming a certain temperature, pressure, initial concentration, exsolved gas, degassing path and decompression rate, the code generated time-dependent profiles of H 2 O and CO 2 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 [START_REF] Humphreys | Magma ascent rates in explosive eruptions: Constraints from H2O diffusion in melt inclusions[END_REF] using the COMSOL multiphysics software to model only H 2 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. [START_REF] Lloyd | NanoSIMS results from olivine-hosted melt embayments: Magma ascent rate during explosive basaltic eruptions[END_REF] developed a model calculating simultaneously, for a range of decompression rates, three volatile element profiles at once: H 2 O, CO 2 and S. In a following study, an improved development of a MATLAB code by [START_REF] Ferguson | Magma decompression rates during explosive eruptions of Kīlauea volcano, Hawaii, recorded by melt embayments[END_REF] took H 2 O, CO 2 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 0 ) for basaltic compositions. The addition of this new parameter: (M 0 ), the preexisting (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 [START_REF] Ferguson | Magma decompression rates during explosive eruptions of Kīlauea volcano, Hawaii, recorded by melt embayments[END_REF]. Another study subsequently build the FORTRAN 77 model from [START_REF] Liu | Melt pockets in phenocrysts and decompression rates of silicic magmas before fragmentation[END_REF] in MATLAB with updated diffusion coefficients for rhyolitic compositions and a best fit search algorithm [START_REF] Myers | Ascent rates of rhyolitic magma at the onset of three caldera-forming eruptions[END_REF]. This code was later updated taking the pressure at which degassing stops as a free parameter [START_REF] Myers | Evolution of magma decompression and discharge during a Plinian event (Late Bronze-Age eruption, Santorini) from multiple eruption-intensity proxies[END_REF]. 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 [START_REF] Moussallam | Fast ascent rate during the 2017-2018 Plinian eruption of Ambae (Aoba) volcano: a petrological investigation[END_REF]. One of the most recent model is tuned to intermediate magma compositions, contains a specific S solubility relation and a general H 2 O diffusion coefficient relation for intermediate magma, and looks for the decompression rate and initial pressure [START_REF] Newcombe | Water-in-olivine magma ascent chronometry: Every crystal is a clock[END_REF]. The latest model to date is written in Matlab and uses grid searches to find decompression rate, initial concentration of H 2 O and S and exsolved gas content [START_REF] Moussallam | How to turn off a lava lake? A petrological investigation of the 2018 intra-caldera and submarine eruptions of Ambrym volcano[END_REF].

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 [START_REF] Moussallam | Fast ascent rate during the 2017-2018 Plinian eruption of Ambae (Aoba) volcano: a petrological investigation[END_REF], 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.

The aim of this article is to provide the community with a user-friendly and crossoperating 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 [START_REF] Girona | DIPRA: A user-friendly program to model multi-element diffusion in olivine with applications to timescales of magmatic processes[END_REF], which can extract timescales from diffusion zoning in olivine crystals, and is also a widely distributed MATLAB program.

Code architecture

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 [START_REF] Ferguson | Magma decompression rates during explosive eruptions of Kīlauea volcano, Hawaii, recorded by melt embayments[END_REF][START_REF] Humphreys | Magma ascent rates in explosive eruptions: Constraints from H2O diffusion in melt inclusions[END_REF][START_REF] Liu | Melt pockets in phenocrysts and decompression rates of silicic magmas before fragmentation[END_REF][START_REF] Lloyd | NanoSIMS results from olivine-hosted melt embayments: Magma ascent rate during explosive basaltic eruptions[END_REF][START_REF] Moussallam | Fast ascent rate during the 2017-2018 Plinian eruption of Ambae (Aoba) volcano: a petrological investigation[END_REF][START_REF] Moussallam | How to turn off a lava lake? A petrological investigation of the 2018 intra-caldera and submarine eruptions of Ambrym volcano[END_REF][START_REF] Myers | Ascent rates of rhyolitic magma at the onset of three caldera-forming eruptions[END_REF][START_REF] Myers | Evolution of magma decompression and discharge during a Plinian event (Late Bronze-Age eruption, Santorini) from multiple eruption-intensity proxies[END_REF][START_REF] Newcombe | Water-in-olivine magma ascent chronometry: Every crystal is a clock[END_REF]. 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 [START_REF] Degraffenried | Modeling diffusion in 1D within melt embayments: correcting for 3D geometry[END_REF]. 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): where 𝐷 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, 𝐶 𝑠𝑎𝑡 , the saturation concentration calculated at the mouth of the embayment by using a solubility model, 𝐶 𝑖 , the initial concentration along the embayment and t, the time. These diffusion equations show that diffusion coefficients depend heavily on H 2 O concentration(𝐶 𝐻 2 𝑂 ), 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.

𝜕𝐶 𝜕𝑡 = 𝜕 𝜕𝑥 (𝐷(𝑥, 𝑡) 𝜕 𝐶 𝜕𝑥 ) (1)
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 𝐶 𝑖𝑛𝑖𝑡 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 [START_REF] Witham | SolEx: A model for mixed COHSCl-volatile solubilities and exsolved gas compositions in basalt[END_REF], and VolatileCalc [START_REF] Newman | VolatileCalc: a silicate melt-H2O-CO2 solution model written in Visual Basic for excel[END_REF] depending on magma types.

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 [START_REF] Newman | VolatileCalc: a silicate melt-H2O-CO2 solution model written in Visual Basic for excel[END_REF], Solex [START_REF] Witham | SolEx: A model for mixed COHSCl-volatile solubilities and exsolved gas compositions in basalt[END_REF] 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.

where 𝐶 𝐻 2 𝑂 is water concentration in wt. %; R, the gas constant in J.mol -1 .K -1 and T, the temperature in Kelvin [START_REF] Freda | Water diffusion in natural potassic melts[END_REF]. For melt with intermediate composition (55 < SiO 2 < 70 wt. %), the water diffusion coefficient (𝐷 𝐻 2 𝑂 𝑡 in 𝑚 2 . 𝑠 -1 ) is calculated using a combination of equations (5a-i), from [START_REF] Ni | A general model of water diffusivity in calc-alkaline silicate melts and glasses[END_REF], as previously used by [START_REF] Newcombe | Water-in-olivine magma ascent chronometry: Every crystal is a clock[END_REF]:

𝐷 𝐻 2 𝑂 𝑡 = 𝐷 𝐻 2 𝑂 𝑚 (1 - 𝑑𝑋 𝑂𝐻 2𝑑𝑋 ) + 𝐷 𝑂𝐻 𝑑𝑋 𝑂𝐻 2𝑑𝑋 , ( 5a
)
𝐷 𝐻 2 𝑂 𝑚 = 𝐷 0 exp (𝑎𝑋) , ( 5b
)
𝐷 𝑂𝐻 = 𝑐𝑜𝑛𝑠𝑡 , ( 5c
) 𝑑𝑋 𝑂𝐻 2𝑑𝑋 = 1-2𝑋 √4𝑋(𝑋-1)(1- 4 𝐾 )+1 , ( 5d
) 𝑙𝑛𝐾 = 2.6𝑌 𝑆𝐼 - 4339𝑌 𝑆𝐼 𝑇 , (5e)
for basalts to rhyolites [START_REF] Nowak | Argon and CO2 on the race track in silicate melts: A tool for the development of a CO2 speciation and diffusion model[END_REF][START_REF] Zhang | Silicate melt properties and volcanic eruptions: SILICATE MELT PROPERTIES[END_REF]. The program also calculates diffusion coefficient of S (𝐷 𝑆 ), in 𝑚 2 . 𝑠 -1 :

𝐷 𝑆 (𝑖) = exp (-8.21 - 27692-651.6 𝐶 𝐻 2 𝑂 𝑇 ) (8)
[ [START_REF] Zhang | Silicate melt properties and volcanic eruptions: SILICATE MELT PROPERTIES[END_REF]. It is to be noted that the sulfur diffusivity determined by [START_REF] Zhang | Silicate melt properties and volcanic eruptions: SILICATE MELT PROPERTIES[END_REF] is applicable only to basaltic melts under reduced conditions [START_REF] Zhang | Silicate melt properties and volcanic eruptions: SILICATE MELT PROPERTIES[END_REF].

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 (𝑁𝐸𝑟𝑟 𝑗 ) for each 𝑗 of N sets of parameters,

𝑁𝐸𝑟𝑟 𝑗 = ∑ 1 𝑢 𝑗 2 𝑖 (𝐶 𝑖𝑗 -𝐶𝑠 𝑖𝑗 ) 2 (9)
where 𝑢 𝑗 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 𝐶𝑠 𝑖𝑗 _𝑚𝑎𝑥), 𝑖 denotes discrete points along the distance to the embayment mouth, 𝐶 𝑖𝑗 is the calculated concentration and 𝐶𝑠 𝑖𝑗 is the measured concentration. Among the generated profiles, the one that gives the minimum 𝑁𝐸𝑟𝑟 𝑗 is the best fit. One may calculate 𝑁𝐸𝑟𝑟 𝑗 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 𝑁𝐸𝑟𝑟 𝑗 values without dimensions thanks to 𝑢 𝑗 . The term 𝑢 𝑗 effectively constrains the "weight" of three diffusion profiles (H 2 O, CO 2 , and S) as measurements with better uncertainty (in relative term) have a stronger influence to the sum total of three 𝑁𝐸𝑟𝑟. EMBER then scales the 𝑁𝐸𝑟𝑟 values between each volatiles before adding them. If one chooses not to weight and scale each volatile species' error, 𝑢 𝑗 = 1.

EMBER assesses statistical variations of 𝑁𝐸𝑟𝑟 by a Monte Carlo simulation accounting for the uncertainty of the measurement C ± 𝜎, 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 𝑁𝐸𝑟𝑟 computed by comparing a calculated profile against m numbers of the randomly generated profiles 𝐶𝑠 𝑖 + 𝜀 (Figure 1, Step 1). A similar approach is taken for the uncertainty of 𝑥. The result of this simulation is represented by a mean (𝑁𝑒𝑟𝑟_𝑚) and its related 1-σ confidence interval for a parameter set j (Figure 1,

Step 2).

Propagation of errors

The uncertainties on the fit parameters 𝐶𝑖, 𝑑𝑃𝑑𝑡 𝑎𝑛𝑑 𝑀 0 (𝑡𝑜 𝑎 𝑙𝑒𝑠𝑠𝑒𝑟 𝑒𝑥𝑡𝑒𝑛𝑡) are assessed by the statistical distributions of the parameters that are extracted from the calculated profiles 𝑦 that satisfy the following condition (Figure 1.

Step 3):

Pr(𝑎 -𝜎 1 ≤ 𝑦 ≤ 𝑎 + 𝜎 2 ) = 68% (10)
Here, 𝑎 is the value for which 𝑁𝑒𝑟𝑟_𝑚 is the smallest, and 𝜎 1 𝑎𝑛𝑑 𝜎 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 𝑁𝐸𝑟𝑟. Equation (10) states that there are cases of calculated profiles for which 𝑁𝐸𝑟𝑟 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.

Examples of model outputs

Organization of the graphical user interface

EMBER runs inside a user-friendly Graphical User Interface (GUI) upon execution of the application file. The GUI is separated in two parts: the input section, dedicated to generating the grid search and decompression conditions (left hand side Part1, panels a, b, c; Fig. 2), and the results section dedicated to displaying the results once calculation is terminated (Part 2, panels a, b; Fig. 2). Panel 2c display a process log during the calculations.

Figure 3 shows the result of a simulation using the volatile profiles measured in an embayment from the study of [START_REF] Ferguson | Magma decompression rates during explosive eruptions of Kīlauea volcano, Hawaii, recorded by melt embayments[END_REF] as an example. Unticking the "Weighting and scaling of error" checkbox will remove the scaling of NErr. The 𝑢 𝑗 parameter from equation ( 8) is set to 1 in this case and NErr is not scaled before being added for the cumulative error calculation. In such case, as all concentration are treated as weight percent internally, the H 2 O profile will usually weight more than the others as it usually records much larger variations in absolute concentration.

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 0 or C i . 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 0 ). Figure 3d shows the influence of the initial concentration on each volatile best fit error calculation.

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 i and M 0 , and a copy of the results for each volatile with the respective exsolved gas content.

Calculation and performance test

A priori requirements

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 i and the starting pressure of ascent (P Start ). In absence of such condition, we recommend the user to determine P Start , C i and T from melt inclusions and geo-thermobarometry studies. Finally, volatile concentration measurements on surrounding glass are required to assess P end , 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.

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 [START_REF] Humphreys | Magma ascent rates in explosive eruptions: Constraints from H2O diffusion in melt inclusions[END_REF], 1974 Fuego [START_REF] Lloyd | NanoSIMS results from olivine-hosted melt embayments: Magma ascent rate during explosive basaltic eruptions[END_REF], 1500, 1650 and 1959

Kīlauea [START_REF] Ferguson | Magma decompression rates during explosive eruptions of Kīlauea volcano, Hawaii, recorded by melt embayments[END_REF], 27 ka Taupo [START_REF] Myers | Ascent rates of rhyolitic magma at the onset of three caldera-forming eruptions[END_REF], 767 ka Long Valley [START_REF] Myers | Ascent rates of rhyolitic magma at the onset of three caldera-forming eruptions[END_REF], 2 Ma Yellowstone [START_REF] Myers | Ascent rates of rhyolitic magma at the onset of three caldera-forming eruptions[END_REF], 2017-2018 Ambae/Aoba [Moussallam et al., 2019], Late Bronze Age (LBA) Santorini [START_REF] Myers | Evolution of magma decompression and discharge during a Plinian event (Late Bronze-Age eruption, Santorini) from multiple eruption-intensity proxies[END_REF] and December 2018

Ambrym [START_REF] Moussallam | How to turn off a lava lake? A petrological investigation of the 2018 intra-caldera and submarine eruptions of Ambrym volcano[END_REF] eruptions. We used the uncertainties reported in each study when provided. Also when provided, we directly used the specific grid search (e.g. [START_REF] Ferguson | Magma decompression rates during explosive eruptions of Kīlauea volcano, Hawaii, recorded by melt embayments[END_REF] provides a range of M 0 ). 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 reanalysed results were calculated using the same set of values for M 0 as the original literature data to which they are compared. In detail this was M 0 = 0.1, 1.6, 1.6, 0.8 wt. % respectively for IkiE1, ReticE1, ReticE2 and KeaE1 [START_REF] Ferguson | Magma decompression rates during explosive eruptions of Kīlauea volcano, Hawaii, recorded by melt embayments[END_REF], M 0 = 0.4, 0, 0, 3.2 wt. % respectively for AF2, AD5, AE38 and PG11 [START_REF] Moussallam | How to turn off a lava lake? A petrological investigation of the 2018 intra-caldera and submarine eruptions of Ambrym volcano[END_REF] and M 0 = 0 wt. % for all embayments from Ambae [START_REF] Moussallam | Fast ascent rate during the 2017-2018 Plinian eruption of Ambae (Aoba) volcano: a petrological investigation[END_REF]. For all other studies of which the exsolved gas content was not specified we used M 0 = 0 wt. %. In some cases [START_REF] Myers | Ascent rates of rhyolitic magma at the onset of three caldera-forming eruptions[END_REF], 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 [START_REF] Humphreys | Magma ascent rates in explosive eruptions: Constraints from H2O diffusion in melt inclusions[END_REF] 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 [START_REF] Humphreys | Magma ascent rates in explosive eruptions: Constraints from H2O diffusion in melt inclusions[END_REF] 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 2 O diffusion coefficient as well as a single episode of ascent in all our calculation (as opposed to two-step ascent; [START_REF] Lloyd | NanoSIMS results from olivine-hosted melt embayments: Magma ascent rate during explosive basaltic eruptions[END_REF]), ( 2) the calculation of diffusion coefficient (e.g., taking into account the change of diffusivity due to H 2 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).

Monte Carlo simulation

We ran tests on three samples with variable number of Monte Carlo iterations (m=3, 11, 101, 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.

How well does EMBER constrain M 0 ?

Calculation of the exsolved gas content has been introduced within the embayment method with the calculations made by [START_REF] Ferguson | Magma decompression rates during explosive eruptions of Kīlauea volcano, Hawaii, recorded by melt embayments[END_REF]. 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 . 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

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 [START_REF] Ferguson | Magma decompression rates during explosive eruptions of Kīlauea volcano, Hawaii, recorded by melt embayments[END_REF] gives equal weight (𝑢 𝑗 = 1) to all three volatile species, except on a few cases. By doing so, the best fit 𝑁𝐸𝑟𝑟 𝑗 calculated from H 2 O profiles is several orders of magnitude higher than the best fit 𝑁𝐸𝑟𝑟 𝑗 calculated for CO 2 or S profiles. That difference is due to the range of concentration variation, as H 2 O usually varies within a few wt. % and CO 2 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 2 O 𝑁𝐸𝑟𝑟 𝑗 value, making the constrains brought by CO 2 or S profiles almost negligible.

Another approach is to weight the 𝑁𝐸𝑟𝑟 𝑗 values either with analytical uncertainty [START_REF] Myers | Ascent rates of rhyolitic magma at the onset of three caldera-forming eruptions[END_REF] or maximum measured concentration in the embayment [START_REF] Newcombe | Water-in-olivine magma ascent chronometry: Every crystal is a clock[END_REF]

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 [START_REF] Myers | Ascent rates of rhyolitic magma at the onset of three caldera-forming eruptions[END_REF], including the ones associated with the highest and lowest decompression rate recorded, and to 5 profiles for [START_REF] Myers | Evolution of magma decompression and discharge during a Plinian event (Late Bronze-Age eruption, Santorini) from multiple eruption-intensity proxies[END_REF].

The impact of calculation strategies

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 [START_REF] Moussallam | Fast ascent rate during the 2017-2018 Plinian eruption of Ambae (Aoba) volcano: a petrological investigation[END_REF] 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.

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 [START_REF] Myers | Evolution of magma decompression and discharge during a Plinian event (Late Bronze-Age eruption, Santorini) from multiple eruption-intensity proxies[END_REF] for the Minoan eruption of Santorini, were a fixed initial concentration of 5.2 wt. % H 2 O (except for one embayment at 5.6 wt. % H 2 O) was used for all model calculations. Such strategy was also used by [START_REF] Humphreys | Magma ascent rates in explosive eruptions: Constraints from H2O diffusion in melt inclusions[END_REF] for Mt St Helens calculations.

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 grid search) and using the plateau values in volatile content recorded in embayment interiors as the model starting conditions (C i ).

This difference in starting assumptions can lead to significant differences in the resulting decompression rate. Recalculating the embayment data from [START_REF] Myers | Evolution of magma decompression and discharge during a Plinian event (Late Bronze-Age eruption, Santorini) from multiple eruption-intensity proxies[END_REF] 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.

Decompression rates vs eruption parameters

A recent literature compilation showed a clear relationship between magma decompression rate and whether an eruption is explosive or effusive in character [START_REF] Cassidy | Controls on explosive-effusive volcanic eruption styles[END_REF]. 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 [START_REF] Newhall | The volcanic explosivity index (VEI) an estimate of explosive magnitude for historical volcanism[END_REF]. 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 [START_REF] Ferguson | Magma decompression rates during explosive eruptions of Kīlauea volcano, Hawaii, recorded by melt embayments[END_REF][START_REF] Moussallam | Fast ascent rate during the 2017-2018 Plinian eruption of Ambae (Aoba) volcano: a petrological investigation[END_REF]. However, our compilation shows no clear correlation between decompression rate and eruption magnitude (Fig. 7a).

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 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., [START_REF] Martel | Magma storage conditions and control of eruption regime in silicic volcanoes: experimental evidence from Mt[END_REF][START_REF] Gonnermann | The Fluid Mechanics Inside a Volcano[END_REF][START_REF] Cassidy | Rapid and slow: Varying magma ascent rates as a mechanism for Vulcanian explosions[END_REF]). 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.

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., [START_REF] Woods | Dimensions and dynamics of co-ignimbrite eruption columns[END_REF][START_REF] Mastin | A multidisciplinary effort to assign realistic source parameters to models of volcanic ash-cloud transport and dispersion during eruptions[END_REF][START_REF] Pyle | Sizes of Volcanic Eruptions[END_REF]). The eruption intensity, estimated from the mass eruption rate, has been shown to be positively correlated with decompression rate [START_REF] Barth | Magma decompression rate correlates with explosivity at basaltic volcanoes -Constraints from water diffusion in olivine[END_REF][START_REF] Ferguson | Magma decompression rates during explosive eruptions of Kīlauea volcano, Hawaii, recorded by melt embayments[END_REF][START_REF] Newcombe | Water-in-olivine magma ascent chronometry: Every crystal is a clock[END_REF].

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.

Conclusion

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

 We found that decompression rate calculations are particularly sensitive to variations of M 0 especially for the H 2 O diffusion profile. Variations of P end 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 [START_REF] Humphreys | Magma ascent rates in explosive eruptions: Constraints from H2O diffusion in melt inclusions[END_REF].

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

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

Computer code availability

 Title of the software: EMBER  Developers: Guillaume Georgeais, Kenneth T. Koga, Yves Moussallam, Estelle F.

Rose-Koga

 Hardware used: EMBER was run on a computer with 4 core (4.5GHz) and 16GB of RAM  Link to the software, tutorial and test using data from [START_REF] Ferguson | Magma decompression rates during explosive eruptions of Kīlauea volcano, Hawaii, recorded by melt embayments[END_REF] DOI: 10.6084/m9.figshare.13663811

 Software required: MATLAB Runtime Environment 2019b and up https://www.mathworks.com/products/compiler/matlab-runtime.html The outputs are a decompression rate and an initial (dissolved) volatile concentration as well as their associated confidence interval. This process of calculation is repeated for each value of exsolved volatile content (M 0 ). Estimation of the M 0 value and its related uncertainty is done in a similar manner but by considering the range of M 0 variation as the 4 th dimension of the starting matrix. ReaE1 studied previously by [START_REF] Ferguson | Magma decompression rates during explosive eruptions of Kīlauea volcano, Hawaii, recorded by melt embayments[END_REF] [START_REF] Lloyd | NanoSIMS results from olivine-hosted melt embayments: Magma ascent rate during explosive basaltic eruptions[END_REF] were made using their original H 2 O diffusion coefficient. All the results displayed are calculated for M 0 =0 excepted when original studies provided a range of M 0 [START_REF] Ferguson | Magma decompression rates during explosive eruptions of Kīlauea volcano, Hawaii, recorded by melt embayments[END_REF]. EMBER calculations were made using the "weighting and scaling of error" option which is a source of difference with previously reported literature decompression rate. [START_REF] Shreve | From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics[END_REF], [START_REF] Johnston | Revised Estimates for the Volume of the Late Bronze Age Minoan Eruption, Santorini, Greece[END_REF], [START_REF] Self | The effects and consequences of very large explosive volcanic eruptions[END_REF] and [START_REF] Mason | The size and frequency of the largest explosive eruptions on Earth[END_REF] 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 [START_REF] Gardner | Eruption dynamics and magma withdrawal during the Plinian Phase of the Bishop Tuff Eruption, Long Valley Caldera[END_REF], [START_REF] Mastin | Modeling ash fall distribution from a Yellowstone supereruption[END_REF] and [START_REF] Van Eaton | Ascent dynamics of large phreatomagmatic eruption clouds: The role of microphysics[END_REF] respectively for Long Valley's Bishop Tuff, Yellowstone's Huckleberry Ridge Tuff and Taupo's Oruanui eruptions.

Figure caption

Table caption

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 i = initial concentration, M 0 = exsolved volatile content, P f = pressure of quench, gs stands for "grid search". Every study that model decompression rates from CO 2 and/or S uses the diffusion coefficients calculated from [START_REF] Zhang | Silicate melt properties and volcanic eruptions: SILICATE MELT PROPERTIES[END_REF]. The only two exceptions are the CO 2 diffusivity from [Liu et al., Tables [START_REF] Zhang | H2O diffusion in rhyolitic melts and glasses[END_REF] c = [START_REF] Zhang | Diffusion Data in Silicate Melts[END_REF] e = [START_REF] Freda | Water diffusion in natural potassic melts[END_REF]] b = [START_REF] Nowak | An experimental investigation on diffusion of water in haplogranitic melts[END_REF] d = [START_REF] Ni | A general model of water diffusivity in calc-alkaline silicate melts and glasses[END_REF] *unconfirmed use of a grid search

  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 0 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 2 O profiles in embayments of basaltic composition (Fig 6 (a), (b), (c))

  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 specie. In EMBER, the 𝑁𝐸𝑟𝑟 𝑗 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 2 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.

Figure 1 :

 1 Figure 1 : Flow chart showing the calculation steps used in EMBER to determine the parameters and associated errors resulting in the best fit of natural volatile diffusion profiles.

Figure 2 :

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

Figure 3 :

 3 Figure 3: Subset of plots generated from a calculation on the 1500CE Kīlauea reticulite

Figure 5 :Figure 6 :

 56 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

  Figure 1.

Figure 2 .

 2 Figure 2.

  Figure 3.

  Figure 4.

  Figure 5.

  Figure 6.

Figure 7 .

 7 Figure 7.

  𝑌 𝑆𝐼 , the mole fraction of Si among all cations, 𝐷 0 ,a diffusion parameter in 𝑚 2 . 𝑠 -1 , 𝐾, the equilibrium constant, 𝑋, the mole fraction of each species, P in GPa, T in Kelvin, 𝑎, a dimensionless parameter, 𝐷 𝑂𝐻 the OH diffusivity in 𝑚 2 . 𝑠 -1 , 𝐷 𝐻 2 𝑂 𝑚 , the molecular H 2 O diffusivity in 𝑚 2 . 𝑠 -1 and 𝐷 𝐻 2 𝑂 𝑡 , the total H 2 O diffusivity in 𝑚 2 . 𝑠 -1 .[START_REF] Zhang | H2O diffusion in rhyolitic melts and glasses[END_REF]. The diffusion coefficient of CO 2 (𝐷 𝐶𝑂 2 ),

		𝑙𝑛𝐷 𝑂𝐻 = -16.78 -37.428𝑌 𝑆𝐼 -	39250-27576𝑌 𝑆𝐼 𝑇	,	(5g)
	ln (	𝐷 𝑂𝐻 𝐷 0	) = -56.09 -115.93𝑌 𝑆𝐼 + 160.54√𝑌 𝑆𝐼 -	3970√𝑌 𝑆𝐼 𝑇	,	(5h)
	𝑙𝑛𝐷 0 = 8.02 -31𝑌 𝑆𝐼 + 2.348𝑃𝑌 𝑆𝐼 +	121824𝑌 𝑆𝐼 -118323√𝑌 𝑆𝐼 -(10016𝑌 𝑆𝐼 -3648)𝑃 𝑇	, (5i)
	with 10.661 𝑇	-	1.772𝑃 𝑡
						44620 T	+	57.3P t T	) -√X (0.091 +	4.77 * 10 6 T 2	)]}
							(6b)
	where 𝑚 = -20.79 -5030 𝑇 ⁄ -1.4 𝑃 𝑡 𝑇 ⁄ , and X = ( 𝐶 𝐻 2 𝑂 18.015 ⁄	)/ [ 𝐶 𝐻 2 𝑂 18.015 ⁄	+
	(100 -𝐶 𝐻 2 𝑂 ) 32.49 ⁄ ] [in 𝑚 2 . 𝑠 -1 , is given by,
	𝐷 𝐶𝑂 2 = exp [-14.34 -	17360-0.6527𝑃 𝑡 𝑇	+ (-0.7171 +	1436.8 𝑇 ) 𝐶 𝐻 2 𝑂 ]
						198508-166674𝑌 𝑆𝐼 𝑇	,	(5f)

𝑎 = -94.07 + 74.112𝑌 𝑆𝐼 + For rhyolitic melt, the water diffusion coefficient (in µ𝑚 2 . 𝑠 -1 ) is given by,

𝐷 𝐻 2 𝑂 = 𝐶 𝐻 2 𝑂 exp (10.49 -𝑇 )

(6a)

for low (≤ 2 wt. %) water melt, where P t is the pressure in MPa at time 𝑡

[START_REF] Zhang | H2O diffusion in rhyolitic melts and glasses[END_REF]

. However, if the water content is high (from 2 to 8 wt. %), diffusion coefficient is calculated by, 𝐷 𝐻 2 𝑂 = 𝑋exp (𝑚) {1 + exp [56 + 𝑚 + X (-34.1 +

Table 2

 2

. For figures

(d)

,

(e)

and (f), P ref is the P end value displayed in Table 2 for ReE2 and BT251. (d) The initial exsolved volatile content needed to generate the best fit (M Bestfit ) increases for ReE2 when changing Pend but not systematically for M413. (e) Variation of relative decompression rate modeled from the H 2 O diffusion profile only with P end . (f) Variation of relative decompression rate modeled from all studied species (H 2 O, CO 2 and S for ReE2, and H 2 O and CO 2 for M413) with P end .

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 EMBERrecalculated 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 orginal article or respectively

Table 2 861

 2

	18R1F	1.70	66.7	4.70					850
	18R1O	1.70	66.7	4.30	0	55	119	26	850
	18R2O	0.43	16.9	4.60
		0.47	18.4	3.35

Acknowledgement

We would like to thank M. Humphreys for providing the raw data of the Mt St Helens embayments. GG was supported by a PhD fellowship from the French Government "Ministère de l'Enseignement Supérieur, de la Recherche et de l'Innovation". ERK acknowledges partial funding from Laboratory of Excellence initiative n°ANR-10-LABX-0006, the Région Auvergne and the European Regional Development Funds. Y.M acknowledges funding from INSU and the Région Auvergne Rhone Alpes. This is Laboratory of Excellence ClerVolc contribution number 463. We thank Madison Myers for sharing her code during the review process of our manuscript. We would like to thank Madison Myers and an Anonymous reviewer for their comments on the original manuscript and Marie Edmonds for editorial handling.