% Created with jtex v.1.0.10
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% AGUJournalTemplate.tex: this template file is for articles formatted with LaTeX
%
% This file includes commands and instructions
% given in the order necessary to produce a final output that will
% satisfy AGU requirements, including customized APA reference formatting.
%
% You may copy this file and give it your
% article name, and enter your text.
%
%
% Step 1: Set the \documentclass
%
%

%% To submit your paper:
% TODO - options - draft
\documentclass{agujournal2019}

\usepackage{lineno}
\usepackage[inline]{trackchanges} %for better track changes. final new option will compile document with changes incorporated.
\usepackage{soul}
\linenumbers

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%  imports  %%%%%%%%%%%%%%%%%%%
\usepackage{amsmath}
\usepackage{booktabs}
\usepackage{hyperref}
\usepackage{minted}
\usepackage{natbib}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\usepackage{glossaries}
\makeglossaries

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%  acronyms  %%%%%%%%%%%%%%%%%%%
\newacronym{vesical}{VESIcal}{Volatile Equilibria and Saturation Identification calculator}
\newacronym{etna}{Etna}{Etna, Italy}
\newacronym{csv}{CSV}{Comma Separated Values}
\newacronym{sfvf}{SFVF}{San Francisco Volcanic Field, AZ, USA}
\newacronym{tas}{TAS}{total alkalis vs silica}
\newacronym{mi}{MI}{melt inclusion}
\newacronym{mcmc}{MCMC}{Markov chain Monte Carlo}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%
% As of 2018 we recommend use of the TrackChanges package to mark revisions.
% The trackchanges package adds five new LaTeX commands:
%
%  \note[editor]{The note}
%  \annote[editor]{Text to annotate}{The note}
%  \add[editor]{Text to add}
%  \remove[editor]{Text to remove}
%  \change[editor]{Text to remove}{Text to add}
%
% complete documentation is here: http://trackchanges.sourceforge.net/
%%%%%%%
\draftfalse

%% Enter journal name below.
%% Choose from this list of Journals:
%
% JGR: Atmospheres
% JGR: Biogeosciences
% JGR: Earth Surface
% JGR: Oceans
% JGR: Planets
% JGR: Solid Earth
% JGR: Space Physics
% Global Biogeochemical Cycles
% Geophysical Research Letters
% Paleoceanography and Paleoclimatology
% Radio Science
% Reviews of Geophysics
% Tectonics
% Space Weather
% Water Resources Research
% Geochemistry, Geophysics, Geosystems
% Journal of Advances in Modeling Earth Systems (JAMES)
% Earth's Future
% Earth and Space Science
% Geohealth
%
% ie, \journalname{Water Resources Research}

\journalname{Notebooks Now!}

\begin{document}

%% ------------------------------------------------------------------------ %%
%  Title
%
% (A title should be specific, informative, and brief. Use
% abbreviations only if they are defined in the abstract. Titles that
% start with general keywords then specific terms are optimized in
% searches)
%
%% ------------------------------------------------------------------------ %%

% Example: \title{This is a test title}

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

%% ------------------------------------------------------------------------ %%
%
%  AUTHORS AND AFFILIATIONS
%
%% ------------------------------------------------------------------------ %%

% Authors are individuals who have significantly contributed to the
% research and preparation of the article. Group authors are allowed, if
% each author in the group is separately identified in an appendix.)

% List authors by first name or initial followed by last name and
% separated by commas. Use \affil{} to number affiliations, and
% \thanks{} for author notes.
% Additional author notes should be indicated with \thanks{} (for
% example, for current addresses).

% Example: \authors{A. B. Author\affil{1}\thanks{Current address, Antartica}, B. C. Author\affil{2,3}, and D. E.
% Author\affil{3,4}\thanks{Also funded by Monsanto.}}

\authors{Kayla Iacovino\affil{1}, Simon Matthews\affil{2,3}, Penny Wieser\affil{4}, Gordon M. Moore\affil{1}, Florence Bégué\affil{5}}

% \affiliation{1}{First Affiliation}
% \affiliation{2}{Second Affiliation}
% \affiliation{3}{Third Affiliation}
% \affiliation{4}{Fourth Affiliation}

\affiliation{1}{Jacobs, NASA Johnson Space Center, Houston, TX 77058, USA}
\affiliation{2}{Johns Hopkins University, Department of Earth and Planetary Sciences, Baltimore, MD 21218, USA}
\affiliation{3}{University of Iceland, Institute of Earth Sciences, Askja, Sturlugata 7, 101 Reykjavik, Iceland}
\affiliation{4}{University of Cambridge, Department of Earth Sciences, Downing Street, Cambridge CB2 3EQ, UK}
\affiliation{5}{University of Geneva, Department of Earth Sciences, Geneva, Switzerland}
%(repeat as many times as is necessary)

%% Corresponding Author:
% Corresponding author mailing address and e-mail address:

% (include name and email addresses of the corresponding author.  More
% than one corresponding author is allowed in this LaTeX file and for
% publication; but only one corresponding author is allowed in our
% editorial system.)

% Example: \correspondingauthor{First and Last Name}{email@address.edu}

\correspondingauthor{Kayla Iacovino}{kayla.iacovino@nasa.gov}

%% Keypoints, final entry on title page.

%  List up to three key points (at least one is required)
%  Key Points summarize the main points and conclusions of the article
%  Each must be 140 characters or fewer with no special characters or punctuation and must be complete sentences

% Example:
% \begin{keypoints}
% \item	List up to three key points (at least one is required)
% \item	Key Points summarize the main points and conclusions of the article
% \item	Each must be 140 characters or fewer with no special characters or punctuation and must be complete sentences
% \end{keypoints}

\begin{keypoints}
\item The first comprehensive volatile solubility tool capable of processing large datasets automatically
\item Seven built-in solubility models, with automatic calculation and plotting functionality
\item Built in python and easily usable by scientists with any level of coding skill
\end{keypoints}

%% ------------------------------------------------------------------------ %%
%
%  ABSTRACT and PLAIN LANGUAGE SUMMARY
%
% A good Abstract will begin with a short description of the problem
% being addressed, briefly describe the new data or analyses, then
% briefly states the main conclusion(s) and how they are supported and
% uncertainties.

% The Plain Language Summary should be written for a broad audience,
% including journalists and the science-interested public, that will not have
% a background in your field.
%
% A Plain Language Summary is required in GRL, JGR: Planets, JGR: Biogeosciences,
% JGR: Oceans, G-Cubed, Reviews of Geophysics, and JAMES.
% see http://sharingscience.agu.org/creating-plain-language-summary/)
%
%% ------------------------------------------------------------------------ %%

%% \begin{abstract} starts the second page

