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

Python Notebook Download

Benchmarking models for fugacity coefficients

This notebook benchmarks the models for fugacity coefficients in VolFe where possible.

Functions from Sulfur_X (Ding et al., 2022) and EVo (Liggins et al., 2022) have been copied into “function_from_other_tools.py” to allow comparison.

Python set-up

[1]:

import pandas as pd
import numpy as np
import VolFe as vf
import math
import matplotlib.pyplot as plt
import function_from_other_tools as of

Models for fugacity coefficient of O2

options = y_O2, function = y_O2

‘ideal’ Treat as ideal gas, y = 1 at all P.

[2]:

my_models = [["y_O2", "ideal"]]
my_models = vf.make_df_and_add_model_defaults(my_models)

PT = {"P":500.} # bar
PT["T"]=500. # 'C
vf.y_O2(PT,models=my_models)

[2]:

1.0

‘Shi92’ Shi & Saxena (1992) AmMin 77(9-10):1038-1049

Comparison with results from EVo (Liggins et al., 2022) - values match perfectly.

Calculated with EVo:

[3]:

# Conditions
PT = {"P":1000.} # bar
PT["T"]=1200. # C

[4]:

# EVo
species_list = ["O2"]
of.find_Y(PT['P'], PT["T"]+273.15, species_list)

[4]:

(0, 1.209857826330598, 0, 0, 0, 0, 0, 0, 0, 0, 0)

[5]:

