VESIcal: An open-source thermodynamic model engine for mixed volatile (H₂O-CO₂) solubility in silicate melts

1Introduction¶

Understanding the solubility and degassing of volatiles in silicate melts is a crucial component of modeling volcanic systems. As dissolved components, volatiles (primarily H₂O and CO₂) affect magma viscosity, rheology, and crystal growth. In addition, due to the strong dependence of volatile solubility on pressure, measured volatile concentrations in preserved high-pressure melts (i.e., melt inclusions: liquid magma trapped within crystals at high pressure, then brough to the surface during an eruption) can be used to determine pre-eruptive magmatic storage pressures, and thus depths. Importantly, volatile exsolution-driven overpressure of a magmatic system is likely the trigger of many explosive volcanic eruptions Tait et al., 1989Blake, 1984Stock et al., 2016. Once triggered, further drops in magmatic pressure caused by ascent of magma within a volcanic conduit result in the continuous exsolution of volatiles from the melt. Volatile elements experience a large positive volume change when moving from a dissolved to exsolved free fluid state. This expansion fuels a dramatic increase in the magma’s buoyancy, which can often lead to a runaway effect in which the ascent and degassing of volatile-bearing magma eventually erupts at the surface in an explosive fashion. Working in concert with seismic and gas monitoring data, pre-eruptive magmatic volatile concentrations as well as solubility and degassing modelling can be used in forensic and sometimes in predictive scenarios, helping us to understand and potentially mitigate the effects of explosive eruptions.

All of these processes depend directly on the solubility – or the capacity of a magma to hold in solution – of volatile elements. Over the last several decades, a veritable explosion of new volatile solubility data has opened the door to a plethora of models describing the solubility of H₂O, CO₂, or mixed H₂O-CO₂ fluid in magmas covering a wide compositional, pressure, and temperature range. Volatile solubility is highly dependent upon the composition of the host magma, making already challenging experiments more onerous to perform to encapsulate the range of magmas seen in nature. The most fundamental models Stolper, 1982Dixon et al., 1995Moore et al., 1998 focus on a specific range of magma bulk compositions (e.g., basalt or rhyolite only). Later studies filled in compositional gaps, some with an increased focus on mixed-volatile (H₂O-CO₂) studies, increasing the natural applicability of our models to more systems Liu et al., 2005Iacono-Marziano et al., 2012Iacovino et al., 2013. To date, there have been only a few significant efforts to create a holistic thermodynamic model calibrated by a wide range of data in the literature. The most popular are MagmaSat (the mixed-volatile solubility model built into the software package MELTS v. 1.2.0; Ghiorso & Gualda (2015)) and the model of Papale et al. (2006). Both of these studies have made their source code available; the Papale et al. (2006) FORTRAN source code (titled Solwcad), web app, and a Linux program can be found at http://www.pi.ingv.it/progetti/eurovolc/, and very recently MagmaSat has been made accessible via the ENKI thermodynamic python framework (http://enki-portal.org/).

Despite this communal wealth of solubility models, quantitative calculations of volatile solubility, and by extension saturation pressures, equilibrium fluid compositions, and degassing paths, remains a time-consuming endeavor. Modeling tools that are available are typically unable to process more than one sample at a time, requiring manual entry of the concentrations of 8-10 major oxides, temperature, as well as CO₂ and H₂O concentrations to calculate saturation pressures, or X$_{H_{2}O}$ to calculate dissolved volatile contents. This is particularly problematic for melt inclusion studies, where saturation pressures are calculated for hundreds of inclusions, each with different entrapment temperatures, CO₂, H₂O, and major element concentrations. For example, the saturation pressures from 105 Gakkel ridge melt inclusions calculated in MagmaSat by Bennett et al. (2019) required the manual entry of 1,365 values! The potential for user error in this data entry stage should not be overlooked.

In many cases, newly published solubility models do not include an accompanying tool, requiring users to correctly combine and interpret the relevant equations e.g., Dixon et al., 1995Dixon, 1997Liu et al., 2005Shishkina et al., 2014. This is problematic from a perspective of reproducibility of the multitude of studies utilizing these models, especially given that some of the equations in the original manuscripts contain typos or formatting errors. For some models, an Excel spreadsheet was provided, or available at request from the authors. For example, Newman & Lowenstern (2002) included a simplified version of the Dixon (1997) model as part of “VolatileCalc”, which was written in Visual Basic for Excel. Due to its simplicity, allowing users to calculate saturation pressures, degassing paths, isobars and isopleths with a few button clicks and pop-up boxes, this tool has proved extremely popular (with 836 citations at the time of writing). However, to calculate saturation pressures using VolatileCalc, the user must individually enter the SiO₂, H₂O, CO₂ content and temperature of every single sample into pop-up boxes. Similarly, the Excel spreadsheet for the Moore et al. (1998) model calculates dissolved H₂O contents based on the concentration of 9 oxides, temperature, and the fraction of X$_{H_{2}O}$ in the vapor, which must be pasted in for every sample. Finally, Allison et al. (2019) provide an Excel spreadsheet that allows users to calculate fugacities, partial pressures, isobars, isopleths and saturation pressures. Again, parameters for each sample must be entered individually, with no way to calculate large numbers of samples automatically.

Some of these published models and tools are at risk of being lost to time, since spreadsheet tools (particularly earlier studies published before journal-provided hosting of data and electronic supplements was commonplace) must be obtained by request to the author. Even if the files are readily available, programs used to open and operate them may not support depreciated file formats. More recently, authors have provided web-hosted interfaces to calculating saturation pressures and dissolved volatile contents (e.g., Iacono-Marziano et al. (2012); http://calcul-isto.cnrs-orleans.fr/, and Ghiorso & Gualda (2015); http://melts.ofm-research.org/CORBA_CTserver/GG-H2O-CO2.html). Ghiorso & Gualda (2015) also provide a Mac application. While more accessible in the present time, this does not negate the issue of the longevity of these models. The link provided in the Iacono-Marziano et al. (2012) manuscript returns an error “this site cannot be reached”, although email contact with the author directed us towards the newer link given above. Similarly, the link to the H₂O-CO₂ equation of state web calculator that Duan & Zhang (2006) provided in their manuscript returns a 404 error.

While we certainly advocate for the continued refinement of solubility models, including the completion of new experiments in poorly studied yet critical compositional spaces such as andesites Wieser et al., 2021, a perhaps more crucial step at this juncture is in the development of a tool that can apply modern computational solutions to making our current knowledge base of volatile solubility in magmas accessible and enduring.

Here we present VESIcal (Volatile Equilibria and Saturation Identification calculator): a python-based thermodynamic volatile solubility model engine that incorporates seven popular volatile solubility models under one proverbial roof. The models included in VESIcal are (also see Table 1):

1. MagmaSat: VESIcal’s default model. The mixed-volatile solubility model within MELTS v. 1.2.0 Ghiorso & Gualda, 2015
2. Dixon: The simplification of the Dixon (1997) model as implemented in VolatileCalc Newman & Lowenstern, 2002
   • DixonWater and DixonCarbon available as pure-fluid models
3. MooreWater: (Moore et al. (1998); water only, but H₂O fluid concentration can be specified)
4. Liu: Liu et al., 2005
   • LiuCarbon and LiuWater available as pure-fluid models
5. IaconoMarziano: Iacono-Marziano et al., 2012
   • IaconoMarzianoWater and IaconoMarzianoCarbon available as pure-fluid models
6. ShishkinaIdealMixing: Shishkina et al., 2014 using pure-H₂O and pure-CO₂ models and assuming ideal mixing. In general, the pure-fluid versions of this model should be used
   • ShishkinaWater and ShishkinaCarbon available as pure-fluid models
7. AllisonCarbon: Allison et al., 2019, carbon only
   • AllisonCarbon_vesuvius (default; phonotephrite from Vesuvius, Italy)
   • AllisonCarbon_sunset (alkali basalt from Sunset Crater, AZ, USA)
   • AllisonCarbon_sfvf (basaltic andesite from San Francisco Volcanic Field, AZ, USA)
   • AllisonCarbon_erebus (phonotephrite from Erebus, Antarctica)
   • AllisonCarbon_etna (trachybasalt from Etna, Italy)
   • AllisonCarbon_stromboli (alkali basalt from Stromboli, Italy)

As any individual model is only valid within its calibrated range (see below), and each model is parameterized and expressed differently (e.g., empirical vs. thermodynamic models), it is impractical to simply combine them into one large model. Instead, VESIcal is a single tool that can access and utilize all built-in models, giving VESIcal an extensive pressure-temperature-composition calibration range Figure 1. VESIcal is capable of performing a wide array of calculations on large datasets automatically, with built-in functionality for extracting data from an Excel or CSV file. In addition, the code is written such that it is flexible (sample, calculation type, and model type can be chosen discreetly) and extensible (VESIcal code can be imported for use in python scripts, and the code is formatted such that new volatile models can be added).

Importantly, VESIcal has been designed for practicality and ease of use. It is designed to be used by anyone, from someone who is completely unfamiliar with coding to an adept programmer. The non-coder user can interact with VESIcal through a webapp (https://vesical.anvil.app) or directly within this manuscript, which utilizes the user-friendly Jupyter Notebook format, allowing them to upload a file with data, execute the various example calculations provided below, and save the results to an Excel or CSV file to work with outside of VESIcal. This notebook also incorporates built-in plotting options for easy visualization of user data and calculated results. More experienced programmers may wish to use the more advanced functionality provided by VESIcal, including the ability to hybridize models (e.g., use one model for H₂O and another for CO₂) or write their own routines and code calling VESIcal methods. VESIcal is an open source tool and as such is far less prone to the preservation issues discussed above. Because the VESIcal code is hosted on GitHub, every change to the code is tracked publicly Perkel, 2016. VESIcal’s current release (version 1.0.1) is also archived on Zenodo, which provides a static citable DOI (10.5281/zenodo.5095382) for the current version of the code, along with a snapshot of the GitHub repository at the time of release.

A detailed history of volatile solubility modeling and the implications of VESIcal are explored in detail in the companion manuscript to this work, Wieser et al. (2021).

2Research Methodology¶

Navigating the array of models implemented in VESIcal can be challenging. How can a user determine which model best suits their needs? MagmaSat (the default model in VESIcal) is the most widely calibrated in P-T-X space, and so we recommend it for the majority of cases. Where a user wishes to use the other implemented models, we provide some tools to help choose the most appropriate model (see Supplement). These tools are described in more detail in Section 3.4 on comparing user data to model calibrations.

A list of model names recognized by VESIcal can be retrieved by executing the command v.get_model_names(), assuming VESIcal has been imported as v as is demonstrated in worked examples below. Note that the model names as listed in the previous section are given in terms of how to call them within VESIcal (e.g., model='MooreWater'). Allison et al. (2019) provides unique model equations for each of the six alkali-rich mafic magmas investigated in their study. The default model in VESIcal is that calibrated for Vesuvius magmas, whose calibration has the widest pressure range of the study Table 1. Setting a model name of 'AllisonCarbon' within VESIcal will thus result in calculations using the AllisonCarbon_vesuvius model equations.

All of the calculations implemented in VESIcal can be performed using any of the models included. The code is structured by calculation rather than by model, which provides an intuitive way for users to interact with the code and compare outputs from multiple models. Each calculation class is instantiated with the model name and any applicable data as arguments. It then performs five key functions: 1) creates the requested model object and performs any necessary pre-processing (e.g., ensuring relevant data are present; normalizing data); 2) takes user input and performs the mathematical calculation; 3) does any necessary processing of the output (e.g., normalizing totals); 4) checks that the model is being used within its calibrated range; and 5) stores calculated outputs in an intuitive and manipulatable format (e.g., a python dictionary, a figure, or a pandas DataFrame). Results of calculations can be saved to one or more Excel or CSV files. To demonstrate that VESIcal returns results which are comparable with pre-existing tools, we have performed a number of tests, which are described in the Supplementary Information (Text S2). For single-sample calculations, the calculation object has the following attributes that can be called by the user: model_name, sample (both provided by the user), model (an instance of the Model class used to run the calculations of interest), result (the result of the calculations), and calib_check (the results of the calibration check).

