Magma Decompression Rate Calculations With EMBER: A User‐Friendly Software to Model Diffusion of H2O, CO2, and S in Melt Embayments

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 EMBayment‐Estimated Rates (EMBER), 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 H2O, CO2, 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 reprocessed “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.

(also called reentrants). 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 & 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., Ferguson et al., 2016;Lloyd et al., 2014;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 endeavor. 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 (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 2 O and CO 2 concentrations until the fragmentation pressure is reached. The process was repeated several GEORGEAIS ET AL. Notes. 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 Zhang et al. (2007). The only two exceptions are the CO 2 diffusivity from Liu et al. (2007) which is calculated from Behrens and Zhang (2001) and the S diffusivity from Ferguson et al. (2016), which is calculated from Freda et al. (2005). * Unconfirmed use of a grid search. a Zhang and Behrens (2000). b Nowak and Behrens (1997). c Zhang et al. (2010). d Ni and Zhang (2018). e Freda et al. (2003).

Table 1 Comparison of Previously Published Codes
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 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. Lloyd et al. (2014) 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 Ferguson et al. (2016) 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 modeled profile and made the grid search more complex and the result better constrained (Ferguson et al., 2016). Another study subsequently build 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 . 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 (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 2 O and S and exsolved gas content (Moussallam et al., 2021).
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.
The aim of this article is to provide the community with a user-friendly and cross-operating system MAT-LAB 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 1980s 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 & Costa, 2013), which can extract timescales from diffusion zoning in olivine crystals, and is also a widely distributed MATLAB program.

Diffusion Model
Volatile element diffusion in embayments can be regarded as a 1D process, because of their elongated, tubelike 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., 2019Moussallam et al., , 2021Myers et al., 2018Myers et al., , 2021Newcombe 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 & 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 (Equation 1): 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, sat C , the saturation concentration calculated at the mouth of the embayment by using a solubility model, i C , the initial concentration along the embayment, and t, the time. These diffusion equations show that diffusion coefficients depend heavily on H 2 O concentration   H O 2 C , 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.
The boundary condition at the interior is defined by an absence of mass flux (Neumann condition) (Equation 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 init C 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 & Lowenstern, 2002) 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 (Newman & 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. 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 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 2 O, CO 2 , 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 solver from MATLAB.
exp 11.836 0.139 ln exp 11.924 1.003 ln * exp RT where H O 2 C is water concentration in wt. %; R, the gas constant in J mol −1 K −1 and T, the temperature in Kelvin (Freda et al., 2003). For melt with intermediate composition (55 < SiO 2 < 70 wt. %), the water diffusion coefficient ( H O 2 t D in  2 1 m s ) is calculated using a combination of Equations 5a-5i, from Ni and Zhang (2018), as previously used by Newcombe et al. (2020): with SI Y , the mole fraction of Si among all cations, 0 D , a diffusion parameter in  2 1 m s , K, the equilibrium constant, X, the mole fraction of each species, P in GPa, T in Kelvin, a, a dimensionless parameter, OH D the OH diffusivity in For rhyolitic melt, the water diffusion coefficient (in  2 1 µm s ) is given by, 10.661 1.772 exp 10.49 for low (≤2 wt. %) water melt, where t P is the pressure in MPa at time t (Zhang & Behrens, 2000). However, if the water content is high (from 2 to 8 wt. %), diffusion coefficient is calculated by, for basalts to rhyolites (Nowak et al., 2004;Zhang et al., 2007). The program also calculates diffusion coef- (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).

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 favors 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, where j u 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 Cs _ max ij ), i denotes discrete points along the distance to the embayment mouth, ij C is the calculated concentration, and Cs 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 j u . The term j u 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 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,  1 j u .
EMBER assesses statistical variations of NErr 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 NErr computed by comparing a calculated profile against m numbers of the randomly generated profiles 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−σ confidence interval for a parameter set j (Figure 1, Step 2).

Propagation of Errors
The uncertainties on the fit parameters   0 Ci, dPdt and to a lesser extent M 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): Here, a is the value for which Nerr_m is the smallest, and   1 2 and respectively the lower and upper boundaries of the confidence interval, calculated using the prctile MATLAB function. For every j calculated profile (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.

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 Part 1, panels a, b, c; Figure 2), and the results section dedicated to displaying the results once calculation is terminated (Part 2, panels a, b; Figure 2). Panel 2c displays a process log during the calculations. Figure 3 shows the result of a simulation using the volatile profiles measured GEORGEAIS ET AL.  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 three 3D-plots, of the evolution of the calculative error for each volatile (1) versus the decompression rate and initial concentration, (2) the evolution of the decompression rate versus exsolved volatile content and initial volatile concentration for each studied volatile specie, and (3)   in an embayment from the study of Ferguson et al. (2016) as an example. Unticking the "Weighting and scaling of error" checkbox will remove the scaling of NErr. The j u 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.

