This page was generated from docs/benchmarking/solubility_functions.ipynb. Interactive online version: Binder badge.

Python Notebook Download

Benchmarking models for solubility functions

This notebook benchmarks the models for solubility functions in VolFe where possible.

Python set-up

[1]:

import pandas as pd
import VolFe as vf
import math
import matplotlib as plt

[2]:

# VolFe version
vf.__version__

[2]:

'1.0.2'

Composition and conditions used for benchmarking

Average of high-SiO2 pillow-rim glasses in HSDP from Mauna Kea volcano from Brounce et al. (2017).

[3]:

my_analysis = {
        "Sample": "Hawaiian basalt",
        "T_C": 1200.0,  # Temperature in 'C
        "SiO2": 51.29,  # wt%
        "TiO2": 2.50,  # wt%
        "Al2O3": 13.70,  # wt%
        "FeOT": 11.04,  # wt%
        "MnO": 0.02,  # wt%
        "MgO": 6.70,  # wt%
        "CaO": 11.03,  # wt%
        "Na2O": 2.27,  # wt%
        "K2O": 0.43,  # wt%
        "P2O5": 0.,  # wt%
        "H2O": 0.,  # wt%
        "CO2ppm": 0.,  # ppm
        "STppm": 0.,  # ppm
        "Xppm": 0.0,  # ppm
        "Fe3FeT": 0.1}

my_analysis = pd.DataFrame(my_analysis, index=[0])

PT = {"P":1000.}
PT["T"]=1200.

melt_wf=vf.melt_comp(0.,my_analysis)
melt_wf['CO2'] = my_analysis.loc[0.,"CO2ppm"]/1000000.
melt_wf["H2OT"] = my_analysis.loc[0,"H2O"]/100.
melt_wf['ST'] = my_analysis.loc[0.,"STppm"]/1000000.
melt_wf['CT'] = (melt_wf['CO2']/vf.species.loc['CO2','M'])*vf.species.loc['C','M']
melt_wf['HT'] = (melt_wf['H2OT']/vf.species.loc['H2O','M'])*(2.*vf.species.loc['H','M'])
melt_wf['XT'] = my_analysis.loc[0.,"Xppm"]/1000000.
melt_wf["Fe3FeT"] = my_analysis.loc[0.,"Fe3FeT"]

Models for solubility function of hydrogen (H2)

option = hydrogen, function = C_H2

For the following models, this is the first time they have been coded up so VolFe represents the benchmark for future use:

  • ‘Basalt_Hughes24’ Basalt in Table S4 from Hughes et al. (2024) based on experimetnal data from Hirschmann et al. (2012)

  • ‘Andesite_Hughes24’ Andesite in Table S4 from Hughes et al. (2024) based on experimental data from Hirschmann et al. (2012)

These are the only available models currently available in VolFe.

Models for solubility function of water (H2O)

option = water, function = C_H2O

‘Basalt_Hughes24’ Fig.S2 from Hughes et al. (2024)

Hughes EC, Liggins P, Saper L, Stolper EM (2024). The efects of oxygen fugacity and sulfur on the pressure of vapor-saturation of magma. American Mineralogist 109(3):422-438 https://doi.org/10.2138/am-2022-8739

This is a constant, with a value of 4.6114e-6:

[4]:

