diff --git a/docs/Makefile b/docs/Makefile new file mode 100644 index 0000000..d4bb2cb --- /dev/null +++ b/docs/Makefile @@ -0,0 +1,20 @@ +# Minimal makefile for Sphinx documentation +# + +# You can set these variables from the command line, and also +# from the environment for the first two. +SPHINXOPTS ?= +SPHINXBUILD ?= sphinx-build +SOURCEDIR = . +BUILDDIR = _build + +# Put it first so that "make" without argument is like "make help". +help: + @$(SPHINXBUILD) -M help "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) + +.PHONY: help Makefile + +# Catch-all target: route all unknown targets to Sphinx using the new +# "make mode" option. $(O) is meant as a shortcut for $(SPHINXOPTS). +%: Makefile + @$(SPHINXBUILD) -M $@ "$(SOURCEDIR)" "$(BUILDDIR)" $(SPHINXOPTS) $(O) diff --git a/docs/_static/custom.css b/docs/_static/custom.css new file mode 100644 index 0000000..c7d1eb1 --- /dev/null +++ b/docs/_static/custom.css @@ -0,0 +1,6 @@ +.math { + text-align: left; +} +.eqno { + float: right; +} \ No newline at end of file diff --git a/docs/conf.py b/docs/conf.py new file mode 100644 index 0000000..48008e0 --- /dev/null +++ b/docs/conf.py @@ -0,0 +1,32 @@ +# Configuration file for the Sphinx documentation builder. +# +# For the full list of built-in configuration values, see the documentation: +# https://www.sphinx-doc.org/en/master/usage/configuration.html + +# -- Project information ----------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#project-information + +project = 'evo-outgas' +copyright = '2022, Philippa Liggins' +author = 'Philippa Liggins' +release = '1.0.1' + +# -- General configuration --------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#general-configuration + +extensions = ['sphinx.ext.todo', 'sphinx.ext.viewcode', 'sphinx.ext.autodoc', 'sphinx.ext.napoleon', 'sphinx.ext.mathjax'] + +templates_path = ['_templates'] +exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] + +numfig = True + + +# -- Options for HTML output ------------------------------------------------- +# https://www.sphinx-doc.org/en/master/usage/configuration.html#options-for-html-output + +html_theme = 'sphinx_rtd_theme' +html_static_path = ['_static'] + +def setup(app): + app.add_css_file('custom.css') diff --git a/docs/evo_doc.rst b/docs/evo_doc.rst new file mode 100644 index 0000000..f827a45 --- /dev/null +++ b/docs/evo_doc.rst @@ -0,0 +1,1072 @@ +.. role:: raw-latex(raw) + :format: latex +.. + +.. _`c:models`: + +EVo Model Description +===================== + +EVo is designed to model the isothermal ascent of magma during eruption, +assuming that magma degassing during ascent is primarily due to the +resultant pressure change. A 3-phase system is used, where the silicate +melt and the gas phase coexist in constant thermochemical equilibrium, +alongside a solid graphite phase under sufficiently reduced conditions. +The physical-chemical separation of gas and melt is implemented as being +either impossible (closed-system), or instantaneous (open system - see +:ref:`section_run_types`). The highly complex nature of volcanic +systems necessitates simplifications in my approach; crystallisation +kinetics and sulfide formation (see :ref:`section_cands`) are +neglected, so in those cases where precipitation of solids during ascent +is important in driving degassing, EVo may underestimate the volume and +composition of gases evolved from the melt. + +.. _evo_schem: +.. figure:: images/evo_overview.* + + A schematic overview of EVo, and it’s place in a physical system. + +:numref:`evo_schem` shows a schematic of how EVo represents a volcanic system, calculating +the composition of the gas and melt phase at stages from the pressure of +volatile saturation to the surface. + +.. _section_overview_of_model: + +Overview of the chemical model +------------------------------ + +EVo computes the gas and melt volatile composition for systems of up to +5 volatile elements (C, O, H, S and N, in combinations of OH, COH, SOH, +COHS and COHSN), using an ‘equilibrium constants and mass balance’ +method first established by Holloway (1987) and used since by +several degassing models, including Iacono-Marziano et al. (2012), +Gaillard and Scaillet (2014), and Burgisser et al. (2015). +As stated by Holloway (1987), this method "requires that enough equilibrium +constant equations are used so that, +when combined with the equation for mass balance, there are an equal +number of equations and unknowns". By combining the equilibrium constant +equations which define homogenous gas-gas equilibria, and solubility +laws which control heterogeneous gas-melt equilibria, at any given point +the entire gas-melt system can be defined. + +EVo considers 10 gas-phase species in the full COHSN system: |H2O|, |H2|, |O2|, |CO2|, CO, |CH4|, |S2|, |SO2|, |H2S| and |N2|. Subsets of these species are used in calculations if systems +with fewer components are modelled, e.g., by definition the COH system +does not include S or N-bearing species and so comprises just 6 +gas-phase species. In a typical terrestrial volcanic gas, :math:`>`\ 99% +(by mol) will be made up of |H2O|, |H2|, |CO2|, |SO2|, |H2S| and |N2| (Fischer and +Chiodini, 2015), so the 10 species EVo considers are sufficient to model the +bulk gas-phase chemistry, although it omits some trace species such as which +have been detected at various volcanoes (e.g., Mori and Notsu, 1997; +Oppenheimer and Kyle, 2008; Sawyer et al., 2008). +OCS in particular has been omitted due to the lack of solubility data in +magmas. This reflects a limitation of this approach: even when the +homogenous gas-phase equilibria are known (which for very many species +they are), unless there is corresponding data for the solubility of a +species in silicate liquids it cannot be included and a self-consistent +solution be obtained. + +A set of five homogenous gas-phase equilibria, described in :eq:`eq_k1` +to :eq:`eq_k5`, define the speciation of the gas phase + +.. math:: \ce{H2 + \frac{1}{2} O2 <=> H2O} + :label: eq_k1 + +.. math:: \ce{CO + \frac{1}{2} O2 <=> CO2} + :label: eq_k2 + +.. math:: \ce{CH4 + 2O2 <=> CO2 + 2|H2O|} + :label: eq_k3 + +.. math:: \ce{\frac{1}{2} S2 + H2O <=> H2S + \frac{1}{2} O2} + :label: eq_k4 + +.. math:: \ce{\frac{1}{2} S2 + O2 <=> SO2}. + :label: eq_k5 + +At thermochemical equilibrium, the abundance of each species in the gas +phase are related to an equilibrium constant, :math:`K_i`, one for each +reaction in :eq:`eq_k1` to :eq:`eq_k5`. The general equation for these +equilibria is given by + +.. math:: + :label: eq_base_eqbm + + a\ce{A} + b\ce{B} \ce{<=>} c\ce{C} + d\ce{D}, + +which has a general equilibrium constant reaction equation of + +.. math:: + :label: eq_k_base_eqbm + + K_{2.6} = \frac{[\ce{C}]^c\,[\ce{D}]^d}{[\ce{A}]^a\,[\ce{B}]^b}, + +where square brackets denote activity. For real gases, :math:`K_i` can +be calculated in terms of fugacity, :math:`f_i`, so that for, e.g., +Eq. :eq:`eq_k1` + +.. math:: + :label: eq_k1_ex + + K_1 = \frac{f_{\ce{|H2O|}}}{f_{\ce{H2}}\,{f_{\ce{O2}}}^{0.5}}. + +The fugacity :math:`f_i` of a real gas is an effective partial pressure, +which replaces the mechanical partial pressure in an accurate +computation of the chemical equilibrium constant. In EVo ideal mixing of +real gases is assumed, so the fugacity of species :math:`i` is defined +as + +.. math:: + :label: eq_fugacity + + f_i = \gamma_i\,P_i = \gamma_i\,X^v_i\,P, + +where :math:`\gamma_i` is the fugacity coefficient representing the +deviation of the gas from ideality, so for an ideal gas :math:`\gamma_i` += 1. :math:`P_i` is the partial pressure of :math:`i`, :math:`x^v_i` is +the mole fraction of i in the gas mixture, and :math:`P` is the total +pressure. :math:`x^v_i` is such that when the gas phase is composed of +:math:`n` species, each with a mole fraction :math:`x^v_i`, + +.. math:: + :label: eq_sumtoone + + \sum^{n}_{i=1}\,X^v_i = 1. + +The Lewis and Randall rule is applied to the calculation of +:math:`\gamma_i`, which states that the fugacity coefficient of species +i in the gas mixture equals that of the pure species at the same +pressure and temperature. The coefficients :math:`\gamma_{\ce{|H2O|}}`, +:math:`\gamma_{\ce{CO2}}` and :math:`\gamma_{\ce{N2}}` are from +Holland and Powell (1991), +:math:`\gamma_{\ce{H2}}` is from Shaw and Wones (1964), and the remaining +coefficients are from Shi and Saxena (1992). Coefficients taken from Shi and Saxena +(1992) are only calibrated down to 1 bar. These coefficients are therefore set as +:math:`\gamma_i` = 1 for :math:`P\,\leq\,1` – as gas behaviour tends to +be increasingly ideal at lower pressures this is a valid assumption. + +Calculating reaction constants +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The equilibrium constants of reactions :eq:`eq_k1` to :eq:`eq_k5` are +calculated using the Gibbs free energy of the reaction as + +.. math:: \ln K = \frac{\Delta_rG^{\circ}}{-R\,T}. + :label: eq_calc_k + +:math:`\Delta_rG^{\circ}`, the Gibbs free energy of reaction, is the +difference between the Gibbs free energies of formation +(:math:`\Delta_fG^{\circ}`) for products and reactants *in their +standard states*, thus: + +.. math:: \Delta_rG^{\circ} = \Delta_fG^{\circ}_{products} - \Delta_fG^{\circ}_{reactants}. + :label: eq_gibbs_reaction + +Therefore to calculate :math:`\Delta_rG^{\circ}_1`, the Gibbs free +energy of reaction :eq:`eq_k1`: + +.. math:: \ce{H2 + \frac{1}{2} O2 <=> |H2O|}, + +.. math:: \Delta_rG^{\circ}_1 = \Delta_fG^{\circ}_{\ce{|H2O|}} - \left(\Delta_fG^{\circ}_{\ce{H2}} + 0.5 \Delta_fG^{\circ}_{\ce{O2}}\right). + :label: eq_gibbs_1 + +:math:`\Delta_fG^{\circ}` values are taken from the JANAF tables +(Chase, 1998), with linear interpolation between temperatures using +the NumPy ‘interp’ function (Harris et al., 2020). As +equilibrium constants are pressure-independent, and EVo assumes +decompression is isothermal, :math:`K` values are only computed once, +during the initial setup of a decompression calculation. + +.. _`section:volatile_solubility`: + +Volatile solubility and oxygen exchange +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The solubility of a species is usually defined as the maximum +concentration of a volatile which can remain in solution while +co-existing with a pure gas phase of the same species. In a +multicomponent system like that of EVo, it is assumed that the +solubility of a species is related to the partial pressure/fugacity of +the individual species :math:`i`. Over the past 30 years a plethora of +solubility laws, particularly for |H2O| and |CO2|, have been developed. Different +models for single species, or occasionally mixed |H2O|-|CO2| volatiles are still +being frequently introduced to the community; sometimes improvements +based on new data, sometimes for a new P/T, :math:`f_{\ce{O2}}`, or silicate compositional +regime. However, all these models are at least partially dependent on +the abundance of the relevant species in the gas phase. In order to be +as flexible as possible to different volcanic scenarios, and robust to +future developments, EVo has been designed such that the user can select +from a number of different solubility laws for various volatile species; +it is also reasonably easy to add in additional solubility laws to a +single file as they are released. In :numref:`table_solubility_laws` below, the +solubility laws for each species currently implemented in EVo are +listed, along with the compositional, temperature, and pressure ranges +they are calibrated for. + +.. _table_solubility_laws: +.. list-table:: : The solubility laws implemented in EVo. The values listed here are strictly those used to calibrate the model – e.g., when using the Volatile-Calc spreadsheet temperatures of 600-1500 °C are allowed despite the solubility laws only being calibrated at 1200 °C. DC = D-Compress; VC-B = VolatileCalc Basalt simplification. a finds the sulfide capacity or dissolved sulfide/sulfate ratio in the melt, as discussed in :ref:`section_cands`. + :class: borderless + + * - .. image:: images/table1.png + +As is also discussed in Iacovino et al. (2021) and Wieser et al. (2022) with respect to +the development of VESIcal (a Python3 engine for using and comparing a +range of and solubility models), implementing the solubility laws which +are published in academic papers is often problematic, particularly if +the original work has not provided accompanying tools to benchmark +against, given that some original manuscripts contain typos or +formatting errors. Each of the solubility laws built into EVo are +therefore provided in the API for the `solubility_laws.py` file, both for finding the abundance of +a volatile in the melt based on gas phase fugacity, and for finding the +reverse – the gas phase fugacity based on the concentration in the melt. + +Both H and N form relatively simple binary volatile systems, present +either in the gas phase, or dissolved in the silicate melt. However, C +and S can both form additional phases: graphite in the case of carbon, +and sulfide or anhydrite phases in the case of sulfur. In the case where +graphite saturation occurs, |CO2| solubility is forced to be calculated with +Eguchi and Dasgupta (2018), which can account for +both graphite and fluid saturation within the same equation. Sulfur +solubility, and the way EVo handles both C and S saturation is covered +in :ref:`section_cands`. + +EVo takes into account the way that iron dissolved in the melt affects +the redox state of the magma by modelling the exchange of oxygen, +according to + +.. math:: \ce{Fe2O3_{(melt)} <=> 2FeO_{(melt)} + 0.5O2_{(gas)}}. + :label: eq_iron_buffer + +This reaction is not calculated using an equilibrium constant, but the +ratio of the molar fractions of and : +:math:`F = X_{\ce{FeO}}/X_{\ce{Fe2O3}}`. This can be calculated using +Kress and Carmichael (1991), or in the +case of iron rich magma (total Fe :math:`>`\ 15 wt%) the model of +Righter et al. (2013) which is calibrated for Martian melts. This exchange of oxygen between the +silicate melt and the gas phase means that redox equilibrium is always +maintained in the system. + +Mass balance +~~~~~~~~~~~~ + +As mentioned at the start of this section, the ‘mass balance and +equilibrium constant’ method requires a set of mass balance equations +along with the equilibrium constants already described. The total mass +of each species in the magmatic system (:math:`W_{Ti}`) is the sum its +exsolved and dissolved components, thus + +.. math:: W_{Ti} = W_{gT}\,w^v_i + w^m_i. + :label: eq_gasplusmelt + +:math:`W_{gT}` is the total weight fraction of gas in the system, +:math:`w^v_i` is the weight fraction of species i in the gas phase, and +:math:`w^m_i` is the weight fraction of species :math:`i` in the melt. +:math:`w^m_i` is calculated according to a subset of the solubility laws +listed in :numref:`table_solubility_laws`, and as +such their dependencies differ according to the law and species in +question; however, they are always dependent on the current pressure, +and the fugacity of the corresponding species. + +Mass conservation of volatile elements throughout the degassing process +is enforced by keeping the total weight fraction of each element +(:math:`W_{\mathrm{TH}}`, :math:`W_{\mathrm{TC}}`, +:math:`W_{\mathrm{TO}}`, :math:`W_{\mathrm{TS}}`, +:math:`W_{\mathrm{TN}}`) constant + +.. math:: + :label: eq_conserve_o + + \frac{W_{TO}}{M_{\ce{O}}} = 2\frac{W_{T\ce{O2}}}{M_{\ce{O2}}} + 2\frac{W_{T\ce{H2O}}}{M_{\ce{H2O}}} +2\frac{W_{T\ce{SO2}}}{M_{\ce{SO2}}} + 2\frac{W_{T\ce{CO2}}}{M_{\ce{CO2}}} + \frac{W_{T\ce{CO}}}{M_{\ce{CO}}} + \frac{W_{O(Fe)}}{M_{\ce{O}}} + +.. math:: + :label: eq_conserve_h + + \frac{W_{TH}}{2M_{\ce{H}}} = \frac{W_{T\ce{H2}}}{M_{\ce{H2}}} + \frac{W_{T\ce{H2O}}}{M_{\ce{H2O}}} + \frac{W_{T\ce{H2S}}}{M_{\ce{H2S}}} + 2\frac{W_{T\ce{CH4}}}{M_{\ce{CH4}}} + +.. math:: + :label: eq_conserve_c + + \frac{W_{TC}}{M_{\ce{C}}} = \frac{W_{T\ce{CO2}}}{M_{\ce{CO2}}} + \frac{W_{TCO}}{M_{CO}} + \frac{W_{T\ce{CH4}}}{M_{\ce{CH4}}} + \frac{W_{C(graph)}}{M_{\ce{C}}} + +.. math:: + :label: eq_conserve_s + + \frac{W_{TS}}{M_{\ce{S}}} = \frac{W_{T\ce{H2S}}}{M_{\ce{H2S}}} + \frac{W_{T\ce{SO2}}}{M_{\ce{SO2}}} + 2\frac{W_{T\ce{S2}}}{M_{\ce{S2}}} + +.. math:: \frac{W_{TN}}{M_{\ce{N}}} = 2\frac{W_{T\ce{N2}}}{M_{\ce{N2}}} + :label: eq_conserve_n + +where :math:`W_Ti` for each volatile species (not element) is calculated +using :eq:`eq_gasplusmelt`, :math:`M_i` is the +molecular mass of species :math:`i`, :math:`W_{O(Fe)}` is the weight +fraction of oxygen held in iron (FeO/Fe\ :sub:`2`\ O\ :sub:`3`\) within the melt, and +:math:`W_{C(graph)}` is the mass of graphite in the system (see :ref:`section_worked_example`). + +In order to solve for the state of the magmatic system at a given +temperature and pressure, equations for chemical equilibrium :eq:`eq_k1` +to :eq:`eq_k5` and :eq:`eq_iron_buffer`, +solubility laws (:numref:`table_solubility_laws`) +and mass conservation :eq:`eq_sumtoone` and :eq:`eq_conserve_o` to :eq:`eq_conserve_n` are used jointly to algebraically reduce the +system down to the smallest number of equations. In the case of the +COHSN system, this is 4, which requires four variables to be solved for +(:math:`X_{CO}`, :math:`X_{\ce{S2}}`, :math:`X_{\ce{O2}}` and +:math:`X_{\ce{N2}}`). An example of how these equations are derived and +then solved is shown for the COHSN system in :ref:`section_worked_example`. It is the case that the number of +equations, and therefore variables, is always :math:`n-1` for :math:`n` +number of volatile elements involved. These reduced equations are solved +simultaneously, using an iterative method from the SciPy python package +(Virtanen et al., 2020), conserving +the mass fraction of each volatile element to a precision of +1 × 10\ :sup:`−5`. + +Once a solution has been found for the speciation of the system at a +single pressure, the volume of the gas phase (:math:`v`) is calculated +using + +.. math:: v = \left[1 + \frac{\mu P(1-W_{gT})}{RT\rho_lW_{gT}}\right]^{-1}, + :label: eq_gas_frac + +:math:`R` is the universal gas constant (8.3144 J/mol K), :math:`\rho_l` +is the volatile-free melt density, calculated according to +Lesher and Spera (2015), and +:math:`\mu` is the mean molecular weight of the gas phase, calculated as + +.. math:: \mu = \sum^{10}_{i=1} X_i M_i. + :label: eq_mean_mol_mass + +.. _section_worked_example: + +Example of solving equations for the COHSN system +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +This section shows a worked example of how the full system of equations +is algebraically reduced down to the four equations iterated to solve +the system. In the COHSN system, the gas phase is made up of 10 species, +in equilibrium at all times according to + +.. math:: + + \begin{aligned} + K_1 &= \frac{f_{\ce{H2O}}}{f_{\ce{H2}} \cdot \sqrt{f_{\ce{O2}}}} & K_2 &= \frac{f_{\ce{CO2}}}{f_{\ce{CO}} \cdot \sqrt{f_{\ce{O2}}}} \\ + K_3 &= \frac{f_{\ce{CO2}} \cdot f_{\ce{H2O}}^2}{f_{\ce{CH4}}\cdot f_{\ce{O2}}^2} & K_4 &= \frac{f_{\ce{H2S}} \cdot \sqrt{f_{\ce{O2}}}}{f_{\ce{H2O}} \cdot \sqrt{f_{\ce{S2}}}}\\ + K_5 &= \frac{f_{\ce{SO2}}}{f_{\ce{O2}}\,\sqrt{f_{\ce{S2}}}}, + \end{aligned} + +where :math:`f_i` is calculated according to :eq:`eq_fugacity`. +The variables to solve for are +:math:`X_{CO}`, :math:`X_{\ce{S2}}`, :math:`X_{\ce{O2}}` and +:math:`X_{\ce{N2}}`, so the first aim is to calculate the mole fractions +of all other species in terms of these 4 variables (written in blue, +with derived species in red): + +.. math:: + + \begin{aligned} + \textcolor{red}{X_{\ce{CO2}}} &= \frac{K_2 \, \gamma_{CO} \, \textcolor{blue}{X_{CO}} \, \sqrt{\gamma_{\ce{O2}}\,\textcolor{blue}{X_{\ce{O2}}}\,P}}{\gamma_{\ce{CO2}}} + & \textcolor{red}{X_{\ce{SO2}}} &= \frac{K_5 \, \gamma_{\ce{O2}} \, \textcolor{blue}{X_{\ce{O2}}} \, \sqrt{\gamma_{\ce{S2}}\,\textcolor{blue}{X_{\ce{S2}}}\,P}}{\gamma_{\ce{SO2}}} \\ + a &= \frac{\gamma_{\ce{CO2}}\,\textcolor{red}{X_{\ce{CO2}}}\,\gamma_{\ce{H2O}}^2}{K_3\,\gamma_{\ce{CH4}}\,\sqrt{\gamma_{\ce{O2}}\textcolor{blue}{X_{\ce{O2}}}}} & \\ + b &= \frac{\gamma_{\ce{H2O}}}{K_1\,\gamma_{\ce{H2}}\,\sqrt{\gamma_{\ce{O2}}\textcolor{blue}{X_{\ce{O2}}}}} + \frac{K_4\,\gamma_{\ce{|H2O|}}\sqrt{\gamma_{\ce{S2}}\,X_{\ce{S2}}}}{\gamma_{\ce{H2S}}\,\sqrt{\gamma_{\ce{O2}}\,\textcolor{blue}{X_{\ce{O2}}}}} + 1 & \\ + c &= -(1 - \textcolor{blue}{X_{\ce{CO}}} - \textcolor{blue}{X_{\ce{S2}}} - \textcolor{blue}{X_{\ce{O2}}} - \textcolor{blue}{X_{\ce{N2}}} - \textcolor{red}{X_{\ce{CO2}}} - \textcolor{red}{X_{\ce{SO2}}}) & \\ + \textcolor{red}{X_{\ce{H2O}}} &= \frac{-b + \sqrt{b^2 - (4ac)}}{2a} & \\ + \end{aligned} + +Note that :math:`X_{\ce{H2O}}` is calculated by rearranging equations +for :math:`K_1`, :math:`K_3` and :math:`K_4` in terms of +:math:`f_{\ce{H2}}`, :math:`f_{\ce{CH4}}` and :math:`f_{\ce{H2S}}` +respectively, and substituting these equations into :eq:`eq_sumtoone` such that + +.. math:: X_{\ce{H2}} + X_{\ce{CH4}} + X_{\ce{H2S}} = 1-(\textcolor{blue}{X_{\ce{O2}}} + \textcolor{blue}{X_{CO}} + \textcolor{blue}{X_{\ce{S2}}} + \textcolor{blue}{X_{\ce{N2}}} + \textcolor{red}{X_{\ce{CO2}}} + \textcolor{red}{X_{\ce{SO2}}}). + +This is re-arranged into quadratic form for :math:`X_{\ce{H2O}}`, and +solved using the quadratic formula. + +.. math:: + + \begin{aligned} + X_{\ce{H2}} &= \frac{\gamma_{\ce{H2O}}\,\textcolor{red}{X_{\ce{H2O}}}}{K_1\,\gamma_{\ce{H2}}\,\sqrt{\gamma_{\ce{O2}}\,\textcolor{blue}{X_{\ce{O2}}}\,P}}\\ + X_{\ce{CH4}} &= \frac{\gamma_{\ce{CO2}}\,\textcolor{red}{X_{\ce{CO2}}}\,\sqrt{\gamma_{\ce{H2O}}\,\textcolor{red}{X_{\ce{H2O}}}}}{K_3\,\gamma_{\ce{CH4}}\,(\gamma_{\ce{O2}}\,\textcolor{blue}{X_{\ce{O2}}})^2}\\ + X_{\ce{H2S}} &= \frac{K_4\,\gamma_{\ce{H2O}}\,\textcolor{red}{X_{\ce{H2O}}}\,\sqrt{\gamma_{\ce{S2}}\,\textcolor{blue}{X_{\ce{S2}}}}}{\gamma_{\ce{H2S}}\,\sqrt{\gamma_{\ce{O2}}\,\textcolor{blue}{X_{\ce{O2}}}}} + \end{aligned} + +At this point, all 10 species in the gas phase can be described in terms +of the 4 variables being solved for. As decompression occurs, mass +conservation requires that the total amounts of each volatile element +(COHSN) in the system remains constant. These mass conservation +equations are set out in :eq:`eq_conserve_o` to :eq:`eq_conserve_n`. In order to use the molar fractions derived +above in :eq:`eq_gasplusmelt`, to replace the total +amounts of each species :math:`W_{Ti}` in :eq:`eq_conserve_o` to :eq:`eq_conserve_n`, the molar fraction must be converted to weight fractions using + +.. math:: w_i = \frac{X_iM_i}{\sum_j X_jM_j}, + :label: eq_mol2wt + +and the mass of the volatile dissolved in the melt must be calculated. +Substituting :eq:`eq_mol2wt` and solubility laws into :eq:`eq_gasplusmelt` gives + +.. math:: W_{Ti} = W_{gT}\,\frac{X_iM_i}{\sum^{}_{j} X_jM_j} + w^m_i(X_i,P,T...), + :label: eq_expanded_sum + +where :math:`w^m_i(\dots)` denotes the solubility law for species +:math:`i`, which is dependent on :math:`X_i`, the pressure and +temperature of the system, and in some cases other variables such as the :math:`f_{\ce{O2}}` +or silicate melt chemistry. :math:`w^m_i(\dots)` returns the weight +fraction of the species in the melt. As the gas weight +fraction :math:`W_{gT}` is the final unknown in the system, one of the +mass conservation equations from eqs. :eq:`eq_conserve_o` to :eq:`eq_conserve_n` is re-arranged in terms of +:math:`W_{gT}/\sum^{}_{j} X_jM_j`, thereby reducing the number of +equations to solve down to 4. In the COHSN system, the conservation +equation for nitrogen (:eq:`eq_conserve_n`) is chosen, +re-arranged to be substituted into the remaining conservation equations +as + +.. math:: \mathcal{N} = W_{gT}/\sum^{}_{j} X_jM_j = \frac{W_{TN}/M_{\ce{N}} - w^m_{\ce{N2}}(\textcolor{blue}{X_{\ce{N2}}}, \textcolor{blue}{X_{\ce{O2}}}, P)/M_{\ce{N2}}}{2\,\textcolor{blue}{X_{\ce{N2}}}}. + :label: eq_N + +The system of 4 equations to be solved simultaneously can now be derived: + +.. math:: + :label: eq_worked_example_o + + \begin{gathered} + \frac{W_{TO}}{M_{\ce{O}}} = \mathcal{N} (2\textcolor{blue}{X_{\ce{O2}}} + 2\textcolor{red}{X_{\ce{SO2}}} + \textcolor{red}{X_{\ce{H2O}}} + \textcolor{blue}{X_{\ce{CO}}} + 2\textcolor{red}{X_{\ce{CO2}}}) + \frac{w^m_{\ce{H2O}}(\textcolor{red}{X_{\ce{H2O}}},P)}{M_{\ce{H2O}}}\\ + + \frac{2w^m_{\ce{CO2}}(\textcolor{red}{X_{\ce{CO2}}}, \textcolor{blue}{X_{\ce{O2}}}, P, T, melt)}{M_{\ce{CO2}}} + \frac{w^m_{\ce{CO}}(\textcolor{blue}{X_{\ce{CO}}}, P)}{M_{\ce{CO}}} + \frac{W_{O(Fe)}}{M_{\ce{O}}}, + \end{gathered} + +.. math:: + :label: eq_worked_example_h + + \begin{gathered} + \frac{W_{TH}}{2M_{\ce{H}}} = \mathcal{N} (\textcolor{red}{X_{\ce{H2O}}} + \textcolor{red}{X_{\ce{H2}}} + \textcolor{red}{X_{\ce{H2S}}} + 2\textcolor{red}{X_{\ce{CH4}}}) + \frac{w^m_{\ce{H2}}(\textcolor{red}{X_{\ce{H2}}}, P)}{M_{\ce{H2}}} + \ce{w^m_{\ce{H2O}}(\textcolor{red}{X_{\ce{H2O}}}, P)}{M_{\ce{H2O}}}\\ + + \frac{2\,w^m_{\ce{CH4}}(\textcolor{red}{X_{\ce{CH4}}}, P)}{M_{\ce{CH4}}}, + \end{gathered} + +.. math:: + :label: eq_worked_example_s + + \begin{gathered} + \frac{M_{TS}}{M_{\ce{S}}} = \mathcal{N} (\textcolor{red}{X_{\ce{SO2}}} + \textcolor{red}{X_{\ce{H2S}}} + 2\textcolor{blue}{X_{\ce{S2}}}) + \frac{w^m_{\ce{S^{2-}}}(\textcolor{blue}{X_{\ce{S2}}}, \textcolor{blue}{X_{\ce{O2}}}, P, T, melt)}{M_{\ce{S}}}\\ + + \frac{w^m_{\ce{S^{6+}}}(\textcolor{blue}{X_{\ce{S2}}}, \textcolor{blue}{X_{\ce{O2}}}, P, T, melt)}{M_{\ce{S}}}, + \end{gathered} + +.. math:: + :label: eq_worked_example_c + + \begin{gathered} + \frac{W_{TC}}{W_C} = \mathcal{N} (\textcolor{blue}{X_{\ce{CO}}} + \textcolor{red}{X_{\ce{CO2}}} + \textcolor{red}{X_{\ce{CH4}}}) + \frac{w^m_{\ce{CO2}}(\textcolor{red}{X_{\ce{CO2}}}, \textcolor{blue}{X_{\ce{O2}}}, T, melt)}{M_{\ce{CO2}}} + \frac{w^m_{\ce{CO}}(\textcolor{blue}{X_{\ce{CO}}}, P)}{M_{\ce{CO}}}\\ + + \frac{w^m_{\ce{CH4}}(\textcolor{red}{X_{\ce{CH4}}}, P)}{M_{\ce{CH4}}}, + \end{gathered} + +where :math:`melt` is the major element composition of the silicate +melt, and :math:`W_{O(Fe)}` is the weight fraction of oxygen tied up in +FeO and Fe\ :sub:`2`\ O\ :sub:`3` within the silicate melt. This is controlled by the :math:`f_{\ce{O2}}` of the system, and is calculated by finding :math:`F`, the ratio of FeO/Fe\ :sub:`2`\ O\ :sub:`3` in the melt, using either Kress and Carmichael (1991) or Righter et al. (2013), then + +.. math:: + :label: eq_oinmelt + + W_{O(Fe)} = M_{\ce{O}} \frac{W_{\ce{Fe}}}{M_{\ce{Fe}}} \frac{1+3F}{1 + 2F}. + +These four equations are solved numerically using the SciPy +optimize.fsolve function, formulated as a vector equality + +.. math:: + :label: eq_main_solve_matrix + + \begin{pmatrix} + [W_O]\\ + [W_H]\\ + [W_S]\\ + [W_C] + \end{pmatrix} + - + \begin{pmatrix} + [W_O]_{predicted}(X_{\ce{O2}}, X_{\ce{CO}}, X_{\ce{S2}}, X_{\ce{N2}})\\ + [W_H]_{predicted}(X_{\ce{O2}}, X_{\ce{CO}}, X_{\ce{S2}}, X_{\ce{N2}})\\ + [W_S]_{predicted}(X_{\ce{O2}}, X_{\ce{CO}}, X_{\ce{S2}}, X_{\ce{N2}})\\ + [W_C]_{predicted}(X_{\ce{O2}}, X_{\ce{CO}}, X_{\ce{S2}}, X_{\ce{N2}})\\ + \end{pmatrix} + = + \begin{pmatrix} + 0\\ + 0\\ + 0\\ + 0 + \end{pmatrix}. + +The fsolve function makes an initial guess of the value for the vector +:math:`\begin{bmatrix}X_{\ce{O2}}\\X_{\ce{CO}}\\X_{\ce{S2}}\\X_{\ce{N2}} +\end{bmatrix}`, and then iterates, refining the guesses as +appropriate until a solution is found. As the pressure is decreased in +steps by EVo, the initial guess for each new pressure is the solution to +the previous step. For the starting pressure, the entire system is +constrained, so iteration is not required (see +:ref:`section_run_types` for details on how this is achieved). + +.. _section_cands: + +Handling additional phases: carbon and sulfur +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The base degassing model described above deals with the case where there +are 2 phases present: the silicate melt, containing a fraction of +dissolved volatiles; and an exsolved gas phase. However, in the case of +both C and S, additional phases can form outside of this 2-phase system. +Here I describe how EVo handles these scenarios. + +Carbon +^^^^^^ + +A carbon-bearing system can become graphite saturated at low :math:`f_{\ce{O2}}` +(or, alternatively, diamond saturated at very high pressures; Holloway et al., 1992; LaTourrette and Holloway, 1994; Frost and Wood, 1997). :numref:`c_system` +shows a simplified C-O magmatic system, illustrating the relationship +between graphite and other phases. Note that EVo includes H which also +interacts with C as |CH4|; however this has been excluded from +:numref:`c_system` to clarify the relationship between graphite and +the melt/volatile phase. At highly reduced conditions C may +also be soluble in the silicate melt as a range of O and H-bearing +species; however, as the species and reactions involved are still highly +debated these have also been left off :numref:`c_system`. + +.. _c_system: +.. figure:: images/c_system.png + :align: center + :width: 66% + + A simple diagram of the C-O system in silicate melts. Each + colour represents a different phase (gas/volatile (v), silicate melt + (m) and a precipitated graphite phase). Phases and species which + dominate under reducing conditions are to the left, while those which + dominate under oxidising conditions are to the right. Two-way arrows + indicate species that can interact within, and between different, + phases. + +In a graphite-saturated melt, the :math:`f_{\ce{CO2}}` of the system is +controlled by + +.. math:: \ce{C_{(graphite)} + O^{2-}_{(melt)} + O2 <=> CO3^{2-}_{(melt)}} + :label: eq_carb2co2 + +.. math:: \ce{CO3^{2-}_{(melt)} <=> CO2_{(gas)} + O^{2-}_{(melt)}} + :label: eq_co32co2mol + +.. math:: \mathbf{Net:} \quad \ce{C_{(graphite)} + O2 <=> CO2_{(gas)}}. + :label: eq_graphite_net + +As graphite is a pure, solid phase, it has a chemical activity of 1. The equilibrium constant equation for eq. :eq:`eq_graphite_net` is therefore written as + +.. math:: K_{C} = \frac{f_{\ce{CO2}}}{f_{\ce{O2}}}, + :label: eq_graphite_eqb + +showing that in a graphite saturated melt, :math:`f_{\ce{CO2}}` is +entirely controlled by the of the melt (Holloway et al., 1992). +:math:`K_C` is calculated according to the equation of Holloway et al. (1992), + +.. math:: \log_{10}K_C = 40.07639 - 2.53932\times10^{-2} \cdot T + 5.27096\times10^{-6} \cdot T^2 + 0.0267 \cdot (P - 1) / T, + :label: eq_graphite_holloway + +with T in Kelvin and P in bar. + +The silicate system in EVo is checked for graphite saturation by +comparing the value of :math:`f_{\ce{CO2}}` as calculated using the main +model as described above (which assumes no graphite saturation), to that +calculated using :eq:`eq_graphite_eqb` (:math:`f_{\ce{CO2}, graphite}`). If +:math:`f_{\ce{CO2}, graphite} < f_{\ce{CO2}}`, then the melt must be +graphite saturated (e.g., see Fig 7b of Eguchi and Dasgupta, 2018). +If graphite saturation is found, EVo then forces the use of the +Eguchi and Dasgupta (2018) model for |CO2| solubility. +The thermodynamics and structure of this model mean that it is +applicable under both graphite and fluid saturated conditions, unlike +the other option implemented in EVo, and can therefore be fed either +:math:`f_{\ce{CO2}, graphite}` or :math:`f_{\ce{CO2}}`. + +Graphite saturation is only implemented in the ‘closed system’ degassing +option; because in a closed system material cannot be removed, any +graphite present in the melt at volatile saturation is assumed to +decompress with the rest of the gas-melt system. The separate graphite +phase decompresses with the gas-melt system, gradually depleting as +carbon is degassed into the gas phase and the graphite phase replenishes +the carbon dissolved in the melt (see :numref:`c_system`). Graphite +exhaustion (where the melt becomes graphite under-saturated) is detected +if, to solve for a graphite saturated system, the graphite mass in the +system must be negative. The system is first checked for graphite +saturation during the setup of the system, and then after every pressure +step if the system is not known to already be graphite-saturated. + +When solving for a graphite-saturated system, the fugacity of all 3 +carbon-bearing species can be determined solely using the :math:`f_{\ce{O2}}` and :eq:`eq_graphite_eqb`. This removes an extra unknown from all systems containing carbon; +in the example shown in :ref:`section_worked_example`, +:eq:`eq_worked_example_c` is no longer used, and :math:`X_{\ce{CO}}` is +no longer a variable being +solved for. The mass of the carbon reservoir is calculated as the +difference between the total mass of carbon in the system and the mass +calculated as being present in the gas and dissolved in the silicate +melt. + +.. _`section2:sulfur`: + +Sulfur +^^^^^^ + +As with carbon, the sulfur system has added complexities in that +additional phases can form alongside the two accounted for in EVo (the +silicate melt and gas phase). However, unlike with graphite saturation, +there are two additional phases which have to be considered, the +abundances of which vary according to :math:`f_{\ce{O2}}`. The simplified +system is shown in :numref:`s_system`, which similarly to +:numref:`c_system` shows an idealised S-O system, excluding H. + +.. _s_system: +.. figure:: images/ery_sulfur_b.png + :width: 66% + :align: center + + A simple diagram of the S-O system in silicate melts. Each colour represents a different phase. Phases and species which dominate under reducing conditions are to the left, while those which dominate under oxidising conditions are to the right. Two-way arrows indicate species that can interact within, and between different, phases. Taken from Hughes et al. (2021). + +.. _s_curve: +.. figure:: images/s_curve.pdf + :width: 66% + :align: center + + Sulfur speciation as a function of oxygen fugacity, after Jugo et al. (2010). + +Under reducing conditions (approximately :math:`<`\ FMQ, see :numref:`s_curve`) sulfur dissolves in silicate melts as S\ :sup:`2-`: + +.. math:: + :label: eq_sulfide_dissolution + + \ce{S^{2-}_{(melt)} + 0.5O2 <=> O^{2-}_{(melt)} + 0.5S2_{(gas)}}. + +By comparison, under oxidised conditions, sulfur dissolves in the melt +as SO\ :sub:`4`\ :sup:`2-` through + +.. math:: \ce{SO4^{2-}_{(melt)} <=> 1.5O2 + O^{2-}_{(melt)} + 0.5S2_{(gas)}}, + :label: eq_sulfate_dissolution + +where sulfur is instead in it’s S\ :sup:`6+` state (Paris et al., 2001). It is well +established from experimental studies that a silicate melt can contain S +dissolved as one, the other, or both S\ :sup:`2-` and S\ :sup:`6+` across the :math:`f_{\ce{O2}}` range of terrestrial magmas (Carroll and Rutherford, 1988; Nilsson and Peach, 1993; Metrich and Clocchiatti, 1996). Magma produced on other planets across the solar system is thought to +have an :math:`f_{\ce{O2}}` which is less than the :math:`f_{\ce{O2}}` of terrestrial mid-ocean ridge basalts (MORB, e.g., Schmidt et al., 2013b; Zolotov et al., 2013), +which is well within the S\ :sup:`2-` stability field and is therefore expected to +contain minimal S\ :sup:`6+`. + +When a magma becomes sulfur-saturated, the phase precipitated is :math:`f_{\ce{O2}}` +dependent. A melt can become sulfide (S\ :sup:`2-`) saturated + +.. math:: \ce{FeO_{(melt)} + 0.5S2 <=> FeS_{(sulfide)} + 0.5O2}, + :label: eq_scss + +or sulfate (SO\ :sub:`4`\ :sup:`2-`) saturated (usually speciated as anhydrite, CaSO4). Silicate +melts can also be multiply saturated with sulfide + sulfate :math:`\pm` +gas, resulting in a sulfur solubility maximum (Jugo, 2009; Hughes et al., 2021), +with a corresponding solubility minimum where both S\ :sup:`2-` and S\ :sup:`6+` +are dissolved in approximately equal concentrations (Hughes et al., 2021). Due to the +complexity of the sulfur system and the fact that the additional +sulfide/sulfate phases are not pure S (as is the case for graphite), EVo +follows the convention of other multi-component degassing models by not +simulating magmas which are sulfur-saturated, in either phase +(e.g., CHOSETTO, D-Compress: Moretti et al., 2003; Burgisser et al., 2015). +Instead, EVo only deals with sulfur in the gas phase (speciated as |S2|, |SO2| or |H2S|), +and dissolved in the melt as either S\ :sup:`2-` or S\ :sup:`6+`. + +The dissolution of |S2| into a silicate melt has been shown experimentally to +follow eq. :eq:`eq_sulfide_dissolution` (O’Neill and Mavrogenes, 2002; O’Neill, 2021). +As the :math:`\ce{O}^{2-}_{\mathrm{(melt)}}` is viewed as a ‘structural +element’ with abundances far greater than that of S\ :sup:`2-`, eq. :eq:`eq_sulfide_dissolution` can define the solubility law for |S2| as + +.. math:: + :label: eq_sulfide_cap + + C_{\ce{S^{2-}}} = w^m_{\ce{S^{2-}}} \cdot \left(\frac{f_{\ce{O2}}}{f_{\ce{S2}}}\right)^{0.5} + +where :math:`C_{\ce{S^{2-}}}` is the ‘sulfide capacity’ of the melt +(Fincham et al., 1954), analogous to +the equilibrium constant :math:`K` of :eq:`eq_sulfide_dissolution`. +:math:`C_{\ce{S^{2-}}}` is highly compositionally dependent and can, for +given pressure and temperature conditions, be reformulated as a constant +for a given melt composition. Those parametrisations for +:math:`C_{\ce{S^{2-}}}` which are implemented as options in EVo are +listed in :numref:`table_solubility_laws`. While :eq:`eq_sulfate_dissolution` can similarly be +re-arranged to generate a sulfate capacity :math:`C_{\ce{S^{6+}}}` +(Fincham et al., 1954) + +.. math:: + :label: eq_sulfate_cap + + C_{\ce{S^{6+}}} = w^m_{\ce{S^{6+}}} \left(f_{\ce{S2}} f_{\ce{O2}}^3\right)^{-0.5}, + +at the time of writing no experimental validation (on magmatically +relevant silicate compositions), nor expressions equivalent to those +provided by, e.g., O’Neill (2021) for +:math:`C_{\ce{S^{2-}}}` have been published – although results mentioned +in abstract (O’Neill and Mavrogenes, 2019) confirm the relationship and suggest an expression for :math:`C_{\ce{S^{6+}}}` may be forthcoming. To determine the amount of +in the melt, EVo instead uses Nash et al. (2019) which +calculates the ratio of S\ :sup:`6+`/S\ :sup:`2-` as + +.. math:: + :label: eq_nash2019 + + \log_{10} \left(\frac{w^m_{\ce{S^{6+}}}}{w^m_{\ce{S^{2-}}}} \right) = 8 \log_{10} \left(\frac{w^m_{\ce{Fe^{3+}}}}{w^m_{\ce{Fe^{2+}}}} \right) + \frac{8.7436\times10^6}{T^2} - \frac{27703}{T} + 20.273, + +for temperatures (T) of 1000 - 2000 K. + +The amount of S\ :sup:`2-` which can be dissolved in a silicate melt at sulfide saturation, controlled by :eq:`eq_scss`, is referred to as the “Sulfide Content at Sulfide Saturation”, or SCSS. The SCSS has been extensively studied over the past few decades (O’Neill and Mavrogenes, 2002; Liu et al., 2007; Fortin et al., 2015; Wykes et al., 2015; Smythe et al., 2017; O’Neill, 2021), and in EVo is calculated using the model of Liu et al. (2007). EVo checks for sulfide saturation by comparing :math:`w^m_{\ce{S^{2-}}}` as calculated +using :eq:`eq_sulfide_cap` to the SCSS. If :math:`w^m_{\ce{S^{2-}}}` :math:`\geq` SCSS, then the melt is sulfide saturated and EVo produces a warning stating that sulfide saturation has +been reached, and the model is no longer valid. + +.. _section_run_types: + +Run types, set-up options and input parameters +---------------------------------------------- + +Two different run types can be performed while calculating a +decompression path: “Closed system” or “Open system”. Closed-system +runs, the default, assume that the gas is in equilibrium with the melt +and there is no physical-chemical separation between the two phases. +This models a system where the gas bubbles are entrained in the melt and +rise at the same speed as the magma. + +In open-system degassing, a fraction of the gas released as a magma +decompresses is assumed to be lost (or chemically isolated) from the +melt as it is produced, mimicking a scenario where gas bubbles rise to +the surface through a magma. EVo simulates this by resetting the gas +fraction after every pressure step + +.. math:: W_{gT}(P-dP) = a W_{gT}, + :label: eq_open_system_gasloss + +where :math:`dP` is the pressure step (bar) and :math:`a` is the +fraction of the gas phase which is removed from the system after each +increment (using :math:`a` rather than fully removing all exsolved +volatiles allows for an exploration of efficient vs inefficient +open-system degassing). Once :math:`W_{gT}` has been reset, the new +masses of the volatile elements left in the system are calculated before +the pressure is incremented. EVo tracks the evolution of the magma as a +function of pressure in this scenario, and does not track the +accumulated composition of the degassed volatiles, although it does +provide the chemistry of the instantaneous exsolved gas phase at each +pressure step. + +There are three different options which can be used to set-up EVo for a +decompression calculation, discussed below. These options can also be +used to find the equilibrium conditions at a single pressure, either at +a set pressure using option (1), or the volatile saturation pressure, +found using either (2) or (3). + +.. _section_standard_setup: + +1) Where :math:`P_{start}` and :math:`W_{gT}` are known +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +Where the gas mass fraction is known at the pressure of interest (either +for the start of a degassing run or to calculate the composition at a +single point), along with the temperature, the state of the system can +be found using a subset of additional parameters as described in :numref:`table_input_options`. + +.. _table_input_options: +.. list-table:: : EVo input options when :math:`P_{start}` is known, dependent on the volatile system being used. (Fe) indicates that each system can be run with or without oxygen exchange with the melt, this does not affect how the run is set up. \*Theoretically possible, but as yet not implemented as this is an unusual way of running a degassing model. + :class: borderless + + * - .. image:: images/table2.png + +Using these input values, all other relevant parameters for the initial +conditions can be calculated. Where a melt weight fraction of a volatile +is provided, the respective fugacity is calculated using the solubility +laws. + +Solubility models usually provide analytical expressions for the +concentration of a volatile species dissolved in a magma as a function +of either partial pressure, or fugacity. Analytical expression for the +inverse (calculating the fugacity of a species given the concentration +dissolved in the melt) are usually not provided, and in some cases +cannot be found by simple re-arrangement; in these cases (both the models of Dixon (1997) and Eguchi and Dasgupta (2018)) numerical methods are +used to obtain the correct fugacity. + +The 1-4 fugacities (depending on the number of elements considered), P, +T, reactions Eqs. :eq:`eq_k1` to :eq:`eq_k5` and :eq:`eq_sumtoone` are then used +algebraically to find all :math:`X_i` and :math:`f_i`. Once all molar +fractions have been calculated, the masses of each element in the system +can be computed (Eqs. :eq:`eq_conserve_o` to :eq:`eq_conserve_n`), along with the mass of oxygen fixed by FeO\* using :eq:`eq_oinmelt`. + +.. _`section:find_psat`: + +2) Calculating :math:`P_{sat}`, the volatile saturation pressure +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +In many scenarios, having to specify a starting pressure and gas +fraction in order to model a decompression path is inconvenient, as +these variables are often unknown. Instead, it is preferable to +calculate the saturation pressure for a given system volatile content +(e.g., the |H2O| and |CO2| concentration of the melt prior to the onset of +degassing), and start the decompression from there. EVo allows this +setup option by asking for the :math:`f_{\ce{O2}}` at the point of volatile saturation +(usually taken to be the :math:`f_{\ce{O2}}` of the magma source), and the undegassed +volatile content of the melt. + +When the saturation pressure of a magma with a given dissolved volatile +concentration is found, the following equality holds: + +.. math:: + :label: eq_saturation_condition + + \sum_{i=1} P_i - P = 0 + +where :math:`P` is the total pressure and :math:`P_i` is the partial +pressure of species :math:`i` calculated according to it’s corresponding +solubility law, for a fixed concentration in the melt. The saturation +pressure is found by numerically solving :eq:`eq_saturation_condition` +for :math:`P`, where + +.. math:: + :label: eq_satp_solver + + P_i = \frac{fi(w_i^m, P \dots)}{\gamma_i}, + +using the Brent method as implemented in SciPy (Virtanen et al., 2020). This method +was chosen as it allows the solution to be bracketed; this prevents the +numerical solver guessing a negative pressure which is invalid within +the inverted solubility laws. Once :math:`P_{sat}` has been found, +:math:`W_{gT}` is set to 1×10\ :sup:`−6` wt% (an arbitrarily small +number sufficiently close to zero to not affect the results but +successfully initialise the calculation), the atomic volatile masses in +the system are calculated as in eqs. :eq:`eq_conserve_o` to :eq:`eq_conserve_n` and decompression steps can commence. + +.. _`section:evo_atomic_saturation`: + +3) Calculating :math:`P_{sat}` and volatile speciation +~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + +The above two setup options work well for use-cases such as modelling +the degassing path of volcanic samples, or single volcanic systems where +the starting conditions (e.g., :math:`f_{\ce{O2}}`, volatile content) are similar. +However, in planetary science a key component of research is comparing the +gas phase/atmospheres produced when the only variable is the starting :math:`f_{\ce{O2}}`. +In set-up options 1 & 2, this poses a problem when dealing with volatile +elements which exist as multiple dissolved species, particularly for H +and C. + +Using the simple case of the |H2O|-|O2|-|H2| system, H in the melt can be dissolved +either as |H2O|, or |H2|, although in oxidised to moderately reduced magmas, as +seen on Earth, the |H2| component is minimal. A fixed magmatic |H2O| content +imposes a constant :math:`f_{\ce{H2O}}`; since the is being lowered, +:math:`f_{\ce{H2}}` must increase to maintain equilibrium. A greater +:math:`f_{\ce{H2}}` also enforces a higher magma content; therefore, as +the :math:`f_{\ce{O2}}` of a magma is decreased with a fixed |H2O| content, the dissolved H2 must increase proportionally (:numref:`h_with_fo2`). + +.. _h_with_fo2: +.. figure:: images/h_with_fo2.pdf + :width: 66.0% + :align: center + + The mass of |H2O|, |H2| and total H in the silicate melt at volatile + saturation, when the enforced saturation conditions are :math:`f_{\ce{O2}}` and the + dissolved |H2O| content. + +When comparing degassing regimes where the :math:`f_{\ce{O2}}` is varied, +:numref:`h_with_fo2` shows it is not sufficient to simply fix the +magma volatile content using a single species. E.g., volcanic gases +released at FMQ compared to FMQ-6 in :numref:`h_with_fo2` cannot be +directly compared, as the more reduced scenario has almost double the +amount of H dissolved in the melt at the point of volatile saturation, +hence the :math:`f_{\ce{O2}}` is not the only factor being varied. Instead, for meaningful +isolation of the effects of :math:`f_{\ce{O2}}` changes alone, the *total mass of each +volatile element* must be enforced, with the speciation at volatile +saturation allowed to vary freely according to the :math:`f_{\ce{O2}}` conditions. + +Similarly to set-up option 2, this is solved numerically. EVo is +initialised with the :math:`f_{\ce{O2}}` at volatile saturation, and the total atomic masses +of each volatile element in the system. The solver then guesses the +proportion of each element which is dissolved as |H2O|, |CO2|, S\ :sup:`2-` and +|N2|, with the remainder dissolved as other species, or in the gas phase (which is +minimal at volatile saturation, but relevant in order for the system of +equations to work). As in the main model, a system of 4 equations (for +the COHSN system) is solved simultaneously: + +.. math:: + :label: eq_atomic_solve_matrix + + \begin{pmatrix} + [W_H]_{set}\\ + [W_C]_{set}\\ + [W_S]_{set}\\ + [W_N]_{set} + \end{pmatrix} + - + \begin{pmatrix} + [W_H]_{predicted}(w^m_{\ce{H2O}}, w^m_{\ce{CO2}}, w^m_{\ce{S^{2-}}}, w^m_{\ce{N2}})\\ + [W_C]_{predicted}(w^m_{\ce{H2O}}, w^m_{\ce{CO2}}, w^m_{\ce{S^{2-}}}, w^m_{\ce{N2}})\\ + [W_S]_{predicted}(w^m_{\ce{H2O}}, w^m_{\ce{CO2}}, w^m_{\ce{S^{2-}}}, w^m_{\ce{N2}})\\ + [W_N]_{predicted}(w^m_{\ce{H2O}}, w^m_{\ce{CO2}}, w^m_{\ce{S^{2-}}}, w^m_{\ce{N2}})\\ + \end{pmatrix} + = + \begin{pmatrix} + 0\\ + 0\\ + 0\\ + 0 + \end{pmatrix}. + +This is similar in form to eq. :eq:`eq_main_solve_matrix`, but in this case the +variables are the mass fractions of dissolved volatiles, rather than the +molar fractions of species in the gas phase. When a guess is made for +:eq:`eq_atomic_solve_matrix` (structured as +the vector :math:`\begin{bmatrix}w^m_{\ce{|H2O|}}\\w^m_{\ce{CO2}}\\w^m_{\ce{S^{2-}}}\\w^m_{\ce{N2}}\end{bmatrix}`), it is first passed to :eq:`eq_saturation_condition` to find the +saturation pressure, :math:`P_{sat}`, before the +:math:`[W_i]_{predicted}` values are calculated. The key difference +between this setup option and options 1 & 2 is that in this case the +mass of each element in the system is known and the speciation is +calculated; in the previous cases a subset of the melt or gas phase +speciation is known and the solution finds the total masses of each +element. + +Benchmarking +------------ + +During the development of EVo, individual elements of the model were each +tested to ensure correct implementation. For example, the Kress and +Carmichael (1991) relationship between :math:`f_{\ce{O2}}` and ferric/ferrous iron was +tested against the spreadsheet of Iacovino (2021), and individual solubility +laws were tested either against the relevant published calculator (e.g., +excel spreadsheets in the case of Newman and Lowenstern, 2002; Eguchi and Dasgupta, +2018), or against figures in the original publication where these were not provided. + +EVo has also been initially tested against DCompress to ensure numerical +accuracy in the thermodynamic calculation of the gas phase chemistry. In +order to maximise similarity between the two models, for this comparison +EVo is run using the |H2O|, |H2| and |CO2| solubility laws implemented by +Burgisser et al. (2015) in DCompress; to distinguish this setup from EVo +under standard run conditions (with |H2| and |CO2| solubility laws from Gaillard et al. +(2003) and Eguchi and Dasgupta (2018), respectively), this +version is referred to as EVo(DC). The two major differences left +between the two models are therefore the treatment of sulfur solubility +(EVo uses the sulfide capacity, as described in +:ref:`section_cands`, while DCompress has a solubility law for |SO2| and |H2S|), and the source of the equilibrium constants. DCompress uses equations +from Ohmoto and Kerrick (1977) +to calculate equilibrium constants K1-K5, while EVo and EVo(DC) +calculate them using more recent thermochemical data from +Chase (1998) as described +in :ref:`section_overview_of_model`. The differences between the two +methods in K1-5 are shown in :numref:`k_comparison`. + +.. _k_comparison: +.. figure:: images/eqb_constants.png + :width: 50% + :align: center + + The difference between the equilibrium constants as calculated in + DCompress and EVo, expressed as K1(EVo)/K1(DC). + +EVo(DC) and DCompress were run at a single pressure, using input method +(1) as described in :ref:`section_standard_setup`, a gas weight +fraction of 0.001%, and other parameters as listed on the individual +panels within :numref:`dc2evo_initial_conditions` to compare the +calculation of the gas phase chemistry. + +.. _dc2evo_initial_conditions: +.. figure:: images/dc_2evo_comp.png + :width: 70.0% + :align: center + + The absolute fractional difference (calculated as + :math:`\lvert X_{i, \mathrm{EVo(DC)}} - X_{i, \mathrm{DC}} \rvert/X_{i, \mathrm{DC}}` + where :math:`X_i` is the mole fraction of each species :math:`i`) + between the gas phase speciation of DCompress and EVo(DC) when setup + with the same input parameters. The magma composition used in the + calculations is the default basalt composition for DCompress + +In the case of the OH and COH systems, the difference between the +results of the two models are :math:`\leq` 5% for each species. These +differences can be largely explained by differences in the value of +equilibrium constants K1, K2 and K3 (:numref:`k_comparison`). These +differences in the OH/C systems can be reduced to :math:`\leq` 1% by +using the equilibrium constants from DCompress in EVo(DC). The +particularly large difference on K5 at 1473 K is responsible for the +largest difference between the two models, in the SOH and COHS setups. +As the gas phase speciation is determined first using this method, and +the melt S content is calculated based on the gas fugacity, the sulfur +solubility law does not affect the gas phase speciation in this example; +however the difference in the amount of sulfur dissolved in the melt is +shown in the SOH and COHS panels of +:numref:`dc2evo_initial_conditions`, as calculated by the different +solubility laws. + +.. _SOH_evodc_fugacities: +.. figure:: images/SOH_evodc_fugacities.png + :width: 50 % + :align: center + + Gas phase fugacities for the sulfur species. Both EVo and DC have + been plotted, differences are within the linewidth. + +.. _SOH_comp_melt: +.. figure:: images/SOH_evodc_melt.png + :width: 50 % + :align: center + + Dissolved sulfur species + +The difference in the treatment of sulfur solubility between DCompress +and EVo is demonstrated in :numref:`SOH_evodc_fugacities` and +:numref:`SOH_comp_melt`. When given the same +gas-phase fugacities for |S2|, |SO2| and |H2S| across an :math:`f_{\ce{O2}}` range, EVo consistently +predicts a higher total dissolved S content, following an experimentally +well-established trend that in reduced melts (where sulfur dissolves as S\ :sup:`2-` +), the solubility of sulfur increases with reducing :math:`f_{\ce{O2}}` +(e.g., Fincham et al., 1954; Katsura and +Nagashima, 1974; Backnaes and Deubener, 2011; Lesne et al., 2015). +In contrast, DCompress predicts a largely constant total sulfur content, +decreasing slightly with reducing :math:`f_{\ce{O2}}`. + +Decompression benchmarking +~~~~~~~~~~~~~~~~~~~~~~~~~~ + +I now compare the results of EVo (now using it’s default solubility laws) to those of 3 different models of volcanic degassing which include sulfur: DCompress (Burgisser et al., 2015), SolEx (Witham et al., 2012) and Chosetto, the newly released implementation of Moretti et al. (2003). Two decompression +runs are shown in :numref:`decompression`, one oxidised with an of NNO+0.5, and one slightly reduced at NNO-2. +SolEx is only shown in the oxidised example, as it cannot be run with an +lower than NNO+0.5. All 4 models were initialised with 1 wt% , 500 ppm +and 3000 ppm S. Only EVo has the option to automatically find the +volatile saturation point; SolEx begins all runs at 4000 bar, while +DCompress was run by adjusting the starting pressure until the total +sulfur content matched the required 3000 ppm and Chosetto’s starting +pressure was adjusted until it could numerically resolve. + +.. _decompression: +.. figure:: images/decompression.png + + The melt volatile content, |H2O|/|CO2| ratio and S/|CO2| ratio in the gas phase for 4 different models during decompression. Each model was initiated with 1 wt%\ |H2O|, 500 ppm |CO2| and 3000 ppm S. EVo is the only model which allows the saturation pressure to be found automatically; SolEx always starts at 4000 bar, Psat was found manually for DCompress and Chosetto. At NNO-2, Chosetto was started at the first point it could numerically resolve. + +:numref:`decompression` shows (a) the volatile content of the melt, (b) the |H2O|/|CO2| ratio in the gas +phase and (c) the S/|CO2| ratio of the gas phase. All 4 models show close +agreement on the dissolved volatile content as the pressure decreases, +aside from the SolEx S content, which starts to rapidly decrease at a +much higher pressure than other models in :numref:`decompression` a. Chosetto also +shows a higher content in the melt in :numref:`decompression` d. All 4 models +also show similar trends in the |H2O|/|CO2| ratio, although both SolEx and +Chosetto both have systematically more |CO2|-rich gases. + +The largest difference between models lies in the behaviour of sulfur +during degassing. SolEx shows sulfur degassing at much higher pressures +than other models, while DCompress produces a more gradual release with +pressure. EVo most closely resembles Chosetto under oxidised conditions, +with sulfur starting to rapidly degas at 200 bar. Under reduced +conditions the differences between models grow greater; the difference +between Chosetto and EVo potentially being due to the presence of CO and |CH4| in +the gas phase of EVo while Chosetto only considers |CO2|. + +The benchmarking comparisons shown in this section indicate that EVo +compares well to other publish models for multi-species volcanic +degassing. Where models differ, for example in the COHS subplot of :numref:`dc2evo_initial_conditions`, the driving factors are well +understood. The lack of a standard method/dataset with which to +benchmark and compare such models, which often only show subsets of the +common volatile species or differ in other substantive ways, is an issue +across Earth Sciences; however, :numref:`decompression` shows that on a +limited dataset, EVo follows similar trends to all three of the other +published models which include sulfur in their parametrisations, with +the exception of the sulfur/ behaviour in SolEx, which appears to be an +outlier. + +.. |H2O| replace:: H\ :sub:`2`\ O +.. |H2| replace:: H\ :sub:`2` +.. |O2| replace:: O\ :sub:`2` +.. |CO2| replace:: CO\ :sub:`2` +.. |CH4| replace:: CH\ :sub:`4` +.. |H2S| replace:: H\ :sub:`2`\ S +.. |SO2| replace:: SO\ :sub:`2`\ +.. |S2| replace:: S\ :sub:`2` +.. |N2| replace:: N\ :sub:`2` \ No newline at end of file diff --git a/docs/images/SOH_evodc_fugacities.png b/docs/images/SOH_evodc_fugacities.png new file mode 100644 index 0000000..42267c9 Binary files /dev/null and b/docs/images/SOH_evodc_fugacities.png differ diff --git a/docs/images/SOH_evodc_melt.png b/docs/images/SOH_evodc_melt.png new file mode 100644 index 0000000..292415b Binary files /dev/null and b/docs/images/SOH_evodc_melt.png differ diff --git a/docs/images/c_system.png b/docs/images/c_system.png new file mode 100644 index 0000000..daab2b1 Binary files /dev/null and b/docs/images/c_system.png differ diff --git a/docs/images/coh_satp.png b/docs/images/coh_satp.png new file mode 100644 index 0000000..3bb5ea6 Binary files /dev/null and b/docs/images/coh_satp.png differ diff --git a/docs/images/dc_2evo_comp.png b/docs/images/dc_2evo_comp.png new file mode 100644 index 0000000..0fb691b Binary files /dev/null and b/docs/images/dc_2evo_comp.png differ diff --git a/docs/images/decompression.png b/docs/images/decompression.png new file mode 100644 index 0000000..ad7e62a Binary files /dev/null and b/docs/images/decompression.png differ diff --git a/docs/images/eqb_constants.png b/docs/images/eqb_constants.png new file mode 100644 index 0000000..1176d2e Binary files /dev/null and b/docs/images/eqb_constants.png differ diff --git a/docs/images/ery_sulfur_b.png b/docs/images/ery_sulfur_b.png new file mode 100644 index 0000000..dee5f88 Binary files /dev/null and b/docs/images/ery_sulfur_b.png differ diff --git a/docs/images/evo_overview.png b/docs/images/evo_overview.png new file mode 100644 index 0000000..ca5c8cd Binary files /dev/null and b/docs/images/evo_overview.png differ diff --git a/docs/images/h_with_fo2.pdf b/docs/images/h_with_fo2.pdf new file mode 100644 index 0000000..f40b534 Binary files /dev/null and b/docs/images/h_with_fo2.pdf differ diff --git a/docs/images/model_schematic.pdf b/docs/images/model_schematic.pdf new file mode 100644 index 0000000..578567e Binary files /dev/null and b/docs/images/model_schematic.pdf differ diff --git a/docs/images/s_curve.pdf b/docs/images/s_curve.pdf new file mode 100644 index 0000000..ed6960a Binary files /dev/null and b/docs/images/s_curve.pdf differ diff --git a/docs/images/table1.png b/docs/images/table1.png new file mode 100644 index 0000000..4b11458 Binary files /dev/null and b/docs/images/table1.png differ diff --git a/docs/images/table2.png b/docs/images/table2.png new file mode 100644 index 0000000..49e6da8 Binary files /dev/null and b/docs/images/table2.png differ diff --git a/docs/index.rst b/docs/index.rst new file mode 100644 index 0000000..5ddc0c9 --- /dev/null +++ b/docs/index.rst @@ -0,0 +1,58 @@ +Welcome to EVo's documentation! +======================================= + +EVo is a volcanic outgassing model for COHSN elements. + +The current implementation has been developed in Python 3 and tested on Windows, Linux and macOS. + +Installation/Usage: +******************* + +EVo can be used through local installation, either through a CLI or a webapp. + +To install locally, EVo must be downloaded from GitHub using +:: + + git clone git@github.com:pipliggins/EVo.git + +into the project directory where you wish to use EVo. EVo must then be locally pip-installed: +:: + + cd EVO + python -m pip install ".[streamlit]" + +From this point, EVo can either be imported into your python scripts as a regular module using +:: + + install evo + +and run using +:: + + evo.main('chem_file', 'env_file', 'output_options_file') + +Or EVo can be run directly from the terminal from inside the `evo` directory: +:: + + cd EVO/evo + python dgs.py input/chem.yaml input/env.yaml --output input/output.yaml + +Alternatively you can interact with EVo using the webapp interface, by running: +:: + + cd EVO/webapp + streamlit run streamlit-app.py + + +.. toctree:: + :maxdepth: 2 + :caption: Contents: + + evo_doc + +Indices and tables +================== + +* :ref:`genindex` +* :ref:`modindex` +* :ref:`search` diff --git a/docs/make.bat b/docs/make.bat new file mode 100644 index 0000000..954237b --- /dev/null +++ b/docs/make.bat @@ -0,0 +1,35 @@ +@ECHO OFF + +pushd %~dp0 + +REM Command file for Sphinx documentation + +if "%SPHINXBUILD%" == "" ( + set SPHINXBUILD=sphinx-build +) +set SOURCEDIR=. +set BUILDDIR=_build + +%SPHINXBUILD% >NUL 2>NUL +if errorlevel 9009 ( + echo. + echo.The 'sphinx-build' command was not found. Make sure you have Sphinx + echo.installed, then set the SPHINXBUILD environment variable to point + echo.to the full path of the 'sphinx-build' executable. Alternatively you + echo.may add the Sphinx directory to PATH. + echo. + echo.If you don't have Sphinx installed, grab it from + echo.https://www.sphinx-doc.org/ + exit /b 1 +) + +if "%1" == "" goto help + +%SPHINXBUILD% -M %1 %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% +goto end + +:help +%SPHINXBUILD% -M help %SOURCEDIR% %BUILDDIR% %SPHINXOPTS% %O% + +:end +popd diff --git a/docs/requirements.txt b/docs/requirements.txt new file mode 100644 index 0000000..cf6bc2d --- /dev/null +++ b/docs/requirements.txt @@ -0,0 +1,2 @@ +sphinx==5.3.0 +sphinx_rtd_theme \ No newline at end of file diff --git a/pyproject.toml b/pyproject.toml index 50c6490..680c50d 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -26,12 +26,13 @@ dependencies = ["numpy>=1.21.0", "scipy>=1.7.0", "matplotlib>=3.4.2", "mpmath>=1.2.1", - "pyyaml", "ruamel.yaml", "pandas>=2.1.0",] [project.optional-dependencies] dev = ["ruff", "pre-commit"] +doc = ["sphinx == 5.3.0"] +web = ["streamlit"] test = ["pytest"] [project.urls] diff --git a/src/evo/__init__.py b/src/evo/__init__.py index a2bfe00..ac5e89c 100644 --- a/src/evo/__init__.py +++ b/src/evo/__init__.py @@ -1,9 +1,14 @@ """ EVo + A Python model for volcanic degassing, using the equilibrium constants and mass balance method. """ -from evo.dgs import run_evo as run_evo +# ----------------- IMPORTS ----------------- # +from evo.dgs import run_evo +from evo.multirun import multirun -__all__ = ["run_evo"] +__version__ = "1.0.2" +__author__ = "Philippa Liggins" +__all__ = ["run_evo", "multirun"] diff --git a/src/evo/multirun.py b/src/evo/multirun.py index dcf549d..13915ec 100644 --- a/src/evo/multirun.py +++ b/src/evo/multirun.py @@ -3,10 +3,12 @@ from shutil import copyfile import numpy as np -import yaml +import ruamel.yaml from evo.dgs import run_evo +ryaml = ruamel.yaml.YAML() + def amend_env(file, **kwargs): """ @@ -25,7 +27,7 @@ def amend_env(file, **kwargs): """ with open(file) as f: - env_doc = yaml.full_load(f) + env_doc = ryaml.load(f) for param, val in kwargs.items(): if isinstance(val, np.float64): @@ -34,7 +36,7 @@ def amend_env(file, **kwargs): env_doc[param] = val with open("multirun.yaml", "w") as f: - yaml.dump(env_doc, f) + ryaml.dump(env_doc, f) def multirun(**kwargs): diff --git a/src/evo/readin.py b/src/evo/readin.py index 32fd13c..f99f703 100644 --- a/src/evo/readin.py +++ b/src/evo/readin.py @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------ import numpy as np -import yaml +import ruamel.yaml from numpy import exp, log10 from numpy import log as ln @@ -18,6 +18,8 @@ # bundled scripts from evo.dgs_classes import Gas, Melt, Molecule, Output, RunDef, ThermoSystem +ryaml = ruamel.yaml.YAML() + # ------------------------------------------------------------------------ # FUNCTION DEFINTIONS # ------------------------------------------------------------------------ @@ -44,7 +46,7 @@ def readin_env(f): """ # creates a dictionary of the environment parameters from env.yaml - x = yaml.full_load(f) + x = ryaml.load(f) # setup run definitions run = RunDef() @@ -115,7 +117,7 @@ def readin_chem(f, run, sys): The instantiated Melt class """ - data = yaml.full_load(f) + data = ryaml.load(f) ele_names = [] chems = [] @@ -194,7 +196,7 @@ def readin_output(f): The instantiated Output class """ - x = yaml.full_load(f) + x = ryaml.load(f) out = Output() out.set_outputs(x) diff --git a/webapp/chem.yaml b/webapp/chem.yaml new file mode 100644 index 0000000..d29f49a --- /dev/null +++ b/webapp/chem.yaml @@ -0,0 +1,10 @@ +SIO2: 47.95 +TIO2: 1.67 +AL2O3: 17.32 +FEO: 10.24 +MNO: 0.17 +MGO: 5.76 +CAO: 10.93 +NA2O: 3.45 +K2O: 1.99 +P2O5: 0.51 diff --git a/webapp/env.yaml b/webapp/env.yaml new file mode 100644 index 0000000..612f7e4 --- /dev/null +++ b/webapp/env.yaml @@ -0,0 +1,103 @@ +#run type +COMPOSITION: basalt + +#RUN_TYPE: choose either closed or open +RUN_TYPE: closed + +# Runs only one pressure step, the starting pressure. +SINGLE_STEP: false +# Search for the saturation point based on melt composition; ignores P_Start and WgT +FIND_SATURATION: false +# Finds saturation pressure based on volatile contents as atomic mass fractions (H, S, N, C) + fO2 +ATOMIC_MASS_SET: false + +GAS_SYS: oh +FE_SYSTEM: true +S_SAT_WARN: false + +# model parameters +T_START: 1473.15 +P_START: 3000.0 +P_STOP: 1.0 +DP_MIN: 0.1 +DP_MAX: 100.0 +MASS: 100.0 +WgT: 1.0e-05 + +# Fraction of gas phase lost at each pressure step during open system degassing +LOSS_FRAC: 0.9999 + +# select models +DENSITY_MODEL: spera2000 +# Relate fo2 to Fe2/Fe3 ratio; pick from kc1991 (Kress & Carmicheal 1991), r2013 (Righter et al., 2013)) +FO2_MODEL: kc1991 +# Rock buffer definitions, options: frost1991; +FMQ_MODEL: frost1991 + +# select solubility laws from options. Some options do not include OH-, CO3, CH4 species. +# OH system models, pick from: burguisser2015 (H2O only); +H2O_MODEL: burguisser2015 +# H2 model, pick from: burguisser2015 (ONLY USE TO COMPARE TO DCOMPRESS; ERROR!); gaillard2003 +H2_MODEL: gaillard2003 +# C models, pick from: burguisser2015 (CO2 & CO3); eguchi2018 (CO2, CO3 & graphite) +C_MODEL: burguisser2015 +# CO model, pick either armstrong2015, or None (preferred for higher fO2 cases above ~ IW+1) +CO_MODEL: armstrong2015 +# CH4 model, pick either ardia2013, or None (preferred for higher fO2 cases above ~ IW+1) +CH4_MODEL: ardia2013 +# S models, pick from: sulfide capacity: oneill2002, oneill2020; sulfate capacity: nash2019 ; SCSS: liu2007 +SULFIDE_CAPACITY: oneill2020 +SULFATE_CAPACITY: nash2019 +SCSS: liu2007 +# N model, pick from: libourel2003 (N2 only) +N_MODEL: libourel2003 + +# fO2 can be set relative to a mineral buffer (NNO, FMQ or IW), as an absolute value, or through the FeO/Fe2O3 ratio in the chem.yaml file. +# Or it can be set implicitly in the OH system by initialising using fH2 rather than fO2. +FO2_buffer_SET: true +FO2_buffer: FMQ +FO2_buffer_START: 0.0 + +# SET TOTAL VOLATILES AS ATOMIC MASS FRACTIONS (ppm) + +ATOMIC_H: 0.0 +ATOMIC_C: 0.0 +ATOMIC_S: 0.0 +ATOMIC_N: 0.0 + +# VOLATILE SETTING OPTIONS: +# OH: Either fO2 (see above) OR fH2 +# COH: Pick 2 from the fO2 (see above), fH2, fH2O or wtH2O. You cannot set up with just fO2 and fH2, or specify both fH2O and WtH2O. +# SOH: Pick 2 from the fO2 (see above), fH2, fH2O or wtH2O. You cannot set up with just fO2 and fH2, or specify both fH2O and WtH2O +# COHS: Pick 3 from the fO2 (see above), fH2, fH2O, wtH2O fCO2 or WtCO2. You cannot set up with just fO2 and fH2, or specify either both fH2O and WtH2O or both fCO2 and WtCO2 +# COHSN: Pick 4 from the fO2 (see above), fH2, fH2O, wtH2O fCO2 or WtCO2, plus WtN. You cannot set up with just fO2 and fH2, or specify either both fH2O and WtH2O or both fCO2 and WtCO2 + +# To select an option, switch the boolean specified by 'set' to 'true' and provide a value in the corresponding 'start' option. +# Switch all other volatile setting options to false to avoid triggering an error. + +# select initial volatile compositions, either as gas fugacities + +FH2_SET: false +FH2_START: 0.0 + +FH2O_SET: false +FH2O_START: 0.0 + +FCO2_SET: false +FCO2_START: 0.0 + +# or melt wt fractions - a wt FRACTION still dissolved in the melt at the starting pressure (NOT WHEN EVERYTHING IS DISSOLVED), so 0.1 wt% = 0.001 +WTH2O_SET: false +WTH2O_START: 0.0 + +WTCO2_SET: false +WTCO2_START: 0.0 + +SULFUR_SET: false +SULFUR_START: 0.0 + +NITROGEN_SET: false +NITROGEN_START: 0.0 + +GRAPHITE_SATURATED: false +GRAPHITE_START: 0.0 diff --git a/webapp/evo_logo.png b/webapp/evo_logo.png new file mode 100644 index 0000000..1efd25c Binary files /dev/null and b/webapp/evo_logo.png differ diff --git a/webapp/output.yaml b/webapp/output.yaml new file mode 100644 index 0000000..e0274dc --- /dev/null +++ b/webapp/output.yaml @@ -0,0 +1,7 @@ +# Graphical output selection + +Plt_melt_species: True +Plt_gas_species_wt: False +Plt_gas_species_mol: True +Plt_gas_fraction: True +Plt_fo2_dFMQ: False \ No newline at end of file diff --git a/webapp/streamlit-app.py b/webapp/streamlit-app.py new file mode 100644 index 0000000..1fc5300 --- /dev/null +++ b/webapp/streamlit-app.py @@ -0,0 +1,617 @@ +"""Contains code to enable a web interface for EVo, using StreamLit.""" + +import evo +import streamlit as st +from pathlib import Path +import numpy as np +import ruamel.yaml +import plotly.express as px +import plotly.graph_objects as go + +ryaml = ruamel.yaml.YAML() + + +def amend_chem(file, chem_dict): + """ + Duplicates the chem.yaml file, then edits with the gridsearch values relevant + to the specific run. + """ + + with open(file, "r") as f: + chem_doc = ryaml.load(f) + + for param, val in chem_dict.items(): + if type(val) is np.float64: + chem_doc[param] = float(val) + elif val == 0.0: + pass + else: + chem_doc[param] = val + + with open(file, "w") as f: + ryaml.dump(chem_doc, f) + + +def amend_env(file, env_dict): + """ + Duplicates the env.yaml file, then edits with the gridsearch values relevant + to the specific run. + """ + + with open(file, "r") as f: + env_doc = ryaml.load(f) + + for param, val in env_dict.items(): + # forces the scientific notation to be output correctly to be read as a float + if type(val) is np.float64: + env_doc[param] = float(val) + else: + env_doc[param] = val + + with open(file, "w") as f: + ryaml.dump(env_doc, f) + + +oxides = {} +env_options = {} +default_laws = { + "FO2_MODEL": "kc1991", + "H2_MODEL": "gaillard2003", + "C_MODEL": "burguisser2015", + "CO_MODEL": "armstrong2015", + "CH4_MODEL": "ardia2013", + "SULFIDE_CAPACITY": "oneill2020", +} +volatile_params = {} +models = {} + +st.set_page_config(layout="wide") + +# # Use this to control which things are visible. +# # i.e. use 'disabled' to grey out options which are not +# # valid based on the run options. +# https://docs.streamlit.io/library/api-reference/widgets/st.radio +# if "startType" not in st.session_state: +# # st.session_state.visibility = "visible" +# # st.session_state.disabled = False +# # st.session_state.horizontal = False + +if "default_sollaws" not in st.session_state: + st.session_state.default_sollaws = True + st.session_state.no_loaded_data = True + +title_image, _, title_col = st.columns([4, 1, 15]) + +with title_image: + st.image("webapp/evo_logo.png") + +with title_col: + st.title("EVo: A volcanic degassing model") + + st.text( + "This is a web interface for the EVo volcanic degassing model.\n" + + "Use the input options below to select setup options and load your " + + "input options by clicking 'Load data'.\n" + + "Then click 'Run EVo' to generate a downloadable dataframe of results " + + "and example plots." + ) + +col1, _, col2, col3 = st.columns([2, 0.5, 3, 3], gap="large") + +with col1: + st.header("Major Oxide Composition") + with st.expander("Major Oxide Composition", expanded=False): + sio2 = st.number_input("SiO2", min_value=0.0, max_value=100.0, value=47.95) + tio2 = st.number_input("TiO2", min_value=0.0, max_value=100.0, value=1.67) + al2o3 = st.number_input("Al2O3", min_value=0.0, max_value=100.0, value=17.32) + feo = st.number_input("FeO", min_value=0.0, max_value=100.0, value=10.24) + # set this to zero and greyed out if fo2 is set + fe2o3 = st.number_input("Fe2O3", min_value=0.0, max_value=100.0, value=0.0) + mno = st.number_input("MnO", min_value=0.0, max_value=100.0, value=0.17) + mgo = st.number_input("MgO", min_value=0.0, max_value=100.0, value=5.76) + cao = st.number_input("CaO", min_value=0.0, max_value=100.0, value=10.93) + na2o = st.number_input("Na2O", min_value=0.0, max_value=100.0, value=3.45) + k2o = st.number_input("K2O", min_value=0.0, max_value=100.0, value=1.99) + p2o5 = st.number_input("P2O5", min_value=0.0, max_value=100.0, value=0.51) + nio = st.number_input("NiO", min_value=0.0, max_value=100.0, value=0.0) + cr2o3 = st.number_input("Cr2O3", min_value=0.0, max_value=100.0, value=0.0) + + oxides = { + "SIO2": sio2, + "TIO2": tio2, + "AL2O3": al2o3, + "FEO": feo, + "FE2O3": fe2o3, + "MNO": mno, + "MGO": mgo, + "CAO": cao, + "NA2O": na2o, + "K2O": k2o, + "P2O5": p2o5, + "NIO": nio, + "CR2O3": cr2o3, + } + + st.subheader("Starting conditions") + with st.expander("Starting conditions", expanded=True): + t_start = st.number_input( + "Temperature (K)", min_value=600.0, max_value=2000.0, value=1473.15 + ) + # work out how to grey this out if you don't need to select a starting pressure + p_start = st.number_input( + "Starting Pressure (bar)", min_value=1.0, max_value=10000.0, value=3000.0 + ) + p_stop = st.number_input( + "Final Pressure (bar)", min_value=1e-5, max_value=10000.0, value=1.0 + ) + dp_max = st.number_input( + "Maximum pressure step size (bar)", + min_value=1.0, + max_value=100.0, + value=100.0, + ) + dp_min = st.number_input( + "Minimum pressure step size (bar)", + min_value=1e-5, + max_value=10.0, + value=0.1, + ) + # grey this one out too + wgt = st.number_input( + "Gas weight fraction at starting pressure", + min_value=0.000000001, + max_value=1.0, + value=0.00001, + ) + mass = st.number_input( + "Mass of the magma packet (g)", min_value=0.0, value=100.0 + ) + + env_options.update( + { + "T_START": t_start, + "P_START": p_start, + "P_STOP": p_stop, + "DP_MAX": dp_max, + "DP_MIN": dp_min, + "WgT": wgt, + "MASS": mass, + } + ) + +with col2: + st.header("Model setup options") + rock_type = st.selectbox("Magma type", ("basalt", "phonolite", "rhyolite")) + run_type = st.selectbox( + "Run Type - Closed System or Open System", ("closed", "open") + ) + if run_type == "Open system": + loss_frac = st.number_input( + "Fraction of gas phase lost at each P step during open system degassing", + min_value=0.0, + max_value=0.9999, + value=0.99, + ) + else: + loss_frac = 0.0 + + single_step = st.checkbox("Run for a single pressure?") + + env_options.update( + {"COMPOSITION": rock_type, "RUN_TYPE": run_type, "SINGLE_STEP": single_step} + ) + + saturation = st.radio( + "Choose how to find the starting pressure:", + ( + "Set manually", + "Find volatile saturation pressure", + "Set volatile element mass, and find volatile saturation pressure", + ), + ) + if saturation == "Set manually": + env_options["FIND_SATURATION"] = False + env_options["ATOMIC_MASS_SET"] = False + elif saturation == "Find volatile saturation pressure": + env_options["FIND_SATURATION"] = True + env_options["ATOMIC_MASS_SET"] = False + else: + env_options["FIND_SATURATION"] = False + env_options["ATOMIC_MASS_SET"] = True + + # Eventually use this to control which volatile input options are shown + gas_sys = st.selectbox( + "Select the volatile element system:", ("OH", "COH", "SOH", "COHS", "COHSN") + ) + fe_equil = st.checkbox("Allow redox equilibration between gas and melt?") + s_sat_warn = st.checkbox( + "Stop the run if sulfide saturation is reached? " + + "(EVo is not suitable for use under sulfide-saturated conditions)" + ) + + env_options.update( + { + "GAS_SYS": gas_sys, + "FE_SYSTEM": fe_equil, + "S_SAT_WARN": s_sat_warn, + } + ) + + sol_laws = st.checkbox( + "Use default solubility laws? (if unchecked, dropdown options will be offered)", + value=True, + key="default_sollaws", + ) + # These will be greyed out if default solubility laws are chosen + + with st.expander("Solubility laws", expanded=not (sol_laws)): + fo2_model = st.selectbox( + "Ferric/Ferrous -> fO2 conversions model", + ("kc1991", "righter2015"), + disabled=st.session_state.default_sollaws, + ) + h2o_model = st.selectbox( + "H2O solubility law", + ("burguisser2015",), + disabled=st.session_state.default_sollaws, + ) + h2_model = st.selectbox( + "H2 solubility law", + ("gaillard2003", "burguisser2015"), + disabled=st.session_state.default_sollaws, + ) + co2_model = st.selectbox( + "CO2 solubility law", + ("burguisser2015", "eguchi2018"), + disabled=st.session_state.default_sollaws, + ) + co_model = st.selectbox( + "CO solubility law", + ("armstrong2015", "None"), + disabled=st.session_state.default_sollaws, + ) + ch4_model = st.selectbox( + "CH4 solubility law", + ("ardia2013", "None"), + disabled=st.session_state.default_sollaws, + ) + sulfide_model = st.selectbox( + "sulfide capacity law", + ("oneill2020", "oneill2002"), + disabled=st.session_state.default_sollaws, + ) + n2_model = st.selectbox( + "N2 solubility law", + ("libourel2003",), + disabled=st.session_state.default_sollaws, + ) + + if sol_laws is True: + models = default_laws + else: + models = { + "FO2_MODEL": fo2_model, + "H2_MODEL": h2_model, + "C_MODEL": co2_model, + "CO_MODEL": co_model, + "CH4_MODEL": ch4_model, + "SULFIDE_CAPACITY": sulfide_model, + } + +with col3: + st.header("Volatile content") + buffer_type = st.selectbox( + "Select the mineral buffer for the oxygen fugacity:", ("FMQ", "IW", "NNO") + ) + fo2 = st.number_input("Starting fO2 relative to the mineral buffer:", value=0.0) + + volatile_params["FO2_buffer"] = buffer_type + volatile_params["FO2_buffer_START"] = fo2 + + vol_setup = st.radio( + "Choose a way to set the system volatile content:", + ("Set gas phase", "Set melt content", "Set atomic weight fractions"), + ) + + if vol_setup == "Set gas phase": + st.subheader("Set gas phase fugacities") + fh2_set = st.checkbox("Set fH2 (bar)", value=False) + fh2_start = st.number_input("fh2", value=0.0, label_visibility="collapsed") + + fh2o_set = st.checkbox("Set fH2O (bar)", value=False) + fh2o_start = st.number_input("fh2o", value=0.0, label_visibility="collapsed") + + fco2_set = st.checkbox("Set fCO2 (bar)", value=False) + fco2_start = st.number_input("fco2", value=0.0, label_visibility="collapsed") + + volatile_params.update( + { + "FH2_SET": fh2_set, + "FH2_START": fh2_start, + "FH2O_SET": fh2_set, + "FH2O_START": fh2o_start, + "FCO2_SET": fco2_set, + "FCO2_START": fco2_start, + } + ) + + for param in [ + "wth2o_set", + "wtco2_set", + "sulfur_set", + "nitrogen_set", + "graphite_saturated", + ]: + volatile_params[param.upper()] = False + + for param in [ + "wth2o_start", + "wtco2_start", + "sulfur_start", + "nitrogen_start", + "graphite_start", + "atomic_h", + "atomic_c", + "atomic_s", + "atomic_n", + ]: + volatile_params[param.upper()] = 0.0 + + elif vol_setup == "Set melt content": + st.subheader("Set melt volatile content") + wth2o_set = st.checkbox("Set melt H2O content (wt fraction)", value=False) + wth2o_start = st.number_input("wth2o", value=0.0, label_visibility="collapsed") + + wtco2_set = st.checkbox("Set melt CO2 content (wt fraction)", value=False) + wtco2_start = st.number_input("wtco2", value=0.0, label_visibility="collapsed") + + sulfur_set = st.checkbox("Set melt sulfur content (wt fraction)", value=False) + sulfur_start = st.number_input( + "sulfur", value=0.0, label_visibility="collapsed" + ) + + nitrogen_set = st.checkbox("Set melt N2 content (wt fraction)", value=False) + nitrogen_start = st.number_input( + "nitrogen", value=0.0, label_visibility="collapsed" + ) + + graphite_saturated = st.checkbox("Is the melt graphite saturated?", value=False) + graphite_start = st.number_input( + "Weight fraction of graphite in the melt", value=0.0 + ) + + volatile_params.update( + { + "WTH2O_SET": wth2o_set, + "WTH2O_START": wth2o_start, + "WTCO2_SET": wtco2_set, + "WTCO2_START": wtco2_start, + "SULFUR_SET": sulfur_set, + "SULFUR_START": sulfur_start, + "NITROGEN_SET": nitrogen_set, + "NITROGEN_START": nitrogen_start, + "GRAPHITE_SATURATED": graphite_saturated, + "GRAPHITE_START": graphite_start, + } + ) + + for param in ["fh2_set", "fh2o_set", "fco2_set"]: + volatile_params[param.upper()] = False + + for param in [ + "fh2_start", + "fh2o_start", + "fco2_start", + "atomic_h", + "atomic_c", + "atomic_s", + "atomic_n", + ]: + volatile_params[param.upper()] = 0.0 + + elif vol_setup == "Set atomic weight fractions": + st.subheader( + "Set the mass fraction of each volatile element in the system (ppm)" + ) + + atomic_h = st.number_input("Weight fraction of H in the system", value=0.0) + atomic_c = st.number_input("Weight fraction of C in the system", value=0.0) + atomic_s = st.number_input("Weight fraction of S in the system", value=0.0) + atomic_n = st.number_input("Weight fraction of N in the system", value=0.0) + + volatile_params.update( + { + "ATOMIC_H": atomic_h, + "ATOMIC_C": atomic_c, + "ATOMIC_S": atomic_s, + "ATOMIC_N": atomic_n, + } + ) + + for param in [ + "fh2_set", + "fh2o_set", + "fco2_set", + "wth2o_set", + "wtco2_set", + "sulfur_set", + "nitrogen_set", + "graphite_saturated", + ]: + volatile_params[param.upper()] = False + + for param in [ + "fh2_start", + "fh2o_start", + "fco2_start", + "wth2o_start", + "wtco2_start", + "sulfur_start", + "nitrogen_start", + "graphite_start", + ]: + volatile_params[param.upper()] = 0.0 + +_, col4, _, col5, _ = st.columns([1, 1, 1, 1, 1]) + +with col4: + load_data = st.button("Load input data", use_container_width=True) + if load_data: + st.session_state.no_loaded_data = False + filepath = filepath = Path(__file__).parent + amend_chem(filepath / "chem.yaml", oxides) + env_params = env_options | models | volatile_params + amend_env(filepath / "env.yaml", env_params) # update to get all options + +with col5: + # Want to grey this out as an option until 'load data' has been pressed... + run_evo = st.button( + "Run EVo", + type="primary", + disabled=st.session_state.no_loaded_data, + use_container_width=True, + ) + if run_evo: + st.session_state.no_loaded_data = True + filepath = filepath = Path(__file__).parent + df = evo.main( + filepath / "chem.yaml", filepath / "env.yaml", filepath / "output.yaml" + ) + +if run_evo: + st.header("Results:") + st.dataframe(df.style.format("{:.3e}")) + st.download_button( + "Download results", df.to_csv(index=False), "results.csv", "text/csv" + ) + + col1, col2 = st.columns(2, gap="large") + + with col1: + colors = px.colors.qualitative.Prism + fig = go.Figure() + fig.add_traces( + go.Scatter( + x=df["mH2O"], + y=df["P"], + mode="lines", + line=dict(color=colors[0]), + name="H2O", + ) + ) + fig.add_traces( + go.Scatter( + x=df["mH2"], + y=df["P"], + mode="lines", + line=dict(color=colors[1]), + name="H2", + ) + ) + fig.add_traces( + go.Scatter( + x=df["mO2"], + y=df["P"], + mode="lines", + line=dict(color=colors[2]), + name="O2", + ) + ) + fig.add_traces( + go.Scatter( + x=df["mCO2"], + y=df["P"], + mode="lines", + line=dict(color=colors[3]), + name="CO2", + ) + ) + fig.add_traces( + go.Scatter( + x=df["mCO"], + y=df["P"], + mode="lines", + line=dict(color=colors[4]), + name="CO", + ) + ) + fig.add_traces( + go.Scatter( + x=df["mCH4"], + y=df["P"], + mode="lines", + line=dict(color=colors[5]), + name="CH4", + ) + ) + fig.add_traces( + go.Scatter( + x=df["mH2S"], + y=df["P"], + mode="lines", + line=dict(color=colors[6]), + name="H2S", + ) + ) + fig.add_traces( + go.Scatter( + x=df["mSO2"], + y=df["P"], + mode="lines", + line=dict(color=colors[7]), + name="SO2", + ) + ) + fig.add_traces( + go.Scatter( + x=df["mS2"], + y=df["P"], + mode="lines", + line=dict(color=colors[8]), + name="S2", + ) + ) + if not (df["mN2"] == 0).all(): + fig.add_traces( + go.Scatter( + x=df["mN2"], + y=df["P"], + mode="lines", + line=dict(color=colors[9]), + name="N2", + ) + ) + fig.update_xaxes(type="log") + fig.update_yaxes(autorange="reversed") + + fig.update_layout( + title="Gas composition vs pressure", + xaxis_title="Gas composition (mol fraction)", + yaxis_title="Pressure (bar)", + legend_title="Species", + ) + + st.plotly_chart(fig) + + with col2: + colors = px.colors.qualitative.Plotly + fig = go.Figure() + fig.add_traces( + go.Scatter( + x=df["Exsol_vol%"], + y=df["P"], + mode="lines", + line=dict(color=colors[0]), + ) + ) + + fig.update_xaxes(type="log") + fig.update_yaxes(autorange="reversed") + + fig.update_layout( + title="Gas volume vs pressure", + xaxis_title="Gas volume %", + yaxis_title="Pressure (bar)", + ) + + st.plotly_chart(fig)