Outputs of the Calculation
Figures generated by EMBER have two purposes: (a) to display the results of the calculation and (b) 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; Figure 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.

A Priori Requirements
As with all computational software, meaningful results in EMBER will only be achieved with an appropriate data set. 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 & 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 the absence of such a 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 reanalyzed 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), 1974Fuego (Lloyd et al., 2014), 1500, 1650and 1959Kī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 , Late Bronze Age 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 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 reanalyzed 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 (Ferguson et al., 2016), M 0 = 0.4, 0, 0, 3.2 wt. % respectively for AF2, AD5, AE38, and PG11 (Moussallam et al., 2021), and M 0 = 0 wt. % for all embayments from Ambae . For all other studies of which the exsolved gas content was not specified we used M 0 = 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 1980s 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 reanalysis and earlier decompression rates estimates come from (a) 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 (Lloyd et al., 2014)), (b) 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 (c) the weighting and scaling of error associated with each volatile diffusion profile in all our calculations (Figure 4).

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 ( Figure 5) and a stability of the GEORGEAIS ET AL.  (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. 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), modeled 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. GEORGEAIS ET AL.  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.

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 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 (Figures 3a and 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 (Figure 6a). This variation is most pronounced on H 2 O profiles in embayments of basaltic composition (Figures 6a-6c). 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 (Figure 6d) (which must be determined from measured data) and the critical influence of the modeled degassing path. GEORGEAIS ET AL.

10.1029/2020GC009542
12 of 20 Figure 6. Evolution of modeled relative decompression rates from H 2 O (a), CO 2 (b), and S (c) as a function of the exsolved volatile content for a selection of embayments. Each case in (a-c) is calculated with a fixed P end , reported in Table 2. For figures (d-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 P end 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 .
The best fit M 0 values are directly affected by the choice of P end (Figure 6d) 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 2 O 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 path associated with at 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 (Figure 6d). Therefore, we advise the user to be cautious when choosing P end for decompression rate modeling, especially in case constrained by one volatile species (Figures 6e and 6f).

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 (  1 j u ) to all three volatile species, except on a few cases. By doing so, the best fit NErr j calculated from H 2 O profiles is several orders of magnitude higher than the best fit NErr j 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 NErr j value, making the constrains brought by CO 2 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 specie. 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 not to 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.

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 modeled profiles of each embayments can be found in the supplementary materials. We restricted the calculations to three profiles per eruption for Myers et al. (2018), including the ones associated with the highest and lowest decompression rate recorded, and to five profiles for Myers et al. (2021).

The Impact of Calculation Strategies
Our hope with developing EMBER is to minimize differences in modeling 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 allows 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  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, 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 Humphreys et al. (2008) for Mt St Helens calculations.
In this study, we recalculated all embayment profiles using the same strategy throughout, but we do not pass judgment on the validity of one strategy over another. Our aim is to present a data set, which is as comparable as possible. The data presented in Figure 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 Myers et al. (2021) with our protocol leads to decompression rates 3-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 Versus Eruption Parameters
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 Volcanic Explosivity Index, as an approximation of eruption explosivity, is related to eruption magnitude and/or plume height (Newhall & 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 (Figure 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 2 = 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., Cassidy et al., 2015;Gonnermann & Manga, 2007;Martel et al., 1998). 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 2 = 0.86). The calculated decompression rates for rhyolitic eruptions on the other hand show 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., Mastin et al., 2009;Pyle, 2015;Woods & Wohletz, 1991). 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 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 2 = 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 EM-BER reproduced existing results using the parameters from the original studies, and second, 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 data set shows no significant correlation between magma decompression rate and eruption magnitude when considering the entire data set 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 2 = 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 2 = 0.86). Additionally, there is no significant correlation between decompression rate and plume height when considering the entire data set. 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 2 = 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. data of the Mt St Helens embayments. Guillaume Georgeais was supported by a PhD fellowship from the French Government "Ministère de l'Enseignement Supérieur, de la Recherche et de l'Innovation." Estelle F. Rose-Koga acknowledges partial funding from Laboratory of Excellence initiative n°ANR-10-LABX-0006, the Région Auvergne and the European Regional Development Funds. Yves Moussallam acknowledges funding from INSU and the Région Auvergne Rhone Alpes. This is Laboratory of Excellence ClerVolc contribution number 463. The authors thank Madison Myers for sharing her code during the review process of our manuscript. The authors would like to thank Madison Myers and an Anonymous reviewer for their comments on the original manuscript and Marie Edmonds for editorial handling.