\begin{abstract}
Thermodynamics has been fundamental to the interpretation of geologic data and modeling of geologic systems for decades. However, more recent advancements in computational capabilities and a marked increase in researchers' accessibility to computing tools has outpaced the functionality and extensibility of currently available modeling tools. Here we present VESIcal (Volatile Equilibria and Saturation Identification calculator): the first comprehensive modeling tool for H\textsubscript{2}O, CO\textsubscript{2}, and mixed (H\textsubscript{2}O-CO\textsubscript{2}) solubility in silicate melts that: \textbf{a)} allows users access to seven of the most popular models, plus easy inter-comparison between models; \textbf{b)} provides universal functionality for all models (e.g., functions for calculating saturation pressures, degassing paths, etc.); \textbf{c)} can process large datasets (1,000's of samples) automatically; \textbf{d)} can output computed data into an Excel spreadsheet or CSV file for simple post-modeling analysis; \textbf{e)} integrates plotting capabilities directly within the tool; and \textbf{f)} provides all of these within the framework of a python library, making the tool extensible by the user and allowing any of the model functions to be incorporated into any other code capable of calling python. The tool is presented within this manuscript, which may be read as a static PDF but is better experienced via the Jupyter notebook version of this manuscript. We present here worked examples accessible to python users with a range of skill levels. The basic functions of VESIcal can also be accessed via a web app (\href{https://vesical.anvil.app}{https://vesical.anvil.app}). The VESIcal python library is open-source and available for download at \href{https://github.com/kaylai/VESIcal}{https://github.com/kaylai/VESIcal} and can be installed via pip.
\end{abstract}

\section*{Plain Language Summary}
Geologists use numerical models to understand and predict how volcanoes behave during storage (pre-eruption), eruption, and the composition and amount of volcanic gas released into the atmosphere of Earth and other planets. Most models are made by performing experiments on a limited dataset and creating a model that applies to that dataset. Some models combine lots of these individual models to make a generalized model that can apply to lots of different volcanoes. Many of these different models exist, and they all have specific uses, limitations, and pitfalls. Here we present the first tool, VESIcal, which acts as a simple interface to seven of the most commonly used models. VESIcal is written in python, so users can use VESIcal as an application or include it in their own models. VESIcal is the first tool that allows geologists to model thousands of data points automatically and provides a simple platform to compare results from different models in a way never before possible.

\subsection{Introduction}

Understanding the solubility and degassing of volatiles in silicate melts is a crucial component of modeling volcanic systems. As dissolved components, volatiles (primarily H\textsubscript{2}O and CO\textsubscript{2}) 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 \cite{tait1989, blake1984, stock2016}. 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\textsubscript{2}O, CO\textsubscript{2}, or mixed H\textsubscript{2}O-CO\textsubscript{2} 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 \cite{stolper1982, dixon1995, moore1998} 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\textsubscript{2}O-CO\textsubscript{2}) studies, increasing the natural applicability of our models to more systems \cite{liu2005, iacono2012, iacovino2013}. 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; \cite{ghiorso2015}) and the model of \cite{papale2006}. Both of these studies have made their source code available; the \cite{papale2006} FORTRAN source code (titled Solwcad), web app, and a Linux program can be found at \href{http://www.pi.ingv.it/progetti/eurovolc/}{http://www.pi.ingv.it/progetti/eurovolc/}, and very recently MagmaSat has been made accessible via the ENKI thermodynamic python framework (\href{http://enki-portal.org/}{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\textsubscript{2} and H\textsubscript{2}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\textsubscript{2}, H\textsubscript{2}O, and major element concentrations. For example, the saturation pressures from 105 Gakkel ridge melt inclusions calculated in MagmaSat by \cite{bennett2019} 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 \cite{dixon1995, dixon1997, liu2005, shishkina2014}. 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, \cite{newman2002} included a simplified version of the \cite{dixon1997} 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\textsubscript{2}, H\textsubscript{2}O, CO\textsubscript{2} content and temperature of every single sample into pop-up boxes. Similarly, the Excel spreadsheet for the \cite{moore1998} model calculates dissolved H\textsubscript{2}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, \cite{allison2019} 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., \cite{iacono2012}; \href{http://calcul-isto.cnrs-orleans.fr/}{http://calcul-isto.cnrs-orleans.fr/}, and \cite{ghiorso2015}; \href{http://melts.ofm-research.org/CORBA\_CTserver/GG-H2O-CO2.html}{http://melts.ofm-research.org/CORBA\_CTserver/GG-H2O-CO2.html}). \cite{ghiorso2015} 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 \cite{iacono2012} 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\textsubscript{2}O-CO\textsubscript{2} equation of state web calculator that \cite{duan2006} 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 \cite{wieser2021}, 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 \acrshort{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 \acrshort{vesical} are (also see Table~\ref{table-calibration-ranges}):

\begin{enumerate}
\item MagmaSat: \acrshort{vesical}'s default model. The mixed-volatile solubility model within MELTS v. 1.2.0 \cite{ghiorso2015}
\item Dixon: The simplification of the \cite{dixon1997} model as implemented in VolatileCalc \cite{newman2002}\begin{itemize}
\item DixonWater and DixonCarbon available as pure-fluid models
\end{itemize}


\item MooreWater: (\cite{moore1998}; water only, but H\textsubscript{2}O fluid concentration can be specified)
\item Liu: \cite{liu2005}\begin{itemize}
\item LiuCarbon and LiuWater available as pure-fluid models
\end{itemize}


\item IaconoMarziano: \cite{iacono2012}\begin{itemize}
\item IaconoMarzianoWater and IaconoMarzianoCarbon available as pure-fluid models
\end{itemize}


\item ShishkinaIdealMixing: \cite{shishkina2014} using pure-H\textsubscript{2}O and pure-CO\textsubscript{2} models and assuming ideal mixing. In general, the pure-fluid versions of this model should be used\begin{itemize}
\item ShishkinaWater and ShishkinaCarbon available as pure-fluid models
\end{itemize}


\item AllisonCarbon: \cite{allison2019}\begin{itemize}
\item AllisonCarbon\_vesuvius (default; phonotephrite from Vesuvius, Italy)
\item AllisonCarbon\_sunset (alkali basalt from Sunset Crater, AZ, USA)
\item AllisonCarbon\_sfvf (basaltic andesite from San Francisco Volcanic Field, AZ, USA)
\item AllisonCarbon\_erebus (phonotephrite from Erebus, Antarctica)
\item AllisonCarbon\_etna (trachybasalt from \acrshort{etna}, Italy)
\item AllisonCarbon\_stromboli (alkali basalt from Stromboli, Italy)
\end{itemize}
\end{enumerate}

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, \acrshort{vesical} is a single tool that can access and utilize all built-in models, giving \acrshort{vesical} an extensive pressure-temperature-composition calibration range \cite{Figure~\ref{fig1}}. \acrshort{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 \acrshort{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 (\acrshort{vesical} code can be imported for use in python scripts, and the code is formatted such that new volatile models can be added).

Importantly, \acrshort{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 \acrshort{vesical} through a webapp (\href{https://vesical.anvil.app}{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 \acrshort{csv} file to work with outside of \acrshort{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 \acrshort{vesical}, including the ability to hybridize models (e.g., use one model for H\textsubscript{2}O and another for CO\textsubscript{2}) or write their own routines and code calling \acrshort{vesical} methods. \acrshort{vesical} is an open source tool and as such is far less prone to the preservation issues discussed above. Because the \acrshort{vesical} code is hosted on GitHub, every change to the code is tracked publicly \cite{perkel2016}. \acrshort{vesical}'s current release (version 1.0.1) is also archived on Zenodo, which provides a static citable DOI \cite{https://doi.org/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 \acrshort{vesical} are explored in detail in the companion manuscript to this work, \cite{wieser2021}.

\subsection{Research Methodology}

Navigating the array of models implemented in \acrshort{vesical} can be challenging. How can a user determine which model best suits their needs? MagmaSat (the default model in \acrshort{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~\ref{comparing-user-data-to-model-calibrations} on comparing user data to model calibrations.

\begin{table}
\centering
\caption[]{Calibration ranges of \acrshort{vesical} models}
\label{table-calibration-ranges}
\begin{tabular}{p{\dimexpr 0.167\linewidth-2\tabcolsep}p{\dimexpr 0.167\linewidth-2\tabcolsep}p{\dimexpr 0.167\linewidth-2\tabcolsep}p{\dimexpr 0.167\linewidth-2\tabcolsep}p{\dimexpr 0.167\linewidth-2\tabcolsep}p{\dimexpr 0.167\linewidth-2\tabcolsep}}
\toprule
Model/Reference & Species & P (bar) & T (°C) & Compositional range & Notes \\
\hline
\textbf{MagmaSat}\newline
 \cite{ghiorso2015} & H\textsubscript{2}O & 0--20,000\textsuperscript{1} & 550-1420\textsuperscript{1} & Very broad compositional range of natural silicate melts: subalkaline picrobasalts to rhyolites, including a variety of mafic and silicic alkaline compositions & \textsuperscript{1}Ranges extracted from Fig. 2d of \cite{ghiorso2015} \\
CO\textsubscript{2} & 0--30,000\textsuperscript{1} & 1139-1400\textsuperscript{1} \\
H\textsubscript{2}O - CO\textsubscript{2} & 0--10,000\textsuperscript{1} & 800-1400\textsuperscript{1} \\
\textbf{Dixon}\newline
 Simplification of \cite{dixon1997} used in VolatileCalc (\cite{newman2002}) & H\textsubscript{2}O - CO\textsubscript{2} & 0-5000\textsuperscript{1}\newline
0-2000\textsuperscript{2}\newline
0-1000\textsuperscript{3} & 600-1500\textsuperscript{1}\newline
(1200)\textsuperscript{4} & Alkali basalts:\newline
40-49 wt\% SiO2 & \textsuperscript{1}Warnings implemented in VolatileCalc (\cite{newman2002}).\newline
 \textsuperscript{2}Calibration range suggested by \cite{lesne2011}\newline
 \textsuperscript{3}Calibration range suggested by \cite{iacono2012}\newline
 \textsuperscript{4}Calibration temperature of \cite{dixon1997}\newline
 \\
\textbf{MooreWater}\newline
\cite{moore1998} & H\textsubscript{2}O & 0--3000\textsuperscript{1} & 700--1200\textsuperscript{1} & Broad compositional range: subalkaline basalts to rhyolites, alkaline trachybasalts-andesites, foidites, phonolites & \textsuperscript{1}Author-suggested calibration range. The calibration dataset spans 190 to 6067 bar, and 800-1200°C \\
\textbf{Liu}\newline
\cite{liu2005} & H\textsubscript{2}O - CO\textsubscript{2} & 0-5000\textsuperscript{1} & 700--1200\textsuperscript{1} & Haplogranites and rhyolites & \textsuperscript{1}Author-suggested calibration range for the mixed fluid model. The calibration dataset covers 750-5510 bar and 800-1150°C for the Carbon model, and 1-5000 bar and 700-1200 °C for the water model \\
\textbf{Iacono-Marziano}\newline
\cite{iacono2012} & H\textsubscript{2}O - CO\textsubscript{2} & 95--10,500 (mostly \textless 5000)\textsuperscript{1} & 1100-1400 (preferably 1200- 1300)\textsuperscript{2} & Predominantly mafic compositions: subalkaline and alkaline basalts-andesites & \textsuperscript{1}Range of calibration dataset, as authors do not specifically state a calibration range. We note that the vast majority of experiments were conducted at \textless 5000 bar.

\textsuperscript{2}Authors state that most experiments were conducted between 1200-1300°C (whole range 1100-1400°C

 \\
\textbf{Shishkina}\newline
\cite{shishkina2014} & H\textsubscript{2}O\textsuperscript{1} & 0--5000\textsuperscript{2} & 1050--1400 (preferably 1150- 1250)\textsuperscript{2,~3} & Mafic and intermediate compositions: Subalkaline basalts- basaltic andesites, alkali basanites-phonolites. SiO2\textless 65 wt\%. & \textsuperscript{1}Although their empirical expressions are for pure fluids, they were mostly calibrated on mixed CO\textsubscript{2}-H\textsubscript{2}O experiments.

\textsuperscript{2}Author-suggested range

\textsuperscript{3}Note, this model contains no temperature term.

 \\
CO\textsubscript{2}\textsuperscript{1} & 500--5000\textsuperscript{2} & 1200--1250\textsuperscript{2,~3} & Predominantly mafic compositions: subalkaline basalts, alkaline basanites, trachybasalts \\
\textbf{AllisonCarbon}\newline
\cite{allison2019} & CO\textsubscript{2} & 0-7000\textsuperscript{1} & 1200\textsuperscript{2}\newline
({\textasciitilde}1000-1400) & Alkali-rich mafic magmas from six volcanic fields. Separate model coefficients for each composition. & \textsuperscript{1}Author-suggested range. The calibration dataset spans: (\acrshort{sfvf}: 4133--6141 bar, Sunset Crater: 4071-6098 bar, Erebus: 4078-6175 bar, Vesuvius: 269-6175 bar, \acrshort{etna}=485-6199, Stromboli=524-6080)

\textsuperscript{2}Note, all calculations performed at 1200°C (the experimental temperature). Authors suggest results generally applicable between 1000-1400°C

 \\
\bottomrule
\end{tabular}
\end{table}

\begin{figure}[!htbp]
\centering
\includegraphics[width=0.7\linewidth]{files/Figure1-906bf392f0d1b3aaeeaa0686319646f1.png}
\caption[]{Illustrations showing the calibrated ranges of \acrshort{vesical} models in pressure-temperature space. Due to difficulty in differentiating between pure-CO\textsubscript{2} and mixed fluid experiments in the literature, plots are subdivided into: experiments performed with pure-CO\textsubscript{2} or mixed (H\textsubscript{2}O-CO\textsubscript{2}) fluid; and pure-H\textsubscript{2}O fluid.}
\label{fig1}
\end{figure}

A list of model names recognized by \acrshort{vesical} can be retrieved by executing the command \texttt{v.get\_model\_names()}, assuming \acrshort{vesical} has been imported as \texttt{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 \acrshort{vesical} (e.g., \texttt{model='MooreWater'}). \cite{allison2019} provides unique model equations for each of the six alkali-rich mafic magmas investigated in their study. The default model in \acrshort{vesical} is that calibrated for Vesuvius magmas, whose calibration has the widest pressure range of the study Table~\ref{table-calibration-ranges}. Setting a model name of \texttt{'AllisonCarbon'} within \acrshort{vesical} will thus result in calculations using the AllisonCarbon\_vesuvius model equations.

All of the calculations implemented in \acrshort{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 \acrshort{csv} files. To demonstrate that \acrshort{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: \texttt{model\_name}, \texttt{sample} (both provided by the user), \texttt{model} (an instance of the Model class used to run the calculations of interest), \texttt{result} (the result of the calculations), and \texttt{calib\_check} (the results of the calibration check).

\subsubsection{Model Calibrations and Benchmarking}

The pressure, temperature, and composition calibration ranges of the seven models implemented in \acrshort{vesical} are shown in Table~\ref{table-calibration-ranges} and Figure~\ref{fig1}. \acrshort{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 \acrshort{vesical}, including information about their calibration range in terms of melt composition, pressure, and temperature, are given in \cite{wieser2021}

Testing was undertaken to ensure that \acrshort{vesical} faithfully reproduces the results of all incorporated models. When possible, all models were benchmarked by testing \acrshort{vesical} outputs against those of a relevant published calculator (e.g., web apps or Excel macros). The models of \cite{shishkina2014} and \cite{liu2005} were published with no such tool and so testing instead compares \acrshort{vesical} outputs to experimental conditions or analyses and, where possible, plots \acrshort{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, \acrshort{vesical} reproduced the results from previous tools (e.g., web apps, Excel spreadsheets) to within $\pm$1\% relative and often on the order of $\pm$0.1\% relative.

MagmaSat, \acrshort{vesical}'s default model, underwent three tests, the results of which are shown in Figure~\ref{fig2}: 1. Comparison of saturation pressures from MORB melt inclusions in \acrshort{vesical} to those published by \cite{bennett2019}, who used the MagmaSat Mac App (R\textsuperscript{2}=0.99998; Figure~\ref{fig2}); 2. Comparison of fluid composition (X$_{H_{2}O}$) calculated with \acrshort{vesical} and the web app (R\textsuperscript{2}=0.999, identical considering the web app returns 2dp; Figure~\ref{fig2}); 3. Comparison of isobars for the Early Bishop Tuff calculated with \acrshort{vesical} (star symbols) and isobars published in Figure 14 of \cite{ghiorso2015} (Figure~\ref{fig2}). \acrshort{vesical} outputs using the model of \cite{dixon1997} were tested against outputs from the VolatileCalc Excel spreadsheet \cite{newman2002} and a widely used Excel macro \cite{tucker2019}.

\begin{figure}[!htbp]
\centering
\includegraphics[width=0.7\linewidth]{files/Figure2-19c54e526b8b3d94d21028935441720c.png}
\caption[]{Benchmarking of \acrshort{vesical} against MagmaSat. a. Comparison of saturation pressures calculated with \acrshort{vesical} against those by \cite{bennett2019} using the MagmaSat app for Mac. Samples are all MORB melt inclusions, and pressures were calculated at a temperature unique to each sample. b. Equilibrium fluid compositions calculated with \acrshort{vesical} against those calculated with the MagmaSat web app. c. Individual points along the 1,000, 2,000, and 3,000 bar isobars for the Early Bishop Tuff rhyolite calculated with \acrshort{vesical} (stars) and plotted atop isobars published in Figure 14 of \cite{ghiorso2015}.}
\label{fig2}
\end{figure}

\subsubsection{Format of the python library}

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

Specific details on how to perform model calculations are discussed in Section~\ref{example-uses} and include worked examples. The \acrshort{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:

\begin{enumerate}
\item \texttt{calculate\_dissolved\_volatiles()}
\item \texttt{calculate\_equilibrium\_fluid\_composition()}
\item \texttt{calculate\_saturation\_pressure()}
\item \texttt{calculate\_isobars\_and\_isopleths()} (plus functionality for plotting; only for mixed volatiles models)
\item \texttt{calculate\_degassing\_path()} (plus functionality for plotting; only for mixed volatiles models).
\end{enumerate}

Figure~\ref{fig3} 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\textsubscript{2}O\textsuperscript{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:

\begin{minted}[breaklines]{python}
saturation_pressure_calculation = calculate_saturation_pressure(sample=mysample, temperature=850.0)
\end{minted}

where \texttt{mysample} is a variable (\acrshort{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:

\begin{minted}[breaklines]{python}
saturation_pressure_calculation.result
saturation_pressure_calculation.calib_check
\end{minted}

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:

\begin{minted}[breaklines]{python}
saturation_pressure = calculate_saturation_pressure(sample=mysample, temperature=850.0).result
\end{minted}

If a different model is desired, for example \cite{dixon1997}, it can be passed as:

\begin{minted}[breaklines]{python}
calculate_saturation_pressure(sample=mysample, temperature=850.0, model='Dixon').result
\end{minted}

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 \texttt{calculate()} method sends the user input to the appropriate model.

\begin{figure}[!htbp]
\centering
\includegraphics[width=0.7\linewidth]{files/Figure3-0c92d500cb6c6ce8d65f8ff749429745.png}
\caption[]{Flowchart illustrating the basic organization of the python library. First, a user chooses a calculation to perform and calls one of the five core calculation classes. Here, any necessary parameters are passed such as sample composition, pressure, and temperature. A check is run to ensure the calculation is being performed within model-specified limits. The \texttt{Calculate()} class then calls on one of the \texttt{Model()} classes. The default model is MagmaSat, but a user may specify a different model when defining the calculation parameters. Standard pre-processing is then performed on the input data, and this pre-processing step is unique to each model. The processed data are then fed into a model-specific method to perform the desired core calculation.}
\label{fig3}
\end{figure}

Users can process individual samples (single-sample calculations) or entire datasets (batch calculations; Figure~\ref{fig4}). If processing more than one sample, the ``simplest'' way to interact with \acrshort{vesical} is via batch calculations. Here, the user provides input data in the form of a Microsoft Excel spreadsheet (.xlsx file) or \acrshort{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 \acrshort{csv} file with a simple command (see Section~\ref{exporting-data}).

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~\ref{code-hybridization}). This also gives the user more flexibility in integrating any \acrshort{vesical} model function into some other python code.

\begin{figure}[!htbp]
\centering
\includegraphics[width=0.7\linewidth]{files/Figure4-65ff494af9e775fa245a775022e6c236.png}
\caption[]{Flowchart illustrating the different operational paths. On top, batch calculation is shown, in which an Excel or \acrshort{csv} file with any amount of samples is fed into the model, calculations are performed, and the original user data plus newly calculated values are returned and can be saved as an Excel or \acrshort{csv} file. Below, single-sample calculation is shown. These methods can run calculations on one sample at a time, but multi-sample calculations can be performed iteratively with code written by the user. Calculated values are returned as a variable. For single-sample calculations, more advanced modeling options can be set, and hybridization of models can be performed.}
\label{fig4}
\end{figure}

\subsubsection{Running the code}

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

\acrshort{vesical} was born from functionality provided by ENKI and so all the files necessary to use \acrshort{vesical} are hosted on the ENKI server (\href{http://enki-portal.org/}{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 \acrshort{vesical} while retaining all of its functionality is within this very manuscript, in the form of a Jupyter notebook. Because this manuscript and \acrshort{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 \acrshort{vesical} by creating a Jupyter notebook on the ENKI server and importing \acrshort{vesical} as is demonstrated in the code below.

Computation time on the ENKI server is limited by the server itself. \acrshort{vesical} may run faster if installed locally. Advanced instructions on installing \acrshort{vesical} on your own computer are provided in the Supplement (Supplementary Text S1). Note that \acrshort{vesical} requires installation of the ENKI thermoengine library to function properly. Thermoengine is written in python but is based on the original MELTS code \cite{ghiorso1995, ghiorso2015}, 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; \href{https://gitlab.com/ENKI-portal/ThermoEngine}{https://gitlab.com/ENKI-portal/ThermoEngine}).

The most limited but simplest method to interacting with \acrshort{vesical} is through the web app (\href{https://vesical.anvil.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 \acrshort{csv} file). Some, but not all, optional parameters can be set.

\textbf{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 \acrshort{csv} file into the same folder containing this notebook and then changing the filename in Section~\ref{loading-viewing-preparing}.

\subsubsection{Documentation}

This manuscript serves as an introduction to the \acrshort{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 (\href{https://vesical.readthedocs.io/}{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., ``\texttt{v.calculate\_saturation\_pressure?}'').

Video tutorials are also available on the \acrshort{vesical} YouTube (\href{https://www.youtube.com/channel/UCpvCCs5KMXzOxXWm0seF8Qw}{https://www.youtube.com/channel/UCpvCCs5KMXzOxXWm0seF8Qw}). The first tutorial covers the basics of \acrshort{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.

\subsubsection{Generic methods for calculating mixed-fluid properties}

\acrshort{vesical} provides a set of methods for calculating the properties of mixed CO\textsubscript{2}-H\textsubscript{2}O fluids, which can be used with any combination of H\textsubscript{2}O and CO\textsubscript{2} solubility model. The use of generic methods allows additional models to be added to \acrshort{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).

\subsection{Workable example uses}\label{example-uses}

In this section we detail how to use the various functions available in \acrshort{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 Section~\ref{initialize-packages} and Section~\ref{comparing-user-data-to-model-calibrations} 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 \texttt{plot()} function \cite{Section~\ref{plotting}} 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):

\begin{itemize}
\item Dataset from \href{/example\_data-9296285330fa856e8b827784ac193203.xlsx}{\texttt{example\_data.xlsx}}  loaded in Section~\ref{batch-processing} (variable name \texttt{myfile})
\item Single composition defined in Section~\ref{defining-single-sample} (variable name \texttt{mysample})
\item Sample 10* extracted from \href{/example\_data-9296285330fa856e8b827784ac193203.xlsx}{\texttt{example\_data.xlsx}} dataset in Section~\ref{extracting-single-melt-composition} (variable name \texttt{sample\_10})
\end{itemize}

Calculations performed on single samples or on a dataset imported from an Excel or \acrshort{csv} file containing many samples are executed in two distinct ways. Note that single sample calculations require that the argument \texttt{sample} be defined. To return the numerical result of the calculation, the \texttt{.result} method must be called, as shown below. Batch calculations are performed on the dataset itself, after that dataset is imported into \acrshort{vesical}. Thus, the \texttt{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:

\begin{description}
\item[\textit{Single sample calculations}] \texttt{myvariable = v.name\_of\_the\_core\_calculation(sample=mysample, argument1=value1, argument2=value2).result}
\item[\textit{Batch calculations}] \texttt{myvariable = myfile.name\_of\_the\_core\_calculation(argument1=value1, argument2=value2)}
\end{description}

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

Workable examples detailed here are:

\begin{enumerate}
\item Section~\ref{loading-viewing-preparing}\begin{itemize}
\item Section~\ref{batch-processing}
\item Section~\ref{defining-single-sample}
\item Section~\ref{extracting-single-melt-composition}
\item Section~\ref{normalizing-transforming-data}
\end{itemize}


\item Section~\ref{calculating-dissolved-volatile-concentrations}
\item Section~\ref{calculating-equilibrium-fluid-compositions}
\item Section~\ref{calculating-saturation-pressures}
\item Section~\ref{calculating-plotting-isobars-isopleths}
\item Section~\ref{calculating-plotting-degassing-paths}
\item Section~\ref{plotting-multiple-calculations}
\item Section~\ref{comparing-results-from-multiple-models}
\item Section~\ref{code-hybridization}
\item Section~\ref{exporting-data}
\end{enumerate}

\subsubsection{Calculation class arguments and their definitions}\label{class-arguments-and-definitions}

Each section below details what arguments are required or optional inputs and gives examples of how to perform the calculations. Table~\ref{table2} 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:

\begin{description}
\item[\texttt{sample}] \textit{Single sample calculations only} The composition of a sample. A \acrshort{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.

\item[\texttt{temperature}, \texttt{pressure}, and \texttt{X\_fluid}] the temperature in °C, the pressure in bars, and the mole fraction of H\textsubscript{2}O in the H\textsubscript{2}O-CO\textsubscript{2} fluid, XH\textsubscript{2}O\textsuperscript{fluid}. In all cases, \texttt{X\_fluid} is optional, with a default value of 1 (pure H\textsubscript{2}O fluid). Note that that \texttt{X\_fluid} argument is only used for calculation of dissolved volatile concentrations.

\begin{description}
\item[\textit{For single sample calculations}] Temperature, pressure, and X fluid should be specified as a numerical value.
\item[\textit{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 \texttt{temperature},  \texttt{pressure}, and/or \texttt{X\_fluid}, respectively. Note for batch calculations that if temperature, pressure, or XH\textsubscript{2}O\textsuperscript{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.
\end{description}

\item[\texttt{verbose}] \textit{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.

\item[\texttt{print\_status}] \textit{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.

\item[\texttt{model}] Always an optional argument referring to the name of the desired solubility model to use. The default is always ``MagmaSat''.

\end{description}

\begin{table}
\centering
\caption[]{Matrix of all arguments used in the five core calculations, the nature of the argument (required or optional) and the input type or default value.

SS = Single-sample. Batch = batch processing. Color of cells corresponds to the type of argument: green=required; orange=optional; gray=argument not used.
Values in cells indicate the unit or type of data to input for required arguments or the default value in the case of optional arguments.}
\label{table2}
\begin{tabular}{p{\dimexpr 0.111\linewidth-2\tabcolsep}p{\dimexpr 0.111\linewidth-2\tabcolsep}p{\dimexpr 0.111\linewidth-2\tabcolsep}p{\dimexpr 0.111\linewidth-2\tabcolsep}p{\dimexpr 0.111\linewidth-2\tabcolsep}p{\dimexpr 0.111\linewidth-2\tabcolsep}p{\dimexpr 0.111\linewidth-2\tabcolsep}p{\dimexpr 0.111\linewidth-2\tabcolsep}p{\dimexpr 0.111\linewidth-2\tabcolsep}}
\toprule
 & \multicolumn{2}{p{\dimexpr 0.222\linewidth-2\tabcolsep}}{dissolved\_volatiles} & \multicolumn{2}{p{\dimexpr 0.222\linewidth-2\tabcolsep}}{equilibrium\_fluid\_comp} & \multicolumn{2}{p{\dimexpr 0.222\linewidth-2\tabcolsep}}{saturation\_pressure} & isobars\_isopleths & degassing\_path \\
\hline
 & SS & Batch & SS & Batch & SS & Batch & SS & SS \\
\hline
sample & wt\% oxides &  & wt\% oxides &  & wt\% oxides &  & wt\% oxides & wt\% oxides \\
\hline
temperature & °C & °C & °C & °C & °C & °C & °C & °C \\
\hline
pressure & bars & bars & bars & None &  &  &  & 'saturation' \\
\hline
pressure\_list &  &  &  &  &  &  & bars &  \\
\hline
X\_fluid & 1 & 1 &  &  &  &  &  &  \\
\hline
isopleth\_list &  &  &  &  &  &  & None &  \\
\hline
verbose & False &  & False &  & False &  &  &  \\
\hline
model & 'MagmaSat' & 'MagmaSat' & 'MagmaSat' & 'MagmaSat' & 'MagmaSat' & 'MagmaSat' & 'MagmaSat' & 'MagmaSat' \\
\hline
print\_status &  & True &  & False &  & True & True &  \\
\hline
smooth\_isobars &  &  &  &  &  &  & True &  \\
\hline
smooth\_isopleths &  &  &  &  &  &  & True &  \\
\hline
fractionate\_vapor &  &  &  &  &  &  &  & 0.0 \\
\hline
init\_vapor &  &  &  &  &  &  &  & 0.0 \\
\hline
\bottomrule
\end{tabular}
\end{table}

\subsubsection{Initialize packages}\label{initialize-packages}

For any code using the \acrshort{vesical} library, the library must be imported for use. Here we import \acrshort{vesical} as v. Any time we wish to initialize a \acrshort{vesical} object, that class name must be preceded by `\texttt{v.}' (e.g., \texttt{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.

\begin{minted}[breaklines]{python}
import VESIcal as v
import pandas as pd

#The following are options for formatting this manuscript
pd.set_option('display.max_colwidth', 0)
from IPython.display import display, HTML
%matplotlib inline
import warnings
warnings.filterwarnings("ignore", category=DeprecationWarning)
\end{minted}

\subsubsection{Loading, viewing, and preparing user data}\label{loading-viewing-preparing}

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 \href{/example\_data-9296285330fa856e8b827784ac193203.xlsx}{\texttt{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.

\paragraph{Batch processing}\label{batch-processing}

Batch calculations are always facilitated via the \texttt{BatchFile()} class, which the user uses to specify the filename corresponding to sample data. Loading in data is as simple as calling \texttt{BatchFile(filename)}. Optionally, \texttt{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.

\textbf{Structure of the input file:} A file containing compositions (and optional pressure, temperature, or XH\textsubscript{2}O\textsuperscript{fluid} information) on one or multiple samples can be loaded into \acrshort{vesical}. The loaded file must be a Microsoft Excel file with the extension .xls or .xlsx or \acrshort{csv} file with the extension .csv. The file must be laid out in the same manner as the example file \href{/example\_data-9296285330fa856e8b827784ac193203.xlsx}{\texttt{example\_data.xlsx}}. The basic structure is also shown in Table~\ref{table3}.

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 \acrshort{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~\ref{extracting-single-melt-composition}). An index column can be specified by the user using the argument \texttt{label} (see documentation below). The following columns must contain compositional information as oxides. The only allowable oxides are: SiO\textsubscript{2}, TiO\textsubscript{2}, Al\textsubscript{2}O\textsubscript{3}, Fe\textsubscript{2}O\textsubscript{3}, FeO, Cr\textsubscript{2}O\textsubscript{3}, MnO, MgO, CaO, NiO, CoO, Na\textsubscript{2}O, K\textsubscript{2}O, P\textsubscript{2}O\textsubscript{5}, H\textsubscript{2}O, and CO\textsubscript{2}. Currently, \acrshort{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 \acrshort{vesical} assumes a particular formatting of column names, we highly recommend that users examine their data after loading into \acrshort{vesical} and before performing calculations. The user data, as it will be used by \acrshort{vesical}, can be viewed at any time with \texttt{myfile.get\_data()} (see generation of Table~\ref{table3} below).

Pressure, temperature, or XH\textsubscript{2}O\textsuperscript{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 \acrshort{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\textsubscript{2}O in an H\textsubscript{2}O-CO\textsubscript{2} fluid, ranging from 0-1). Sample compositions may be translated between wt\%, mol fraction, and mol cations if necessary.

\begin{description}
\item[\textbf{Class structure}] \texttt{BatchFile(filename, sheet\_name=0, file\_type='excel', units='wtpt\_oxides', label='Label', default\_normalization='none', default\_units='wtpt\_oxides', dataframe=None)}

\item[\textbf{Required inputs}] \texttt{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.

\item[\textbf{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.] \texttt{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.

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

\texttt{units}: The units in which data are input. The default value is \texttt{'wtpt\_oxides'} for data as wt\% oxides. The user can pass \texttt{'mol\_oxides'} for data in mol fraction oxides or \texttt{'mol\_cations'} for data in mol fraction cations.

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

\texttt{default\_units}: The type of composition to return by default, one of: \texttt{'wtpt\_oxides'} (wt\% oxides, default), \texttt{'mol\_oxides'} (mol fraction oxides), or \texttt{'mol\_cations'} (mol fraction cations).

\texttt{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 \acrshort{vesical}'s recognized oxides will be set as the index column. The index column will be used to select samples by name.

\texttt{dataframe}: This argument is used for transforming a pandas DataFrame object into a \acrshort{vesical} BatchFile object. For convenience, this functionality is also defined as a separate function \texttt{BatchFile\_from\_DataFrame(dataframe, units='wtpt\_oxides', label='Label')}.

\item[\textbf{Outputs}] A special type of python object defined in the \acrshort{vesical} code known as an BatchFile object.

\end{description}

\begin{minted}[breaklines]{python}
myfile = v.BatchFile('Supplement/Datasets/example_data.xlsx')
\end{minted}

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

Usage of \texttt{get\_data()} allows the user to retrieve the data as originally entered or in any units and with any normalization supported by \acrshort{vesical}.

\begin{description}
\item[\textbf{Method structure}] \texttt{get\_data(self, normalization=None, units=None, asBatchFile=False)}

\item[\textbf{Optional inputs}] \texttt{normalization} or \texttt{units} may be passed, with options as defined in the description of BatchFile above.

\texttt{asBatchFile} Default is False. If True, will return a \acrshort{vesical} BatchFile object.

\item[\textbf{Outputs}] A pandas dataframe or BatchFile object with all user data.

\end{description}

\begin{code}[H]
\centering
\caption[]{User input data: Compositions, pressures, and temperatures for several silicate melts as supplied in the file \href{/example\_data-9296285330fa856e8b827784ac193203.xlsx}{\texttt{example\_data.xlsx}}}
\label{table3}\begin{minted}[breaklines]{python}
data = myfile.get_data()
data
\end{minted}

\begin{minted}[breaklines]{text}
SiO2      TiO2      Al2O3  Fe2O3  Cr2O3  \
Kil3 -6_1a                 48.249207  2.222114  11.692194  0.00   0.0     
Kil3 -6_3a                 48.295691  2.165357  11.755584  0.00   0.0     
Kil3 -6_4a                 49.124079  2.360984  12.172833  0.00   0.0     
10*                       47.960000  0.780000  18.770000  0.00   0.0     
19*                       49.640000  0.710000  18.050000  0.00   0.0     
25                        50.320000  0.720000  18.030000  0.00   0.0     
SAT -M12 -1                 62.600000  0.630000  17.300000  2.01   0.0     
SAT -M12 -2                 62.600000  0.630000  17.300000  2.01   0.0     
SAT -M12 -4                 62.600000  0.630000  17.300000  2.01   0.0     
samp. P1968a              76.974880  0.085516  3.110636   0.00   0.0     
samp. P1968b              76.943845  0.133125  3.169657   0.00   0.0     
samp. P1968c              77.187205  0.119506  3.167827   0.00   0.0     
samp. HPR3 -1_XL -3         75.413966  0.095164  14.077692  0.00   0.0     
samp. HPR3 -1_XL -4_INCL -1  76.613586  0.095843  13.476762  0.00   0.0     
AW -6                      48.030000  2.840000  18.120000  0.00   0.0     
AW -46                     52.980000  2.180000  20.490000  0.00   0.0     
KI -07                     44.610000  4.370000  14.410000  0.00   0.0     

                                FeO       MnO        MgO  NiO  CoO        CaO  \
Kil3 -6_1a                 0.000000   0.079999  14.183817  0.0  0.0  9.892732    
Kil3 -6_3a                 0.000000   0.084045  13.403980  0.0  0.0  10.052578   
Kil3 -6_4a                 0.000000   0.098809  11.997699  0.0  0.0  10.308188   
10*                       10.920000  0.150000  6.860000   0.0  0.0  12.230000   
19*                       10.540000  0.190000  6.430000   0.0  0.0  12.090000   
25                        10.110000  0.140000  5.650000   0.0  0.0  12.780000   
SAT -M12 -1                 2.010000   0.060000  2.650000   0.0  0.0  5.640000    
SAT -M12 -2                 2.010000   0.060000  2.650000   0.0  0.0  5.640000    
SAT -M12 -4                 2.010000   0.060000  2.650000   0.0  0.0  5.640000    
samp. P1968a              4.788883   0.000000  12.549439  0.0  0.0  1.207910    
samp. P1968b              4.763435   0.000000  12.446403  0.0  0.0  1.231728    
samp. P1968c              4.814076   0.000000  12.229534  0.0  0.0  1.184773    
samp. HPR3 -1_XL -3         0.654992   0.125882  0.012003   0.0  0.0  0.636124    
samp. HPR3 -1_XL -4_INCL -1  0.620769   0.113495  0.032069   0.0  0.0  0.624350    
AW -6                      9.600000   0.230000  3.080000   0.0  0.0  7.570000    
AW -46                     5.540000   0.200000  2.000000   0.0  0.0  7.100000    
KI -07                     10.600000  0.170000  7.690000   0.0  0.0  11.550000   

                              Na2O       K2O      P2O5       H2O       CO2  \
Kil3 -6_1a                 1.810522  0.352014  0.210479  0.424695  0.002873   
Kil3 -6_3a                 2.268198  0.373328  0.204452  0.425984  0.006786   
Kil3 -6_4a                 2.001863  0.396512  0.238996  0.437758  0.004984   
10*                       1.950000  0.210000  0.170000  4.500000  0.047900   
19*                       1.990000  0.200000  0.170000  5.100000  0.111300   
25                        1.800000  0.240000  0.230000  5.200000  0.043700   
SAT -M12 -1                 4.050000  1.610000  0.240000  2.620000  0.000000   
SAT -M12 -2                 4.050000  1.610000  0.240000  5.030000  0.000000   
SAT -M12 -4                 4.050000  1.610000  0.240000  6.760000  0.000000   
samp. P1968a              0.138963  1.133084  0.000000  4.340000  0.007000   
samp. P1968b              0.140993  1.170806  0.000000  5.850000  0.012300   
samp. P1968c              0.138201  1.158924  0.000000  5.754571  0.010663   
samp. HPR3 -1_XL -3         3.703110  5.128392  0.000000  5.943750  0.010000   
samp. HPR3 -1_XL -4_INCL -1  3.676972  4.579799  0.000000  5.340000  0.008000   
AW -6                      6.040000  3.080000  1.410000  1.420000  0.129800   
AW -46                     5.680000  3.160000  0.660000  4.760000  0.343900   
KI -07                     3.930000  1.740000  0.920000  2.900000  0.113100   

                                        CITATION                ROCK TYPE  \
Kil3 -6_1a                 Tucker et al. (2019)    Basalt                    
Kil3 -6_3a                 Tucker et al. (2019)    Basalt                    
Kil3 -6_4a                 Tucker et al. (2019)    Basalt                    
10*                       Roggensack (2001)       Basalt                    
19*                       Roggensack (2001)       Basalt                    
25                        Roggensack (2001)       Basalt                    
SAT -M12 -1                 Moore et al. (1998)     Andesite                  
SAT -M12 -2                 Moore et al. (1998)     Andesite                  
SAT -M12 -4                 Moore et al. (1998)     Andesite                  
samp. P1968a              Myers et al. (2019)     Rhyolite                  
samp. P1968b              Myers et al. (2019)     Rhyolite                  
samp. P1968c              Myers et al. (2019)     Rhyolite                  
samp. HPR3 -1_XL -3         Mercer et al. (2015)    Rhyolite                  
samp. HPR3 -1_XL -4_INCL -1  Mercer et al. (2015)    Rhyolite                  
AW -6                      Iacovino et al. (2016)  Phonotephrite             
AW -46                     Iacovino et al. (2016)  Basaltic -Trachyandesite   
KI -07                     Iacovino et al. (2016)  Basanite                  

                           Press         Temp  
Kil3 -6_1a                 62.5    1299.094712  
Kil3 -6_3a                 128.0   1283.419991  
Kil3 -6_4a                 100.0   1255.153759  
10*                       2000.0  1200.000000  
19*                       2000.0  1200.000000  
25                        2000.0  1200.000000  
SAT -M12 -1                 703.0   1100.000000  
SAT -M12 -2                 1865.0  1100.000000  
SAT -M12 -4                 2985.0  1050.000000  
samp. P1968a              300.0   900.000000   
samp. P1968b              300.0   900.000000   
samp. P1968c              300.0   900.000000   
samp. HPR3 -1_XL -3         300.0   0.000000     
samp. HPR3 -1_XL -4_INCL -1  0.0     900.000000   
AW -6                      1500.0  1050.000000  
AW -46                     4000.0  1000.000000  
KI -07                     1500.0  1100.000000
\end{minted}
\end{code}

For the rest of this manuscript, data will be pulled from the \href{/example\_data-9296285330fa856e8b827784ac193203.xlsx}{\texttt{example\_data.xlsx}} file (Supplemental Dataset S1), which contains compositional information for basalts \cite{tucker2019, roggensack2001}, andesites \cite{moore1998}, rhyolites \cite{mercer2015, myers2019}, and alkaline melts (phototephrite, basaltic-trachyandesite, and basanite from \cite{iacovino2016}). Several additional example datasets from the literature are available in the Supplement (Supplementary Datasets S2-S5; Table~\ref{table4}). These include experimentally produced alkaline magmas from \cite{iacovino2016}, basaltic melt inclusions from Kilauea \cite{tucker2019} and Gakkel Ridge \cite{bennett2019}, basaltic melt inclusions from Cerro Negro volcano, Nicaragua \cite{roggensack2001}, and rhyolite melt inclusions from the Taupo Volcanic Center, New Zealand \cite{myers2019} and a topaz rhyolite from the Rio Grande Rift \cite{mercer2015}. Where available, the calibration datasets for \acrshort{vesical} models are also provided (Supplementary Datasets S6-S7).

\begin{code}[H]
\centering
\caption[]{Example datasets included with \acrshort{vesical}}
\label{table4}\begin{minted}[breaklines]{python}
pd.read_excel("tables/Table_Example_Data.xlsx", index_col="Filename")
\end{minted}

\begin{minted}[breaklines]{text}
Explanation  \
Filename                                                  
example_data.xlsx  Example data used in this manuscript   
alkaline.xlsx      Experimental glasses                   
basalts.xlsx       Melt inclusion glasses                 
cerro_negro.xlsx   Melt inclusion glasses                 
rhyolites.xlsx     Melt inclusion glasses                 

                                  Compositions  \
Filename                                         
example_data.xlsx  Wide comp. range              
alkaline.xlsx      Basanite to Tephriphonolite   
basalts.xlsx       Basaltic                      
cerro_negro.xlsx   Basaltic                      
rhyolites.xlsx     Rhyolitic                     

                                                                                                                    Citations  
Filename                                                                                                                       
example_data.xlsx  Iacovino et al. (2016); Mercer et al. (2015); Myers et al. (2019); Roggensack (2001); Tucker et al. (2019)  
alkaline.xlsx      Iacovino et al. (2016)                                                                                      
basalts.xlsx       Tucker et al. (2019); Bennett et al. (2019)                                                                 
cerro_negro.xlsx   Roggensack (2001)                                                                                           
rhyolites.xlsx     Mercer et al. (2015); Myers et al. (2019)
\end{minted}
\end{code}

\paragraph{Defining a single sample}\label{defining-single-sample}

More advanced functionality of \acrshort{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 \texttt{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 \acrshort{vesical} Sample object to represent our sample composition.

\begin{description}
\item[\textbf{Class structure}] \texttt{Sample(composition, units='wtpt\_oxides', default\_normalization='none', default\_units='wtpt\_oxides')}
\item[\textbf{Required inputs}] \texttt{composition}: The composition of the sample in the format specified by the units parameter. The default is oxides in wt\%.
\item[\textbf{Optional inputs}] \texttt{units}, \texttt{default\_normalization}, and \texttt{default\_units} have the same meaning here as in the BatchFile class desribed above.
\item[\textbf{Outputs}] A special type of python object defined in the \acrshort{vesical} code known as a Sample object.
\end{description}

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

\begin{minted}[breaklines]{python}
mysample = v.Sample({'SiO2':  77.3,
             'TiO2':   0.08,
             'Al2O3': 12.6,
             'Fe2O3':  0.207,
             'Cr2O3':  0.0,
             'FeO':    0.473,
             'MnO':    0.0,
             'MgO':    0.03,
             'NiO':    0.0,
             'CoO':    0.0,
             'CaO':    0.43,
             'Na2O':   3.98,
             'K2O':    4.88,
             'P2O5':   0.0,
             'H2O':    6.5,
             'CO2':    0.05})
\end{minted}

To see the composition of mysample, use the \texttt{get\_composition(species=None, normalization=None, units=None, exclude\_volatiles=False, asSampleClass=False)} method. By default, the composition is returned exactly as input above. \texttt{species} can be set as an element or oxide (e.g., ``Si'' or ``SiO\textsubscript{2}'')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 \texttt{BatchFile.get\_data()} function, a sample composition can be returned as a dictionary (default) or as a \acrshort{vesical} Sample object (if \texttt{asSampleClass} is set to True).

\begin{minted}[breaklines]{python}
mysample.get_composition()
\end{minted}

\begin{minted}[breaklines]{text}
SiO2     77.300
TiO2     0.080 
Al2O3    12.600
Fe2O3    0.207 
Cr2O3    0.000 
FeO      0.473 
MnO      0.000 
MgO      0.030 
NiO      0.000 
CoO      0.000 
CaO      0.430 
Na2O     3.980 
K2O      4.880 
P2O5     0.000 
H2O      6.500 
CO2      0.050 
dtype: float64
\end{minted}

The oxides considered by \acrshort{vesical} are:

\begin{minted}[breaklines]{python}
print(v.oxides)
\end{minted}

\begin{minted}[breaklines]{text}
['SiO2', 'TiO2', 'Al2O3', 'Fe2O3', 'Cr2O3', 'FeO', 'MnO', 'MgO', 'NiO', 'CoO', 'CaO', 'Na2O', 'K2O', 'P2O5', 'H2O', 'CO2']
\end{minted}

\paragraph{Extracting a single sample from a batch file}\label{extracting-single-melt-composition}

Defined within the \texttt{BatchFile()} class, the method \texttt{get\_sample\_composition()} allows for the extraction of a melt composition from a loaded Excel or \acrshort{csv} file.

\begin{description}
\item[\textbf{Method structure}] \texttt{myfile.get\_sample\_composition(samplename, species=None, normalization=None, units=None, asSampleClass=False)}

\item[\textbf{Required inputs}] \texttt{samplename}: The name of the sample, as a string, as defined in the `Label' column of the input file.

\item[\textbf{Optional inputs}] \texttt{species}: This is used if only the concentration of a single species (either oxide or element) is desired.

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

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

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

\item[\textbf{Outputs}] The bulk composition stored in a dictionary or Sample object.

\end{description}

\begin{minted}[breaklines]{python}
"""To get composition from a specific sample in the input data:"""
sample_10 = myfile.get_sample_composition('10*', asSampleClass=True)

"""To see the extracted sample composition, uncomment the line below by removing the # and execute this code cell"""
#sample_10.get_composition()
\end{minted}

\begin{minted}[breaklines]{text}
'To see the extracted sample composition, uncomment the line below by removing the # and execute this code cell'
\end{minted}

\paragraph{Normalizing and transforming data}\label{normalizing-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 \cite{hughes2019}. As an example, when normalized, a volatile-free basalt with a measured SiO\textsubscript{2} content of 46 wt\% and an analytical total of 97 wt\% actually contains 47.4 wt\% SiO\textsubscript{2} (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, \acrshort{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\textsubscript{2}O or CO\textsubscript{2} are passed to the function, their un-normalized values will be retained in addition to the normalized non-volatile oxides, summing to \textgreater 100\%.

Normalization can be performed on a Sample object or on all samples within a BatchFile object using the \texttt{get\_composition()} or \texttt{get\_data()} methods (e.g., \texttt{myfile.get\_composition(normalization='standard')} or \texttt{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), \texttt{BatchFile.get\_composition()} returns only compositional data, where as \texttt{BatchFile.get\_data()} returns all data stored in the BatchFile object. The \texttt{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 \texttt{asSampleClass} or \texttt{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.

\begin{minted}[breaklines]{python}
"""Retrieve the standard normalization for one sample"""
mysample_normalized = mysample.get_composition(normalization="standard", asSampleClass=True)
#print(mysample_normalized.get_composition())

"""Retrieve the standard normalization for all samples in a BatchFile"""
myfile_normalized = myfile.get_data(normalization="standard", asBatchFile=True)
#print(myfile_normalized.get_data())
\end{minted}

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 \texttt{additionalvolatiles} and \texttt{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\textsubscript{2}O and CO\textsubscript{2} is expressed solely in terms of the absolute melt SiO\textsubscript{2} 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 \acrshort{vesical}'s additionalvolatiles normalization routine. For consistency with the web app, \acrshort{vesical} automatically uses the additionalvolatiles normaliztion during calculations with this model.

The implementation of MagmaSat in \acrshort{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 \acrshort{vesical}'s fixedvolatiles routine, and so this normalization is forced on all samples for such calculations with MagmaSat in \acrshort{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\textsubscript{2} content of 47.4 wt\%, 1000 ppm dissolved CO\textsubscript{2}, and an anhydrous (volatile-free) total of 96.77 wt\%:

\begin{minted}[breaklines]{python}
mybasalt = v.Sample({'SiO2': 47,
                    'TiO2': 1.01,
                    'Al2O3': 17.46,
                    'Fe2O3': 0.89,
                    'FeO': 7.18,
                    'MgO': 7.63,
                    'CaO': 12.44,
                    'Na2O': 2.65,
                    'K2O': 0.03,
                    'P2O5': 0.08,
                    'CO2': 0.1})
\end{minted}

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

\begin{minted}[breaklines]{python}
"""Normalize three ways"""
mybasalt_std = mybasalt.get_composition(normalization="standard", asSampleClass=True)
mybasalt_add = mybasalt.get_composition(normalization="additionalvolatiles", asSampleClass=True)
mybasalt_fix = mybasalt.get_composition(normalization="fixedvolatiles", asSampleClass=True)

"""Choose a model to test"""
mymodel = "IaconoMarziano"

for basalt, normtype in zip([mybasalt, mybasalt_std, mybasalt_add, mybasalt_fix],
                        ["Raw", "standard", "additionalvolatiles", "fixedvolatiles"]):
    print(str(normtype) +
          " Saturation Pressure = " +
          str(v.calculate_saturation_pressure(sample=basalt, temperature=1200, model=mymodel).result) +
          " bars")
\end{minted}

\begin{minted}[breaklines]{text}
Raw Saturation Pressure = 1848.0344238069026 bars
standard Saturation Pressure = 1906.5480423917331 bars
additionalvolatiles Saturation Pressure = 1848.2699609182075 bars
fixedvolatiles Saturation Pressure = 1848.263700904032 bars
\end{minted}

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

\subsubsection{Comparing User Data to Model Calibrations: Which Model Should I Use?}\label{comparing-user-data-to-model-calibrations}

MagmaSat is the most thermodynamically robust model implemented in \acrshort{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 \acrshort{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 \cite{wieser2021}.

For the remainder of this section, all example calculations are carried out with MagmaSat, the default model of \acrshort{vesical}. To use any other \acrshort{vesical} model, simply add \texttt{model=} and the name of the desired model in quotes to any calculation (e.g., \texttt{v.calculate\_dissolved\_volatiles(temperature=900, pressure=1000, model="Dixon")}. The model names recognized by \acrshort{vesical} are: MagmaSat, ShishkinaIdealMixing, Dixon, IaconoMarziano, Liu, AllisonCarbon, MooreWater. For more advanced use cases such as hybridizing models (see Section~\ref{code-hybridization}), pure-H\textsubscript{2}O and pure-CO\textsubscript{2} 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. \acrshort{vesical} includes some built-in functionality for comparing melt compositions from user loaded data to those in the datasets upon which each of the \acrshort{vesical} models is calibrated using the method \texttt{calib\_plot}. This can be visualized as a total alkalis vs silica (\acrshort{tas}) diagram (with fields and labels via the python tasplot library by J. Stevenson; \href{https://bitbucket.org/jsteven5/tasplot/src/master/}{https://bitbucket.org/jsteven5/tasplot/src/master/}; Figure~\ref{fig5} a) or as any x-y plot in which x and y are oxides (Figure~\ref{fig5} b).

\begin{description}
\item[\textbf{Method structure}] \texttt{calib\_plot(user\_data=None, model='all', plot\_type='TAS', zoom=None, save\_fig=False)}

\item[\textbf{Optional inputs}] \texttt{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.

\texttt{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'])

\texttt{plot\_type}: The default value is `\acrshort{tas}', which returns a total alkalis vs silica (\acrshort{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'.

\texttt{zoom}: The default is None in which case axes will be set to the default of 35$\leq$x$\leq$100 wt\% and 0$\leq$y$\leq$25 wt\% for \acrshort{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)].

\texttt{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...''.

\item[\textbf{Outputs}] The function returns fig and axes matplotlib objects defining a \acrshort{tas} or x-y plot of user data and model calibration data.

\end{description}

\begin{code}[H]
\centering
\begin{minted}[breaklines]{python}
v.calib_plot(user_data=myfile)
v.calib_plot(user_data=myfile, model='IaconoMarziano', plot_type='xy', x='SiO2', y='K2O', save_fig=False)
v.show()
\end{minted}

\includegraphics[width=0.7\linewidth]{files/fdc01f287b397f28e613546610e87968.png}

\includegraphics[width=0.7\linewidth]{files/f5b2871753ba03246b81384ed134b455.png}
\caption[]{Example calibration plots. a. The default plot with user\_data defined as myfile and no other options set. This produces a \acrshort{tas} diagram with the user data plotted atop data from calibration datasets for all models. b. A plot with all options specified. This example produces an x-y plot for user\_data (myfile) and the Iacono-Marziano calibration dataset where x and y are SiO\textsubscript{2} and K\textsubscript{2}O concentration in wt\%. Symbol shapes correspond to the volatile composition of experiments used to calibrate the model.}
\label{fig5}
\end{code}

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 (\href{https://matplotlib.org/3.2.1/index.html}{https://matplotlib.org/3.2.1/index.html}) if more complex plotting is desired.

\subsubsection{Calculating dissolved volatile concentrations}\label{calculating-dissolved-volatile-concentrations}


\bigskip
\centerline{\rule{13cm}{0.4pt}}
\bigskip

The \texttt{calculate\_dissolved\_volatiles()} function calculates the concentration of dissolved H\textsubscript{2}O and CO\textsubscript{2} in the melt at a given pressure-temperature condition and with a given H\textsubscript{2}O-CO\textsubscript{2} fluid composition, defined as the mole fraction of H\textsubscript{2}O in an H\textsubscript{2}O-CO\textsubscript{2} fluid (XH\textsubscript{2}O\textsuperscript{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\textsubscript{2}O/(H\textsubscript{2}O+CO\textsubscript{2}) in the system is equal to the user input mole fraction of H\textsubscript{2}O/(H\textsubscript{2}O+CO\textsubscript{2}) in the fluid. This is done by increasing the H\textsubscript{2}O and CO\textsubscript{2} 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\textsubscript{2}O or CO\textsubscript{2} in the system if the mole fraction of H\textsubscript{2}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\textsubscript{2}O\textsuperscript{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\textsubscript{2}O\textsuperscript{fluid} within 0.0001 of the user defined value.

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

\begin{description}
\item[\textbf{Method structure}] Single sample: \texttt{calculate\_dissolved\_volatiles(sample, temperature, pressure, X\_fluid=1, verbose=False, model='MagmaSat').result}

BatchFile batch process: \texttt{myfile.calculate\_dissolved\_volatiles(temperature, pressure, X\_fluid=1, print\_status=True, model='MagmaSat')}

\item[\textbf{Standard inputs}] \texttt{sample}, \texttt{temperature}, \texttt{pressure}, \texttt{X\_fluid}, \texttt{model} (see Section~\ref{class-arguments-and-definitions}).

\item[\textbf{Unique optional inputs}] \texttt{verbose}: \textit{Only for single sample calculations.} Default value is False, in which case H\textsubscript{2}O and CO\textsubscript{2} concentrations are returned. If set to True, additional parameters are returned in a dictionary: H\textsubscript{2}O and CO\textsubscript{2} concentrations in the fluid in mole fraction, temperature, pressure, and proportion of the fluid in the system in wt\%.

\texttt{print\_status}: \textit{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.

\item[\textbf{Calculated outputs}] If the single-sample method is used, a dictionary with keys `H2O' and `CO2' corresponding to the calculated dissolved H\textsubscript{2}O and CO\textsubscript{2} 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 \texttt{verbose} is set to True).

If the BatchFile method is used, a pandas DataFrame is returned with sample information plus calculated dissolved H\textsubscript{2}O and CO\textsubscript{2} 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.

\end{description}

\begin{minted}[breaklines]{python}
"""Calculate dissolved volatiles for sample 10*"""
v.calculate_dissolved_volatiles(sample=sample_10, temperature=900.0, pressure=2000.0,
                                X_fluid=0.5, verbose=True).result
\end{minted}

\begin{minted}[breaklines]{text}
{'H2O_liq': 2.69352739399806,
 'CO2_liq': 0.0638439414375309,
 'XH2O_fl': 0.500092686493868,
 'XCO2_fl': 0.499907313506132,
 'FluidProportion_wt': 0.18407321260435108}
\end{minted}

\begin{minted}[breaklines]{python}
"""Calculate dissolved for all samples in an BatchFile object"""
dissolved = myfile.calculate_dissolved_volatiles(temperature=900.0, pressure=2000.0, X_fluid=1, print_status=True)
\end{minted}

\begin{code}[H]
\centering
\caption[]{Modeled dissolved volatile concentrations}
\label{table5}\begin{minted}[breaklines]{python}
dissolved
\end{minted}

\begin{minted}[breaklines]{text}
SiO2      TiO2      Al2O3  Fe2O3  Cr2O3  \
Kil3 -6_1a                 48.249207  2.222114  11.692194  0.00   0.0     
Kil3 -6_3a                 48.295691  2.165357  11.755584  0.00   0.0     
Kil3 -6_4a                 49.124079  2.360984  12.172833  0.00   0.0     
10*                       47.960000  0.780000  18.770000  0.00   0.0     
19*                       49.640000  0.710000  18.050000  0.00   0.0     
25                        50.320000  0.720000  18.030000  0.00   0.0     
SAT -M12 -1                 62.600000  0.630000  17.300000  2.01   0.0     
SAT -M12 -2                 62.600000  0.630000  17.300000  2.01   0.0     
SAT -M12 -4                 62.600000  0.630000  17.300000  2.01   0.0     
samp. P1968a              76.974880  0.085516  3.110636   0.00   0.0     
samp. P1968b              76.943845  0.133125  3.169657   0.00   0.0     
samp. P1968c              77.187205  0.119506  3.167827   0.00   0.0     
samp. HPR3 -1_XL -3         75.413966  0.095164  14.077692  0.00   0.0     
samp. HPR3 -1_XL -4_INCL -1  76.613586  0.095843  13.476762  0.00   0.0     
AW -6                      48.030000  2.840000  18.120000  0.00   0.0     
AW -46                     52.980000  2.180000  20.490000  0.00   0.0     
KI -07                     44.610000  4.370000  14.410000  0.00   0.0     

                                FeO       MnO        MgO  NiO  CoO  ...  \
Kil3 -6_1a                 0.000000   0.079999  14.183817  0.0  0.0  ...   
Kil3 -6_3a                 0.000000   0.084045  13.403980  0.0  0.0  ...   
Kil3 -6_4a                 0.000000   0.098809  11.997699  0.0  0.0  ...   
10*                       10.920000  0.150000  6.860000   0.0  0.0  ...   
19*                       10.540000  0.190000  6.430000   0.0  0.0  ...   
25                        10.110000  0.140000  5.650000   0.0  0.0  ...   
SAT -M12 -1                 2.010000   0.060000  2.650000   0.0  0.0  ...   
SAT -M12 -2                 2.010000   0.060000  2.650000   0.0  0.0  ...   
SAT -M12 -4                 2.010000   0.060000  2.650000   0.0  0.0  ...   
samp. P1968a              4.788883   0.000000  12.549439  0.0  0.0  ...   
samp. P1968b              4.763435   0.000000  12.446403  0.0  0.0  ...   
samp. P1968c              4.814076   0.000000  12.229534  0.0  0.0  ...   
samp. HPR3 -1_XL -3         0.654992   0.125882  0.012003   0.0  0.0  ...   
samp. HPR3 -1_XL -4_INCL -1  0.620769   0.113495  0.032069   0.0  0.0  ...   
AW -6                      9.600000   0.230000  3.080000   0.0  0.0  ...   
AW -46                     5.540000   0.200000  2.000000   0.0  0.0  ...   
KI -07                     10.600000  0.170000  7.690000   0.0  0.0  ...   

                                        ROCK TYPE   Press         Temp  \
Kil3 -6_1a                 Basalt                   62.5    1299.094712   
Kil3 -6_3a                 Basalt                   128.0   1283.419991   
Kil3 -6_4a                 Basalt                   100.0   1255.153759   
10*                       Basalt                   2000.0  1200.000000   
19*                       Basalt                   2000.0  1200.000000   
25                        Basalt                   2000.0  1200.000000   
SAT -M12 -1                 Andesite                 703.0   1100.000000   
SAT -M12 -2                 Andesite                 1865.0  1100.000000   
SAT -M12 -4                 Andesite                 2985.0  1050.000000   
samp. P1968a              Rhyolite                 300.0   900.000000    
samp. P1968b              Rhyolite                 300.0   900.000000    
samp. P1968c              Rhyolite                 300.0   900.000000    
samp. HPR3 -1_XL -3         Rhyolite                 300.0   0.000000      
samp. HPR3 -1_XL -4_INCL -1  Rhyolite                 0.0     900.000000    
AW -6                      Phonotephrite            1500.0  1050.000000   
AW -46                     Basaltic -Trachyandesite  4000.0  1000.000000   
KI -07                     Basanite                 1500.0  1100.000000   

                          H2O_liq_VESIcal  CO2_liq_VESIcal  \
Kil3 -6_1a                 5.256561         0.0               
Kil3 -6_3a                 5.417720         0.0               
Kil3 -6_4a                 5.353421         0.0               
10*                       4.984021         0.0               
19*                       5.134419         0.0               
25                        5.189068         0.0               
SAT -M12 -1                 5.810439         0.0               
SAT -M12 -2                 5.810439         0.0               
SAT -M12 -4                 5.810439         0.0               
samp. P1968a              6.484749         0.0               
samp. P1968b              6.473813         0.0               
samp. P1968c              6.482109         0.0               
samp. HPR3 -1_XL -3         6.097630         0.0               
samp. HPR3 -1_XL -4_INCL -1  6.138658         0.0               
AW -6                      5.856636         0.0               
AW -46                     5.879457         0.0               
KI -07                     4.918430         0.0               

                          Temperature_C_VESIcal Pressure_bars_VESIcal  \
Kil3 -6_1a                 900.0                  2000.0                 
Kil3 -6_3a                 900.0                  2000.0                 
Kil3 -6_4a                 900.0                  2000.0                 
10*                       900.0                  2000.0                 
19*                       900.0                  2000.0                 
25                        900.0                  2000.0                 
SAT -M12 -1                 900.0                  2000.0                 
SAT -M12 -2                 900.0                  2000.0                 
SAT -M12 -4                 900.0                  2000.0                 
samp. P1968a              900.0                  2000.0                 
samp. P1968b              900.0                  2000.0                 
samp. P1968c              900.0                  2000.0                 
samp. HPR3 -1_XL -3         900.0                  2000.0                 
samp. HPR3 -1_XL -4_INCL -1  900.0                  2000.0                 
AW -6                      900.0                  2000.0                 
AW -46                     900.0                  2000.0                 
KI -07                     900.0                  2000.0                 

                         X_fluid_input_VESIcal     Model  Warnings  
Kil3 -6_1a                 1                     MagmaSat            
Kil3 -6_3a                 1                     MagmaSat            
Kil3 -6_4a                 1                     MagmaSat            
10*                       1                     MagmaSat            
19*                       1                     MagmaSat            
25                        1                     MagmaSat            
SAT -M12 -1                 1                     MagmaSat            
SAT -M12 -2                 1                     MagmaSat            
SAT -M12 -4                 1                     MagmaSat            
samp. P1968a              1                     MagmaSat            
samp. P1968b              1                     MagmaSat            
samp. P1968c              1                     MagmaSat            
samp. HPR3 -1_XL -3         1                     MagmaSat            
samp. HPR3 -1_XL -4_INCL -1  1                     MagmaSat            
AW -6                      1                     MagmaSat            
AW -46                     1                     MagmaSat            
KI -07                     1                     MagmaSat            

[17 rows x 27 columns]
\end{minted}
\end{code}

\subsubsection{Calculating equilibrium fluid compositions}\label{calculating-equilibrium-fluid-compositions}


\bigskip
\centerline{\rule{13cm}{0.4pt}}
\bigskip

The \texttt{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\textsubscript{2}O and CO\textsubscript{2} 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\textsubscript{2}O and CO\textsubscript{2} concentrations in the fluid.

\begin{description}
\item[\textbf{Method structure}] Single sample: \texttt{calculate\_equilibrium\_fluid\_comp(sample, temperature, pressure, verbose=False, model='MagmaSat').result}

BatchFile batch process: \texttt{myfile.calculate\_equilibrium\_fluid\_comp(temperature, pressure=None, print\_status=False, model='MagmaSat')}

\item[\textbf{Standard inputs}] \texttt{sample}, \texttt{temperature}, \texttt{pressure}, \texttt{model} (see Section~\ref{class-arguments-and-definitions}).

\item[\textbf{Unique optional inputs}] \texttt{verbose}: \textit{Only for single sample calculations.} Default value is False, in which case H\textsubscript{2}O and CO\textsubscript{2} concentrations in the fluid in mol fraction are returned. If set to True, additional parameters are returned in a dictionary: H\textsubscript{2}O and CO\textsubscript{2} concentrations in the fluid, mass of the fluid in grams, and proportion of the fluid in the system in wt\%.

\texttt{print\_status}: \textit{Only for batch calculations.} The default value is False. If True is passed, the progress of the calculation will be printed to the terminal.

\item[\textbf{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 \texttt{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.

\end{description}

\begin{minted}[breaklines]{python}
"""Calculate fluid composition for the extracted sample"""
v.calculate_equilibrium_fluid_comp(sample=sample_10, temperature=900.0, pressure=100.0).result
\end{minted}

\begin{minted}[breaklines]{text}
{'CO2': 0.00528661429366132, 'H2O': 0.994713385706339}
\end{minted}

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 \cite{tucker2019} 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\textsubscript{2}O and CO\textsubscript{2}.

\begin{code}[H]
\centering
\caption[]{Isothermally modeled equilibrium fluid compositions}
\label{table6}\begin{minted}[breaklines]{python}
"""Calculate fluid composition for all samples in an BatchFile object"""
eqfluid = myfile.calculate_equilibrium_fluid_comp(temperature=900.0, pressure=1000.0)
eqfluid
\end{minted}

\begin{minted}[breaklines]{text}
SiO2      TiO2      Al2O3  Fe2O3  Cr2O3  \
Kil3 -6_1a                 48.249207  2.222114  11.692194  0.00   0.0     
Kil3 -6_3a                 48.295691  2.165357  11.755584  0.00   0.0     
Kil3 -6_4a                 49.124079  2.360984  12.172833  0.00   0.0     
10*                       47.960000  0.780000  18.770000  0.00   0.0     
19*                       49.640000  0.710000  18.050000  0.00   0.0     
25                        50.320000  0.720000  18.030000  0.00   0.0     
SAT -M12 -1                 62.600000  0.630000  17.300000  2.01   0.0     
SAT -M12 -2                 62.600000  0.630000  17.300000  2.01   0.0     
SAT -M12 -4                 62.600000  0.630000  17.300000  2.01   0.0     
samp. P1968a              76.974880  0.085516  3.110636   0.00   0.0     
samp. P1968b              76.943845  0.133125  3.169657   0.00   0.0     
samp. P1968c              77.187205  0.119506  3.167827   0.00   0.0     
samp. HPR3 -1_XL -3         75.413966  0.095164  14.077692  0.00   0.0     
samp. HPR3 -1_XL -4_INCL -1  76.613586  0.095843  13.476762  0.00   0.0     
AW -6                      48.030000  2.840000  18.120000  0.00   0.0     
AW -46                     52.980000  2.180000  20.490000  0.00   0.0     
KI -07                     44.610000  4.370000  14.410000  0.00   0.0     

                                FeO       MnO        MgO  NiO  CoO  ...  \
Kil3 -6_1a                 0.000000   0.079999  14.183817  0.0  0.0  ...   
Kil3 -6_3a                 0.000000   0.084045  13.403980  0.0  0.0  ...   
Kil3 -6_4a                 0.000000   0.098809  11.997699  0.0  0.0  ...   
10*                       10.920000  0.150000  6.860000   0.0  0.0  ...   
19*                       10.540000  0.190000  6.430000   0.0  0.0  ...   
25                        10.110000  0.140000  5.650000   0.0  0.0  ...   
SAT -M12 -1                 2.010000   0.060000  2.650000   0.0  0.0  ...   
SAT -M12 -2                 2.010000   0.060000  2.650000   0.0  0.0  ...   
SAT -M12 -4                 2.010000   0.060000  2.650000   0.0  0.0  ...   
samp. P1968a              4.788883   0.000000  12.549439  0.0  0.0  ...   
samp. P1968b              4.763435   0.000000  12.446403  0.0  0.0  ...   
samp. P1968c              4.814076   0.000000  12.229534  0.0  0.0  ...   
samp. HPR3 -1_XL -3         0.654992   0.125882  0.012003   0.0  0.0  ...   
samp. HPR3 -1_XL -4_INCL -1  0.620769   0.113495  0.032069   0.0  0.0  ...   
AW -6                      9.600000   0.230000  3.080000   0.0  0.0  ...   
AW -46                     5.540000   0.200000  2.000000   0.0  0.0  ...   
KI -07                     10.600000  0.170000  7.690000   0.0  0.0  ...   

                                        CITATION                ROCK TYPE  \
Kil3 -6_1a                 Tucker et al. (2019)    Basalt                    
Kil3 -6_3a                 Tucker et al. (2019)    Basalt                    
Kil3 -6_4a                 Tucker et al. (2019)    Basalt                    
10*                       Roggensack (2001)       Basalt                    
19*                       Roggensack (2001)       Basalt                    
25                        Roggensack (2001)       Basalt                    
SAT -M12 -1                 Moore et al. (1998)     Andesite                  
SAT -M12 -2                 Moore et al. (1998)     Andesite                  
SAT -M12 -4                 Moore et al. (1998)     Andesite                  
samp. P1968a              Myers et al. (2019)     Rhyolite                  
samp. P1968b              Myers et al. (2019)     Rhyolite                  
samp. P1968c              Myers et al. (2019)     Rhyolite                  
samp. HPR3 -1_XL -3         Mercer et al. (2015)    Rhyolite                  
samp. HPR3 -1_XL -4_INCL -1  Mercer et al. (2015)    Rhyolite                  
AW -6                      Iacovino et al. (2016)  Phonotephrite             
AW -46                     Iacovino et al. (2016)  Basaltic -Trachyandesite   
KI -07                     Iacovino et al. (2016)  Basanite                  

                           Press         Temp  XH2O_fl_VESIcal  \
Kil3 -6_1a                 62.5    1299.094712  0.000000          
Kil3 -6_3a                 128.0   1283.419991  0.000000          
Kil3 -6_4a                 100.0   1255.153759  0.000000          
10*                       2000.0  1200.000000  0.984531          
19*                       2000.0  1200.000000  0.974997          
25                        2000.0  1200.000000  0.990107          
SAT -M12 -1                 703.0   1100.000000  1.000000          
SAT -M12 -2                 1865.0  1100.000000  1.000000          
SAT -M12 -4                 2985.0  1050.000000  1.000000          
samp. P1968a              300.0   900.000000   0.977773          
samp. P1968b              300.0   900.000000   0.996799          
samp. P1968c              300.0   900.000000   0.997028          
samp. HPR3 -1_XL -3         300.0   0.000000     0.997770          
samp. HPR3 -1_XL -4_INCL -1  0.0     900.000000   0.997273          
AW -6                      1500.0  1050.000000  0.261572          
AW -46                     4000.0  1000.000000  0.897441          
KI -07                     1500.0  1100.000000  0.826014          

                          XCO2_fl_VESIcal Temperature_C_VESIcal  \
Kil3 -6_1a                 0.000000         900.0                  
Kil3 -6_3a                 0.000000         900.0                  
Kil3 -6_4a                 0.000000         900.0                  
10*                       0.015469         900.0                  
19*                       0.025003         900.0                  
25                        0.009893         900.0                  
SAT -M12 -1                 0.000000         900.0                  
SAT -M12 -2                 0.000000         900.0                  
SAT -M12 -4                 0.000000         900.0                  
samp. P1968a              0.022227         900.0                  
samp. P1968b              0.003201         900.0                  
samp. P1968c              0.002972         900.0                  
samp. HPR3 -1_XL -3         0.002230         900.0                  
samp. HPR3 -1_XL -4_INCL -1  0.002727         900.0                  
AW -6                      0.738428         900.0                  
AW -46                     0.102559         900.0                  
KI -07                     0.173986         900.0                  

                         Pressure_bars_VESIcal     Model  \
Kil3 -6_1a                 1000.0                MagmaSat   
Kil3 -6_3a                 1000.0                MagmaSat   
Kil3 -6_4a                 1000.0                MagmaSat   
10*                       1000.0                MagmaSat   
19*                       1000.0                MagmaSat   
25                        1000.0                MagmaSat   
SAT -M12 -1                 1000.0                MagmaSat   
SAT -M12 -2                 1000.0                MagmaSat   
SAT -M12 -4                 1000.0                MagmaSat   
samp. P1968a              1000.0                MagmaSat   
samp. P1968b              1000.0                MagmaSat   
samp. P1968c              1000.0                MagmaSat   
samp. HPR3 -1_XL -3         1000.0                MagmaSat   
samp. HPR3 -1_XL -4_INCL -1  1000.0                MagmaSat   
AW -6                      1000.0                MagmaSat   
AW -46                     1000.0                MagmaSat   
KI -07                     1000.0                MagmaSat   

                                                          Warnings  
Kil3 -6_1a                 Sample not saturated at these conditions  
Kil3 -6_3a                 Sample not saturated at these conditions  
Kil3 -6_4a                 Sample not saturated at these conditions  
10*                                                                 
19*                                                                 
25                                                                  
SAT -M12 -1                                                           
SAT -M12 -2                                                           
SAT -M12 -4                                                           
samp. P1968a                                                        
samp. P1968b                                                        
samp. P1968c                                                        
samp. HPR3 -1_XL -3                                                   
samp. HPR3 -1_XL -4_INCL -1                                            
AW -6                                                                
AW -46                                                               
KI -07                                                               

[17 rows x 26 columns]
\end{minted}
\end{code}

Below, we calculate equilibrium fluid compositions for the same dataset using temperatures and pressures as defined in the input data (Table~\ref{table3}). 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. \acrshort{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.

\begin{code}[H]
\centering
\caption[]{\textbf{Modeled equilibrium fluid compositions with unique temperatures.} Warnings ``Bad temperature'' and ``Bad pressure'' indicate that no data (or 0.0 value data) was given for the temperature or pressure of that sample, in which case the calculation of that sample is skipped.}
\label{table7}\begin{minted}[breaklines]{python}
"""Calculate fluid composition for all samples with unique pressure and temperature values for each sample.
Pressure and temperature values are taken from columns named "Press" and "Temp" in the example BatchFile"""

eqfluid_wtemps = myfile.calculate_equilibrium_fluid_comp(temperature='Temp', pressure='Press')
eqfluid_wtemps
\end{minted}

\begin{minted}[breaklines]{text}
SiO2      TiO2      Al2O3  Fe2O3  Cr2O3  \
Kil3 -6_1a                 48.249207  2.222114  11.692194  0.00   0.0     
Kil3 -6_3a                 48.295691  2.165357  11.755584  0.00   0.0     
Kil3 -6_4a                 49.124079  2.360984  12.172833  0.00   0.0     
10*                       47.960000  0.780000  18.770000  0.00   0.0     
19*                       49.640000  0.710000  18.050000  0.00   0.0     
25                        50.320000  0.720000  18.030000  0.00   0.0     
SAT -M12 -1                 62.600000  0.630000  17.300000  2.01   0.0     
SAT -M12 -2                 62.600000  0.630000  17.300000  2.01   0.0     
SAT -M12 -4                 62.600000  0.630000  17.300000  2.01   0.0     
samp. P1968a              76.974880  0.085516  3.110636   0.00   0.0     
samp. P1968b              76.943845  0.133125  3.169657   0.00   0.0     
samp. P1968c              77.187205  0.119506  3.167827   0.00   0.0     
samp. HPR3 -1_XL -3         75.413966  0.095164  14.077692  0.00   0.0     
samp. HPR3 -1_XL -4_INCL -1  76.613586  0.095843  13.476762  0.00   0.0     
AW -6                      48.030000  2.840000  18.120000  0.00   0.0     
AW -46                     52.980000  2.180000  20.490000  0.00   0.0     
KI -07                     44.610000  4.370000  14.410000  0.00   0.0     

                                FeO       MnO        MgO  NiO  CoO  ...  \
Kil3 -6_1a                 0.000000   0.079999  14.183817  0.0  0.0  ...   
Kil3 -6_3a                 0.000000   0.084045  13.403980  0.0  0.0  ...   
Kil3 -6_4a                 0.000000   0.098809  11.997699  0.0  0.0  ...   
10*                       10.920000  0.150000  6.860000   0.0  0.0  ...   
19*                       10.540000  0.190000  6.430000   0.0  0.0  ...   
25                        10.110000  0.140000  5.650000   0.0  0.0  ...   
SAT -M12 -1                 2.010000   0.060000  2.650000   0.0  0.0  ...   
SAT -M12 -2                 2.010000   0.060000  2.650000   0.0  0.0  ...   
SAT -M12 -4                 2.010000   0.060000  2.650000   0.0  0.0  ...   
samp. P1968a              4.788883   0.000000  12.549439  0.0  0.0  ...   
samp. P1968b              4.763435   0.000000  12.446403  0.0  0.0  ...   
samp. P1968c              4.814076   0.000000  12.229534  0.0  0.0  ...   
samp. HPR3 -1_XL -3         0.654992   0.125882  0.012003   0.0  0.0  ...   
samp. HPR3 -1_XL -4_INCL -1  0.620769   0.113495  0.032069   0.0  0.0  ...   
AW -6                      9.600000   0.230000  3.080000   0.0  0.0  ...   
AW -46                     5.540000   0.200000  2.000000   0.0  0.0  ...   
KI -07                     10.600000  0.170000  7.690000   0.0  0.0  ...   

                               H2O       CO2                CITATION  \
Kil3 -6_1a                 0.424695  0.002873  Tucker et al. (2019)     
Kil3 -6_3a                 0.425984  0.006786  Tucker et al. (2019)     
Kil3 -6_4a                 0.437758  0.004984  Tucker et al. (2019)     
10*                       4.500000  0.047900  Roggensack (2001)        
19*                       5.100000  0.111300  Roggensack (2001)        
25                        5.200000  0.043700  Roggensack (2001)        
SAT -M12 -1                 2.620000  0.000000  Moore et al. (1998)      
SAT -M12 -2                 5.030000  0.000000  Moore et al. (1998)      
SAT -M12 -4                 6.760000  0.000000  Moore et al. (1998)      
samp. P1968a              4.340000  0.007000  Myers et al. (2019)      
samp. P1968b              5.850000  0.012300  Myers et al. (2019)      
samp. P1968c              5.754571  0.010663  Myers et al. (2019)      
samp. HPR3 -1_XL -3         5.943750  0.010000  Mercer et al. (2015)     
samp. HPR3 -1_XL -4_INCL -1  5.340000  0.008000  Mercer et al. (2015)     
AW -6                      1.420000  0.129800  Iacovino et al. (2016)   
AW -46                     4.760000  0.343900  Iacovino et al. (2016)   
KI -07                     2.900000  0.113100  Iacovino et al. (2016)   

                                        ROCK TYPE   Press         Temp  \
Kil3 -6_1a                 Basalt                   62.5    1299.094712   
Kil3 -6_3a                 Basalt                   128.0   1283.419991   
Kil3 -6_4a                 Basalt                   100.0   1255.153759   
10*                       Basalt                   2000.0  1200.000000   
19*                       Basalt                   2000.0  1200.000000   
25                        Basalt                   2000.0  1200.000000   
SAT -M12 -1                 Andesite                 703.0   1100.000000   
SAT -M12 -2                 Andesite                 1865.0  1100.000000   
SAT -M12 -4                 Andesite                 2985.0  1050.000000   
samp. P1968a              Rhyolite                 300.0   900.000000    
samp. P1968b              Rhyolite                 300.0   900.000000    
samp. P1968c              Rhyolite                 300.0   900.000000    
samp. HPR3 -1_XL -3         Rhyolite                 300.0   0.000000      
samp. HPR3 -1_XL -4_INCL -1  Rhyolite                 0.0     900.000000    
AW -6                      Phonotephrite            1500.0  1050.000000   
AW -46                     Basaltic -Trachyandesite  4000.0  1000.000000   
KI -07                     Basanite                 1500.0  1100.000000   

                         XH2O_fl_VESIcal XCO2_fl_VESIcal     Model  \
Kil3 -6_1a                 0.586164        0.413836        MagmaSat   
Kil3 -6_3a                 0.286160        0.713840        MagmaSat   
Kil3 -6_4a                 0.377439        0.622561        MagmaSat   
10*                       0.892371        0.107629        MagmaSat   
19*                       0.918888        0.081112        MagmaSat   
25                        0.955803        0.044197        MagmaSat   
SAT -M12 -1                 1.000000        0.000000        MagmaSat   
SAT -M12 -2                 1.000000        0.000000        MagmaSat   
SAT -M12 -4                 1.000000        0.000000        MagmaSat   
samp. P1968a              0.998764        0.001236        MagmaSat   
samp. P1968b              0.998686        0.001314        MagmaSat   
samp. P1968c              0.998831        0.001169        MagmaSat   
samp. HPR3 -1_XL -3        NaN             NaN              MagmaSat   
samp. HPR3 -1_XL -4_INCL -1 NaN             NaN              MagmaSat   
AW -6                      0.000000        0.000000        MagmaSat   
AW -46                     0.492213        0.507787        MagmaSat   
KI -07                     0.681758        0.318242        MagmaSat   

                                                          Warnings  
Kil3 -6_1a                                                           
Kil3 -6_3a                                                           
Kil3 -6_4a                                                           
10*                                                                 
19*                                                                 
25                                                                  
SAT -M12 -1                                                           
SAT -M12 -2                                                           
SAT -M12 -4                                                           
samp. P1968a                                                        
samp. P1968b                                                        
samp. P1968c                                                        
samp. HPR3 -1_XL -3         Calculation skipped. Bad temperature.     
samp. HPR3 -1_XL -4_INCL -1  Calculation skipped. Bad pressure.        
AW -6                      Sample not saturated at these conditions  
AW -46                                                               
KI -07                                                               

[17 rows x 24 columns]
\end{minted}
\end{code}

\paragraph{Converting fluid composition units}

The fluid composition is always returned in units of mol fraction. Two functions exist to transform only the H\textsubscript{2}O-CO\textsubscript{2} 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\textsubscript{2}O and CO\textsubscript{2} concentrations in the fluid. The default values for column names are set to those that may be returned by \acrshort{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).

\begin{description}
\item[\textbf{Method structure}] Mol fraction to wt\%: \texttt{fluid\_molfrac\_to\_wt(data, H2O\_colname='XH2O\_fl\_VESIcal', CO2\_colname='XCO2\_fl\_VESIcal')}

Wt\% to mol fraction: \texttt{fluid\_wt\_to\_molfrac(data, H2O\_colname='H2O\_fl\_wt', CO2\_colname='CO2\_fl\_wt')}

\item[\textbf{Required inputs}] \texttt{data}: A pandas DataFrame containing columns for H\textsubscript{2}O and CO\textsubscript{2} concentrations in the fluid.

\item[\textbf{Optional inputs}] \texttt{H2O\_colname} and \texttt{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\textsubscript{2}O and CO\textsubscript{2} concentrations in the fluid.

\item[\textbf{Calculated outputs}] The original data passed plus newly calculated values are returned in a DataFrame.

\end{description}

\begin{code}[H]
\centering
\caption[]{Equilibrium fluid compositions converted from mol fraction to wt\%}
\label{table8}\begin{minted}[breaklines]{python}
"""Converting from mol fraction to wt%"""
eqfluid_wt = v.fluid_molfrac_to_wt(eqfluid)
eqfluid_wt
\end{minted}

\begin{minted}[breaklines]{text}
SiO2      TiO2      Al2O3  Fe2O3  Cr2O3  \
Kil3 -6_1a                 48.249207  2.222114  11.692194  0.00   0.0     
Kil3 -6_3a                 48.295691  2.165357  11.755584  0.00   0.0     
Kil3 -6_4a                 49.124079  2.360984  12.172833  0.00   0.0     
10*                       47.960000  0.780000  18.770000  0.00   0.0     
19*                       49.640000  0.710000  18.050000  0.00   0.0     
25                        50.320000  0.720000  18.030000  0.00   0.0     
SAT -M12 -1                 62.600000  0.630000  17.300000  2.01   0.0     
SAT -M12 -2                 62.600000  0.630000  17.300000  2.01   0.0     
SAT -M12 -4                 62.600000  0.630000  17.300000  2.01   0.0     
samp. P1968a              76.974880  0.085516  3.110636   0.00   0.0     
samp. P1968b              76.943845  0.133125  3.169657   0.00   0.0     
samp. P1968c              77.187205  0.119506  3.167827   0.00   0.0     
samp. HPR3 -1_XL -3         75.413966  0.095164  14.077692  0.00   0.0     
samp. HPR3 -1_XL -4_INCL -1  76.613586  0.095843  13.476762  0.00   0.0     
AW -6                      48.030000  2.840000  18.120000  0.00   0.0     
AW -46                     52.980000  2.180000  20.490000  0.00   0.0     
KI -07                     44.610000  4.370000  14.410000  0.00   0.0     

                                FeO       MnO        MgO  NiO  CoO  ...  \
Kil3 -6_1a                 0.000000   0.079999  14.183817  0.0  0.0  ...   
Kil3 -6_3a                 0.000000   0.084045  13.403980  0.0  0.0  ...   
Kil3 -6_4a                 0.000000   0.098809  11.997699  0.0  0.0  ...   
10*                       10.920000  0.150000  6.860000   0.0  0.0  ...   
19*                       10.540000  0.190000  6.430000   0.0  0.0  ...   
25                        10.110000  0.140000  5.650000   0.0  0.0  ...   
SAT -M12 -1                 2.010000   0.060000  2.650000   0.0  0.0  ...   
SAT -M12 -2                 2.010000   0.060000  2.650000   0.0  0.0  ...   
SAT -M12 -4                 2.010000   0.060000  2.650000   0.0  0.0  ...   
samp. P1968a              4.788883   0.000000  12.549439  0.0  0.0  ...   
samp. P1968b              4.763435   0.000000  12.446403  0.0  0.0  ...   
samp. P1968c              4.814076   0.000000  12.229534  0.0  0.0  ...   
samp. HPR3 -1_XL -3         0.654992   0.125882  0.012003   0.0  0.0  ...   
samp. HPR3 -1_XL -4_INCL -1  0.620769   0.113495  0.032069   0.0  0.0  ...   
AW -6                      9.600000   0.230000  3.080000   0.0  0.0  ...   
AW -46                     5.540000   0.200000  2.000000   0.0  0.0  ...   
KI -07                     10.600000  0.170000  7.690000   0.0  0.0  ...   

                           Press         Temp  XH2O_fl_VESIcal  \
Kil3 -6_1a                 62.5    1299.094712  0.000000          
Kil3 -6_3a                 128.0   1283.419991  0.000000          
Kil3 -6_4a                 100.0   1255.153759  0.000000          
10*                       2000.0  1200.000000  0.984531          
19*                       2000.0  1200.000000  0.974997          
25                        2000.0  1200.000000  0.990107          
SAT -M12 -1                 703.0   1100.000000  1.000000          
SAT -M12 -2                 1865.0  1100.000000  1.000000          
SAT -M12 -4                 2985.0  1050.000000  1.000000          
samp. P1968a              300.0   900.000000   0.977773          
samp. P1968b              300.0   900.000000   0.996799          
samp. P1968c              300.0   900.000000   0.997028          
samp. HPR3 -1_XL -3         300.0   0.000000     0.997770          
samp. HPR3 -1_XL -4_INCL -1  0.0     900.000000   0.997273          
AW -6                      1500.0  1050.000000  0.261572          
AW -46                     4000.0  1000.000000  0.897441          
KI -07                     1500.0  1100.000000  0.826014          

                          XCO2_fl_VESIcal  Temperature_C_VESIcal  \
Kil3 -6_1a                 0.000000         900.0                   
Kil3 -6_3a                 0.000000         900.0                   
Kil3 -6_4a                 0.000000         900.0                   
10*                       0.015469         900.0                   
19*                       0.025003         900.0                   
25                        0.009893         900.0                   
SAT -M12 -1                 0.000000         900.0                   
SAT -M12 -2                 0.000000         900.0                   
SAT -M12 -4                 0.000000         900.0                   
samp. P1968a              0.022227         900.0                   
samp. P1968b              0.003201         900.0                   
samp. P1968c              0.002972         900.0                   
samp. HPR3 -1_XL -3         0.002230         900.0                   
samp. HPR3 -1_XL -4_INCL -1  0.002727         900.0                   
AW -6                      0.738428         900.0                   
AW -46                     0.102559         900.0                   
KI -07                     0.173986         900.0                   

                          Pressure_bars_VESIcal     Model  \
Kil3 -6_1a                 1000.0                 MagmaSat   
Kil3 -6_3a                 1000.0                 MagmaSat   
Kil3 -6_4a                 1000.0                 MagmaSat   
10*                       1000.0                 MagmaSat   
19*                       1000.0                 MagmaSat   
25                        1000.0                 MagmaSat   
SAT -M12 -1                 1000.0                 MagmaSat   
SAT -M12 -2                 1000.0                 MagmaSat   
SAT -M12 -4                 1000.0                 MagmaSat   
samp. P1968a              1000.0                 MagmaSat   
samp. P1968b              1000.0                 MagmaSat   
samp. P1968c              1000.0                 MagmaSat   
samp. HPR3 -1_XL -3         1000.0                 MagmaSat   
samp. HPR3 -1_XL -4_INCL -1  1000.0                 MagmaSat   
AW -6                      1000.0                 MagmaSat   
AW -46                     1000.0                 MagmaSat   
KI -07                     1000.0                 MagmaSat   

                                                          Warnings  \
Kil3 -6_1a                 Sample not saturated at these conditions   
Kil3 -6_3a                 Sample not saturated at these conditions   
Kil3 -6_4a                 Sample not saturated at these conditions   
10*                                                                  
19*                                                                  
25                                                                   
SAT -M12 -1                                                            
SAT -M12 -2                                                            
SAT -M12 -4                                                            
samp. P1968a                                                         
samp. P1968b                                                         
samp. P1968c                                                         
samp. HPR3 -1_XL -3                                                    
samp. HPR3 -1_XL -4_INCL -1                                             
AW -6                                                                 
AW -46                                                                
KI -07                                                                

                           H2O_fl_wt  CO2_fl_wt  
Kil3 -6_1a                NaN         NaN         
Kil3 -6_3a                NaN         NaN         
Kil3 -6_4a                NaN         NaN         
10*                       96.304445   3.695555   
19*                       94.106168   5.893832   
25                        97.617908   2.382092   
SAT -M12 -1                 100.000000  0.000000   
SAT -M12 -2                 100.000000  0.000000   
SAT -M12 -4                 100.000000  0.000000   
samp. P1968a              94.740209   5.259791   
samp. P1968b              99.221744   0.778256   
samp. P1968c              99.277291   0.722709   
samp. HPR3 -1_XL -3         99.457027   0.542973   
samp. HPR3 -1_XL -4_INCL -1  99.336700   0.663300   
AW -6                      12.666745   87.333255  
AW -46                     78.179789   21.820211  
KI -07                     66.031544   33.968456  

[17 rows x 28 columns]
\end{minted}
\end{code}

\begin{code}[H]
\centering
\caption[]{Equilibrium fluid compositions converted from wt\% to mol fraction}
\label{table9}\begin{minted}[breaklines]{python}
"""Converting from wt% to mol fraction"""
eqfluid_mol = v.fluid_wt_to_molfrac(eqfluid_wt)
eqfluid_mol
\end{minted}

\begin{minted}[breaklines]{text}
SiO2      TiO2      Al2O3  Fe2O3  Cr2O3  \
Kil3 -6_1a                 48.249207  2.222114  11.692194  0.00   0.0     
Kil3 -6_3a                 48.295691  2.165357  11.755584  0.00   0.0     
Kil3 -6_4a                 49.124079  2.360984  12.172833  0.00   0.0     
10*                       47.960000  0.780000  18.770000  0.00   0.0     
19*                       49.640000  0.710000  18.050000  0.00   0.0     
25                        50.320000  0.720000  18.030000  0.00   0.0     
SAT -M12 -1                 62.600000  0.630000  17.300000  2.01   0.0     
SAT -M12 -2                 62.600000  0.630000  17.300000  2.01   0.0     
SAT -M12 -4                 62.600000  0.630000  17.300000  2.01   0.0     
samp. P1968a              76.974880  0.085516  3.110636   0.00   0.0     
samp. P1968b              76.943845  0.133125  3.169657   0.00   0.0     
samp. P1968c              77.187205  0.119506  3.167827   0.00   0.0     
samp. HPR3 -1_XL -3         75.413966  0.095164  14.077692  0.00   0.0     
samp. HPR3 -1_XL -4_INCL -1  76.613586  0.095843  13.476762  0.00   0.0     
AW -6                      48.030000  2.840000  18.120000  0.00   0.0     
AW -46                     52.980000  2.180000  20.490000  0.00   0.0     
KI -07                     44.610000  4.370000  14.410000  0.00   0.0     

                                FeO       MnO        MgO  NiO  CoO  ...  \
Kil3 -6_1a                 0.000000   0.079999  14.183817  0.0  0.0  ...   
Kil3 -6_3a                 0.000000   0.084045  13.403980  0.0  0.0  ...   
Kil3 -6_4a                 0.000000   0.098809  11.997699  0.0  0.0  ...   
10*                       10.920000  0.150000  6.860000   0.0  0.0  ...   
19*                       10.540000  0.190000  6.430000   0.0  0.0  ...   
25                        10.110000  0.140000  5.650000   0.0  0.0  ...   
SAT -M12 -1                 2.010000   0.060000  2.650000   0.0  0.0  ...   
SAT -M12 -2                 2.010000   0.060000  2.650000   0.0  0.0  ...   
SAT -M12 -4                 2.010000   0.060000  2.650000   0.0  0.0  ...   
samp. P1968a              4.788883   0.000000  12.549439  0.0  0.0  ...   
samp. P1968b              4.763435   0.000000  12.446403  0.0  0.0  ...   
samp. P1968c              4.814076   0.000000  12.229534  0.0  0.0  ...   
samp. HPR3 -1_XL -3         0.654992   0.125882  0.012003   0.0  0.0  ...   
samp. HPR3 -1_XL -4_INCL -1  0.620769   0.113495  0.032069   0.0  0.0  ...   
AW -6                      9.600000   0.230000  3.080000   0.0  0.0  ...   
AW -46                     5.540000   0.200000  2.000000   0.0  0.0  ...   
KI -07                     10.600000  0.170000  7.690000   0.0  0.0  ...   

                          XH2O_fl_VESIcal  XCO2_fl_VESIcal  \
Kil3 -6_1a                 0.000000         0.000000          
Kil3 -6_3a                 0.000000         0.000000          
Kil3 -6_4a                 0.000000         0.000000          
10*                       0.984531         0.015469          
19*                       0.974997         0.025003          
25                        0.990107         0.009893          
SAT -M12 -1                 1.000000         0.000000          
SAT -M12 -2                 1.000000         0.000000          
SAT -M12 -4                 1.000000         0.000000          
samp. P1968a              0.977773         0.022227          
samp. P1968b              0.996799         0.003201          
samp. P1968c              0.997028         0.002972          
samp. HPR3 -1_XL -3         0.997770         0.002230          
samp. HPR3 -1_XL -4_INCL -1  0.997273         0.002727          
AW -6                      0.261572         0.738428          
AW -46                     0.897441         0.102559          
KI -07                     0.826014         0.173986          

                          Temperature_C_VESIcal  Pressure_bars_VESIcal  \
Kil3 -6_1a                 900.0                  1000.0                  
Kil3 -6_3a                 900.0                  1000.0                  
Kil3 -6_4a                 900.0                  1000.0                  
10*                       900.0                  1000.0                  
19*                       900.0                  1000.0                  
25                        900.0                  1000.0                  
SAT -M12 -1                 900.0                  1000.0                  
SAT -M12 -2                 900.0                  1000.0                  
SAT -M12 -4                 900.0                  1000.0                  
samp. P1968a              900.0                  1000.0                  
samp. P1968b              900.0                  1000.0                  
samp. P1968c              900.0                  1000.0                  
samp. HPR3 -1_XL -3         900.0                  1000.0                  
samp. HPR3 -1_XL -4_INCL -1  900.0                  1000.0                  
AW -6                      900.0                  1000.0                  
AW -46                     900.0                  1000.0                  
KI -07                     900.0                  1000.0                  

                             Model                                  Warnings  \
Kil3 -6_1a                 MagmaSat  Sample not saturated at these conditions   
Kil3 -6_3a                 MagmaSat  Sample not saturated at these conditions   
Kil3 -6_4a                 MagmaSat  Sample not saturated at these conditions   
10*                       MagmaSat                                             
19*                       MagmaSat                                             
25                        MagmaSat                                             
SAT -M12 -1                 MagmaSat                                             
SAT -M12 -2                 MagmaSat                                             
SAT -M12 -4                 MagmaSat                                             
samp. P1968a              MagmaSat                                             
samp. P1968b              MagmaSat                                             
samp. P1968c              MagmaSat                                             
samp. HPR3 -1_XL -3         MagmaSat                                             
samp. HPR3 -1_XL -4_INCL -1  MagmaSat                                             
AW -6                      MagmaSat                                             
AW -46                     MagmaSat                                             
KI -07                     MagmaSat                                             

                           H2O_fl_wt  CO2_fl_wt   XH2O_fl   XCO2_fl  
Kil3 -6_1a                NaN         NaN        NaN       NaN        
Kil3 -6_3a                NaN         NaN        NaN       NaN        
Kil3 -6_4a                NaN         NaN        NaN       NaN        
10*                       96.304445   3.695555   0.984531  0.015469  
19*                       94.106168   5.893832   0.974997  0.025003  
25                        97.617908   2.382092   0.990107  0.009893  
SAT -M12 -1                 100.000000  0.000000   1.000000  0.000000  
SAT -M12 -2                 100.000000  0.000000   1.000000  0.000000  
SAT -M12 -4                 100.000000  0.000000   1.000000  0.000000  
samp. P1968a              94.740209   5.259791   0.977773  0.022227  
samp. P1968b              99.221744   0.778256   0.996799  0.003201  
samp. P1968c              99.277291   0.722709   0.997028  0.002972  
samp. HPR3 -1_XL -3         99.457027   0.542973   0.997770  0.002230  
samp. HPR3 -1_XL -4_INCL -1  99.336700   0.663300   0.997273  0.002727  
AW -6                      12.666745   87.333255  0.261572  0.738428  
AW -46                     78.179789   21.820211  0.897441  0.102559  
KI -07                     66.031544   33.968456  0.826014  0.173986  

[17 rows x 30 columns]
\end{minted}
\end{code}

\subsubsection{Calculating saturation pressures}\label{calculating-saturation-pressures}


\bigskip
\centerline{\rule{13cm}{0.4pt}}
\bigskip

The \texttt{calculate\_saturation\_pressure()} function calculates the minimum pressure at which a given silicate melt with known temperature and H\textsubscript{2}O and CO\textsubscript{2} 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 \textgreater 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.

\begin{description}
\item[\textbf{Method structure}] Single sample: \texttt{calculate\_saturation\_pressure(sample, temperature, verbose=False, model='MagmaSat').result}

BatchFile batch process: \texttt{myfile.calculate\_saturation\_pressure(temperature, print\_status=True, model='MagmaSat')}

\item[\textbf{Standard inputs}] \texttt{sample}, \texttt{temperature}, \texttt{model} (see Section~\ref{class-arguments-and-definitions}).

\item[\textbf{Unique optional inputs}] \texttt{verbose}: \textit{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\textsubscript{2}O and CO\textsubscript{2} concentrations in the fluid, mass of the fluid in grams, and proportion of the fluid in the system in wt\%.

\texttt{print\_status}: \textit{Only for batch calculations.} The default value is True, in which case the progress of the calculation will be printed to the terminal.

\item[\textbf{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 \texttt{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.

\end{description}

\begin{minted}[breaklines]{python}
"""Calculate the saturation pressure of the single sample we defined in Section 3.3.1 at 925 degrees C"""
v.calculate_saturation_pressure(sample=mysample, temperature=925.0, verbose=True).result
\end{minted}

\begin{minted}[breaklines]{text}
{'SaturationP_bars': 2960,
 'FluidMass_grams': 0.0018160337487088,
 'FluidProportion_wt': 0.0018160337487087978,
 'XH2O_fl': 0.838064480487942,
 'XCO2_fl': 0.161935519512058}
\end{minted}

\begin{code}[H]
\centering
\caption[]{Isothermally modeled saturation pressures}
\label{table10}\begin{minted}[breaklines]{python}
"""Calculate the saturation pressure for all samples in an BatchFile object at 925 degrees C"""
satPs = myfile.calculate_saturation_pressure(temperature=925.0)
satPs
\end{minted}

\begin{minted}[breaklines]{text}
SiO2      TiO2      Al2O3  Fe2O3  Cr2O3  \
Kil3 -6_1a                 48.249207  2.222114  11.692194  0.00   0.0     
Kil3 -6_3a                 48.295691  2.165357  11.755584  0.00   0.0     
Kil3 -6_4a                 49.124079  2.360984  12.172833  0.00   0.0     
10*                       47.960000  0.780000  18.770000  0.00   0.0     
19*                       49.640000  0.710000  18.050000  0.00   0.0     
25                        50.320000  0.720000  18.030000  0.00   0.0     
SAT -M12 -1                 62.600000  0.630000  17.300000  2.01   0.0     
SAT -M12 -2                 62.600000  0.630000  17.300000  2.01   0.0     
SAT -M12 -4                 62.600000  0.630000  17.300000  2.01   0.0     
samp. P1968a              76.974880  0.085516  3.110636   0.00   0.0     
samp. P1968b              76.943845  0.133125  3.169657   0.00   0.0     
samp. P1968c              77.187205  0.119506  3.167827   0.00   0.0     
samp. HPR3 -1_XL -3         75.413966  0.095164  14.077692  0.00   0.0     
samp. HPR3 -1_XL -4_INCL -1  76.613586  0.095843  13.476762  0.00   0.0     
AW -6                      48.030000  2.840000  18.120000  0.00   0.0     
AW -46                     52.980000  2.180000  20.490000  0.00   0.0     
KI -07                     44.610000  4.370000  14.410000  0.00   0.0     

                                FeO       MnO        MgO  NiO  CoO  ...  \
Kil3 -6_1a                 0.000000   0.079999  14.183817  0.0  0.0  ...   
Kil3 -6_3a                 0.000000   0.084045  13.403980  0.0  0.0  ...   
Kil3 -6_4a                 0.000000   0.098809  11.997699  0.0  0.0  ...   
10*                       10.920000  0.150000  6.860000   0.0  0.0  ...   
19*                       10.540000  0.190000  6.430000   0.0  0.0  ...   
25                        10.110000  0.140000  5.650000   0.0  0.0  ...   
SAT -M12 -1                 2.010000   0.060000  2.650000   0.0  0.0  ...   
SAT -M12 -2                 2.010000   0.060000  2.650000   0.0  0.0  ...   
SAT -M12 -4                 2.010000   0.060000  2.650000   0.0  0.0  ...   
samp. P1968a              4.788883   0.000000  12.549439  0.0  0.0  ...   
samp. P1968b              4.763435   0.000000  12.446403  0.0  0.0  ...   
samp. P1968c              4.814076   0.000000  12.229534  0.0  0.0  ...   
samp. HPR3 -1_XL -3         0.654992   0.125882  0.012003   0.0  0.0  ...   
samp. HPR3 -1_XL -4_INCL -1  0.620769   0.113495  0.032069   0.0  0.0  ...   
AW -6                      9.600000   0.230000  3.080000   0.0  0.0  ...   
AW -46                     5.540000   0.200000  2.000000   0.0  0.0  ...   
KI -07                     10.600000  0.170000  7.690000   0.0  0.0  ...   

                           Press         Temp  SaturationP_bars_VESIcal  \
Kil3 -6_1a                 62.5    1299.094712  60                         
Kil3 -6_3a                 128.0   1283.419991  130                        
Kil3 -6_4a                 100.0   1255.153759  100                        
10*                       2000.0  1200.000000  2500                       
19*                       2000.0  1200.000000  3600                       
25                        2000.0  1200.000000  2750                       
SAT -M12 -1                 703.0   1100.000000  550                        
SAT -M12 -2                 1865.0  1100.000000  1590                       
SAT -M12 -4                 2985.0  1050.000000  2540                       
samp. P1968a              300.0   900.000000   1100                       
samp. P1968b              300.0   900.000000   1790                       
samp. P1968c              300.0   900.000000   1730                       
samp. HPR3 -1_XL -3         300.0   0.000000     2090                       
samp. HPR3 -1_XL -4_INCL -1  0.0     900.000000   1730                       
AW -6                      1500.0  1050.000000  1220                       
AW -46                     4000.0  1000.000000  4800                       
KI -07                     1500.0  1100.000000  1510                       

                          Temperature_C_VESIcal  XH2O_fl_VESIcal  \
Kil3 -6_1a                 925.0                  0.469913          
Kil3 -6_3a                 925.0                  0.215529          
Kil3 -6_4a                 925.0                  0.292354          
10*                       925.0                  0.796514          
19*                       925.0                  0.702654          
25                        925.0                  0.836895          
SAT -M12 -1                 925.0                  1.000000          
SAT -M12 -2                 925.0                  1.000000          
SAT -M12 -4                 925.0                  1.000000          
samp. P1968a              925.0                  0.972472          
samp. P1968b              925.0                  0.972875          
samp. P1968c              925.0                  0.975614          
samp. HPR3 -1_XL -3         925.0                  0.951891          
samp. HPR3 -1_XL -4_INCL -1  925.0                  0.950741          
AW -6                      925.0                  0.231708          
AW -46                     925.0                  0.456750          
KI -07                     925.0                  0.684729          

                          XCO2_fl_VESIcal FluidMass_grams_VESIcal  \
Kil3 -6_1a                 0.530087         0.000836                 
Kil3 -6_3a                 0.784471         0.000038                 
Kil3 -6_4a                 0.707646         0.000635                 
10*                       0.203486         0.001232                 
19*                       0.297346         0.000797                 
25                        0.163105         0.000226                 
SAT -M12 -1                 0.000000         0.012903                 
SAT -M12 -2                 0.000000         0.001052                 
SAT -M12 -4                 0.000000         0.016093                 
samp. P1968a              0.027528         0.007924                 
samp. P1968b              0.027125         0.006671                 
samp. P1968c              0.024386         0.008637                 
samp. HPR3 -1_XL -3         0.048109         0.002941                 
samp. HPR3 -1_XL -4_INCL -1  0.049259         0.002864                 
AW -6                      0.768292         0.000093                 
AW -46                     0.543250         0.000938                 
KI -07                     0.315271         0.000431                 

                         FluidSystem_wt_VESIcal     Model  Warnings  
Kil3 -6_1a                 0.000836               MagmaSat            
Kil3 -6_3a                 0.000038               MagmaSat            
Kil3 -6_4a                 0.000635               MagmaSat            
10*                       0.001232               MagmaSat            
19*                       0.000797               MagmaSat            
25                        0.000226               MagmaSat            
SAT -M12 -1                 0.012903               MagmaSat            
SAT -M12 -2                 0.001052               MagmaSat            
SAT -M12 -4                 0.016093               MagmaSat            
samp. P1968a              0.007924               MagmaSat            
samp. P1968b              0.006671               MagmaSat            
samp. P1968c              0.008637               MagmaSat            
samp. HPR3 -1_XL -3         0.002941               MagmaSat            
samp. HPR3 -1_XL -4_INCL -1  0.002864               MagmaSat            
AW -6                      0.000093               MagmaSat            
AW -46                     0.000938               MagmaSat            
KI -07                     0.000431               MagmaSat            

[17 rows x 28 columns]
\end{minted}
\end{code}

\begin{code}[H]
\centering
\caption[]{\textbf{Modeled saturation pressures with unique temperatures} The warning ``Bad temperature'' indicates that no data (or 0.0 value data) was given for the temperature of that sample, in which case the calculation of that sample is skipped.}
\label{table11}\begin{minted}[breaklines]{python}
"""Calculate the saturation pressure for all samples in an BatchFile object, taking temperature
values from a column named "Temp" in the BatchFile"""
satPs_wtemps = myfile.calculate_saturation_pressure(temperature="Temp")
satPs_wtemps
\end{minted}

\begin{minted}[breaklines]{text}
SiO2      TiO2      Al2O3  Fe2O3  Cr2O3  \
Kil3 -6_1a                 48.249207  2.222114  11.692194  0.00   0.0     
Kil3 -6_3a                 48.295691  2.165357  11.755584  0.00   0.0     
Kil3 -6_4a                 49.124079  2.360984  12.172833  0.00   0.0     
10*                       47.960000  0.780000  18.770000  0.00   0.0     
19*                       49.640000  0.710000  18.050000  0.00   0.0     
25                        50.320000  0.720000  18.030000  0.00   0.0     
SAT -M12 -1                 62.600000  0.630000  17.300000  2.01   0.0     
SAT -M12 -2                 62.600000  0.630000  17.300000  2.01   0.0     
SAT -M12 -4                 62.600000  0.630000  17.300000  2.01   0.0     
samp. P1968a              76.974880  0.085516  3.110636   0.00   0.0     
samp. P1968b              76.943845  0.133125  3.169657   0.00   0.0     
samp. P1968c              77.187205  0.119506  3.167827   0.00   0.0     
samp. HPR3 -1_XL -3         75.413966  0.095164  14.077692  0.00   0.0     
samp. HPR3 -1_XL -4_INCL -1  76.613586  0.095843  13.476762  0.00   0.0     
AW -6                      48.030000  2.840000  18.120000  0.00   0.0     
AW -46                     52.980000  2.180000  20.490000  0.00   0.0     
KI -07                     44.610000  4.370000  14.410000  0.00   0.0     

                                FeO       MnO        MgO  NiO  CoO  ...  \
Kil3 -6_1a                 0.000000   0.079999  14.183817  0.0  0.0  ...   
Kil3 -6_3a                 0.000000   0.084045  13.403980  0.0  0.0  ...   
Kil3 -6_4a                 0.000000   0.098809  11.997699  0.0  0.0  ...   
10*                       10.920000  0.150000  6.860000   0.0  0.0  ...   
19*                       10.540000  0.190000  6.430000   0.0  0.0  ...   
25                        10.110000  0.140000  5.650000   0.0  0.0  ...   
SAT -M12 -1                 2.010000   0.060000  2.650000   0.0  0.0  ...   
SAT -M12 -2                 2.010000   0.060000  2.650000   0.0  0.0  ...   
SAT -M12 -4                 2.010000   0.060000  2.650000   0.0  0.0  ...   
samp. P1968a              4.788883   0.000000  12.549439  0.0  0.0  ...   
samp. P1968b              4.763435   0.000000  12.446403  0.0  0.0  ...   
samp. P1968c              4.814076   0.000000  12.229534  0.0  0.0  ...   
samp. HPR3 -1_XL -3         0.654992   0.125882  0.012003   0.0  0.0  ...   
samp. HPR3 -1_XL -4_INCL -1  0.620769   0.113495  0.032069   0.0  0.0  ...   
AW -6                      9.600000   0.230000  3.080000   0.0  0.0  ...   
AW -46                     5.540000   0.200000  2.000000   0.0  0.0  ...   
KI -07                     10.600000  0.170000  7.690000   0.0  0.0  ...   

                                        ROCK TYPE   Press         Temp  \
Kil3 -6_1a                 Basalt                   62.5    1299.094712   
Kil3 -6_3a                 Basalt                   128.0   1283.419991   
Kil3 -6_4a                 Basalt                   100.0   1255.153759   
10*                       Basalt                   2000.0  1200.000000   
19*                       Basalt                   2000.0  1200.000000   
25                        Basalt                   2000.0  1200.000000   
SAT -M12 -1                 Andesite                 703.0   1100.000000   
SAT -M12 -2                 Andesite                 1865.0  1100.000000   
SAT -M12 -4                 Andesite                 2985.0  1050.000000   
samp. P1968a              Rhyolite                 300.0   900.000000    
samp. P1968b              Rhyolite                 300.0   900.000000    
samp. P1968c              Rhyolite                 300.0   900.000000    
samp. HPR3 -1_XL -3         Rhyolite                 300.0   0.000000      
samp. HPR3 -1_XL -4_INCL -1  Rhyolite                 0.0     900.000000    
AW -6                      Phonotephrite            1500.0  1050.000000   
AW -46                     Basaltic -Trachyandesite  4000.0  1000.000000   
KI -07                     Basanite                 1500.0  1100.000000   

                          SaturationP_bars_VESIcal  XH2O_fl_VESIcal  \
Kil3 -6_1a                 60.0                      0.493184          
Kil3 -6_3a                 110.0                     0.266595          
Kil3 -6_4a                 90.0                      0.337738          
10*                       2540.0                    0.817548          
19*                       3650.0                    0.725724          
25                        2850.0                    0.855214          
SAT -M12 -1                 580.0                     1.000000          
SAT -M12 -2                 1650.0                    1.000000          
SAT -M12 -4                 2610.0                    1.000000          
samp. P1968a              1090.0                    0.972916          
samp. P1968b              1780.0                    0.973133          
samp. P1968c              1720.0                    0.975860          
samp. HPR3 -1_XL -3        NaN                       NaN                
samp. HPR3 -1_XL -4_INCL -1  1730.0                    0.951017          
AW -6                      1280.0                    0.228644          
AW -46                     4910.0                    0.458904          
KI -07                     1590.0                    0.679643          

                          XCO2_fl_VESIcal FluidMass_grams_VESIcal  \
Kil3 -6_1a                 0.506816         0.000610                 
Kil3 -6_3a                 0.733405         0.000700                 
Kil3 -6_4a                 0.662262         0.000807                 
10*                       0.182452         0.001532                 
19*                       0.274276         0.000669                 
25                        0.144786         0.000849                 
SAT -M12 -1                 0.000000         0.003442                 
SAT -M12 -2                 0.000000         0.015280                 
SAT -M12 -4                 0.000000         0.008153                 
samp. P1968a              0.027084         0.008855                 
samp. P1968b              0.026867         0.005916                 
samp. P1968c              0.024140         0.008088                 
samp. HPR3 -1_XL -3        NaN              NaN                       
samp. HPR3 -1_XL -4_INCL -1  0.048983         0.003350                 
AW -6                      0.771356         0.001475                 
AW -46                     0.541096         0.001767                 
KI -07                     0.320357         0.001914                 

                         FluidSystem_wt_VESIcal     Model  \
Kil3 -6_1a                 0.000610               MagmaSat   
Kil3 -6_3a                 0.000700               MagmaSat   
Kil3 -6_4a                 0.000807               MagmaSat   
10*                       0.001532               MagmaSat   
19*                       0.000669               MagmaSat   
25                        0.000849               MagmaSat   
SAT -M12 -1                 0.003442               MagmaSat   
SAT -M12 -2                 0.015280               MagmaSat   
SAT -M12 -4                 0.008153               MagmaSat   
samp. P1968a              0.008855               MagmaSat   
samp. P1968b              0.005916               MagmaSat   
samp. P1968c              0.008088               MagmaSat   
samp. HPR3 -1_XL -3        NaN                     MagmaSat   
samp. HPR3 -1_XL -4_INCL -1  0.003350               MagmaSat   
AW -6                      0.001475               MagmaSat   
AW -46                     0.001767               MagmaSat   
KI -07                     0.001914               MagmaSat   

                                                       Warnings  
Kil3 -6_1a                                                        
Kil3 -6_3a                                                        
Kil3 -6_4a                                                        
10*                                                              
19*                                                              
25                                                               
SAT -M12 -1                                                        
SAT -M12 -2                                                        
SAT -M12 -4                                                        
samp. P1968a                                                     
samp. P1968b                                                     
samp. P1968c                                                     
samp. HPR3 -1_XL -3         Calculation skipped. Bad temperature.  
samp. HPR3 -1_XL -4_INCL -1                                         
AW -6                                                             
AW -46                                                            
KI -07                                                            

[17 rows x 27 columns]
\end{minted}
\end{code}

\subsubsection{Calculating isobars and isopleths}\label{calculating-plotting-isobars-isopleths}


\bigskip
\centerline{\rule{13cm}{0.4pt}}
\bigskip

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\textsubscript{2}O\textsuperscript{fluid} values, respectively, passed as lists.

The calculation is performed by iterating through possible concentrations of H\textsubscript{2}O and CO\textsubscript{2} and calculating the equilibrium state for the system. The iteration begins at a fixed H\textsubscript{2}O concentration, increasing the CO\textsubscript{2} concentration in steps of 0.1 wt\% until a fluid phase is stable. The H\textsubscript{2}O concentration is then increased by 0.5 wt\% and CO\textsubscript{2} is again increased from 0 until a fluid phase is stable. This process is repeated for H\textsubscript{2}O values ranging from 0--15 wt\%. The H\textsubscript{2}O and CO\textsubscript{2} 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\textsubscript{2}O-CO\textsubscript{2} 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 \texttt{smooth\_isobars} or \texttt{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\textsubscript{2}O\textsuperscript{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.

\begin{description}
\item[\textbf{Method structure}] \textit{Only single sample calculations.} \texttt{calculate\_isobars\_and\_isopleths(sample, temperature, pressure\_list, isopleth\_list=None, smooth\_isobars=True, smooth\_isopleths=True, print\_status=True, model="MagmaSat").result}

\item[\textbf{Standard inputs}] \texttt{sample}, \texttt{temperature}, \texttt{model} (see Section~\ref{class-arguments-and-definitions}).

\item[\textbf{Unique required inputs}] \texttt{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.

\item[\textbf{Unique optional inputs}] \texttt{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\textsubscript{2}O (XH\textsubscript{2}O\textsuperscript{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.

\texttt{smooth\_isobars} and \texttt{smooth\_isopleths}: The default value for both of these arguments is True, in which case polynomials will be fit to the computed data points.

\texttt{print\_status}: The default value is True. If True, the progress of the calculations will be printed to the terminal.

\item[\textbf{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\textsubscript{2}O and CO\textsubscript{2} in the melt in wt\%. Columns in the isopleth dataframe are `XH2O\_fl', `H2O\_liq', and `CO2\_liq', corresponding to XH\textsubscript{2}O\textsuperscript{fluid} and dissolved H\textsubscript{2}O and CO\textsubscript{2} in the melt in wt\%.

\end{description}

\begin{minted}[breaklines]{python}
"""Define all variables to be passed to the function for calculating isobars and isopleths"""
"""Define the temperature in degrees C"""
temperature = 1200.0

"""Define a list of pressures in bars:"""
pressures = [1000.0, 2000.0, 3000.0]
\end{minted}

Next, the H\textsubscript{2}O and CO\textsubscript{2} 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.

\begin{minted}[breaklines]{python}
isobars, isopleths = v.calculate_isobars_and_isopleths(sample=sample_10,
                                            temperature=temperature,
                                            pressure_list=pressures,
                                            isopleth_list=[0.25,0.5,0.75]).result
\end{minted}

\begin{minted}[breaklines]{text}
Calculating isobar at 1000.0 bars
 done.                                                                                                                           
Calculating isobar at 2000.0 bars
 done.                                                                                                                           
Calculating isobar at 3000.0 bars
 done.                                                                                                                           
Done!
\end{minted}

\subsubsection{Calculating degassing paths}\label{calculating-plotting-degassing-paths}


\bigskip
\centerline{\rule{13cm}{0.4pt}}
\bigskip

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 \texttt{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.

\begin{description}
\item[\textbf{Method structure}] \textit{Only single-sample calculations.} \texttt{calculate\_degassing\_path(sample, temperature, pressure="saturation", fractionate\_vapor=0.0, init\_vapor=0.0, steps=50, model="MagmaSat").result}

\item[\textbf{Standard inputs}] \texttt{sample}, \texttt{temperature}, \texttt{model} (see Section~\ref{class-arguments-and-definitions}).

\item[\textbf{Unique optional inputs}] \texttt{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.

\texttt{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 \textless 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.

\texttt{init\_vapor}: Default value is 0.0. Specifies the amount of vapor (in wt\%) coexisting with the melt before degassing.

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

\item[\textbf{Calculated outputs}] The function returns a pandas DataFrame with columns as: `Pressure\_bars', `H2O\_liq' and `CO2\_liq' (the concentration of H\textsubscript{2}O and CO\textsubscript{2} in the melt, in wt\%), `XH2O\_fl' and `XCO2\_fl' (the composition of the H\textsubscript{2}O-CO\textsubscript{2} fluid, in mol fraction), and `FluidProportion\_wt' (the proportion of fluid in the fluid-melt system, in wt\%).

\end{description}

\textbf{Note}: The following two cells can take up to a minute to execute.

\begin{minted}[breaklines]{python}
temp = 1200 #temperature in °C

"""Calculate open, closed, and closed + 2 wt% initial vapor"""
closed_df = v.calculate_degassing_path(sample=sample_10, temperature=temp).result
open_df = v.calculate_degassing_path(sample=sample_10, temperature=temp, fractionate_vapor=1.0).result
half_df = v.calculate_degassing_path(sample=sample_10, temperature=temp, fractionate_vapor=0.5).result
exsolved_df = v.calculate_degassing_path(sample=sample_10, temperature=temp, init_vapor=2.0).result
"""Calculate closed -system degassing starting from a pressure of 2000 bars"""
start2000_df = v.calculate_degassing_path(sample=sample_10, temperature=temp, pressure=2000.0).result
\end{minted}

\begin{minted}[breaklines]{text}
[====================] 100% Calculating degassing path...
[====================] 100% Calculating degassing path...
[====================] 100% Calculating degassing path...
[====================] 100% Calculating degassing path...
[====================] 100% Calculating degassing path...
\end{minted}

\subsubsection{Plotting}\label{plotting}

After calculating isobars, isopleths, and degassing paths, any or all of these may be plotted in an H\textsubscript{2}O versus CO\textsubscript{2} 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. \acrshort{vesical}'s \texttt{plot} function takes in lists of pandas DataFrames with calculated isobar, isopleth, and degassing path information (e.g., output from \texttt{calculate\_isobars\_and\_isopleths()} or \texttt{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\textsubscript{2}O and CO\textsubscript{2} 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 \texttt{plot} function. Multiple sets of plottable data can be passed. For example, isobars calculated with two different models can be passed to the \texttt{isobars} argument as a list.

\acrshort{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 (\href{https://matplotlib.org/3.2.1/index.html}{https://matplotlib.org/3.2.1/index.html}) if more complex plotting is desired. If preferred, \acrshort{vesical} outputs can be saved to an Excel or \acrshort{csv} file (see Section~\ref{exporting-data}), 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 \texttt{plot()} should be saved to objects such as fig, ax as:

\texttt{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:

\texttt{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 \texttt{v.show()}.

\begin{description}
\item[\textbf{Method structure}] \texttt{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)}

\item[\textbf{Optional inputs}] \texttt{isobars}: DataFrame object containing isobar information as calculated by \texttt{calculate\_isobars\_and\_isopleths()}. Or a list
of DataFrame objects.

\texttt{isopleths}: DataFrame object containing isopleth information as calculated by \texttt{calculate\_isobars\_and\_isopleths()}. Or a list of DataFrame objects.

\texttt{degassing\_paths}: List of DataFrames with degassing information as generated by \texttt{calculate\_degassing\_path()}.

\texttt{custom\_H2O}: List of floats or array-like shapes of H\textsubscript{2}O concentration values to plot as points. For example \texttt{myfile.data['H2O']} is one array-like shape (here, pandas.Series) of H\textsubscript{2}O values. Must be passed with \texttt{custom\_CO2} and must be same length as \texttt{custom\_CO2}.

\texttt{custom\_CO2}: List of floats or array-like shapes of CO\textsubscript{2} values to plot as points. For example \texttt{myfile.data['CO2']} is one array-like shape of CO\textsubscript{2} values. Must be passed with \texttt{custom\_H2O} and must be same length as \texttt{custom\_H2O}.

\texttt{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.

\texttt{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\textsubscript{2}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\textsubscript{2}O value.

\texttt{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.

\texttt{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.

\texttt{custom\_colors} and \texttt{custom\_symbols}: Custom colors and symbol shapes can be specified for (\texttt{custom\_H2O}, \texttt{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 \texttt{custom\_H2O} and \texttt{custom\_CO2}. If nothing is specified for custom\_colors, \acrshort{vesical}'s default colors will be used. If nothing is specified for \texttt{custom\_symbols}, all points will be plotted as filled circles.

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

\texttt{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.

\item[\textbf{Advanced inputs:} Most users will not need to use these inputs.] \texttt{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.

\texttt{smooth\_isobars} and \texttt{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 \texttt{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 \texttt{calculate\_isobars\_and\_isopleths()} (see documentation on this method).

\item[\textbf{Calculated outputs}] The function returns fig and axes matplotlib objects defining a plot with x-axis as H\textsubscript{2}O wt\% in the melt and y-axis as CO\textsubscript{2} 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\textsubscript{2}O and CO\textsubscript{2} in a melt equilibrated along a path of decreasing pressure, is plotted if passed.

\end{description}

\paragraph{A simple example: Isobars and isopleths}\label{isobars-and-isopleths-example}

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\textsubscript{2}O\textsuperscript{fluid} calculated for sample `10*' at 1,200 °C in Section~\ref{calculating-plotting-isobars-isopleths} onto one plot.

\begin{code}[H]
\centering
\begin{minted}[breaklines]{python}
fig, ax = v.plot(isobars=isobars, isopleths=isopleths)
v.show()
\end{minted}

\includegraphics[width=0.7\linewidth]{files/7f44e41a58e430814cbaf716cf371bf2.png}
\caption[]{Isobars and isopleths calculated for the sample, temperature, pressures, XH\textsubscript{2}O\textsuperscript{fluid} values, and with the model as defined in Section~\ref{calculating-plotting-isobars-isopleths}. Manuscript default values are sample \texttt{10*} at a 1,200 °C with isobars at 1,000, 2,000, and 3,000 bars, isopleths at XH\textsubscript{2}O\textsuperscript{fluid} = 0, 0.25, 0.5, 0.75, and 1 calculated with MagmaSat}
\label{fig6}
\end{code}

When plotting isobars and isopleths via MagmaSat, the values calculated by \texttt{calculate\_isobars\_and\_isopleths()} are used to calculate polynomial fits using Numpy's \texttt{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  \cite{newman2002, iacono2012, iacovino2013} 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 \texttt{smooth\_isobars=True} and \texttt{smooth\_isopleths=True} to \texttt{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 \texttt{smooth\_isobars\_and\_isopleths()}. For this advanced case, we refer the reader to the documentation.

\paragraph{A simple example: Degassing paths}\label{degassing-paths-example}

Here we plot all four degassing paths calculated for sample `10*' at 1,200 °C in Section~\ref{calculating-plotting-degassing-paths} onto one plot. We designate labels of ``Open'', ``Half'', ``Closed'', and ``Exsolved'' for the legend.

\begin{code}[H]
\centering
\begin{minted}[breaklines]{python}
fig, ax = v.plot(degassing_paths=[open_df, half_df, closed_df, exsolved_df],
            degassing_path_labels=["Open", "Half", "Closed", "Exsolved"])
fig, ax = v.plot(degassing_paths=[exsolved_df, start2000_df],
            degassing_path_labels=["Exsolved", "2000 bars"])
v.show()
\end{minted}

\includegraphics[width=0.7\linewidth]{files/60e09b4352190c978578174217f7e4b8.png}

\includegraphics[width=0.7\linewidth]{files/735b6bf806d598d137cdd17b4ac2f77f.png}
\caption[]{Degassing paths calculated for the sample, temperature, degassing style, initial exsolved fluid wt\%, starting pressure, and model as designated in Section~\ref{calculating-plotting-degassing-paths}. Default manuscript values are sample \texttt{10*} at 1,200 °C. ``Open'', ``Half'', and ``Closed'' curves in (a) represent open-system, partially open-system (50\% fractionated fluid), and closed-system degassing paths, respectively, starting at the saturation pressure. The ``Exsolved'' curve in (b) represents closed-system degassing with an initial exsolved fluid wt\% = 2.0. The ``2000'' curve in (b) represents closed-system degassing calculated starting at a pressure of 2,000 bars.}
\label{fig7}
\end{code}

\paragraph{Plotting multiple calculations}\label{plotting-multiple-calculations}

One of the major advantages to \acrshort{vesical} over any other modeling tool is the ability to quickly calculate and plot multiple calculations. \acrshort{vesical}'s \texttt{plot()} function is built on top of the popular Matplotlib python library and is designed to work with any \acrshort{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\textsubscript{2}O on the x-axis and CO\textsubscript{2} on the y-axis. \texttt{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 \acrshort{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.

\subparagraph{Isobars, isopleths, and degassing paths}

In this example we will use data imported in Section~\ref{comparing-user-data-to-model-calibrations} and calculations performed in Sections Section~\ref{calculating-saturation-pressures} and Section~\ref{calculating-plotting-isobars-isopleths}. Of course, all of the data calculated with \acrshort{vesical} can be exported to an Excel or \acrshort{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:

\begin{itemize}
\item Isobars calculated at 1200 °C and pressures of 1,000, 2,000, and 3,000 bars for sample 10*
\item Isopleths calculated at 1200 °C and XH\textsubscript{2}O\textsuperscript{fluid} values of 0, 0.25, 0.5, 0.75, and 1 for sample 10*
\item An open-system degassing path for sample 10*
\item A closed-system degassing path for sample 10*
\end{itemize}

\begin{code}[H]
\centering
\begin{minted}[breaklines]{python}
fig, ax = v.plot(isobars=isobars,
                   isopleths=isopleths,
                   degassing_paths=[open_df, closed_df],
                   degassing_path_labels=["Open System", "Closed System"])
v.show()
\end{minted}

\includegraphics[width=0.7\linewidth]{files/e2400da73ecb125e8d29c6f21be3bd5f.png}
\caption[]{Example of plotting multiple calculations on one plot. Isobars and isopleths as defined in Section~\ref{calculating-plotting-isobars-isopleths} and shown in Section~\ref{isobars-and-isopleths-example} and degassing curves as defined in Section~\ref{calculating-plotting-degassing-paths} and shown in Section~\ref{degassing-paths-example}. Default manuscript values are for sample \texttt{10*} at 1,200 °C with isobars at 1,000, 2,000, and 3,000 bars, isopleths at XH\textsubscript{2}O\textsuperscript{fluid} values of 0, 0.25, 0.5, 0.75, and 1 with an open-system and a closed-system degassing path.}
\label{fig8}
\end{code}

\subparagraph{Isobars, isopleths, and degassing paths for multiple samples}\label{comparing-results-from-multiple-models}

First, we will calculate some new data for two different samples: a basanite \cite{iacovino2016} and a rhyolite \cite{myers2019}. For both samples we will calculate and then plot:

\begin{itemize}
\item Isobars and isopleths at 1100 °C, pressures of 1,000 and 2,000 bars and fluid compositions of XH\textsubscript{2}O\textsuperscript{fluid} of 0.25, 0.5, and 0.75
\item Closed-system degassing paths at 1100 °C
\end{itemize}

\begin{minted}[breaklines]{python}
basanite_sample = myfile.get_sample_composition('KI -07', asSampleClass=True)
rhyolite_sample = myfile.get_sample_composition('samp. P1968a', asSampleClass=True)

basanite_isobars, basanite_isopleths = v.calculate_isobars_and_isopleths(sample=basanite_sample,
                                                                         temperature=1100,
                                                                         pressure_list=[1000, 2000],
                                                                         isopleth_list=[0.25,0.75]).result

rhyolite_isobars, rhyolite_isopleths = v.calculate_isobars_and_isopleths(sample=rhyolite_sample,
                                                                         temperature=1100,
                                                                         pressure_list=[1000, 2000],
                                                                         isopleth_list=[0.25,0.75]).result

basanite_degassing_path = v.calculate_degassing_path(sample=basanite_sample,
                                                     temperature=1100).result

rhyolite_degassing_path = v.calculate_degassing_path(sample=rhyolite_sample,
                                                    temperature=1100).result
\end{minted}

\begin{minted}[breaklines]{text}
Calculating isobar at 1000 bars
 done.                                                                                                                           
Calculating isobar at 2000 bars
 done.                                                                                                                           
Done!
Calculating isobar at 1000 bars
 done.                                                                                                                           
Calculating isobar at 2000 bars
 done.                                                                                                                           
Done!
[====================] 100% Calculating degassing path...
[====================] 100% Calculating degassing path...
\end{minted}

\begin{code}[H]
\centering
\begin{minted}[breaklines]{python}
fig, ax = v.plot(isobars=[basanite_isobars, rhyolite_isobars],
                  isopleths=[basanite_isopleths, rhyolite_isopleths],
                  degassing_paths=[basanite_degassing_path, rhyolite_degassing_path],
                  isobar_labels=["Basanite", "Rhyolite"],
                  isopleth_labels=["Basanite", "Rhyolite"],
                  degassing_path_labels=["Basanite", "Rhyolite"])
v.show()
\end{minted}

\includegraphics[width=0.7\linewidth]{files/1b03bbb3e1578df84c274fde29f2dc70.png}
\caption[]{Example of plotting multiple calculations from multiple samples on the same plot. Note that the colors are automatically set to correspond to each sample for all plotted items (here, isobars, isopleths, and degassing paths). Samples, pressures, temperatures, XH\textsubscript{2}O\textsuperscript{fluid} values, and degassing path styles are defined above in this section. Manuscript default values are for a basanite (sample \texttt{KI-07}) and a rhyolite (sample \texttt{samp. P1968a}) at 1,100 °C, 1,000 and 2,000 bars, and XH\textsubscript{2}O\textsuperscript{fluid} = 0.25 and 0.75 and closed-system degassing.}
\label{fig9}
\end{code}

\subsubsection{Model hybridization (Advanced)}\label{code-hybridization}

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 \href{/s10-model-hybridization}{Supplemental Jupyter notebook S10}.

\subsubsection{Exporting data}\label{exporting-data}

Once batch calculations have been performed, they can be exported to an Excel or \acrshort{csv} file with the \texttt{save\_excel()} and \texttt{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 \texttt{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 \texttt{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 \acrshort{csv} file, each calculation will be saved to its own \acrshort{csv} file, and a file name for each of these is required.

Advanced users note that the \texttt{calculations} argument takes in any pandas DataFrame object, meaning this functionality is not limited to \acrshort{vesical}'s prescribed outputs. The \texttt{save\_excel()} and \texttt{save\_csv()} methods use the pandas \texttt{to\_excel} and \texttt{to\_csv} methods, however not all options are implemented here. If saving to a \acrshort{csv} file, any arguments that can be passed to pandas \texttt{to\_csv} method may be passed to \acrshort{vesical}'s \texttt{save\_csv()}.

\begin{description}
\item[\textbf{Method structures}] \texttt{save\_excel(filename, calculations, sheet\_name=None)}\texttt{save\_csv(filenames, calculations)}
\item[\texttt{save\_excel()} \textbf{Required inputs}] \texttt{filename} (Excel): Name of the file to create. The extension (.xlsx) should be included along with the name itself, all in quotes (e.g., \texttt{filename='myfile.xlsx'}).\texttt{calculations}: A list of variables containing calculated outputs from any of the core BatchFile functions: \texttt{calculate\_dissolved\_volatiles()}, \texttt{calculate\_equilibrium\_fluid\_comp()}, and \texttt{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., \texttt{calculations=[my\_calculation]}).
\item[\texttt{save\_excel()} \textbf{Optional inputs}] \texttt{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 \texttt{calculations}. Otherwise, a list of names for the sheets can be passed, with the names in quotes (e.g. \texttt{sheet\_name=['SaturationPressures'])}. `Original\_User\_data' will always be saved as the first sheet.
\item[\texttt{save\_csv()} \textbf{Required inputs}] \texttt{filenames} (\acrshort{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., \texttt{filenames=["file1.csv", "file2.csv"]}).\texttt{calculations}: same as for \texttt{save\_excel()}. Must be same length as \texttt{filenames}.
\item[Calculated outputs] An Excel or \acrshort{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).
\end{description}

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 \texttt{save\_excel()} to save the original data and a list of calculations given by the calculations argument to an Excel file.

\begin{minted}[breaklines]{python}
myfile.save_excel(filename='testsave.xlsx',
                      calculations=[dissolved, eqfluid, eqfluid_wtemps, satPs, satPs_wtemps],
                      sheet_names=['dissolved', 'eqfluid', 'eqfluid_wtemps', 'SaturationPs', 'SatPs_wtemps'])
\end{minted}

\begin{minted}[breaklines]{text}
Saved testsave.xlsx
\end{minted}

\paragraph{Saving data for re-import into \acrshort{vesical}}

In many cases, it may be preferable to compute large amounts of data using \acrshort{vesical} and then reimport them, either to perform more analysis or to plot the data. Likewise, a user may wish to compute data in \acrshort{vesical} and then send the results to a colleague, who can then re-import that data into \acrshort{vesical} directly. For this case, we suggest using python's pickle package (\href{https://wiki.python.org/moin/UsingPickle}{https://wiki.python.org/moin/UsingPickle}). Any python object, such as the results of a \acrshort{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 \texttt{dissolved} to a pickle file that we name ``dissolved.p''.

\begin{minted}[breaklines]{python}
import pickle

pickle.dump(dissolved, open("dissolved.p", "wb"))
\end{minted}

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

\begin{minted}[breaklines]{python}
import pickle

dissolved = pickle.load(open("dissolved.p", "rb"))
\end{minted}

\subsection{Discussion and Applications}

\subsubsection{Compositional 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\textsubscript{2}O and CO\textsubscript{2} in magmas \cite{papale2006, moore2008, ghiorso2015, wieser2021}. Thus, compositional variance must be accounted for in any study examining solubility in multiple samples. A key use case where \acrshort{vesical} can facilitate the adoption of this practice is in melt inclusion (\acrshort{mi})  studies; specifically, where a single suite of \acrshort{mi} with multiple melt compositions is examined using solubility models to interrogate magmatic degassing processes. Prior to the availability of \acrshort{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 \acrshort{mi}.

Using \acrshort{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 \cite{roggensack2001}. The compositional variation of these \acrshort{mi} \cite{Figure~\ref{fig10}}, 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\textsubscript{2}O-CO\textsubscript{2} solubilities, isobars were computed at 1200 °C and 3,000 bars for all samples using the MagmaSat model in \acrshort{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 Figure~\ref{fig10} and Figure~\ref{fig11} in blue and green, respectively (sample 41b*, IMS=0.81 and 36a*, IMS=0.66 wt\%\textsuperscript{2} at 3,000 bars). A composition representing the average of all \acrshort{mi} in the dataset is shown in orange (``Average Sample'', IMS=0.70 wt\%\textsuperscript{2} at 3,000 bars). A Jupyter notebook to reproduce these calculations is provided in the supplement (\href{/s8-cerronegro-isobar-comparison}{Supplementary Jupyter Notebook S8}).

\begin{figure}[!htbp]
\centering
\includegraphics[width=0.7\linewidth]{files/Figure10-9e087b45e3538e005f5b77c45ff5a72e.png}
\caption[]{Harker style diagrams illustrating the compositional range of \acrshort{mi}s from Cerro Negro volcano from \cite{roggensack2001}. The ``Average Sample'' plotted as an orange dot represents a fictitious sample, calculated as the average of all \acrshort{mi}s in the dataset. Sample 41b* and 36a* are the names of samples that produced isobars with maximum and minimum area under the curve, respectively (see text). Gray diamonds are all other data in the dataset.}
\label{fig10}
\end{figure}

\begin{figure}[!htbp]
\centering
\includegraphics[width=0.7\linewidth]{files/Figure11-22e97fe1c6c5b41d526100546673f693.png}
\caption[]{H\textsubscript{2}O-CO\textsubscript{2} diagram with isobars for \acrshort{mi} from Cerro Negro volcano \cite{roggensack2001} computed by \acrshort{vesical} using MagmaSat at 1200 °C, pressures of 500, 1000, 2000, 3000, and 4,000 bars. Curves shown are polynomials fitted to data computed by \acrshort{vesical}. Blue and green curves correspond to samples 41b* and 36a*, which produced isobars with maximum and minimum area under the curve, respectively. Orange isobars were those computed for a fictitious sample representing the average composition of the \acrshort{mi} dataset. Gray diamonds are all other data in the dataset.}
\label{fig11}
\end{figure}

At all pressures, the integrated mixed-volatile solubility across the Cerro Negro dataset varies as much as 10\% relative \cite{Figure~\ref{fig11}}. For these \acrshort{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 \acrshort{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 \acrshort{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~\ref{fig11}). Alternative plots to the commonly used H\textsubscript{2}O-CO\textsubscript{2} diagram are shown in Figure~\ref{fig12}, in which the same dataset is plotted in terms of computed saturation pressure (at 1200 °C calculated with \acrshort{vesical} using MagmaSat) versus dissolved H\textsubscript{2}O, dissolved CO\textsubscript{2}, and fluid composition (as XH\textsubscript{2}O\textsuperscript{fluid} calculated with \acrshort{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\textsubscript{2}O and CO\textsubscript{2} concentrations lie along expected degassing trends with much less scatter than the H\textsubscript{2}O-CO\textsubscript{2} plot. We can also see from this figure that the fluid composition during this eruption at Cerro Negro remained relatively constant at XH\textsubscript{2}O\textsuperscript{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 \cite{wieser2021}.

\begin{figure}[!htbp]
\centering
\includegraphics[width=0.7\linewidth]{files/Figure12-1bec7dd3b65dcaeffac31f4013fcea62.png}
\caption[]{Saturation pressure at 1200 °C calculated using \acrshort{vesical} with MagmaSat versus measured dissolved H\textsubscript{2}O and CO\textsubscript{2} concentrations and calculated fluid composition in Cerro Negro melt inclusions. These plots meaningfully illustrate degassing processes while avoiding issues associated with commonly used H\textsubscript{2}O-CO\textsubscript{2} diagrams, which occur with even minor compositional variation within a given dataset.}
\label{fig12}
\end{figure}

\subsubsection{Model Comparisons}

One of the possible workflows enabled through \acrshort{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~\ref{table12}. The use of \acrshort{vesical}'s \texttt{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:

\begin{code}[H]
\centering
\caption[]{Melt compositions used for modeling}
\label{table12}\begin{minted}[breaklines]{python}
model_comps = v.BatchFile("tables/Table_Model_Comps.xlsx")
model_comps.data
\end{minted}

\begin{minted}[breaklines]{text}
SiO2  TiO2  Al2O3  Fe2O3   FeO   MnO   MgO   CaO  Na2O   K2O  \
Alkali Basalt  49.00  1.27  19.7   3.74   5.33  0.17  4.82  8.85  4.23  1.00   
Rhyolite       77.19  0.06  12.8   0.00   0.94  0.00  0.03  0.53  3.98  4.65   

               P2O5   H2O   CO2  Cr2O3  NiO  CoO  
Alkali Basalt  0.37  4.51  0.25  0.0    0.0  0.0  
Rhyolite       0.00  0.26  0.05  0.0    0.0  0.0
\end{minted}
\end{code}

\begin{minted}[breaklines]{python}
alkbasalt = model_comps.get_sample_composition("Alkali Basalt", asSampleClass=True)
rhyolite = model_comps.get_sample_composition("Rhyolite", asSampleClass=True)

alkbasalt_isobars, alkbasalt_isopleths = v.calculate_isobars_and_isopleths(sample=alkbasalt,
                                                                           temperature=1200, pressure_list=[500, 1000, 2000],
                                                                           isopleth_list=[0.5],
                                                                           print_status=True).result

rhyolite_isobars, rhyolite_isopleths = v.calculate_isobars_and_isopleths(sample=rhyolite,
                                                                         temperature=800,
                                                                         pressure_list=[500, 1000, 2000],
                                                                         isopleth_list=[0.5]).result

Iac_alkbasalt_isobars, Iac_alkbasalt_isopleths = v.calculate_isobars_and_isopleths(sample=alkbasalt,
                                                                                   temperature=1200,
                                                                                   pressure_list=[500, 1000, 2000],
                                                                                   isopleth_list=[0.5],
                                                                                   model="IaconoMarziano").result

Dixon_alkbasalt_isobars, Dixon_alkbasalt_isopleths = v.calculate_isobars_and_isopleths(sample=alkbasalt,
                                                                                       temperature=1200,
                                                                                       pressure_list=[500, 1000, 2000],
                                                                                       isopleth_list=[0.5],
                                                                                       model="Dixon").result

Shish_alkbasalt_isobars, Shish_alkbasalt_isopleths = v.calculate_isobars_and_isopleths(sample=alkbasalt,
                                                                                       temperature=1200,
                                                                                       pressure_list=[500, 1000, 2000],
                                                                                       isopleth_list=[0.5],
                                                                                       model="ShishkinaIdealMixing").result

Liu_rhyolite_isobars, Liu_rhyolite_isopleths = v.calculate_isobars_and_isopleths(sample=rhyolite,
                                                                                 temperature=800,
                                                                                 pressure_list=[500, 1000, 2000],
                                                                                 isopleth_list=[0.5],
                                                                                 model="Liu").result
\end{minted}

\begin{minted}[breaklines]{text}
Calculating isobar at 500 bars
 done.                                                                                                                           
Calculating isobar at 1000 bars
 done.                                                                                                                           
Calculating isobar at 2000 bars
 done.                                                                                                                           
Done!
Calculating isobar at 500 bars
 done.                                                                                                                           
Calculating isobar at 1000 bars
 done.                                                                                                                           
Calculating isobar at 2000 bars
 done.                                                                                                                           
Done!
\end{minted}

\begin{code}[H]
\centering
\begin{minted}[breaklines]{python}
fig, ax = v.plot(isobars=[alkbasalt_isobars, Iac_alkbasalt_isobars, Dixon_alkbasalt_isobars, Shish_alkbasalt_isobars],
                 isobar_labels=["MagmaSat", "Iacono -Marziano", "Dixon", "Shishkina"])
fig, ax = v.plot(isobars=[rhyolite_isobars, Liu_rhyolite_isobars], isobar_labels=["MagmaSat", "Liu"])
v.show()
\end{minted}

\includegraphics[width=0.7\linewidth]{files/bf565d855d9d59e1370ad7d675446a48.png}

\includegraphics[width=0.7\linewidth]{files/39f644bfdb8f8a1f6df53a9ae256b055.png}
\caption[]{Isobars plotted for an alkali basalt (a) and rhyolite (b) with \acrshort{vesical} for five mixed-fluid solubility models. For alkali basalt, MagmaSat, Iacono-Marziano, Dixon, and ShishkinaIdealMixing were used to create isobars at 1,200 °C. For rhyolite, MagmaSat and Liu were used to create isobars at 800 °C.}
\label{fig13}
\end{code}

It is immediately clear from Figure~\ref{fig10} that major disagreement exists between these models. For the alkali basalt, MagmaSat and Dixon show the best agreement, particularly at pressures \textless 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\textsubscript{2} solubilities at XH\textsubscript{2}O\textsuperscript{fluid} values approaching 0.

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

The models of MagmaSat and Liu show a similar level of disagreement for H\textsubscript{2}O-CO\textsubscript{2} solubility in the rhyolite, with Liu predicting much higher dissolved CO\textsubscript{2} concentrations at low XH\textsubscript{2}O\textsuperscript{fluid} (\textless 20 relative\% or {\textasciitilde}220 ppm at XH\textsubscript{2}O\textsuperscript{fluid}=0.1).

\subsubsection{Sensitivity 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, \acrshort{vesical} allows users to import an Excel or \acrshort{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, \acrshort{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\textsubscript{2} and H\textsubscript{2}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 \acrshort{vesical} also allows the user to combine the code with other Python3 modules. For example, users could utilize Markov chain Monte Carlo (\acrshort{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 \acrshort{vesical}, building on existing tools within the python community.

\subsubsection{Future development}

\acrshort{vesical} represents the first comprehensive volatile solubility modeling tool of its kind, including the feature that \acrshort{vesical} is extensible. \acrshort{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 \acrshort{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 \acrshort{vesical} are very welcome. Users can contribute to \acrshort{vesical}'s code, implementing new models and new features, via github (\href{https://github.com/kaylai/VESIcal}{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 \acrshort{vesical} by creating a ``Pull Request'' inside of github. Changes and enhancements to \acrshort{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: \cite{https://doi.org/10.5281/zenodo.5095382}). Planned features not implemented in this release include: 1. Models to calculate sample oxygen fugacity from Fe\textsuperscript{2+}/$\Sigma$Fe and vice versa; 2. Additional volatiles such as sulfur; 3. More thermodynamic solubility models such as that of \cite{papale2006}; 4. Sensitivity and error analysis functions.

\subsubsection{How to cite \acrshort{vesical} and its models}

To cite computations done using \acrshort{vesical}, please cite this manuscript, the \acrshort{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 \cite{ghiorso2015}''. For example: ``Calculations were performed using \acrshort{vesical} \cite{iacovino2021} with the models of \cite{shishkina2014} and \cite{dixon1997}.'' The web-app always runs on the most up-to-date version of the \acrshort{vesical} code, but it is best practice to note if the web-app was used (``Calculations were performed using the \acrshort{vesical} web-app \cite{iacovino2021}...''). 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, \acrshort{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.

\subsection{Conclusions}

\acrshort{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 \acrshort{vesical}, which employs the most diversely calibrated (chemically and in P-T space) of the group, MagmaSat \cite{ghiorso2015}, as the default model. \acrshort{vesical} can perform five core calculations with any mixed-fluid model and three core calculations with any model (mixed-fluid, CO\textsubscript{2}-only, H\textsubscript{2}O-only). \acrshort{vesical} allows for automatic calculation of large datasets and robust built-in plotting capability.

Alongside model frameworks such as ENKI, \acrshort{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 \acrshort{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.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%  acronyms & glossary  %%%%%%%%%%%%%
\printglossaries
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Open Research
% AGU requires an Availability Statement for the underlying data needed to understand, evaluate, and build upon the reported research at the time of peer review and publication.

%Additionally, authors should include an Availability Statement for the software that has a significant impact on the research. Details and templates are in the Availability Statement section of the Data & Software for Authors Guidance:
% https://www.agu.org/Publish-with-AGU/Publish/Author-Resources/Data-and-Software-for-Authors#availability

%For physical samples, use the IGSN persistent identifier, see the International Geo Sample Numbers section:
% https://www.agu.org/Publish-with-AGU/Publish/Author-Resources/Data-and-Software-for-Authors#IGSN

\section*{Open Research}
The VESIcal software is open source and is hosted on github (\href{https://github.com/kaylai/VESIcal}{https://github.com/kaylai/VESIcal}). The version of VESIcal used in this manuscript is version 1.0.1 and is archived on Zenodo (\cite{https://doi.org/10.5281/zenodo.5095382}). VESIcal runs on top of thermoengine, a python package that is a part of the ENKI framework (\href{http://enki-portal.org/}{http://enki-portal.org/}). The thermoengine library is open source and is available on GitLab (\href{https://gitlab.com/ENKI-portal/ThermoEngine}{https://gitlab.com/ENKI-portal/ThermoEngine}). VESIcal was written in Python3 and should be stable up to at least Python version 3.7.6. In addition to thermoengine, VESIcal requires the following standard libraries (with versions used for testing indicated in brackets): pandas (1.0.1), numpy (1.18.1), matplotlib (3.1.2), cycler (0.10.0), scipy (1.4.1), and sympy (1.5.1). The VESIcal webapp interface runs through Anvil (anvil.works), which executes VESIcal code on a cloud server. The code that facilitates the link between the anvil interface and the VESIcal code is available on the VESIcal github. VESIcal can also be used within a Jupyter Notebook and is hosted on the ENKI JupyterHub (\href{https://server.enki-portal.org/hub/login}{https://server.enki-portal.org/hub/login}) such that the code can be accessed without installation on the user's local machine.

All data sets used in this manuscript are available on the VESIcal github as well as in the Supplementary Material of this manuscript. The example dataset used for worked examples in Section~\ref{example-uses} (\href{/example\_data-9296285330fa856e8b827784ac193203.xlsx}{\texttt{example\_data.xlsx}} file; Supplemental Dataset S1) contains compositional information for basalts \cite{tucker2019, roggensack2001}, andesites \cite{moore1998}, rhyolites \cite{mercer2015, myers2019}, and alkaline melts (phototephrite, basaltic-trachyandesite, and basanite from \cite{iacovino2016}). Several additional example datasets from the literature are available in the Supplement (Supplementary Datasets S2-S5; Table~\ref{table4}). These include experimentally produced alkaline magmas from \cite{iacovino2016} (\href{/alkaline-ab8dc69843f9db52d73ba16705e49ae7.xlsx}{\texttt{alkaline.xlsx}}), basaltic melt inclusions from Kilauea \cite{tucker2019} and Gakkel Ridge \cite{bennett2019} (\href{/basalts-8053fbc358b9099eeda22efa1ad28e1f.xlsx}{\texttt{basalts.xlsx}}), basaltic melt inclusions from Cerro Negro volcano, Nicaragua \cite{roggensack2001} (\href{/cerro\_negro-7895112c5be9cd0244bc4c972a08b834.xlsx}{\texttt{cerro\_negro.xlsx}}), and rhyolite melt inclusions from the Taupo Volcanic Center, New Zealand \cite{myers2019} and a topaz rhyolite from the Rio Grande Rift \cite{mercer2015} (\href{/rhyolites-ab31a74de01c072d655c95793ca0991e.xlsx}{\texttt{rhyolites.xlsx}}). Where available, the calibration datasets for VESIcal models are also provided (Supplementary Datasets S6-S7).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



\acknowledgments
% This section is optional. Include any Acknowledgments here.
This manuscript is dedicated to the memory of Dr. Peter Fox without whom none of this work would have been possible. We thank Peter for his encouragement of this work, his editorial handling of the manuscript, and for blazing a path for bringing executable manuscripts to AGU journals. KI and GMM were supported by the NASA Jacobs JETS Contract \#NNJ13HA01C. PEW acknowledges support from a NERC DTP studentship (NE/L002507/1). The authors thank Jackie Dixon and Bob Myhill for reviews, which greatly helped strengthen this manuscript and the VESIcal code. The authors would also like to thank Mark Ghiorso, Aaron Wolf, and the ENKI team for pushing thermodynamic modeling into the future and for making this publication possible; Chelsea Allison and Giada Iacono-Marziano for discussions on their published models and how to properly implement them in VESIcal; Christy B. Till for support of KI during early coding work with MagmaSat; and presentationgo for style elements used in flowchart figures. Permission for the use of the VESIcal fox logo was graciously provided by DeviantArt user Twai.

%% ------------------------------------------------------------------------ %%
%% References and Citations

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% \bibliography{<name of your .bib file>} don't specify the file extension
%
% don't specify bibliographystyle
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliography{main.bib}

%Reference citation instructions and examples:
%
% Please use ONLY \cite and \citeA for reference citations.
% \cite for parenthetical references
% ...as shown in recent studies (Simpson et al., 2019)
% \citeA for in-text citations
% ...Simpson et al. (2019) have shown...
%
%
%...as shown by \citeA{jskilby}.
%...as shown by \citeA{lewin76}, \citeA{carson86}, \citeA{bartoldy02}, and \citeA{rinaldi03}.
%...has been shown \cite{jskilbye}.
%...has been shown \cite{lewin76,carson86,bartoldy02,rinaldi03}.
%... \cite <i.e.>[]{lewin76,carson86,bartoldy02,rinaldi03}.
%...has been shown by \cite <e.g.,>[and others]{lewin76}.
%
% apacite uses < > for prenotes and [ ] for postnotes
% DO NOT use other cite commands (e.g., \citet, \citep, \citeyear, \citealp, etc.).
% \nocite is okay to use to add references from your Supporting Information
%

\end{document}
