VapoRock: Equilibrium Outgassing

Demonstrates how to run an equilibrium vaporization routine using VapoRock.

Open this code in an executable MyBinder instance (MyBinder links may be slow to load– please be patient!):

Initialization

Import necessary packages:

import vaporock
import numpy as np

Set Melt Composition, Temp., redox

a simplified BSE composition (Schaefer & Fegley, 2009)

melt_comp_wts = {
    'SiO2':45.97,
    'MgO':36.66,
    'Al2O3':4.77,
    'FeO':8.24,
    'CaO':3.78,
    'Na2O':0.35,
    'K2O':0.04}

T = np.linspace(1500,3000,101)
buffer = 'IW'
dlogfO2 = +1.5
logfO2 = vaporock.redox_buffer(T, buffer=buffer, dlogfO2=dlogfO2)

Initialize system & predict gas abundances

system = vaporock.System()
system.set_melt_comp(melt_comp_wts)
system.print_melt_comp()
logP, lnK, lna =  system.eval_gas_abundances(T, logfO2, full_output=True)

MELTS liquid components:
SiO2         1.790509e-01
TiO2         0.000000e+00
Al2O3        4.635782e-02
Fe2O3        2.220446e-16
MgCr2O4      0.000000e+00
Fe2SiO4      5.734456e-02
MnSi1d2O2    0.000000e+00
Mg2SiO4      4.547891e-01
NiSi1d2O2    0.000000e+00
CoSi1d2O2    0.000000e+00
CaSiO3       6.740443e-02
Na2SiO3      5.647079e-03
KAlSiO4      8.492984e-04
Ca3P2O8      0.000000e+00
H2O          0.000000e+00
dtype: float64

Plot vapor species in different ways

vaporock.plot.partial_pressure(T, logP)

<Axes: xlabel='Temperature [K]', ylabel='Partial Pressures [bar]'>

vaporock.plot.species_mol_frac(T, logP)

<Axes: xlabel='Temperature [K]', ylabel='Species Abundance [mol fraction]'>

mol_elem_frac, Ptotal = system.calc_gas_props(T, logP)
vaporock.plot.element_mol_frac(T, mol_elem_frac)

/home/vscode/.local/lib/python3.10/site-packages/pandas/core/internals/blocks.py:395: RuntimeWarning: divide by zero encountered in log10
  result = func(self.values, **kwargs)

<Axes: xlabel='Temperature [K]', ylabel='Elemental Abundance [mol fraction]'>