# VolFe
my_models = [["y_O2", "Shi92"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.y_O2(PT,models=my_models)

[5]:

1.209857826330598

Models for fugacity coefficient of H2

options = y_H2, function = y_H2

‘Shaw64’ Eq. (4) from Shaw & Wones (1964)

Note: used a value of -0.011901 instead of -0.11901 as reported for calculation of C3 to match data in Table 2.

Table 3 in Shaw & Wone’s contains values for the fugacity coefficient for various P and T.

We spot check ten values here (using values for bars not atmospheres) - they all match to three decimal places.

y = 1.133 (P = 500 bar, T = 500)

y = 1.085 (P = 500 bar, T = 850)

y = 1.073 (P = 500 bar, T = 1000)

y = 1.481 (P = 1600 bar, T = 500)

y = 1.296 (P = 1600 bar, T = 850)

y = 1.252 (P = 1600 bar, T = 1000)

y = 1.749 (P = 2300 bar, T = 500)

y = 1.378 (P = 2300 bar, T = 1000)

y = 1.615 (P = 3000 bar, T = 850)

y = 1.515 (P = 3000 bar, T = 1000)

[6]:

my_models = [["y_H2", "Shaw64"]]
my_models = vf.make_df_and_add_model_defaults(my_models)

PT = {"P":500.} # bar
PT["T"]=500. # 'C
y = vf.y_H2(PT,models=my_models)
print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')
PT["T"]=850. # 'C
y = vf.y_H2(PT,models=my_models)
print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')
PT["T"]=1000. # 'C
y = vf.y_H2(PT,models=my_models)
print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')
PT = {"P":1600.} # bar
PT["T"]=500. # 'C
y = vf.y_H2(PT,models=my_models)
print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')
PT["T"]=850. # 'C
y = vf.y_H2(PT,models=my_models)
print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')
PT["T"]=1000. # 'C
y = vf.y_H2(PT,models=my_models)
print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')
PT = {"P":2300.} # bar
PT["T"]=500. # 'C
y = vf.y_H2(PT,models=my_models)
print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')
PT["T"]=1000. # 'C
y = vf.y_H2(PT,models=my_models)
print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')
PT = {"P":3000.} # bar
PT["T"]=850. # 'C
y = vf.y_H2(PT,models=my_models)
print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')
PT["T"]=1000. # 'C
y = vf.y_H2(PT,models=my_models)
print('y = ', y, '(P =',PT['P'],'bars, T =', PT['T'],')')

y =  1.1327128811340401 (P = 500.0 bars, T = 500.0 )
y =  1.0853513685035343 (P = 500.0 bars, T = 850.0 )
y =  1.0733936356394882 (P = 500.0 bars, T = 1000.0 )
y =  1.4812384375778227 (P = 1600.0 bars, T = 500.0 )
y =  1.295939045430643 (P = 1600.0 bars, T = 850.0 )
y =  1.2516574282929442 (P = 1600.0 bars, T = 1000.0 )
y =  1.7496048740312717 (P = 2300.0 bars, T = 500.0 )
y =  1.3780736491913488 (P = 2300.0 bars, T = 1000.0 )
y =  1.6148638426007402 (P = 3000.0 bars, T = 850.0 )
y =  1.5154237931859629 (P = 3000.0 bars, T = 1000.0 )

‘ideal’ Treat as ideal gas, y = 1 at all P.

[7]:

my_models = [["y_H2", "ideal"]]
my_models = vf.make_df_and_add_model_defaults(my_models)

PT = {"P":500.} # bar
PT["T"]=500. # 'C
vf.y_H2(PT,models=my_models)

[7]:

1.0

Models for fugacity coefficient of S2

options = y_S2, function = y_S2

‘ideal’ Treat as ideal gas, y = 1 at all P.

[8]:

my_models = [["y_S2", "ideal"]]
my_models = vf.make_df_and_add_model_defaults(my_models)

PT = {"P":500.} # bar
PT["T"]=500. # 'C
vf.y_S2(PT,models=my_models)

[8]:

1.0

‘Shi92’ Shi & Saxena (1992) AmMin 77(9-10):1038-1049

Comparison with results from EVo (Liggins et al., 2022) - values match perfectly.

[9]:

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

[10]:

# EVo
species_list = ["S2"]
of.find_Y(PT['P'], PT["T"]+273.15, species_list)

[10]:

(0, 0, 0, 0, 0, 0, 1.1913268980917349, 0, 0, 0, 0)

[11]:

# VolFe
my_models = [["y_S2", "Shi92"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.y_S2(PT,models=my_models)

[11]:

1.1913268980917349

Models for fugacity coefficient of CO

options = y_CO, function = y_CO

‘ideal’ Treat as ideal gas, y = 1 at all P.

[12]:

my_models = [["y_CO", "ideal"]]
my_models = vf.make_df_and_add_model_defaults(my_models)

PT = {"P":500.} # bar
PT["T"]=500. # 'C
vf.y_CO(PT,models=my_models)

[12]:

1.0

‘Shi92’ Shi & Saxena (1992) AmMin 77(9-10):1038-1049

Comparison with results from EVo (Liggins et al., 2022) - values match perfectly.

[13]:

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

[14]:

# EVo
species_list = ["CO"]
of.find_Y(PT['P'], PT["T"]+273.15, species_list)

[14]:

(0, 0, 0, 1.2675079568495566, 0, 0, 0, 0, 0, 0, 0)

[15]:

# VolFe
my_models = [["y_CO", "Shi92"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.y_CO(PT,models=my_models)

[15]:

1.2675079568495566

Models for fugacity coefficient of H2O

options = y_H2O, function = y_H2O

‘ideal’ Treat as ideal gas, y = 1 at all P.

[16]:

my_models = [["y_H2O", "ideal"]]
my_models = vf.make_df_and_add_model_defaults(my_models)

PT = {"P":500.} # bar
PT["T"]=500. # 'C
vf.y_H2O(PT,models=my_models)

[16]:

1.0

‘Holland91’ Holland & Powell (1991) CMP 109:265-273

Comparison with results from EVo (Liggins et al., 2022) - values match perfectly.

[17]:

# Conditions
PT = {"P":1000.} # bar
PT["T"]=1200. # C

[18]:

# EVo
species_list = ["H2O"]
of.find_Y(PT['P'], PT["T"]+273.15, species_list)

[18]:

(0.9988859662849531, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)

[19]:

# VolFe
my_models = [["y_H2O", "Holland91"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.y_H2O(PT,models=my_models)

[19]:

0.9988859662849531

[20]:

# Conditions
PT = {"P":5000.} # bar
PT["T"]=1000. # C
my_models = [["y_H2O", "Holland91"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.gas_molar_volume('H2O', PT, my_models)

[20]:

2.9499443446696176

Models for fugacity coefficient of CO2

options = y_CO2, function = y_CO2

‘Shi92’ Shi & Saxena (1992)

Nothing available in original paper to benchmark too.

‘Holland91’ Holland & Powell (1991) CMP 109:265-273 10.1007/BF00306484

Nothing available for Eq. (8,9) and Table 2 or Eq. (4,A1-3) and Table 2.

Comparison with results from EVo (Liggins et al., 2022) - values match perfectly for Eq. (8) and Table 1.

[21]:

# Conditions
PT = {"P":1000.} # bar
PT["T"]=1200. # C

[22]:

# EVo
species_list = ["CO2"]
of.find_Y(PT['P'], PT["T"]+273.15, species_list)

[22]:

(0, 0, 0, 0, 1.301912410138897, 0, 0, 0, 0, 0, 0)

[23]:

# VolFe
my_models = [["y_CO2", "Holland91_eq8_tab1"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.y_CO2(PT,models=my_models)

[23]:

1.301912410138897

‘ideal’ Treat as ideal gas, y = 1 at all P.

[24]:

my_models = [["y_CO2", "ideal"]]
my_models = vf.make_df_and_add_model_defaults(my_models)

PT = {"P":500.} # bar
PT["T"]=500. # 'C
vf.y_CO2(PT,models=my_models)

[24]:

1.0

Models for fugacity coefficient of SO2

options = y_SO2, function = y_SO2

‘ideal’ Treat as ideal gas, y = 1 at all P.

[25]:

my_models = [["y_SO2", "ideal"]]
my_models = vf.make_df_and_add_model_defaults(my_models)

PT = {"P":500.} # bar
PT["T"]=500. # 'C
vf.y_SO2(PT,models=my_models)

[25]:

1.0

‘Shi92’ Shi & Saxena (1992) AmMin 77(9-10):1038-1049

Comparison with results from Sulfur_X (Ding et al., 2023) - values match perfectly.

Calculated using Sulfur_X function:

[26]:

PT = {"P":1000.} # bar
PT["T"]=1200.+273.15 # K
of.phiso2(PT)

[26]:

1.2064562376797612

Calculated using VolFe:

[27]:

my_models = [["y_SO2", "Shi92"]]
my_models = vf.make_df_and_add_model_defaults(my_models)

PT = {"P":1000.} # bar
PT["T"]=1200. # 'C
vf.y_SO2(PT,models=my_models)

[27]:

1.206456237679761

‘Shi92_Hughes23’ Fig.S1 modified from Shi & Saxena (1992) from Hughes et al. (2023) JGSL 180(3)

This is the first time this model has been coded up so VolFe represents the benchmark for future use.

Models for fugacity coefficient of H2S

options = y_H2S, function = y_H2S

‘ideal’ Treat as ideal gas, y = 1 at all P.

[28]:

my_models = [["y_H2S", "ideal"]]
my_models = vf.make_df_and_add_model_defaults(my_models)

PT = {"P":500.} # bar
PT["T"]=500. # 'C
vf.y_H2S(PT,models=my_models)

[28]:

1.0

‘Shi92’ Shi & Saxena (1992) AmMin 77(9-10):1038-1049

Comparison with results from EVo (Liggins et al., 2022) - values match perfectly.

[29]:

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

[30]:

# EVo
species_list = ["H2S"]
of.find_Y(PT['P'], PT["T"]+273.15, species_list)

[30]:

(0, 0, 0, 0, 0, 0, 0, 0, 1.489941993252676, 0, 0)

[31]:

# VolFe
my_models = [["y_H2S", "Shi92"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.y_H2S(PT,models=my_models)

[31]:

1.489941993252676

‘Shi92_Hughes23’ Fig.S1 modified from Shi & Saxena (1992) Hughes et al. (2024) AmMin 109(3):422-438 doi:10.2138/am-2023-8739

This is the first time this model has been coded up so VolFe represents the benchmark for future use.

Models for fugacity coefficient of CH4

options = y_CH4, function = y_CH4

‘ideal’ Treat as ideal gas, y = 1 at all P.

[32]:

my_models = [["y_CH4", "ideal"]]
my_models = vf.make_df_and_add_model_defaults(my_models)

PT = {"P":500.} # bar
PT["T"]=500. # 'C
vf.y_CH4(PT,models=my_models)

[32]:

1.0

‘Shi92’ Shi & Saxena (1992) AmMin 77(9-10):1038-1049

Comparison with results from EVo (Liggins et al., 2022) - values match perfectly.

[33]:

# conditions
PT = {"P":1000.} # bar
PT["T"]=1200. # C

[34]:

# EVo
species_list = ["CH4"]
of.find_Y(PT['P'], PT["T"]+273.15, species_list)

[34]:

(0, 0, 0, 0, 0, 1.2858217235364475, 0, 0, 0, 0, 0)

[35]:

my_models = [["y_CH4", "Shi92"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.y_CH4(PT,models=my_models)

[35]:

1.2858217235364475

Models for fugacity coefficient of OCS

options = y_OCS, function = y_OCS

‘ideal’ Treat as ideal gas, y = 1 at all P.

[36]:

my_models = [["y_OCS", "ideal"]]
my_models = vf.make_df_and_add_model_defaults(my_models)

PT = {"P":500.} # bar
PT["T"]=500. # 'C
vf.y_OCS(PT,models=my_models)

[36]:

1.0

‘Shi92’ Shi & Saxena (1992) AmMin 77(9-10):1038-1049

Comparison with results from EVo (Liggins et al., 2022) - values match perfectly.

[37]:

# conditions
PT = {"P":1000.} # bar
PT["T"]=1200 # K

[38]:

# EVo
species_list = ["OCS"]
of.find_Y(PT['P'], PT["T"]+273.15, species_list)

[38]:

(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1.3560000299582256)

[39]:

# VolFe
my_models = [["y_OCS", "Shi92"]]
my_models = vf.make_df_and_add_model_defaults(my_models)
vf.y_OCS(PT,models=my_models)

[39]:

1.3560000299582256