2.1Model Calibrations and Benchmarking¶

The pressure, temperature, and composition calibration ranges of the seven models implemented in VESIcal are shown in Table 1 and Figure 1. VESIcal abides by statements of caution made by the authors of these models regarding their extrapolation by informing the user if a calculation is being performed outside of a model’s calibrated range. In this case, the code returns a warning message, which is as specific as possible, along with the requested output. We also provide several Jupyter notebooks in the supplementary material (Supplementary Text S3-S4 and Supplementary files S1-S7), allowing users to plot their data amongst the calibrations of the different models to assess their suitability for less objective measures. Detailed descriptions of the seven solubility models implemented in VESIcal, including information about their calibration range in terms of melt composition, pressure, and temperature, are given in Wieser et al. (2021)

Testing was undertaken to ensure that VESIcal faithfully reproduces the results of all incorporated models. When possible, all models were benchmarked by testing VESIcal outputs against those of a relevant published calculator (e.g., web apps or Excel macros). The models of Shishkina et al. (2014) and Liu et al. (2005) were published with no such tool and so testing instead compares VESIcal outputs to experimental conditions or analyses and, where possible, plots VESIcal results against published figures. All models underwent multiple tests, the results of which are shown in the supplement (Supplementary Text S3-S4 and Supplemental Jupyter Notebooks S1-S7). For all models, VESIcal reproduced the results from previous tools (e.g., web apps, Excel spreadsheets) to within ±1% relative and often on the order of ±0.1% relative.

MagmaSat, VESIcal’s default model, underwent three tests, the results of which are shown in Figure 2: 1. Comparison of saturation pressures from MORB melt inclusions in VESIcal to those published by Bennett et al. (2019), who used the MagmaSat Mac App (R²=0.99998; Figure 2a); 2. Comparison of fluid composition (X$_{H_{2}O}$) calculated with VESIcal and the web app (R²=0.999, identical considering the web app returns 2dp; Figure 2b); 3. Comparison of isobars for the Early Bishop Tuff calculated with VESIcal (star symbols) and isobars published in Figure 14 of Ghiorso & Gualda (2015) (Figure 2c). VESIcal outputs using the model of Dixon (1997) were tested against outputs from the VolatileCalc Excel spreadsheet Newman & Lowenstern, 2002 and a widely used Excel macro e.g., Tucker et al., 2019.

2.2Format of the python library¶

In this section, the basic organization and use cases of VESIcal are discussed. VESIcal relies heavily on python pandas, a python package designed for working with tabulated data. Knowledge of pandas is not required to use VESIcal, and we refer the user to the pandas documentation for an overview of the package (https://pandas.pydata.org/pandas-docs/stable/userguide/index.html).

Specific details on how to perform model calculations are discussed in Section 3 and include worked examples. The VESIcal library is written so that users can interact first and foremost with the calculation they want to perform. Five standard calculations can be performed with any model in the library:

1. calculate_dissolved_volatiles()
2. calculate_equilibrium_fluid_composition()
3. calculate_saturation_pressure()
4. calculate_isobars_and_isopleths() (plus functionality for plotting; only for mixed volatiles models)
5. calculate_degassing_path() (plus functionality for plotting; only for mixed volatiles models).

Figure 3 illustrates the basic organization of the code. First, the user determines which calculation they wish to perform by accessing one of the five core calculation classes (listed above). In this step, the user specifies any input parameters needed for the calculation (e.g., sample composition in wt% oxides, pressure in bars, temperature in °C, and fluid composition “X_fluid” in terms of XH₂O^fluid) as well as the model they wish to use. The default model is MagmaSat, but the user may specify any model in the library. As an example, the code to calculate the saturation pressure of some sample using the MagmaSat model would be written as:

saturation_pressure_calculation = calculate_saturation_pressure(sample=mysample, temperature=850.0)

where mysample is a variable (VESIcal Sample object) containing the composition of the sample, and the temperature is given in °C. Here, this line of code creates a Calculate object, which is something that can be given a variable name and stored so that the user can call upon this object for viewing or manipulation later. In this example, we name the object “saturation_pressure_calculation”, but this can be any variable name desired by the user. The Calculate object stores important information about the calculation, including the result. The result of the calculation or calibration check can be accessed as:

saturation_pressure_calculation.result
saturation_pressure_calculation.calib_check

In python, the object creation and attribute access can be combined into a single line, with the understanding that the Calculate object will not be accessible to the user. This usage is used in the remaining examples throughout the manuscript and would be written as:

saturation_pressure = calculate_saturation_pressure(sample=mysample, temperature=850.0).result

If a different model is desired, for example Dixon (1997), it can be passed as:

calculate_saturation_pressure(sample=mysample, temperature=850.0, model='Dixon').result

The core calculation classes each perform two functions: 1) a check is performed to ensure that the user input is within the model’s recommended calibration range; 2) the calculate() method sends the user input to the appropriate model.

Users can process individual samples (single-sample calculations) or entire datasets (batch calculations; Figure 4). If processing more than one sample, the “simplest” way to interact with VESIcal is via batch calculations. Here, the user provides input data in the form of a Microsoft Excel spreadsheet (.xlsx file) or CSV file and instructs the model to perform whatever calculation is desired. The model is run on all samples and returns data formatted like a spreadsheet (using the python pandas package), which contains the user’s original input data plus whatever model outputs were calculated. The user can continue to work with returned data by saving the result to a variable (as is shown in all examples in this manuscript). Data can then be exported to an Excel or CSV file with a simple command (see Section 3.12).

The syntax for processing a single sample is very similar to that for batch calculations but provides the user direct access to more advanced features that cannot be accessed via batch calculations (e.g., specifying fugacity or activity model, hybridizing models; see Section 3.11). This also gives the user more flexibility in integrating any VESIcal model function into some other python code.

2.3Running the code¶

VESIcal can be used in a number of ways: via this Jupyter notebook, via the VESIcal web app, or by directly importing VESIcal into any python script.

