VESIcal Part II: A critical approach to volatile solubility

17 Accurate models of H2O and CO2 solubility in silicate melts are vital for understanding volcanic 18 plumbing systems. These models are used to estimate the depths of magma storage regions from melt 19 inclusion volatile contents, investigate the role of volatile exsolution as a driver of volcanic eruptions, 20 and track the degassing paths followed by magma ascending to the surface. However, despite the large 21 increase in the number of experimental constraints over the last two decades, many recent studies still 22 utilize the earlier generation of models, which were calibrated on experimental datasets with restricted 23 compositional ranges. This may be because many of the available tools for more recent models re24 quire large numbers of input parameters to be hand-typed (e.g., temperature, concentrations of H2O, 25 CO2, and 8–14 oxides), making them difficult to implement on large datasets. Here, we use a new 26 open-source Python3 tool, VESIcal, to critically evaluate the behaviours and sensitivities of different 27 solubility models for a range of melt compositions. Using literature datasets of andesitic-dacitic ex28 perimental products and melt inclusions as case studies, we illustrate the importance of evaluating the 29 calibration dataset of each model. Finally, we highlight the limitations of particular data presentation 30 methods such as isobar diagrams, and provide suggestions for alternatives, and best practices regard31 ing the presentation and archiving of data. This review will aid the selection of the most applicable 32 solubility model for different melt compositions, and identifies areas where additional experimental 33 constraints are required. (242/250 words) 34 Plain Language Summary 35 Being able to accurately model the solubility of H2O and CO2 in magmas is very important for 36 understanding a wide variety of volcanic processes, such as the depths at which magma is stored in 37 the crust, the driving force behind volcanic eruptions, and the release of volatile elements into the 38 atmosphere. However, there has been no easy way for volcanologists to perform calculations on large 39 datasets, or to compare different models. This review uses a new, open-source tool called VESIcal 40 written in the popular programming language Python3. This allows us to compare different models for 41 a wide variety of melt compositions, temperatures and pressures, helping researchers to identify the 42 most suitable model for their study. We also suggest areas where further experimental constraints are 43 required. Finally, we highlight the limitations of particular data presentation methods such as isobar 44 diagrams, and provide suggestions for alternatives, and best practices regarding the presentation and 45 archiving of data. 46