my_models = [["water", "Basalt_Hughes24"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
(vf.C_H2O(PT,melt_wf,models=my_models))

[4]:

4.6114e-06

‘Rhyolite_Hughes25’ Eq. (S1) from Hughes et al. (2025)

Hughes EC, Liggins P, Wieser P, and Stolper EM (2025). VolFe: an open-source Python package for calculating melt-vapor equilibria including silicate melt, carbon, hydrogen, sulfur, and noble gases. Volcanica 8(2):457-481 https://doi.org/10.30909/vol/imvc1781

Based on data in Fig. 3 of Blank et al. (1993)

This is a constant, with a value of 5.3851e-6:

[5]:

my_models = [["water", "Rhyolite_Hughes25"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
(vf.C_H2O(PT,melt_wf,models=my_models))

[5]:

5.3851e-06

Models for solubility function of carbon monoxide (CO)

option = carbon monoxide, function = C_CO

For the following model, this is the first time they have been coded up so VolFe represents the benchmark for future use:

  • ‘Basalt_Hughes24’ CO in Table S4 from Hughes et al. (2024) based on data from Armstrong et al. (2015), Stanley et al., (2014), and Wetzel et al., (2013)

This is the only available model currently available in VolFe.

Models for solubility function of carbon doxide (CO2)

option = carbon dioxide, function = C_CO3

The following models do not have material available in the original papers for benchmarking:

  • ‘Basalt_Dixon97’ Eq. (7) from Dixon et al. (1997) AmMin 82(3-4)368-378 doi:10.2138/am-1997-3-415

  • ‘NorthArchBasalt_Dixon97’ Eq. (8) from Dixon et al. (1997) AmMin 82(3-4)368-378 doi:10.2138/am-1997-3-415

  • ‘Basalt_Lesne11’ Eq. (25,26) from Lesne et al. (2011) CMP 162:153-168 doi:10.1007/s00410-010-0585-0

  • ‘VesuviusAlkaliBasalt_Lesne11’ VES-9 in Table 4 from Lesne et al. (2011) CMP 162:153-168 doi:10.1007/s00410-010-0585-0

  • ‘EtnaAlkaliBasalt_Lesne11’ ETN-1 in Table 4 from Lesne et al. (2011) CMP 162:153-168 doi:10.1007/s00410-010-0585-0

  • ‘StromboliAlkaliBasalt_Lense11’ PST-9 in Table 4 from Lesne et al. (2011) CMP 162:153-168 doi:10.1007/s00410-010-0585-0

  • ‘Basanite_Holloway94’ Basanite in Table 5 from Holloway and Blank (1994) RiMG 30:187-230 doi:10.1515/9781501509674-012

  • ‘Leucitite_Thibault94’ Leucitite from Thibault & Holloway (1994) CMP 116:216-224 doi:10.1007/BF00310701

‘MORB_Dixon95’ Bullet (5) of summary from Dixon et al. (1995)

Dixonet al. (1995) Journal of Petrology 36(6):1607-1631

In Table 5, Dixon et al. list P, CO2 concentration, and fCO2, from which we can calculate C_CO3 at 1200 ‘C.

P = 1000 bar, fCO2 = 1240 bar, and CO2 = 471 ppm

First we convert ppm CO2 to xmCO2 by rearranging the equation in superscript h of Table 5.

[6]:

def xm(y,M):
    value = ((y*M)/10**4.)/(4400.+((44.*y)/10**4.))
    return value

xm = xm(471,36.596)

Then calculate C_CO3 from xmCO2/fCO2

[7]:

xm/1240.

[7]:

3.1577348476492196e-07

And compare to VolFe

[8]:

my_models = [["carbon dioxide", "MORB_Dixon95"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.C_CO3(PT,melt_wf,models=my_models)

[8]:

3.150061941762063e-07

They agree to two significant figures.

‘Rhyolite_Blank93’ Fig.2 caption from Blank et al. (1993)

Blank et al. (1993) EPSL 119:27-36

The value of xmCO2 (mole fraction) at fCO2 (bar) of 1000 bar at P = 750 bar and T = 850 ‘C using Fig. 2 from Blank et al. (1993) is 0.0004173 (File Blank93_Fig2.png digitised using https://plotdigitizer.com/app).

The C_CO3 calcualted from xmCO2/fCO2 and VolFe agree to one significant figure.

[9]:

# xmCO2/fCO2
0.0004173/1000.

[9]:

4.173e-07

[10]:

my_models = [["carbon dioxide", "Rhyolite_Blank93"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
PT_ = {'P':750.}
PT_['T'] = 850.
vf.C_CO3(PT_,melt_wf,models=my_models)

[10]:

4.2355673190981494e-07

Models for solubility function of methane (CH4)

option = methane, function = C_CH4

The following models do not have material available in the original papers for benchmarking:

  • ‘Basalt_Ardia13’ Eq. (7a) from Ardia et al. (2013)

There are no other models for this solubility function available in VolFe currently.

Models for solubility function of sulfide (*S2-)

option = sulfide, function = C_S

The following models do not have material available in the original papers for benchmarking:

  • ‘Boulliung23_eq6’ Eq. (6) from Boulliung & Wood (2023) CMP 178:56 doi:10.1007/s00410-023-02033-9

  • ‘Boulliung23_eq7_12’ Eq. (7) and (12) [P-effect] from Boulliung & Wood (2023) CMP 178:56 https://doi.org10.1007/s00410-023-02033-9

  • ‘Gorojovsky26_eq8’ Eq. (8) from Gorojovsky, L.R. and Wood, B.J., (2026). Solubility and speciation of sulfur in silicate melts under crustal conditions. EPSL 687:120088 https://doi.org/10.1016/j.epsl.2026.120088

  • ‘Gorojovsky26_eq8_13’ Eq. (8) and (13) from Gorojovsky, L.R. and Wood, B.J., (2026). Solubility and speciation of sulfur in silicate melts under crustal conditions. EPSL 687:120088 https://doi.org/10.1016/j.epsl.2026.120088

  • ‘Gorojovsky26_eq9_13’ Eq. (9) and (13) from Gorojovsky, L.R. and Wood, B.J., (2026). Solubility and speciation of sulfur in silicate melts under crustal conditions. EPSL 687:120088 https://doi.org/10.1016/j.epsl.2026.120088

‘ONeill21’ Eq. (10.34) from O’Neill (2021)

O’Neill HSC (2021). The thermodynamic controls on sulfide saturation in silicate melts with application to ocean floor basalts. Magma redox geochemistry p177-213 https://doi.org/10.1002/9781119473206.ch10

Supplementary spreadsheet (ONeill21.xlsx)

Cell AN14: -2.419

Matches to 3 decimal places: Note spreadsheet uses +273 to convert to K, rather than 273.15 used in VolFe so T in spreadsheet = 1200.15 ‘C

[11]:

my_models = [["sulfide", "ONeill21"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
math.log(vf.C_S(PT,melt_wf,models=my_models))

[11]:

-2.4188118891936226

‘Boulliung23_eq7’ Eq. (7) from Boulliung & Wood (2023)

Boulliung & Wood (2023) CMP 178:56 https://doi.org10.1007/s00410-023-02033-9

Supplementary Information - Supplementary file1 used for benchmarking (Boulliung23.xlsx)

Cell K22: -5.1286877

Matches to 3 decimal places

[12]:

my_models = [["sulfide", "Boulliung23_eq7"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
math.log10(vf.C_S(PT,melt_wf,models=my_models))

[12]:

-5.1289519698860575

‘Thomas26_eq15’ Eq. (15) from Thomas & Wood (2026)

Thomas, R.W. and Wood, B.J., 2026. Sulfur speciation in silicate melts at high pressure. GCA 417:37-51 https://doi.org/10.1016/j.gca.2026.02.003

Supplementary Material - Supplementary Data 2 used for benchmarking (Thomas26.xlsx)

Cell IG14: -5.11

Matches to precision (2 decimal places) given on spreadsheet.

[4]:

my_models = [["sulfide", "Thomas26_eq15"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
math.log10(vf.C_S(PT,melt_wf,models=my_models))

[4]:

-5.1137053079338735

‘Gorojovsky26_eq9’ Eq. (9) from Gorojovsky & Wood (2026)

Gorojovsky, L.R. and Wood, B.J., (2026). Solubility and speciation of sulfur in silicate melts under crustal conditions. EPSL 687:120088 https://doi.org/10.1016/j.epsl.2026.120088

Supplementary Material - S3_Sulfur_fO2_calculator used for benchmarking (Gorojovsky26.xlsx)

Note spreadsheet uses +273 to convert to K, rather than 273.15 used in VolFe so T in spreadsheet = 1200.15 ‘C

Cell GC9: -5.36

Matches to precision (2 decimal places) given on spreadsheet.

[4]:

my_models = [["sulfide", "Gorojovsky26_eq9"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
math.log10(vf.C_S(PT,melt_wf,models=my_models))

[4]:

-5.364628528792665

Models for solubility function of sulfate (S6+)

option = sulfate, function = C_SO4

The following models do not have material available in the original papers for benchmarking:

  • ‘Boulliung22nP’ Eq. (5) from Boulliung & Wood (2023) GCA and Eq. (8) for P from Boulliung & Wood (2022)

  • ‘Boulliung22wP’ Eq. (5) from Boulliung & Wood (2023) GCA and Eq. (8) for P from Boulliung & Wood (2022)

  • ‘Boulliung23_eq9_12’ Eq. (9) and (12) [P-effect] from Boulliung & Wood (2023) CMP 178:56 https://doi.org/10.1007/s00410-023-02033-9

  • ‘Thomas26_eq21’ Eq. (21) from Thomas, R.W. and Wood, B.J., 2026. Sulfur speciation in silicate melts at high pressure. GCA 417:37-51 https://doi.org/10.1016/j.gca.2026.02.003

  • ‘Gorojovsky26_eq10’ Eq. (10) from Gorojovsky, L.R. and Wood, B.J., (2026). Solubility and speciation of sulfur in silicate melts under crustal conditions. EPSL 687:120088 https://doi.org/10.1016/j.epsl.2026.120088

  • ‘Gorojovsky26_eq10_14’ Eq. (10) and (14) from Gorojovsky, L.R. and Wood, B.J., (2026). Solubility and speciation of sulfur in silicate melts under crustal conditions. EPSL 687:120088 https://doi.org/10.1016/j.epsl.2026.120088

  • ‘Gorojovsky26_eq11_14’ Eq. (11) and (14) from Gorojovsky, L.R. and Wood, B.J., (2026). Solubility and speciation of sulfur in silicate melts under crustal conditions. EPSL 687:120088 https://doi.org/10.1016/j.epsl.2026.120088

‘ONeill22’ Eq. (12a) without H2O dilution from O’Neill & Mavrogenes (2022)

Appendix A. Supplementary material - Supplementary data 2. Tab = Table S6 S redox calculator (ONeill22.xlsx),

Cell AH12: 12.95

Matches to 2 decimal places. Note spreadsheet uses +273 to convert to K, rather than 273.15 used in VolFe so T in spreadsheet = 1200.15 ‘C

[13]:

my_models = [["sulfate", "ONeill22"],["KOSg2", "ONeill22"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
math.log(vf.C_SO4(PT,melt_wf,models=my_models)/(vf.KOSg2(PT,models=my_models)))

[13]:

12.954840211020695

‘Boulliung23_eq9’ Eq. (9) from Boulliung et al. (2023)

Supplementary Information - Supplementary file1 used for benchmarking (Boulliung23.xlsx)

Cell K24: 8.73571472

Matches to 2 decimal places

[14]:

my_models = [["sulfate", "Boulliung23_eq9"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
math.log10(vf.C_SO4(PT,melt_wf,models=my_models))

[14]:

8.736100842146033

‘Boulliung23_eq11’ Eq. (11) from Boulliung et al. (2023)

Supplementary Information - Supplementary file1 used for benchmarking (Boulliung23.xlsx)

Cell M24: 8.82627819

Matches to 3 decimal places

[16]:

my_models = [["sulfate", "Boulliung23_eq11"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
math.log10(vf.C_SO4(PT,melt_wf,models=my_models))

[16]:

8.826711528634405

‘Thomas26_eq22’ Eq. (22) from Thomas & Wood (2026)

Thomas, R.W. and Wood, B.J., 2026. Sulfur speciation in silicate melts at high pressure. GCA 417:37-51 https://doi.org/10.1016/j.gca.2026.02.003

Supplementary Material - Supplementary Data 2 used for benchmarking (Thomas26.xlsx)

Cell IG14: 8.68

Matches to precision (2 decimal places) given on spreadsheet.

[4]:

my_models = [["sulfate", "Thomas26_eq22"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
math.log10(vf.C_SO4(PT,melt_wf,models=my_models))

[4]:

8.6842501689359

‘Gorojovsky26_eq11’ Eq. (11) from Gorojovsky & Wood (2026)

Gorojovsky, L.R. and Wood, B.J., (2026). Solubility and speciation of sulfur in silicate melts under crustal conditions. EPSL 687:120088 https://doi.org/10.1016/j.epsl.2026.120088

Supplementary Material - S3_Sulfur_fO2_calculator used for benchmarking (Gorojovsky26.xlsx)

Note spreadsheet uses +273 to convert to K, rather than 273.15 used in VolFe so T in spreadsheet = 1200.15 ‘C

Cell GB9: 8.57

Matches to precision (2 decimal places) given on spreadsheet.

[5]:

my_models = [["sulfate", "Gorojovsky26_eq11"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
math.log10(vf.C_SO4(PT,melt_wf,models=my_models))

[5]:

8.571669080232766

Models for solubility function of hydrogen sulfide (H2S)

option = water, function = C_H2S

‘Basalt_Hughes24’ Fig.S6 from Hughes et al. (2024)

Based on experimental data Moune et al. (2009) and calculations in Lesne et al. (2011) (see file Hughes24_hydrogensulfide.xlxs).

This is a constant with a value of 10.231:

[17]:

my_models = [["hydrogen sulfide", "Basalt_Hughes24"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.C_H2S(PT,melt_wf,models=my_models)

[17]:

10.23

‘BasalticAndesite_Hughes24’ Fig.S6 from Hughes et al. (2024)

Based on experimental data Moune et al. (2009) and calculations in Lesne et al. (2011) (see file Hughes24_hydrogensulfide.xlxs).

This is a constant with a value of 6.8193:

[18]:

my_models = [["hydrogen sulfide", "BasalticAndesite_Hughes24"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.C_H2S(PT,melt_wf,models=my_models)

[18]:

6.82

Models for solubility function of “X”

option = species X solubility, function = C_X

‘Ar_Basalt_Hughes25’ Hughes et al. (2025)

Hughes EC, Liggins P, Wieser P, and Stolper EM (2025). VolFe: an open-source Python package for calculating melt-vapor equilibria including silicate melt, carbon, hydrogen, sulfur, and noble gases. Volcanica 8(2):457-481 https://doi.org/10.30909/vol/imvc1781

Based on data from Iacono-Marziano et al. (2010)

This is a constant with a value of 0.0799:

[19]:

my_models = [["species X solubility", "Ar_Basalt_Hughes25"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.C_X(PT,melt_wf,models=my_models)

[19]:

0.0799

‘Ar_Rhyolite_Hughes25’ Hughes et al. (2025)

Hughes EC, Liggins P, Wieser P, and Stolper EM (2025). VolFe: an open-source Python package for calculating melt-vapor equilibria including silicate melt, carbon, hydrogen, sulfur, and noble gases. Volcanica 8(2):457-481 https://doi.org/10.30909/vol/imvc1781

Based on data from Iacono-Marziano et al. (2010)

This is a constant with a value of 0.4400:

[20]:

my_models = [["species X solubility", "Ar_Rhyolite_Hughes25"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.C_X(PT,melt_wf,models=my_models)

[20]:

0.44

‘Ne_Basalt_Hughes25’ Hughes et al. (2025)

Hughes EC, Liggins P, Wieser P, and Stolper EM (2025). VolFe: an open-source Python package for calculating melt-vapor equilibria including silicate melt, carbon, hydrogen, sulfur, and noble gases. Volcanica 8(2):457-481 https://doi.org/10.30909/vol/imvc1781

Based on data from Iacono-Marziano et al. (2010)

This is a constant with a value of 0.1504:

[21]:

my_models = [["species X solubility", "Ne_Basalt_Hughes25"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.C_X(PT,melt_wf,models=my_models)

[21]:

0.1504

‘Ne_Rhyolite_Hughes25’ Hughes et al. (2025)

Hughes EC, Liggins P, Wieser P, and Stolper EM (2025). VolFe: an open-source Python package for calculating melt-vapor equilibria including silicate melt, carbon, hydrogen, sulfur, and noble gases. Volcanica 8(2):457-481 https://doi.org/10.30909/vol/imvc1781

Based on data from Iacono-Marziano et al. (2010)

This is a constant with a value of 0.8464:

[22]:

my_models = [["species X solubility", "Ne_Rhyolite_Hughes25"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.C_X(PT,melt_wf,models=my_models)

[22]:

0.8464

[user specified number]

User can type a number that will be used instead (i.e., a constant value)

This is user specified - for instance with a value of 10:

[23]:

my_models = [["species X solubility", 10.]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.C_X(PT,melt_wf,models=my_models)

[23]:

10.0

Models for equilibrium functions of CO2 speciation

option = Cspeccomp, function = KCOm

‘Andesite_Botcharnikov06’ Eq. (8) from Botcharnikov et al. (2006)

‘Dacite_Botcharnikov06’ Eq. in the text from Botcharnikov et al. (2006), based on data from Behrens et al. (2004)

Fig. 12 from Botcharnikov et al. (2006) shows the fit to the data used to derive their models for both Andesite and Dacite (see file Botcharnikov06_Fig12.jpg).

As they are straight lines, we read the data at the intersection of the axes

Andesite:

Intersection with y-axis: 1000/T = 0.60 and lnK2 = 0.2

Intersection with x-axis: 1000/T = 1.41 and lnK2 = 7.0

Dacite:

Intersection with y-axis: 1000/T = 0.60 and lnK2 = -2.0

Intersection with x-axis: 1000/T = 1.50 and lnK2 = 7.0

And then calculate at the same temperatures using VolFe - there is good agreement for both models.

[24]:

my_models = [["Cspeccomp", "Andesite_Botcharnikov06"]]
my_models = vf.make_df_and_add_model_defaults(my_models)

PT = {"T":(1000./0.60)-273.15}
T1 = math.log(vf.KCOm(PT,melt_wf,models=my_models))
PT = {"T":(1000./1.41-273.15)}
T2 = math.log(vf.KCOm(PT,melt_wf,models=my_models))

print(T1, T2)

0.0889999999999986 7.107649999999999

[25]:

my_models = [["Cspeccomp", "Dacite_Botcharnikov06"]]
my_models = vf.make_df_and_add_model_defaults(my_models)

PT = {"T":(1000./0.6)-273.15}
T1 = math.log(vf.KCOm(PT,melt_wf,models=my_models))
PT = {"T":(1000./1.5)-273.15}
T2 = math.log(vf.KCOm(PT,melt_wf,models=my_models))

print(T1, T2)

-1.8178000000000019 6.9905

‘Basalt’ Assume all oxidised carbon in the melt is present as carbonate ions.

Value is “infinite”.

[26]:

my_models = [["Cspeccomp", "Basalt"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.KCOm(PT,melt_wf,models=my_models)

[26]:

'infinite'

‘Rhyolite’ Assume all oxidised carbon in the melt is present as molecular CO2.

Value is 0.

[27]:

my_models = [["Cspeccomp", "Rhyolite"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.KCOm(PT,melt_wf,models=my_models)

[27]:

0.0