VESIcal was born from functionality provided by ENKI and so all the files necessary to use VESIcal are hosted on the ENKI server (http://enki-portal.org/). A unique personal coding environment can be initiated by logging into the ENKI production server using a GitLab username and password (which is free to obtain; see directions on the ENKI website for specifics). The simplest way to use VESIcal while retaining all of its functionality is within this very manuscript, in the form of a Jupyter notebook. Because this manuscript and VESIcal python library files are hosted on the ENKI server, code can be manipulated and executed in the code cells below. Making changes won’t affect the public version of this manuscript. Likewise, any user can write their own python code using VESIcal by creating a Jupyter notebook on the ENKI server and importing VESIcal as is demonstrated in the code below.

Computation time on the ENKI server is limited by the server itself. VESIcal may run faster if installed locally. Advanced instructions on installing VESIcal on your own computer are provided in the Supplement (Supplementary Text S1). Note that VESIcal requires installation of the ENKI thermoengine library to function properly. Thermoengine is written in python but is based on the original MELTS code Ghiorso & Sack, 1995Ghiorso & Gualda, 2015, which contains MacOS-specific header files. The result is that thermoengine is most easily installed on MacOS but can be installed on Windows and Linux operating systems via Docker (see thermoengine documentation for installation instructions; https://gitlab.com/ENKI-portal/ThermoEngine).

The most limited but simplest method to interacting with VESIcal is through the web app (https://vesical.anvil.app). The web app can currently perform three of the five core calculations in batch process mode (via upload of an Excel or CSV file). Some, but not all, optional parameters can be set.

To run the code in this notebook, nothing needs to be installed. Simply execute the code cells below, changing parameters as desired. Custom data may be processed by uploading an Excel or CSV file into the same folder containing this notebook and then changing the filename in Section 3.3.

2.4Documentation¶

This manuscript serves as an introduction to the VESIcal library aimed at python users of all levels. However, the code itself is documented with explanations of each method, its input parameters, and its returned values. This documentation can be accessed at our readthedocs website (https://vesical.readthedocs.io/). The documentation for any function can be viewed in a Jupyter Notebook by typing the function followed by a question mark and executing the cell (e.g., “v.calculate_saturation_pressure?”).

Video tutorials are also available on the VESIcal YouTube (https://www.youtube.com/channel/UCpvCCs5KMXzOxXWm0seF8Qw). The first tutorial covers the basics of VESIcal, and subsequent videos cover the calculation of saturation pressures, dissolved volatile contents, and degassing paths. More videos for specific features and uses are planned.

2.5Generic methods for calculating mixed-fluid properties¶

VESIcal provides a set of methods for calculating the properties of mixed CO₂-H₂O fluids, which can be used with any combination of H₂O and CO₂ solubility model. The use of generic methods allows additional models to be added to VESIcal by defining only the (simpler) expressions describing pure fluid solubility. Non-ideality of mixing in the fluid or magma phases can be incorporated by specifying activity and fugacity models. A complete description of these methods, including all relevant equations, can be found in the Supplement (Supplementary Text S2).

3Workable example uses¶

In this section we detail how to use the various functions available in VESIcal through worked examples. The python code presented below may be copied and pasted into a script or can be edited and executed directly within the Jupyter notebook version of this manuscript. For all examples, code in Sections 3.2 and 3.4 must be executed to initialize the model and import data from the provided companion Excel file. The following sections then may be executed on their own and do not need to be executed in order.

In each example below, a generic “method structure” is given along with definitions of unique, required, and optional user inputs. The method structure is simply for illustrative purposes and gives default values for every argument (input). In some cases, executing the method structure as shown will not produce a sensible result. For example, the default values for the plot() function Section 3.10 contain no data, and so no plot would be produced. Users should replace the default values shown with values corresponding to the samples or conditions of interest.

All examples will use the following sample data by default (but this can be changed by the user):

• Dataset from example_data.xlsx loaded in Section 3.3.1 (variable name myfile)
• Single composition defined in Section 3.3.2 (variable name mysample)
• Sample 10* extracted from example_data.xlsx dataset in Section 3.3.3 (variable name sample_10)

Calculations performed on single samples or on a dataset imported from an Excel or CSV file containing many samples are executed in two distinct ways. Note that single sample calculations require that the argument sample be defined. To return the numerical result of the calculation, the .result method must be called, as shown below. Batch calculations are performed on the dataset itself, after that dataset is imported into VESIcal. Thus, the sample argument does not need to be defined discretely, since sample compositional information is stored within the dataset object. The two basic formats for performing calculations are:

Single sample calculations

myvariable = v.name_of_the_core_calculation(sample=mysample, argument1=value1, argument2=value2).result

Batch calculations

myvariable = myfile.name_of_the_core_calculation(argument1=value1, argument2=value2)

where VESIcal has been imported as v, myvariable is some arbitrary variable name to which the user wishes to save the calculated output, name_of_the_core_calculation is one of the five core calculations, mysample is a variable containing compositional information in wt% oxides, myfile is a variable containing an BatchFile object created by importing an Excel or CSV file, and argument1, argument2, value1, and argument2 are two required or optional arguments and their user-assigned values, respectively.

Workable examples detailed here are:

1. Loading, viewing, and preparing user data
   • Loading a Batch file
   • Defining a single sample composition
   • Extracting a single sample from a Batch file
   • Normalizing and transforming data
2. Calculating dissolved volatile concentrations
3. Calculating equilibrium fluid compositions
4. Calculating saturation pressures
5. Calculating and plotting isobars and isopleths
6. Calculating and plotting degassing paths
7. Plotting multiple calculations
8. Comparing results from multiple models
9. Code hybridization (Advanced)
10. Exporting data

3.1Calculation class arguments and their definitions¶

Each section below details what arguments are required or optional inputs and gives examples of how to perform the calculations. Table 2 lists all arguments, both required and optional, used in the five core calculations. Many of the function arguments have identical form and use across all calculations, and so we list these here. Any special cases are noted in the section describing that calculation.

The most commonly used arguments are:

sample

Single sample calculations only The composition of a sample. A VESIcal Sample object is created to hold compositional information about sample. A Sample object can be created from a dictionary or pandas Series containing values, with compositions of oxides in wt%, oxides in mol fraction, or cations in mol fraction. This argument is not needed for batch calculations since they are performed on BatchFile objects, which already contain sample information. See examples for details.

temperature, pressure, and X_fluid

the temperature in °C, the pressure in bars, and the mole fraction of H₂O in the H₂O-CO₂ fluid, XH₂O^fluid. In all cases, X_fluid is optional, with a default value of 1 (pure H₂O fluid). Note that that X_fluid argument is only used for calculation of dissolved volatile concentrations.

For single sample calculations

Temperature, pressure, and X fluid should be specified as a numerical value.

For batch calculations

Temperature, pressure, and X_fluid can either be specified as a numerical value or as strings referring to the names of columns within the file containing temperature, pressure, or X_fluid values for each sample. If a numerical value is passed for either temperature, pressure, or X_fluid, that will be the value used for one or all samples. If, alternatively, the user wishes to use temperature, pressure, and/or X_fluid information in their BatchFile object, the title of the column containing temperature, pressure, or X_fluid data should be passed in quotes (as a string) to temperature, pressure, and/or X_fluid, respectively. Note for batch calculations that if temperature, pressure, or XH₂O^fluid information exists in the BatchFile but a single numerical value is defined for one or both of these variables, both the original information plus the values used for the calculations will be returned.

verbose

Only for single sample calculations Always an optional argument with a default value of False. If set to True, additional values of interest, which were calculated during the main calculation, are returned in addition to the results of the calculation.

print_status

Only for batch calculations Always an optional argument, which sometimes defaults to True and other times defaults to False (see specific calculation section for details). If set to True, the progress of the calculation will be printed to the terminal. The user may desire to see the status of the calculation, as some calculations using MagmaSat can be somewhat slow, particularly for large datasets.

model

Always an optional argument referring to the name of the desired solubility model to use. The default is always “MagmaSat”.

3.2Initialize packages¶

For any code using the VESIcal library, the library must be imported for use. Here we import VESIcal as v. Any time we wish to initialize a VESIcal object, that class name must be preceded by ‘v.’ (e.g., v.calculate_saturation_pressure). Specific examples of this usage follow. Here we also import some other python libraries that we will be using in the worked examples below.

3.3Loading, viewing, and preparing user data¶

All of the following examples will use data loaded in the code cells in this section. Both batch processing of data loaded from a file and single-sample processing are shown. An example file called example_data.xlsx is included with this manuscript. You can load in your own data by first ensuring that your file is in the same folder as this notebook and then by replacing the filename in the code cell below with the name of your file. The code cell below must be executed for the examples in the rest of this section to function properly.

3.3.1Batch processing¶

Batch calculations are always facilitated via the BatchFile() class, which the user uses to specify the filename corresponding to sample data. Loading in data is as simple as calling BatchFile(filename). Optionally, units can be used to specify whether the data are in wt% oxides, mol fraction oxides, or mol fraction cations. Calculations will always be performed and returned with melt composition in the default units (wt% oxides unless changed by the user) and fluid composition in mol fraction.

Structure of the input file: A file containing compositions (and optional pressure, temperature, or XH₂O^fluid information) on one or multiple samples can be loaded into VESIcal. The loaded file must be a Microsoft Excel file with the extension .xls or .xlsx or CSV file with the extension .csv. The file must be laid out in the same manner as the example file example_data.xlsx. The basic structure is also shown in Table 3.

Any extraneous columns that are not labeled as oxides or input parameters will be ignored during calculations. The first column titled ‘Label’ contains sample names. Note that the default assumption on the part of VESIcal is that this column will be titled ‘Label’. If no ‘Label’ column is found, the first non-oxide column name will be set as the index column, meaning this is how samples can be accessed by name (see Section 3.3.3). An index column can be specified by the user using the argument label (see documentation below). The following columns must contain compositional information as oxides. The only allowable oxides are: SiO₂, TiO₂, Al₂O₃, Fe₂O₃, FeO, Cr₂O₃, MnO, MgO, CaO, NiO, CoO, Na₂O, K₂O, P₂O₅, H₂O, and CO₂. Currently, VESIcal can only read these oxide names exactly as written (e.g., with no leading or trailing spaces and with correct capitalization), but functionality to interpret variations in how these oxides are entered is planned (e.g., such that "sio2. " would be understood as “SiO2”). All of these oxides need not be included; if for example your samples contain no NiO concentration information, you can omit the NiO column. Omitted oxide data will be set to 0 wt% concentration. If other oxide columns not listed here are included in your file, they will be ignored during calculations. Notably, the order of the columns does not matter, as they are indexed by name rather than by position. Compositions can be entered either in wt% (the default), mol%, or mole fraction. If mol% or mole fraction data are loaded, this must be specified when importing the tile.

Because VESIcal assumes a particular formatting of column names, we highly recommend that users examine their data after loading into VESIcal and before performing calculations. The user data, as it will be used by VESIcal, can be viewed at any time with myfile.get_data() (see generation of Table 3 below).

Pressure, temperature, or XH₂O^fluid data may optionally be included, if they are known. Column names for these data do not matter, as they can be specified by the user as will be shown in following examples.

The standard units used by VESIcal are always pressure in bars, temperature in °C, melt composition as oxides in wt%, and fluid composition as mol fraction (typically specified as X_fluid, the mol fraction of H₂O in an H₂O-CO₂ fluid, ranging from 0-1). Sample compositions may be translated between wt%, mol fraction, and mol cations if necessary.

Class structure

BatchFile(filename, sheet_name=0, file_type='excel', units='wtpt_oxides', label='Label', default_normalization='none', default_units='wtpt_oxides', dataframe=None)

Required inputs

filename: A file name must be passed in quotes. This file must be in the same folder as the notebook or script that is calling it. This imports the data from the file name given and saves it to a variable of your choosing.

Optional inputs By default, the BatchFile class assumes that loaded data is in units of wt%; alternatively, data in mole fraction oxides or cations may be loaded.

sheet_name: If importing data from an Excel file, this argument is used to specify which sheet to import. Only one sheet can be imported to a single BatchFile object. The default is ‘0’, which imports the first sheet in the file, regardless of its name.

file_type: Specifies whether the file being imported is an Excel or CSV file. This argument is never strictly necessary, as BatchFile() will automatically detect whether an imported file is Excel or CSV if the file extension is one of .xls or .xslx (Excel) or .csv (CSV).

units: The units in which data are input. The default value is 'wtpt_oxides' for data as wt% oxides. The user can pass 'mol_oxides' for data in mol fraction oxides or 'mol_cations' for data in mol fraction cations.

default_normalization: The type of normalization to apply to the data by default. One of: None, 'standard', 'fixedvolatiles', or 'additionalvolatiles'. These normalization types are described in the section on normalization below.

default_units: The type of composition to return by default, one of: 'wtpt_oxides' (wt% oxides, default), 'mol_oxides' (mol fraction oxides), or 'mol_cations' (mol fraction cations).

label: This is optional but can be specified if the column title referring to sample names is anything other than “Label”. The default value is “Label”. The default value is “Label”. If no “Label” column is present and the label argument is not specified, the first column whose first row is not one of VESIcal’s recognized oxides will be set as the index column. The index column will be used to select samples by name.

dataframe: This argument is used for transforming a pandas DataFrame object into a VESIcal BatchFile object. For convenience, this functionality is also defined as a separate function BatchFile_from_DataFrame(dataframe, units='wtpt_oxides', label='Label').

Outputs

A special type of python object defined in the VESIcal code known as an BatchFile object.

Once the BatchFile object is created and assigned to a variable, the user can then access the data loaded from their file as variable.get_data(). In this example, the variable corresponding to the BatchFile object is named myfile and so the data in that file can be accessed with myfile.get_data(). Below, myfile.get_data() is saved to a variable we name data. The variable data is a pandas DataFrame object, which makes displaying the data itself quite simple and aesthetically pleasing, since pandas DataFrames mimic spreadsheets.

Usage of get_data() allows the user to retrieve the data as originally entered or in any units and with any normalization supported by VESIcal.

Method structure

get_data(self, normalization=None, units=None, asBatchFile=False)

Optional inputs

normalization or units may be passed, with options as defined in the description of BatchFile above.

asBatchFile Default is False. If True, will return a VESIcal BatchFile object.

Outputs

A pandas dataframe or BatchFile object with all user data.

For the rest of this manuscript, data will be pulled from the example_data.xlsx file (Supplemental Dataset S1), which contains compositional information for basalts Tucker et al., 2019Roggensack, 2001, andesites Moore et al., 1998, rhyolites Mercer et al., 2015Myers et al., 2019, and alkaline melts (phototephrite, basaltic-trachyandesite, and basanite from Iacovino et al. (2016)). Several additional example datasets from the literature are available in the Supplement (Supplementary Datasets S2-S5; Table 4). These include experimentally produced alkaline magmas from Iacovino et al. (2016, `alkaline.xlsx`), basaltic melt inclusions from Kilauea Tucker et al., 2019 and Gakkel Ridge Bennett et al., 2019, `basalts.xlsx`, basaltic melt inclusions from Cerro Negro volcano, Nicaragua Roggensack, 2001, `cerro_negro.xlsx`, and rhyolite melt inclusions from the Taupo Volcanic Center, New Zealand Myers et al., 2019 and a topaz rhyolite from the Rio Grande Rift Mercer et al., 2015, `rhyolites.xlsx`. Where available, the calibration datasets for VESIcal models are also provided (Supplementary Datasets S6-S7).

3.3.2Defining a single sample¶

More advanced functionality of VESIcal is facilitated directly through the five core calculation classes. Each calculation requires its own unique inputs, but all calculations require that a sample composition be passed. We can pass in a sample either as a python dictionary or pandas Series. Below, we define a sample and name it mysample. Oxides are given in wt%. Only the oxides shown here can be used, but not all oxides are required. Any extra oxides (or other information not in the oxide list) the user defines will be ignored during calculations.

Much like is done to create a BatchFile object, we can create a VESIcal Sample object to represent our sample composition.

Class structure

Sample(composition, units='wtpt_oxides', default_normalization='none', default_units='wtpt_oxides')

Required inputs

composition: The composition of the sample in the format specified by the units parameter. The default is oxides in wt%.

Optional inputs

units, default_normalization, and default_units have the same meaning here as in the BatchFile class desribed above.

Outputs

A special type of python object defined in the VESIcal code known as a Sample object.

To manually input a bulk composition, fill in the oxides in wt% below:

To see the composition of mysample, use the get_composition(species=None, normalization=None, units=None, exclude_volatiles=False, asSampleClass=False) method. By default, the composition is returned exactly as input above. species can be set as an element or oxide (e.g., “Si” or “SiO₂”)to return the float value for only that species. The composition can automatically be normalized using any of the standard normalization functions listed above and can be returned in any of the units discussed above. As with the BatchFile.get_data() function, a sample composition can be returned as a dictionary (default) or as a VESIcal Sample object (if asSampleClass is set to True).

The oxides considered by VESIcal are:

3.3.3Extracting a single sample from a batch file¶

Defined within the BatchFile() class, the method get_sample_composition() allows for the extraction of a melt composition from a loaded Excel or CSV file.

Method structure

myfile.get_sample_composition(samplename, species=None, normalization=None, units=None, asSampleClass=False)

Required inputs

samplename: The name of the sample, as a string, as defined in the ‘Label’ column of the input file.

Optional inputs

species: This is used if only the concentration of a single species (either oxide or element) is desired.

normalization: This is optional and determines the style of normalization performed on a sample. The default value is None, which returns the value-for-value un-normalized composition. Other normalization options are described in the BatchFile class description above.

units: The default is wt% oxides. Other options are described in the BatchFile class description above.

asSampleClass: Can be True or False (default). If set to False, this will return a dictionary with compositional values. If set to True, this will return a Sample object with compositional data stored within.

Outputs

The bulk composition stored in a dictionary or Sample object.

3.3.4Normalizing and transforming data¶

Before performing model calculations on your data, it may be desired to normalize the input composition to a total of 100 wt%. For a user to decide whether normalization is prudent, is important to understand the influence any normalization, or lack thereof, to a composition will have on modeling results. Electron microprobe analyses of major elements in silicate glasses combined with volatile element analyses by SIMS and FTIR often sum to less than 100 wt%. This deficiency is normally attributed to subsurface charging, matrix corrections, and unknown redox states of Fe and S during analyses by electron microprobe see Hughes et al., 2019. As an example, when normalized, a volatile-free basalt with a measured SiO₂ content of 46 wt% and an analytical total of 97 wt% actually contains 47.4 wt% SiO₂ (46/0.97; a 3% relative change in silica content). Many studies report major element data normalized to 100% with volatiles listed separately. The result is that, value for value, literature datasets can have totals several wt% less than 100 (if raw data are reported) or several wt% higher than 100 (if major elements are normalized anhydrous).

To deal with this variation, VESIcal provides users with four options for normalization. Normalization types are: - None (no normalization) - ‘standard’: Normalizes an input composition to 100%. - ‘fixedvolatiles’: Normalizes major element oxides to 100 wt%, including volatiles. The volatile wt% will remain fixed, whilst the other major element oxides are reduced proportionally so that the total is 100 wt%. - ‘additionalvolatiles’: Normalizes major element oxide wt% to 100%, assuming it is volatile-free. If H₂O or CO₂ are passed to the function, their un-normalized values will be retained in addition to the normalized non-volatile oxides, summing to >100%.

Normalization can be performed on a Sample object or on all samples within a BatchFile object using the get_composition() or get_data() methods (e.g., myfile.get_composition(normalization='standard') or mysample.get_composition(normalization='additionalvolatiles'). Note that, since a BatchFile object may have other data in addition to sample compositions (e.g., information on pressure, temperature, other user notes), BatchFile.get_composition() returns only compositional data, where as BatchFile.get_data() returns all data stored in the BatchFile object. The normalization argument can be passed to either. In the example below, we obtain the standard normalization of mysample and myfile and save these to new Sample and BatchFile objects called mysample_normalized and myfile_normalized. Note that asSampleClass or asBatchFile must be set to True in order to return a Sample or BatchFile object. Without this argument, a dictionary or pandas DataFrame will be returned and new Sample or BatchFile objects will need to be constructed from those in order to perform calculations on the normalized datasets.

The Liu and all six AllisonCarbon models are not sensitive to normalization because they contain no compositional terms. Similarly, the expressions for Shishkina and MooreWater contain compositional terms expressed solely in terms of anhydrous cation fractions; the additionalvolatiles and fixedvolatiles normalization routines do not affect the relative abundances of major elements (and therefor anhydrous cation fractions). Thus, Shishkina and MooreWater are only affected by the standard normalization routine. In contrast, the Dixon model is highly sensitive to the choice of normalization because its compositional term for both H₂O and CO₂ is expressed solely in terms of the absolute melt SiO₂ content.

The expressions of Iacono-Marziano are parameterized in terms of hydrous cation fractions and NBO/O, and so this model is sensitive to additionalvolatiles or fixedvolatiles normalization routines, which will change the relative proportions of volatiles to major elements. Even so, the effect of normalization on volatile solubility calculations is relatively small and of similar magnitude to the discrepancy between the hydrous total and 100 for the hydrous model. Thus, the choice of normalization is only important when data has hydrous totals that differ significantly from 100%. The Iacono-Marziano web app normalizes input data a la VESIcal’s additionalvolatiles normalization routine. For consistency with the web app, VESIcal automatically uses the additionalvolatiles normaliztion during calculations with this model.

The implementation of MagmaSat in VESIcal is sensitive to the relative proportion of major and volatile element components rather than the absolute concentrations entered (as with the whole MELTS family of models). Thus, calculations using raw, fixed- and additionalvolatile routines yield different results. If the hydrous total of an input composition is less than 100%, the fixedvolatile routine effectively reduces the relative proportion of volatiles to major elements, so calculated saturation pressures go down. Conversely, if inputs have high hydrous totals, the fixedvolatile routine increases the relative proportion of volatiles in the system, so the saturation pressure goes up. As with Iacono-Marziano, the percent discrepancy between calculations for different normalization routines is similar to the difference between the total and 100%. For saturation pressure calculations, the MagmaSat app automatically normalizes input data in a manner identical to VESIcal’s fixedvolatiles routine, and so this normalization is forced on all samples for such calculations with MagmaSat in VESIcal. Further discussion on the effect of normalization in MagmaSat is provided in Supplementary Text S5 (and Supplementary Figs S22-S26).

For example, consider a basalt with a measured SiO₂ content of 47.4 wt%, 1000 ppm dissolved CO₂, and an anhydrous (volatile-free) total of 96.77 wt%:

We can apply each normalization routine to this sample and examine how this will affect the saturation pressure predicted by each model:

Because the compositional effect on H₂O solubility is smaller, so are the changes in calculated saturation pressures for a pure-H₂O system, but they can still be significant for H₂O-rich liquids (where high H₂O contents can change totals and therefor SiO₂ contents more dramatically).

3.4Comparing User Data to Model Calibrations: Which Model Should I Use?¶

MagmaSat is the most thermodynamically robust model implemented in VESIcal, and thus it is the most generally appropriate model to use (n.b. that it is also the most computationally expensive). However, one of the strengths of VESIcal is its ability to utilize up to seven different solubility models. Each of these models is based on its own calibration dataset, meaning the pressure-temperature-composition space over which models are calibrated is quite variable from model to model. The individual model calibrations are discussed in detail in this manuscript’s companion paper VESIcal Part II, Wieser et al., 2021.

For the remainder of this section, all example calculations are carried out with MagmaSat, the default model of VESIcal. To use any other VESIcal model, simply add model= and the name of the desired model in quotes to any calculation (e.g., v.calculate_dissolved_volatiles(temperature=900, pressure=1000, model="Dixon"). The model names recognized by VESIcal are: MagmaSat, ShishkinaIdealMixing, Dixon, IaconoMarziano, Liu, AllisonCarbon, MooreWater. For more advanced use cases such as hybridizing models (see Section 3.11), pure-H₂O and pure-CO₂ models from within a mixed-fluid model can be used by adding ‘Water’ or ‘Carbon’ to the model name (e.g., DixonCarbon; note that MagmaSat does not have this functionality).

Determination of the appropriate model to use with any sample is crucial to the correct application of these models, and so we stress the importance of understanding how a model’s calibration space relates to the sample at hand. VESIcal includes some built-in functionality for comparing melt compositions from user loaded data to those in the datasets upon which each of the VESIcal models is calibrated using the method calib_plot. This can be visualized as a total alkalis vs silica (TAS) diagram (with fields and labels via the python tasplot library by J. Stevenson; https://bitbucket.org/jsteven5/tasplot/src/master/; Figure 5 a) or as any x-y plot in which x and y are oxides (Figure 5 b).

Method structure

calib_plot(user_data=None, model='all', plot_type='TAS', zoom=None, save_fig=False)

Optional inputs

user_data: The default value is None, in which case only the model calibration set is plotted. User provided sample data describing the oxide composition of one or more samples. Multiple samples can be passed as an BatchFile object or pandas DataFrame. A single sample can be passed as a pandas Series.

model: The default value is ‘all’, in which case all model calibration datasets will be plotted. Otherwise, any model can be plotted by passing the name of the model desired (e.g., ‘Liu’). Multiple models can be plotted by passing them as strings within a list (e.g., [‘Liu’, ‘Dixon’])

plot_type: The default value is ‘TAS’, which returns a total alkalis vs silica (TAS) diagram. Any two oxides can be plotted as an x-y plot by setting plot_type=‘xy’ and specifying x- and y-axis oxides, e.g., x=‘SiO2’, y=‘Al2O3’.

zoom: The default is None in which case axes will be set to the default of 35≤x≤100 wt% and 0≤y≤25 wt% for TAS type plots and the best values to show the data for xy type plots. The user can pass “user_data” to plot the figure where the x and y axes are scaled down to zoom in and only show the region surrounding the user_data. A list of tuples may be passed to manually specify x and y limits. Pass in data as [(x_min, x_max), (y_min, y_max)]. For example, the default limits here would be passed in as [(35,100), (0,25)].

save_fig: The default value is False, in which case the plot will be generated and displayed but not saved. If the user wishes to save the figure, the desired filename (including the file extension, e.g., .png) can be passed here. Note that all plots in this Jupyter notebook can be saved by right clicking the plot and choosing “Save Image As...”.

Outputs

The function returns fig and axes matplotlib objects defining a TAS or x-y plot of user data and model calibration data.

Using the functionality built into python and the matplotlib library, user data can be plotted on its own at any time, including before any calculations are performed. Almost any plot type imaginable can be produced, and users should refer to the maptlotlib documentation (https://matplotlib.org/3.2.1/index.html) if more complex plotting is desired.

3.5Calculating dissolved volatile concentrations¶

The calculate_dissolved_volatiles() function calculates the concentration of dissolved H₂O and CO₂ in the melt at a given pressure-temperature condition and with a given H₂O-CO₂ fluid composition, defined as the mole fraction of H₂O in an H₂O-CO₂ fluid (XH₂O^fluid). The default MagmaSat model relies on the underlying functionality of MELTS, whose basic function is to calculate the equilibrium phase assemblage given the bulk composition of the system and pressure-temperature conditions. To calculate dissolved volatile concentrations thus requires computing the equilibrium state of a system at fixed pressure and temperature over a range of bulk volatile concentrations until a solution is found that satisfies the user defined fluid composition.

First, the function makes an initial guess at the appropriate bulk volatile concentrations by finding the minimum dissolved volatile concentrations in the melt at saturation, while asserting that the weight fraction of H₂O/(H₂O+CO₂) in the system is equal to the user input mole fraction of H₂O/(H₂O+CO₂) in the fluid. This is done by increasing the H₂O and CO₂ concentrations appropriately until a fluid phase is stable. Once fluid saturation is determined, the code then performs directional, iterative, and progressively more refined searches, increasing the proportion of H₂O or CO₂ in the system if the mole fraction of H₂O calculated in the fluid is greater than or less than that defined by the user, respectively. Four iterative searches are performed; the precision of the match between the calculated and defined XH₂O^fluid increases from 0.1 in the first iteration to 0.01, 0.001, and finally to 0.0001. Thus, the calculated dissolved volatile concentrations correspond to a system with XH₂O^fluid within 0.0001 of the user defined value.

For non-MagmaSat models, dissolved volatile concentrations are calculated directly from model equations.

Method structure

Single sample: calculate_dissolved_volatiles(sample, temperature, pressure, X_fluid=1, verbose=False, model='MagmaSat').result

BatchFile batch process: myfile.calculate_dissolved_volatiles(temperature, pressure, X_fluid=1, print_status=True, model='MagmaSat')

Standard inputs

sample, temperature, pressure, X_fluid, model (see Section 3.1).

Unique optional inputs

verbose: Only for single sample calculations. Default value is False, in which case H₂O and CO₂ concentrations are returned. If set to True, additional parameters are returned in a dictionary: H₂O and CO₂ concentrations in the fluid in mole fraction, temperature, pressure, and proportion of the fluid in the system in wt%.

print_status: Only for batch calculations. The default value is True, in which case the progress of the calculation will be printed to the terminal. The user may desire to see the status of the calculation, as this particular function can be quite slow, averaging between 3-5 seconds per sample.

Calculated outputs

If the single-sample method is used, a dictionary with keys ‘H2O’ and ‘CO2’ corresponding to the calculated dissolved H₂O and CO₂ concentrations in the melt is returned (plus additional variables ‘temperature’ in °C, ‘pressure’ in bars, ‘XH2O_fl’, ‘XCO2_fl’, and ‘FluidProportion_wtper’ (the proportion of the fluid in the system in wt%) if verbose is set to True).

If the BatchFile method is used, a pandas DataFrame is returned with sample information plus calculated dissolved H₂O and CO₂ concentrations in the melt, the fluid composition in mole fraction, and the proportion of the fluid in the system in wt%. Pressure (in bars) and Temperature (in °C) columns are always returned.

3.6Calculating equilibrium fluid compositions¶

The calculate_equilibrium_fluid_comp() function calculates the composition of a fluid phase in equilibrium with a given silicate melt with known pressure, temperature, and dissolved H₂O and CO₂ concentrations. The calculation is performed simply by calculating the equilibrium state of the given sample at the given conditions and determining if that melt is fluid saturated. If the melt is saturated, fluid composition and mass are reported back. If the calculation finds that the melt is not saturated at the given pressure and temperature, values of 0.0 will be returned for the H₂O and CO₂ concentrations in the fluid.

Method structure

Single sample: calculate_equilibrium_fluid_comp(sample, temperature, pressure, verbose=False, model='MagmaSat').result

BatchFile batch process: myfile.calculate_equilibrium_fluid_comp(temperature, pressure=None, print_status=False, model='MagmaSat')

Standard inputs

sample, temperature, pressure, model (see Section 3.1).

Unique optional inputs

verbose: Only for single sample calculations. Default value is False, in which case H₂O and CO₂ concentrations in the fluid in mol fraction are returned. If set to True, additional parameters are returned in a dictionary: H₂O and CO₂ concentrations in the fluid, mass of the fluid in grams, and proportion of the fluid in the system in wt%.

print_status: Only for batch calculations. The default value is False. If True is passed, the progress of the calculation will be printed to the terminal.

Calculated outputs

If the single-sample method is used, a dictionary with keys ‘H2O’ and ‘CO2’ is returned (plus additional variables ‘FluidMass_grams’ and ‘FluidProportion_wtper’ if verbose is set to True).

If the BatchFile method is used, a pandas DataFrame is returned with sample information plus calculated equilibrium fluid compositions, mass of the fluid in grams, and proportion of the fluid in the system in wt%. Pressure (in bars) and Temperature (in °C) columns are always returned.

Below we calculate equilibrium fluid compositions for all samples at a single temperature of 900 °C and a single pressure of 1,000 bars. Note that some samples in this dataset have quite low volatile concentrations (e.g., the Tucker et al. (2019) basalts from Kilauea), and so are below saturation at this P-T condition. The fluid composition for undersaturated samples is returned as values of 0 for both H₂O and CO₂.

Below, we calculate equilibrium fluid compositions for the same dataset using temperatures and pressures as defined in the input data (Table 3). Note that Samples “samp. HPR3-1_XL-3” and “samp. HPR3-1_XL-4_INCL-1” have a user-defined value of 0.0 for temperature and pressure, respectively. VESIcal automatically skips the calculation of equilibrium fluids for these samples and returns a warning to the user, which are both printed to the terminal below and appended to the “Warnings” column in the returned data.

3.6.1Converting fluid composition units¶

The fluid composition is always returned in units of mol fraction. Two functions exist to transform only the H₂O-CO₂ fluid composition between mol fraction and wt% and can be applied to returned data sets from calculations. Both functions require that the user provide the dataframe containing fluid composition information plus the names of the columns corresponding to the H₂O and CO₂ concentrations in the fluid. The default values for column names are set to those that may be returned by VESIcal core calculations, such that they need not be specified unless the user has changed them or is supplying their own data (e.g., imported data not processed through a core calculation).

Method structure

Mol fraction to wt%: fluid_molfrac_to_wt(data, H2O_colname='XH2O_fl_VESIcal', CO2_colname='XCO2_fl_VESIcal')

Wt% to mol fraction: fluid_wt_to_molfrac(data, H2O_colname='H2O_fl_wt', CO2_colname='CO2_fl_wt')

Required inputs

data: A pandas DataFrame containing columns for H₂O and CO₂ concentrations in the fluid.

Optional inputs

H2O_colname and CO2_colname: The default values are ‘XH2O_fl’ and ‘XCO2_fl’ if input data are in mol fraction or ‘H2O_fl_wt’ and ‘CO2_fl_wt’ if the data are in wt%. Strings containing the name of the columns corresponding to the H₂O and CO₂ concentrations in the fluid.

Calculated outputs

The original data passed plus newly calculated values are returned in a DataFrame.

3.7Calculating saturation pressures¶

The calculate_saturation_pressure() function calculates the minimum pressure at which a given silicate melt with known temperature and H₂O and CO₂ concentrations would be saturated with fluid. For MagmaSat, this is calculated by finding the pressure at which the smallest amount of vapor is present. This function also calculates the composition of the vapor in equilibrium with the melt at those conditions.

The function works by calculating the equilibrium state of the given melt at very high pressure (20,000 bars). If no fluid is present at this pressure, the melt is undersaturated, and pressure is decreased in steps of 1,000 bars until the mass of vapor is >0 grams. If fluid is present, the saturation limit is found by increasing the pressure iteratively until the point at which no fluid is present. At this point, the pressure space is narrowed and searched in steps of 100 bars and then in steps of 10 bars until the saturation pressure is found. Thus, these calculations are accurate to 10 bars.

For non-MagmaSat models, we use Brent’s minimization method (via scipy’s root_scalar optimization function) to find the pressure that satisfies the computational constraints. This is achieved by iterative calculation of the dissolved volatile concentration over a range of pressures and minimizing the difference between computed and given concentrations. This is only practical for non-MagmaSat models, where the dissolved volatiles calculation is extremely fast.

Method structure

Single sample: calculate_saturation_pressure(sample, temperature, verbose=False, model='MagmaSat').result

BatchFile batch process: myfile.calculate_saturation_pressure(temperature, print_status=True, model='MagmaSat')

Standard inputs

sample, temperature, model (see Section 3.1).

Unique optional inputs

verbose: Only for single sample calculations. Default value is False, in which case the saturation pressure in bars is returned. If set to True, additional parameters are returned in a dictionary: saturation pressure in bars, H₂O and CO₂ concentrations in the fluid, mass of the fluid in grams, and proportion of the fluid in the system in wt%.

print_status: Only for batch calculations. The default value is True, in which case the progress of the calculation will be printed to the terminal.

Calculated outputs

If the single-sample method is used, the saturation pressure in bars is returned as a numerical value (float) (plus additional variables ‘XH2O_fl’, ‘XCO2_fl’, ‘FluidMass_grams’, and ‘FluidProportion_wtper’ if verbose is set to True).

If the BatchFile method is used, a pandas DataFrame is returned with sample information plus calculated saturation pressures, equilibrium fluid compositions, mass of the fluid in grams, and proportion of the fluid in the system in wt%. Temperature (in °C) is always returned.

3.8Calculating isobars and isopleths¶

In this example, we demonstrate how isobars (lines of constant pressure) and isopleths (lines of constant fluid composition) can be calculated for any one composition. A single melt composition can be extracted from a loaded batch file, or a composition can be entered by hand and stored within a dictionary. Due to computational intensity, isobars and isopleths can only be computed for one sample composition at a time.

Once a single composition is defined, conditions over which to calculate isobars and isopleths must be specified. The generated plot is isothermal, so only one temperature can be chosen. Isobars and isopleths can be calculated for any number of pressures or XH₂O^fluid values, respectively, passed as lists.

The calculation is performed by iterating through possible concentrations of H₂O and CO₂ and calculating the equilibrium state for the system. The iteration begins at a fixed H₂O concentration, increasing the CO₂ concentration in steps of 0.1 wt% until a fluid phase is stable. The H₂O concentration is then increased by 0.5 wt% and CO₂ is again increased from 0 until a fluid phase is stable. This process is repeated for H₂O values ranging from 0–15 wt%. The H₂O and CO₂ concentrations from each system for which a fluid phase was found to be stable are saved and written to a pandas DataFrame, which is returned upon completion of the calculation.

Isobars and isopleths are computed at fixed H₂O-CO₂ points for any given pressure. To generate curves using the MagmaSat model, polynomials are fit to computed points using numpy’s polyfit method. This can be optionally disabled by setting smooth_isobars or smooth_isopleths to False. The curvature of the isobars depends strongly on the number of points used to fit a polynomial, deemed “control points”, with curve fits becoming more accurate to the model as the number of control points increases. We found that above five control points, changes to the shape of the curve fits becomes negligible. Thus, as a compromise between accuracy and computation time, and to maintain consistency, MagmaSat isobars are always computed with 5 control points at XH₂O^fluid values of 0, 0.25, 0.5, 0.75, and 1. Because non-MagmaSat models compute extremely quickly, all non-MagmaSat models use 51 control points per isobar and do not utilize polynomial fits to the data by default.

Method structure

Only single sample calculations. calculate_isobars_and_isopleths(sample, temperature, pressure_list, isopleth_list=None, smooth_isobars=True, smooth_isopleths=True, print_status=True, model="MagmaSat").result

Standard inputs

sample, temperature, model (see Section 3.1).

Unique required inputs

pressure_list: A list of all pressures in bars at which to calculate isobars. If only one value is passed it can be as float instead of list.

Unique optional inputs

isopleth_list: The default value is None in which case only isobars will be calculated. A list of all fluid composition values, in mole fraction H₂O (XH₂O^fluid), at which to calculate isopleths. Values can range from 0–1. If only one value is passed it can be as float instead of list. N.b. that, due to the method of isobar smoothing using control points as outlined above, each isopleth value passed here not equal to one of the five standard control point values (0, 0.25, 0.5, 0.75, or 1) will result in an an additional control point being used to smooth the isobars. Thus, entering additional isopleth values results not only in more isopleth outputs but also in “smoother” (i.e., more well constrained) isobars.

smooth_isobars and smooth_isopleths: The default value for both of these arguments is True, in which case polynomials will be fit to the computed data points.

print_status: The default value is True. If True, the progress of the calculations will be printed to the terminal.

Calculated outputs

The function returns two pandas DataFrames: the first has isobar data, and the second has isopleth data. Columns in the isobar dataframe are ‘Pressure’, ‘H2Omelt’, and ‘CO2melt’, corresponding to pressure in bars and dissolved H₂O and CO₂ in the melt in wt%. Columns in the isopleth dataframe are ‘XH2O_fl’, ‘H2O_liq’, and ‘CO2_liq’, corresponding to XH₂O^fluid and dissolved H₂O and CO₂ in the melt in wt%.

Next, the H₂O and CO₂ dissolved in the melt at saturation is calculated at the specified temperature and over the range of specified pressures. Note that, because this function calculates two things (isobars and isopleths), two variable names must be given (below, “isobars, isopleths”). This calculation can be quite slow, and so it is recommended to set print_status to True.

3.9Calculating degassing paths¶

A degassing path is a series of volatile concentrations both in the melt and fluid that a magma will follow during decompression. In the calculation, the saturation pressure is computed, and then the system is equilibrated along a trajectory of decreasing pressure values at discrete steps. The default number of steps to calculate is 50, but this can be defined by the user by setting the argument steps to any integer value. A detailed explanation of how non-MagmaSat models handle the calculation of mixed-fluid composition can be found in the supplement (Supplementary Text S2). If so desired, this calculation can be performed for any initial pressure, but the default is the saturation pressure. If a pressure is specified that is above the saturation pressure, the calculation will simply proceed from the saturation pressure, since the magma cannot degas until it reaches saturation.

Completely open-system, completely closed-system or partially open-system degassing paths can be calculated by specifying what proportion of the fluid to fractionate. The fluid fractionation value can range between 0 (closed-system: no fluid is removed, all is retained at each pressure step) and 1 (open-system: all fluid is removed, none is retained at each pressure step). Closed and partially open-system runs allow the user to specify the initial presence of exsolved fluid that is in equilibrium with the melt at the starting pressure.

Method structure

Only single-sample calculations. calculate_degassing_path(sample, temperature, pressure="saturation", fractionate_vapor=0.0, init_vapor=0.0, steps=50, model="MagmaSat").result

Standard inputs

sample, temperature, model (see Section 3.1).

Unique optional inputs

pressure: The pressure at which to begin the degassing calculations, in bars. Default value is ‘saturation’, which runs the calculation with the initial pressure at the saturation pressure. If a pressure greater than the saturation pressure is input, the calculation will start at saturation, since this is the first pressure at which any degassing will occur.

fractionate_vapor: Proportion of vapor removed at each pressure step. Default value is 0.0 (completely closed-system degassing). Specifies the type of calculation performed, either closed system (0.0) or open system (1.0) degassing. If any value between <1.0 is chosen, user can also specify the ‘init_vapor’ argument (see below). A value in between 0 and 1 will remove that proportion of vapor at each step. For example, for a value of 0.2, the calculation will remove 20% of the vapor and retain 80% of the vapor at each pressure step.

init_vapor: Default value is 0.0. Specifies the amount of vapor (in wt%) coexisting with the melt before degassing.

steps: Default value is 50. Specifies the number of steps in pressure space at which to calculate dissolved volatile concentrations.

Calculated outputs

The function returns a pandas DataFrame with columns as: ‘Pressure_bars’, ‘H2O_liq’ and ‘CO2_liq’ (the concentration of H₂O and CO₂ in the melt, in wt%), ‘XH2O_fl’ and ‘XCO2_fl’ (the composition of the H₂O-CO₂ fluid, in mol fraction), and ‘FluidProportion_wt’ (the proportion of fluid in the fluid-melt system, in wt%).

Note: The following two cells can take up to a minute to execute.

3.10Plotting¶

After calculating isobars, isopleths, and degassing paths, any or all of these may be plotted in an H₂O versus CO₂ plot with one simple function call. The plot will be printed directly in the notebook or, if the code is run as script in a command line, the plot will appear it its own window, at which point it can be saved as an image file. VESIcal’s plot function takes in lists of pandas DataFrames with calculated isobar, isopleth, and degassing path information (e.g., output from calculate_isobars_and_isopleths() or calculate_degassing_path()) and plots data as isobars (lines of constant pressure), isopleths (lines of constant fluid composition), and degassing paths (lines indicating the concentrations of H₂O and CO₂ in a melt equilibrated along a path of decreasing pressure).

Labels can be assigned to isobars, isopleths, and/or degassing paths separately. Any or all of these data can be passed to the plot function. Multiple sets of plottable data can be passed. For example, isobars calculated with two different models can be passed to the isobars argument as a list.

VESIcal’s plotting function is entirely based on python’s matplotlib library, which comes standard with many installations of python. With matplotlib, users can create a large variety of plots (note that direct matplotlib functionality is used to create custom plots in several of this manuscript’s supplementary Jupyter notebooks), and users should refer to the maptlotlib documentation (https://matplotlib.org/3.2.1/index.html) if more complex plotting is desired. If preferred, VESIcal outputs can be saved to an Excel or CSV file (see Section 3.12), and plotting can be done in any plotting program desired (e.g., MS Excel).

The function returns both fig and axes matplotlib objects, which can be further edited by the user or plotted directly. Following matplotlib convention, the results of plot() should be saved to objects such as fig, ax as:

fig, ax = v.plot([options])

where [options] represents any optional inputs as defined here. Variables fig and ax can then be edited further using matplotlib tools. For example, the user might wish to set the minimum x-axis value to 0.5 as:

ax.set_xlim(left=0.5)

In Jupyter Notebook, a plot is automatically shown, but in the command line, the plot will only display after executing v.show().

Method structure

plot(isobars=None, isopleths=None, degassing_paths=None, custom_H2O=None, custom_CO2=None, isobar_labels=None, isopleth_labels=None, degassing_path_labels=None, custom_labels=None, custom_colors="VESIcal", custom_symbols=None, markersize=10, save_fig=False, extend_isobars_to_zero=True, smooth_isobars=False, smooth_isopleths=False)

Optional inputs

isobars: DataFrame object containing isobar information as calculated by calculate_isobars_and_isopleths(). Or a list of DataFrame objects.

isopleths: DataFrame object containing isopleth information as calculated by calculate_isobars_and_isopleths(). Or a list of DataFrame objects.

degassing_paths: List of DataFrames with degassing information as generated by calculate_degassing_path().

custom_H2O: List of floats or array-like shapes of H₂O concentration values to plot as points. For example myfile.data['H2O'] is one array-like shape (here, pandas.Series) of H₂O values. Must be passed with custom_CO2 and must be same length as custom_CO2.

custom_CO2: List of floats or array-like shapes of CO₂ values to plot as points. For example myfile.data['CO2'] is one array-like shape of CO₂ values. Must be passed with custom_H2O and must be same length as custom_H2O.

isobar_labels: Labels for the plot legend. Default is None, in which case each plotted line will be given the generic legend name of “Isobars n”, with n referring to the nth isobars passed. Isobar pressure is given in parentheses. The user can pass their own labels as a list of strings. If more than one set of isobars is passed, the labels should refer to each set of isobars, not each pressure.

isopleth_labels: Labels for the plot legend. Default is None, in which case each plotted isopleth will be given the generic legend name of “Isopleth n”, with n referring to the nth isopleths passed. Isopleth XH₂O values are given in parentheses. The user can pass their own labels as a list of strings. If more than one set of isopleths is passed, the labels should refer to each set of isopleths, not each XH₂O value.

degassing_path_labels: Labels for the plot legend. Default is None, in which case each plotted line will be given the generic legend name of “Pathn”, with n referring to the nth degassing path passed. The user can pass their own labels as a list of strings.

custom_labels: Labels for the plot legend. Default is None, in which case each group of custom points will be given the generic legend name of “Customn”, with n referring to the nth degassing path passed. The user can pass their own labels as a list of strings.

custom_colors and custom_symbols: Custom colors and symbol shapes can be specified for (custom_H2O, custom_CO2) points. A list of color values or symbol types readable by Matplotlib (see Matplotlib documentation) can be entered. The length of this list must be equal to the lengths of custom_H2O and custom_CO2. If nothing is specified for custom_colors, VESIcal’s default colors will be used. If nothing is specified for custom_symbols, all points will be plotted as filled circles.

markersize: The size of the symbols can be specified here. If not specified, the default value is marker size 10.

save_fig: Default value is False, in which case the figure will not be saved. If a string is passed, the figure will be saved with the string as the filename. The string must include the file extension.

Advanced inputs: Most users will not need to use these inputs.

extend_isobars_to_zero: If set to True (the default), isobars will be extended to the plot axes, which are at x=0 and y=0, even if there is a finite solubility at zero partial pressure.

smooth_isobars and smooth_isopleths: If set to True, isobar or isopleth data will be fit to a polynomial and plotted. If set to False (the default), the raw input data will be plotted. Note that MagmaSat calculate_isobars_and_isopleths() calculations return already “smoothed” data (that is, the raw data are fit to polynomials before being returned). Raw “unsmoothed” data can be returned by MagmaSat calculate_isobars_and_isopleths() (see documentation on this method).

Calculated outputs

The function returns fig and axes matplotlib objects defining a plot with x-axis as H₂O wt% in the melt and y-axis as CO₂ wt% in the melt. Isobars, or lines of constant pressure at which the sample magma composition is saturated, and isopleths, or lines of constant fluid composition at which the sample magma composition is saturated, are plotted if passed. Degassing paths, or the concentration of dissolved H₂O and CO₂ in a melt equilibrated along a path of decreasing pressure, is plotted if passed.

3.10.1A simple example: Isobars and isopleths¶

Here we plot the isobars at 1,000, 2,000, and 3,000 bars and isopleths at 0.25, 0.5, and 0.75 XH₂O^fluid calculated for sample ‘10*’ at 1,200 °C in Section 3.8 onto one plot.

When plotting isobars and isopleths via MagmaSat, the values calculated by calculate_isobars_and_isopleths() are used to calculate polynomial fits using Numpy’s polyfit. These polynomial fits, not the raw calculated data, are what have been plotted above. This method of fitting polynomial curves to these data is common in the literature e.g., Newman & Lowenstern, 2002Iacono-Marziano et al., 2012Iacovino et al., 2013 and is likely a very close approximation of the true saturation surface. Non-MagmaSat models do not calculate polynomial fits by default, but this can be done by passing smooth_isobars=True and smooth_isopleths=True to plot().

A user may wish to apply custom formatting to the plot, in which case the polynomial fits can be calculated and returned as a pandas DataFrame, which the user can then plot up manually using Matplotlib, MS Excel, or some other preferred method. To calculate polynomial fits to isobar and isopleth data, isobars and isopleths can be passed to smooth_isobars_and_isopleths(). For this advanced case, we refer the reader to the documentation.

3.10.2A simple example: Degassing paths¶

Here we plot all four degassing paths calculated for sample ‘10*’ at 1,200 °C in Section 3.9 onto one plot. We designate labels of “Open”, “Half”, “Closed”, and “Exsolved” for the legend.

3.10.3Plotting multiple calculations¶

One of the major advantages to VESIcal over any other modeling tool is the ability to quickly calculate and plot multiple calculations. VESIcal’s plot() function is built on top of the popular Matplotlib python library and is designed to work with any VESIcal generated data. It can automatically plot and label one or multiple calculations. In addition, it can plot, as a scatter plot, any x-y points. The plot function always generates plots with H₂O on the x-axis and CO₂ on the y-axis. scatterplot() will take in and plot any x-y data with custom x- and y-axis labels. Generating other commonly used petrologic plots (e.g. Harker style diagrams) is already possible with Matplotlib, and so VESIcal does not duplicate this functionality, however this may be added in future updates.

It may be tempting to plot multiple calculations on multiple samples and compare them, however we strongly caution against plotting data that do not correspond. For example, isobars and isopleths are calculated isothermally. If degassing paths are also plotted, the user should ensure that the degassing paths were calculated at the same temperature as the isobars and isopleths.

3.10.3.1Isobars, isopleths, and degassing paths¶

In this example we will use data imported in Section 3.4 and calculations performed in Sections 3.7 and 3.8. Of course, all of the data calculated with VESIcal can be exported to an Excel or CSV file for manipulation and plotting as desired. However, some examples of plotting that can be done within this notebook or in a python script are shown below. Here we plot:

• Isobars calculated at 1200 °C and pressures of 1,000, 2,000, and 3,000 bars for sample 10*
• Isopleths calculated at 1200 °C and XH₂O^fluid values of 0, 0.25, 0.5, 0.75, and 1 for sample 10*
• An open-system degassing path for sample 10*
• A closed-system degassing path for sample 10*

3.10.3.2Isobars, isopleths, and degassing paths for multiple samples¶

First, we will calculate some new data for two different samples: a basanite sample KI-07 from Iacovino et al., 2016 and a rhyolite sample samp. P1968a from Myers et al., 2019. For both samples we will calculate and then plot:

• Isobars and isopleths at 1100 °C, pressures of 1,000 and 2,000 bars and fluid compositions of XH₂O^fluid of 0.25, 0.5, and 0.75
• Closed-system degassing paths at 1100 °C

3.11Model hybridization (Advanced)¶

One of the advantages of implementing the solubility models in a generic python module is the flexibility this affords the user in changing the way solubility models are defined and used. In particular, the structure allows any combination of pure-fluid models to be used together in modeling mixed-fluids, and fugacity or activity models can be quickly changed without modifying code. This allows advanced users to see how changing a fugacity or activity model implemented in any particular solubility model would affect model results. Instructions for hybridizing models can be found in Supplemental Jupyter notebook S10.

3.12Exporting data¶

Once batch calculations have been performed, they can be exported to an Excel or CSV file with the save_excel() and save_csv() commands. These operations require that the user define a filename (what to name your new file) and a list of the calculation results to save to this file or files.

Note that this requires that calculations have been assigned to variable names, which has been done in all of the given examples. For example, to calculate saturation pressures of an imported file saved to the variable ‘myfile’ and simply print the output, the user can type myfile.calculate_saturation_pressures([options]), where ‘[options]’ are the required and optional inputs. However, to save this result to a variable (e.g., called ‘my_satPs’) so that it can be accessed later, the correct python syntax would be my_satPs = myfile.calculate_saturation_pressures([options]).

Multiple calculations can be saved at once. If saving to an Excel file, each calculation is saved as its own sheet within a single file. If desired, the user can define the names of each of these sheets. If not specified, the sheets will be named ‘Original_User_Data’, which contains the original input data, and then ‘CalcN’ where N is the nth calculation in a list of calculations. If saving multiple calculations to a CSV file, each calculation will be saved to its own CSV file, and a file name for each of these is required.

Advanced users note that the calculations argument takes in any pandas DataFrame object, meaning this functionality is not limited to VESIcal’s prescribed outputs. The save_excel() and save_csv() methods use the pandas to_excel and to_csv methods, however not all options are implemented here. If saving to a CSV file, any arguments that can be passed to pandas to_csv method may be passed to VESIcal’s save_csv().

Method structures

save_excel(filename, calculations, sheet_name=None)

save_csv(filenames, calculations)

save_excel() Required inputs

filename (Excel): Name of the file to create. The extension (.xlsx) should be included along with the name itself, all in quotes (e.g., filename='myfile.xlsx').

calculations: A list of variables containing calculated outputs from any of the core BatchFile functions: calculate_dissolved_volatiles(), calculate_equilibrium_fluid_comp(), and calculate_saturation_pressure(). This must be passed as a list type variable, even if only one calculation is given. This is done by enclosing the variable in square brackets (e.g., calculations=[my_calculation]).

save_excel() Optional inputs

sheet_name: The default value is None, in which case sheets will be saved as ‘Original_User_data’ (the data input by the user) followed by ‘CalcN’ where N is the nth calculation in calculations. Otherwise, a list of names for the sheets can be passed, with the names in quotes (e.g. sheet_name=['SaturationPressures']). ‘Original_User_data’ will always be saved as the first sheet.

save_csv() Required inputs

filenames (CSV): Name of the file or files to create. The extension (.csv) should be included. If more than one filename is passed, it should be passed as a list. This is done by enclosing the filenames in square brackets (e.g., filenames=["file1.csv", "file2.csv"]).

calculations: same as for save_excel(). Must be same length as filenames.

Calculated outputs

An Excel or CSV file or files will be saved to the active directory (i.e., the same folder as this manuscript notebook or wherever the code is being used).

Here we save five of the calculations performed on an imported data file earlier in this manuscript. The original user-input data are stored in the BatchFile object ‘myfile’. In the following line we use the method save_excel() to save the original data and a list of calculations given by the calculations argument to an Excel file.

3.12.1Saving data for re-import into VESIcal¶

In many cases, it may be preferable to compute large amounts of data using VESIcal and then reimport them, either to perform more analysis or to plot the data. Likewise, a user may wish to compute data in VESIcal and then send the results to a colleague, who can then re-import that data into VESIcal directly. For this case, we suggest using python’s pickle package (https://wiki.python.org/moin/UsingPickle). Any python object, such as the results of a VESIcal calculation, can be “pickled” or saved as a python-readable file. To use pickle, users must first import the pickle module, then “dump” the desired contents to a pickle file. The pickled data can be accessed by “loading” the pickled file.

Below we pickle our computed dissolved volatile concentrations by dumping our variable dissolved to a pickle file that we name “dissolved.p”.

In another python file or terminal session, dissolved can be loaded back in via:

4Discussion and Applications¶

4.1Compositional Variation Within Datasets and Best Practices¶

It has been clearly shown that the composition of a melt plays a strong role in determining the solubility of H₂O and CO₂ in magmas Papale et al., 2006Moore, 2008Ghiorso & Gualda, 2015Wieser et al., 2021. Thus, compositional variance must be accounted for in any study examining solubility in multiple samples. A key use case where VESIcal can facilitate the adoption of this practice is in melt inclusion (MI) studies; specifically, where a single suite of MI with multiple melt compositions is examined using solubility models to interrogate magmatic degassing processes. Prior to the availability of VESIcal, the difficulty associated with performing multiple model calculations on multiple samples resulted in very few studies accounting for any compositional variance within their datasets. Indeed, until now, it has been difficult to even assess whether the potentially minimal compositional variance within a suite of melt inclusions from a single volcanic eruption would have any measurable effect on solubilities calculated for different MI.

Using VESIcal, we can address the question: what is the quantitative effect of compositional variation within a single suite of melt inclusions upon calculated melt inclusion saturation pressures? And, how does this affect conclusions that might be drawn regarding volcanic degassing and eruptive processes? To investigate this, we use a dataset of basaltic melt inclusions from Cerro Negro volcano, Nicaragua Roggensack, 2001. The compositional variation of these MI Figure 10, while relatively restricted, results in quite variable mixed-fluid solubilities from sample to sample. To determine the end-member compositions within the dataset corresponding to the samples with the maximum and minimum combined H₂O-CO₂ solubilities, isobars were computed at 1200 °C and 3,000 bars for all samples using the MagmaSat model in VESIcal. Maximum and minimum samples were taken as the isobar curves with the smallest and largest integral (area under the curve). We refer to this value as the “integrated mixed-volatile solubility” value, IMS, in units of concentration squared. The samples that produced maximum and minimum integrated solubilities are shown in Figures 10 and 11 in blue and green, respectively (sample 41b*, IMS=0.81 and 36a*, IMS=0.66 wt%² at 3,000 bars). A composition representing the average of all MI in the dataset is shown in orange (“Average Sample”, IMS=0.70 wt%² at 3,000 bars). A Jupyter notebook to reproduce these calculations is provided in the supplement (Supplementary Jupyter Notebook S8).

At all pressures, the integrated mixed-volatile solubility across the Cerro Negro dataset varies as much as 10% relative Figure 11. For these MI, this results in as much as 11.5% relative error in the calculation of saturation pressures (average error for the entire dataset of 6.8% relative). It is noteworthy that this error is not systematic either in terms of absolute value or sign. For example, when calculated using their own compositions, saturation pressures for maximum and minimum samples 41b* and 36a* are 3050 and 3090 bars, respectively. But, saturation pressures calculated for both of these MI using the dataset’s average composition are 3020 and 3250 bars, respectively. That is an error of -30 and +160 bars or -1% and +5% respectively. Errors in these calculations, thus, may be quite small. But, in any case, removing this error completely is a simple task using VESIcal, and so we recommend that studies adopt the practice of calculating volatile solubilities (and associated values) in melts using the composition unique to each melt investigated.

Even in cases where solubility values (e.g., saturation pressures) are not calculated, the error highlighted above plagues any isobar diagram over which multiple melt compositions are plotted (e.g., Figure 11). Alternative plots to the commonly used H₂O-CO₂ diagram are shown in Figure 12, in which the same dataset is plotted in terms of computed saturation pressure (at 1200 °C calculated with VESIcal using MagmaSat) versus dissolved H₂O, dissolved CO₂, and fluid composition (as XH₂O^fluid calculated with VESIcal using MagmaSat). These plots avoid the issues discussed above as they are compositionally independent, since the saturation pressure is calculated individually for each sample composition. Degassing trends are more accurately represented; H₂O and CO₂ concentrations lie along expected degassing trends with much less scatter than the H₂O-CO₂ plot. We can also see from this figure that the fluid composition during this eruption at Cerro Negro remained relatively constant at XH₂O^fluid $\sim$0.8 from reservoir to surface, suggesting a scenario approaching closed-system degassing (i.e., melt volatile concentrations are buffered by the co-existing fluid composition). This is discussed in more detail in the companion paper Wieser et al., 2021.

4.2Model Comparisons¶

One of the possible workflows enabled through VESIcal is the ability to compute and compare (numerically and graphically) results from several models at once. To illustrate this point, we will take two single samples within the calibrated compositional range of several models, calculate isobars at multiple pressures, and plot the results. This is a common way to compare the solubility surface computed by different models for a single melt composition, and it is particularly useful since it quickly highlights the significant variation that exists between published models. The results of this exercise are shown here, and a Jupyter notebook to reproduce the code and calibration checks is available in the Supplement (Supplementary Jupyter Notebook S9).

We use a fictitious alkali basalt that we name “alkbasalt” and a fictitious rhyolite whose compositions are given in Table 12. The use of VESIcal’s calib_plot() function (see supplement) illustrates that the composition of the alkali basalt is within the compositional calibration ranges of four mixed-fluid solubility models: MagmaSat, Iacono-Marziano, Dixon, and ShishkinaIdealMixing. The rhyolite is within the ranges of MagmaSat and Liu. Isobars were calculated with these models at 1200 °C for alkbasalt and 800 °C for rhyolite and pressures of 500, 1,000, and 2,000 bars, using the below code:

It is immediately clear from Figure 10 that major disagreement exists between these models. For the alkali basalt, MagmaSat and Dixon show the best agreement, particularly at pressures <2000 bars. However, the mismatch between these models (and, indeed, between all models) increases with pressure. The Iacono-Marziano model is calibrated for highly depolymerized alkali basalts resulting in an increased capacity of the melt to dissolve CO$_3^{2-}$. That may explain why this model predicts significantly higher CO₂ solubilities at XH₂O^fluid values approaching 0.

The ShishkinaIdealMixing model displays nearly linear isobars, with finite solubility below ~1 wt% dissolved H₂O. This is a consequence of the model calibration; the pure-H₂O solubility expression of ShishkinaIdealMixing is not calibrated with any experiments at low PH₂O. This results in a finite solubility at low dissolved H₂O concentrations, such that the zero-pressure solubility is not zero. This produces significant model error at low but non-zero values of XH₂O^fluid. Thus, we caution the user against using the ShishkinaIdealMixing model at low but non-zero XH₂O^fluid or when fluids deviate far from pure H₂O or pure CO₂. In general, the Shishkina model should only be used for modeling pure-H₂O or pure-CO₂ fluids. This is discussed in more detail in Wieser et al. (2021).

The models of MagmaSat and Liu show a similar level of disagreement for H₂O-CO₂ solubility in the rhyolite, with Liu predicting much higher dissolved CO₂ concentrations at low XH₂O^fluid (<20 relative% or ~220 ppm at XH₂O^fluid=0.1).

4.3Sensitivity and error analysis¶

To date, very few studies have compared the sensitivity of their pressure estimates to the choice of solubility model, or propagated errors inherent to measurements of volatile concentrations in melts using SIMS, FTIR and Raman Spectroscopy into an error bar in terms of saturation pressure. In contrast, VESIcal allows users to import an Excel or CSV spreadsheet with each row containing the major element and volatile contents of each inclusion, as well as a temperature at which to evaluate solubility. Using the batch calculation functions, VESIcal will automatically calculate the saturation pressure for each row, using a user-specified model. Thus, users can more easily compare results from different solubility models, to robustly assess their applicability for the system of interest. Additionally, users could load a different spreadsheet, where the CO₂ and H₂O concentrations are adjusted to reflect the analytical uncertainty on the instrument used, allowing error bars on the saturation pressure to be calculated for every single inclusion. The modular and open-source nature of VESIcal also allows the user to combine the code with other Python3 modules. For example, users could utilize Markov chain Monte Carlo (MCMC; e.g., the python library emcee) methods to robustly calculate error distributions for each sample. In future releases, automatic sensitivity and error analysis on datasets and calculated results may be implemented directly within VESIcal, building on existing tools within the python community.

4.4Future development¶

VESIcal represents the first comprehensive volatile solubility modeling tool of its kind, including the feature that VESIcal is extensible. VESIcal is written so that implementing new or yet-to-be-implemented solubility models is as simple as possible. To implement a new model, python code describing the model equations needs to be written, and this model name needs to be added to a list of model names within the code. To make this as simple as possible such that the original authors of VESIcal are not the only people who can develop the code, planned future work includes the creation of detailed instructions (including instructional videos) illustrating this process.

Likewise, new features can be added at any time, and enthusiastic members of the community who wish to help bring such features to VESIcal are very welcome. Users can contribute to VESIcal’s code, implementing new models and new features, via github (https://github.com/kaylai/VESIcal). The repository is public, but we encourage users who wish to contribute to the code to fork the repository into their private workspace on github. Once edits to the code are complete, the new code can be added to VESIcal by creating a “Pull Request” inside of github. Changes and enhancements to VESIcal will correspond to a change in the code’s version number. The published version of the code documented in this manuscript and archived on Zenodo is version 1.0.1 (DOI: Iacovino et al. (2021)). Planned features not implemented in this release include: 1. Models to calculate sample oxygen fugacity from Fe²⁺/ΣFe and vice versa; 2. Additional volatiles such as sulfur; 3. More thermodynamic solubility models such as that of Papale et al. (2006); 4. Sensitivity and error analysis functions.

4.5How to cite VESIcal and its models¶

To cite computations done using VESIcal, please cite this manuscript, the VESIcal version number, as well as the model(s) used. Note that if a model was not specified during calculations, the default model of MagmaSat was used and should be cited as “MagmaSat Ghiorso & Gualda, 2015”. For example: “Calculations were performed using VESIcal v. 1.0.1, Iacovino et al., 2021 with the models of Shishkina et al. (2014) and Dixon (1997, "VolatileCalc").” The web-app always runs on the most up-to-date version of the VESIcal code, but it is best practice to note if the web-app was used (“Calculations were performed using the VESIcal web-app v. 1.0.1, Iacovino et al., 2021...”). We also encourage users to be as explicit as possible as to the conditions used for modelling. This includes stating the pressure, temperature, volatile concentration, and bulk magma composition used in modelling. In the best case, VESIcal users will provide their code (e.g., as a Jupyter Notebook or .py file) along with their publication such that it can be easily replicated.

5Conclusions¶

VESIcal is a thermodynamic mixed-volatile solubility engine designed to meet the growing computational needs of the igneous petrology community. Seven commonly used volatile solubility models are built into VESIcal, which employs the most diversely calibrated (chemically and in P-T space) of the group, MagmaSat Ghiorso & Gualda, 2015, as the default model. VESIcal can perform five core calculations with any mixed-fluid model and three core calculations with any model (mixed-fluid, CO₂-only, H₂O-only). VESIcal allows for automatic calculation of large datasets and robust built-in plotting capability.

Alongside model frameworks such as ENKI, VESIcal represents an early step forward toward creating a generalized thermodynamic framework to model whole scale magmatic processes. Such a framework builds upon the key tenets of VESIcal; namely: fundamental thermodynamic underpinning; inclusion of existing modeling strategies; python powered, open-source, and extensible code base; high usability at all levels; benchmarking and testing; and power as a responsive and predictive tool.