47
The most abundant volatile components found in terrestrial magmatic systems are H 2 O and 48 CO 2 . It has been known for nearly a century (Bowen, 1928;Tuttle & Bowen, 1958) that these volatile 49 species have profound effects on the chemical and material properties of magmas (e.g., phase equilib-50 ria, melting temperatures, magma viscosity and density; Burnham, 1979  The solubility of a volatile species is defined at a given pressure and temperature as the maxi-57 mum concentration that can be dissolved within a silicate melt of a specified composition. Ignoring 58 disequilibrium effects, if the volatile content of the system exceeds this solubility limit, a separate 59 fluid/vapour phase will exsolve from the magma. In this review, we favour the term fluid because of 60 the supercritical nature of exsolved volatile phases at magmatic temperatures. In general terms, a 61 magma is described as volatile undersaturated when there is no fluid phase, and volatile saturated 62 once a fluid phase is present (also referred to as vapour undersaturated/saturated, or fluid undersatu-63 rated/saturated). In detail, different volatile species do not act as independent entities, but influence  Calibration temperature of Dixon (1997) *if normalized (not recommended), different proportions of FeO and Fe2O3 will slightly change the normalized SiO2 content  Range of calibration dataset, as authors do not state a range. We note that the vast majority of experiments were conducted at <5000 bar. 2 Authors state that most experiments were conducted between 1200-1300˚C (whole range 1100-1400˚C.  Although their empirical expressions are for pure fluids, they were mostly calibrated on mixed CO2-H2O experiments. 2 Author-suggested range 3 Note, this model contains no temperature term.  Although this model is for pure CO2, it was calibrated on mixed CO2-H2O experiments. 2 Author-suggested range. The calibration dataset spans: (SFVF:4133-6141 bar, Sunset Crater: 4071-6098 bar, Erebus: 4078-6175 bar, Vesuvius: 269-6175 bar, Etna=485-6199 bar, Stromboli=524-6080 bar). 3 Note, all calculations and experiments were performed at 1200˚C. Authors suggest applicable between 1000-1400˚C  Dixon, 1997). Interestingly, this is the only 231 model discussed here which considers more than one species for dissolved H 2 O in the melt.

232
For CO 2 solubility, Dixon (1997) adapted the model of Dixon et al. (1995) to account for the effect of melt composition, based on observations from experiments that CO 2 solubility increases from tholeiitic (49 wt% SiO 2 ) to basanitic (46 wt% SiO 2 ) to leucitic (44.1 wt% SiO 2 ) melts at 1200 • C, 1 kbar. A linear regression with CO 2 solubility was achieved using a composition parameter (Π) expressed in terms of the cation fractions, X i (Dixon, 1997): Π = −6.50(X Si 4+ + X Al 3+ ) + 20.17(X Ca 2+ + 0.8X K + + 0.7X Na Redlich-Kwong equation of state (Holloway, 1977), with the correction of Flowers (1979). 235 This simplified expression was designed to aid the investigations of volatile solubility in the suite 236 of lavas from the North Arch, where it effectively captures the observed 5× decrease in CO 2 solubility 237 from 40 to 49 wt% SiO 2 . However, this simplified parameterization became very widely used in a wide 238 variety of tectonic settings following its implementation in the Excel-based tool VolatileCalc (Newman 239 and Lowenstern, 2002). Here, we refer to this model as VolatileCalc-Basalt, to differentiate it from the 240 full Π parameterization of Dixon (1997). 241 The advantage of the Π-SiO 2 simplification is that users only have to input the concentration of T is temperature in Kelvin, b 1 -b 4 are fit parameters, and the P W and P CO2 terms account for the partial pressures of each volatile species in the co-existing fluid, with: Similarly, they provide the following expression for H 2 O: [H 2 O t ] wt% = a 1 P 0.5 w + a 2 P w + a 3 P 1.5 w T + a 4 P 1.5 w + P CO2 (a 5 P 0.5 w + a 6 P w ) ity, which is a recalibration of the earlier models of Papale (1999) and Papale (1997 The NBO/O term represents the number of non-bridging oxygens divided by oxygen, expressing the availability of oxygen to form carbonate groups within the melt. NBO/O can be calculated from mol fraction of different oxides, X i , on an anhydrous or hydrous basis: In both cases, mole fractions are calculated on a hydrous basis (Iacono-Marziano, written comms). We note for completeness that in the original publication, equation  Fig. 2).

396
The empirical nature of the fitting terms incorporating melt composition, pressure and temper-    Where P is the pressure in MPa, and Π * is a compositional parameter expressed in terms of the cation fractions of 7 species: We note for completeness that the expression provided in Shishkina   positions of all the models discussed thus far (Fig. 2 to reduce the number of interaction parameters for volatile-melt species, and seems to be a justified   To aid comparisons between models, a number of silicate melt compositions (Table 1) Dixon et al. (1995). MORB2 is the MORB composition given in Table 3     to all the comparisons thus far, the w 0    the co-existing fluid ranging from 0 (interception with the y axis) to 1 (interception with the x axis).

645
The treatment of mixed fluids differs quite considerably in each model. tents merging into a concave-down shape (Fig. 6).

736
There is also significant deviation between different models for Etna melts (Fig. 7c), which is far 737 greater than that observed for H 2 O (Fig. 5). The A-2019 model, developed specifically for the com-  Interestingly, MagmaSat also underpredicts CO 2 concentrations at a given pressure relative to     Differences in the treatment of H 2 O-CO 2 mixing for rhyolitic melts are more subtle than for 779 basaltic compositions (Fig. 9). Unlike for basalts, the differences in isobar positions mostly result from 780 large differences between the pure CO 2 solubility predicted by different models rather than treatment to MagmaSat at <4 kbar, but rapidly rises to higher CO 2 contents at higher pressures, predicting al-most as much dissolved CO 2 as VolatileCalc-Rhyolite at ∼12 kbar (Fig. 8c). For Aluto, the curvature 803 of the P-2006 model is even more prominent, predicting drastically lower CO 2 contents than all other 804 models at <6 kbar, and then rapidly rising, predicting higher CO 2 solubility than even VolatileCalc-

805
Rhyolite at >9 kbar (Fig. 8d). These large deviations between models, as well as the large errors on 806 the interaction terms for CO 2 solubility in MagmaSat ( When all solubility models are compared (4 applicable to rhyolites, 6 to basalts), there is sub-815 stantial overlap between curves calculated for MORB1 at 1200 • C and Mono Craters at 800 • C (com-816 pare Fig. 11a vs. Fig. 8a). To get around this problem of large differences between models, we com-817 pare the predictions from the three models which can be applied to both Rhyolites and Basalts: Mag-818 maSat ( Fig. 10a-b), P-2006 ( Fig. 10c-d) and VolatileCalc-Basalt and -Rhyolite (Fig. 10e- Overall, these comparisons demonstrate that at <5 kbar, the difference in solubility between 830 basalts and rhyolites is relatively subtle and easily overwhelmed by differences in predictions from 831 different solubility models, particularly given some models predict that solubility increases with tem-  for H 2 O contents between 0-1.5 wt% (Fig. 11a- (Table 1) at 1000 • C, and H2O contents between 0-6 wt%. Note that the y scale for parts a-c is significantly smaller than parts d-f.
To investigate the effect of H 2 O re-equilibration on melt inclusion saturation pressures in arcs, 907 we repeat the sensitivity test described above, using the major element composition of a Estimating the initial CO 2 contents of melt inclusions is also challenging. While the total CO 2 920 content of the inclusion is not affected by diffusive re-equilibration, CO 2 may be partitioned from 921 the melt phase into a vapour bubble. Cooling following melt inclusion entrapment is accompanied by 922 the formation of a denser mineral phase from a less dense silicate melt, and differential thermal con-  Saturation pressures in rhyolitic magmas are also very sensitive to melt CO 2 contents (Fig. 13).  The lack of consensus as to whether increasing temperature increases or decreases the solubility 1007 of H 2 O and CO 2 indicates that this effect is relatively subtle, and overwhelmed by analytical errors 1008 associated with measuring experimental products (and other sources of experimental scatter; e.g., Fig.   1009 16a-b). This makes it very difficult for empirical models to fully constrain the temperature sensitiv-  (Fig. 14). As 1020 the magnitude of these coefficients is small, the temperature effect on H 2 O solubility is small, and 1021 only visible at higher pressures (because of the P part of these terms; Fig. 14a  temperature. In both the hydrous and anhydrous models, C CO2 is positive (0.12±0.02 and 0.14±0.02 1025 respectively) and larger in magnitude than C H2O , so CO 2 solubility decreases with increasing tempera-1026 ture (see Fig. 14c). Similarly, temperature sensitivity in rhyolitic melts was evaluated by calculating isobars at 0.5 1053 and 2 kbar for 700 and 900 • C using the Mono Craters rhyolite composition. As for the basaltic exam-1054 ple, the directionality and magnitude of effect of temperature on saturation pressures for melts with 1055 volatile contents indicated by the colored stars is shown in Fig. 15c-  It is noteworthy that the temperature sensitivity of CO 2 solubility predicted by L-2005 and 1064 VolatileCalc-Rhyolite is much greater than that shown by any of the basaltic models (Fig. 14d-e   1065 vs. Fig 15d-e), and significant considering other sources of error associated with saturation pressure  (Fig. 16c vs.e). These discrepancies likely reflect this model being extrapolated towards 1159 the limits of its calibration dataset in terms of both pressure (most experiments were conducted at <5 1160 kbar) and melt composition (Fig. 17, see the next section for more discussion).

1161
-53-manuscript submitted to Earth and Space Science Error bars on all plots shows the 2σ uncertainties from measurements of volatile contents in experimental products.
Fe 3+ /FeT ratios were calculated from author-stated buffers using MELTS for Excel . shown in Fig. 18. Anhydrous molar fractions are used to calculate compositional parameters in parts c-f, because when accounting for discrepancies between isobars (e.g., on Fig. 18, the H2O content and therefore hydrous cation fraction varies as a function of the pressure).
-55- The extreme sensitivity to the Fe 3+ /Fe T ratio makes it very difficult to assess the accuracy of values than the calibration dataset ( Fig. 17d-f). While the effect of NBO/O is more convoluted be-1242 cause it also affects the solubility of H 2 O (which feeds back into the expression for CO 2 ), it is readily 1243 apparent that the positive coefficient attached to the AI term combined with the negative coefficient attached to the MgO+FeO term causes this model to predict higher CO 2 solubilities than the calibra-1245 tion dataset for the andesitic-dacitic melt inclusions considered here.

1246
The discrepancy between isobars for S-2014 and IM-2012 relative to MagmaSat are relatively 1247 similar for the Popocatépetl and Soufriére Hills melt compositions, while discrepancies for saturation 1248 pressures differ markedly (Fig. 18a-b vs. c-d). This is because the volatile contents of Popocatépetl Ethiopian Rift from Iddon and Edmonds (2020). We calculate isobars for a representative inclusion 1292 composition (BJ08 7; Fig. 19a), and then we compare these to the isobars calculated for each individ-1293 ual melt inclusion composition at 1 and 3 kbar (Fig. 19a-b). 3 kbar isobars calculated from the com-   sitions and different solubility models using a workflow similar to that used here for andesites (e.g., 1335 isobar diagrams as in Fig. 16, plots of melt composition vs. calibration datasets) will help users select 1336 a suitable model. As well as examining melt compositions, users should also evaluate whether they 1337 are extrapolating temperature-sensitive models beyond the calibration range (as discussed here for 1338 VolatileCalc-Rhyolite).

1339
In general, if a natural silicate melt composition is poorly represented by experimental data, 1340 MagmaSat is probably the best model to use, as its thermodynamic nature is more suitable to extrap-    Figure 16 shows that it is currently 1400 impossible to differentiate a potential failure in any given solubility model from anomalies in any given 1401 set of experiments (e.g., the differential effect of addition of H 2 O on CO 2 solubility in different experi-1402 ments; Fig. 16b. vs f).

1403
One of the challenges when assessing CO 2 solubility in andesitic-dacitic melts is the fact that 1404 CO 2 is present as both carbonate and molecular CO 2 . Carbon species do appear separately in FTIR 1405 spectra, but the accuracy of FTIR-derived volatile concentrations can be affected by peak over- (only yields total carbon), but may help to resolve issues with FTIR as a result of increased under-1409 standing of the optimal analysis conditions for volatiles in silicate glasses of the last few decades. on having a suite of standards with similar major element compositions and a range of volatile con-1412 tents (and these standards are often characterized by FTIR, so are subject to the caveats mentioned 1413 above). 1414 Second, the effect of redox on volatile solubility across the range encountered in terrestrial mag-1415 mas is still poorly constrained (section 5). This discrepancy largely reflects the fact that the redox 1416 conditions at which many experiments in the literature were conducted are uncertain and/or highly 1417 variable (e.g., Botcharnikov et al., 2006). Because of this uncertainty, many calibration datasets are 1418 built without being able to constrain the quantities of Fe 2 O 3 and FeO for each experimental run. and FeO proportions are accurately measured are needed to be certain that this behaviour is not real.

1423
It is also noteworthy that almost all the andesitic experiments were performed at higher oxygen fugac- provide an avenue to better constrain this parameter in future (and past) experimental products.

1429
It is also worth noting that all the models discussed here only consider the effect of redox 1430 through terms for Fe 2+ and Fe 3+ in the melt, constraining their applicability to melts more oxidising 1431 than the IW buffer. In more reducing conditions, the co-existing CO 2 -rich phase may be graphite or 1432 diamond rather than a CO 2 -rich vapour phase (Eguchi & Dasgupta, 2018), and the dissolved volatile 1433 species may be CO, CH 4 and H 2 (Mysen et al., 2009). This means that extreme caution is required 1434 when applying these solubility models to highly reducing conditions such as those found on other  to VolatileCalc (Newman & Lowenstern, 2002) and Solwcad (Papale et al., 2006) to draw extensive 1467 comparisons between the behaviour of 9 different solubility models for a range of melt compositions. 1468 We show that these models predict surprisingly different volatile solubilities, particularly for pure indicates that re-evaluation of published magma storage depths calculated using older models may be 1482 warranted. 1483 We also investigate the sensitivity of different models to variation in parameters such as H 2 O 1484 content (with relevance to diffusive re-equilibration), CO 2 content (with relevance to melt inclusion 1485 vapour bubble growth), temperature and oxygen fugacity. We suggest that by performing similar sen-1486 sitivity tests in the future, the uncertainties affecting calculations of volatile solubility in magmatic 1487 systems (and therefore the limitations of each study) can be quantified. We also demonstrate that iso-  Finally, we identify that further experimental constraints are required to accurately estimate 1496 volatile solubility in andesitic-dacitic melts, and that further work is needed to understand the effect of 1497 temperature, redox, and non-ideal mixing between H 2 O-CO 2 on volatile solubility. Graph of Dixon, 1997 Raw North Arch Fixed Volatiles North Arch Figure S1. Π vs. SiO 2 graph shown in Fig. 2C of Dixon (1997). Red dots show the raw data presented in Table 1