6a. Solubility constant
You can add your own parameterisation to any of the model dependent variables.
Here, we show an example of adding a new parameterisation for the H2S solubility constant (CH2S, function name: C_H2S).
This is what the function for the current solubility constant for H2S looks like in model_dependent_variables.py:
###################################
# solubility constant for H2S #####
###################################
# C_H2S = wmH2S/fH2S (ppm H2S, fH2S bar)
def C_H2S(PT, melt_wf, models=default_models):
"""
Solubility constant for disolving H2S in the melt: C_H2S = wmH2S/fH2S (ppmw/bar)
Parameters
----------
PT: dict
Dictionary of pressure-temperature conditions
pressure (bars) as "P"
temperature ('C) as "T"
melt_wf: dict
Dictionary of melt composition (SiO2, TiO2, etc.)
Not normally used unless model option requires melt composition.
models: pandas.DataFrame
Minimum requirement is dataframe with index of "hydrogen sulfide" and column
label of "option"
Returns
-------
float
Solubility constant for H2S as <class 'mpfr'>
Model options for hydrogen sulfide
-------------
- 'Basalt_Hughes24' [default] Fig.S6 from Hughes et al. (2024) based on experimental
data Moune et al. (2009) and calculations in Lesne et al. (2011)
- 'BasalticAndesite_Hughes24' Fig.S6 from Hughes et al. (2024) based on experimental
data Moune et al. (2009) and calculations in Lesne et al. (2011)
"""
model = models.loc["hydrogen sulfide", "option"]
# fitted to basalt data from Moune+ 2009 CMP 157:691–707 and Lesne+ 2015 ChemGeol
# 418:104–116
if model == "Basalt_Hughes24":
K = 10.23
# fitted to basaltic andesite data from Moune+ 2009 CMP 157:691–707 and Lesne+ 2015
# ChemGeol 418:104–116
elif model == "BasalticAndesite_Hughes24":
K = 6.82
return K
In the current function, there are two options for H:sub:`2`S: solubility basalt and basaltic andesite, using data from Moune et al. (2009) and Lesne et al. (2015) fitted in Hughes et al. (2024).
Let us say you have a new parameterisation that adheres to having wmH2S in ppm but fH2S is in kbar (this example is completely made up). It depends on pressure, temperature, and melt composition using the form:
lnC = lnC0 – ∆V(P - P0)/RT + ∆H/R(1/T0 – 1/T)
where ∆V = 10 cm3/mol, P0 = 1 bar, ∆H = -13 kJ/mol, T0 = 1400 K, and C0 depends on the composition of the melt through the anhydrous cation mole fractions according to:
C0 = AxSiO2 + BxMgO + DxK2O
where A = 2.5, B = 6.8, and D = -9.0.
The anhydrous cation fraction is calculated using the function melt_cation_proportion(run,melt_wf,setup,species,"no","no") in the melt_gas.py file: the first “no” means it does not include any volatiles and the second “no” means it ignores iron speciation. Similarly (i.e., with the same arguments and in melt_gas.py), melt_normalise_wf() calculates the normalised weight fraction, melt_mole_fraction() the mole fractions, and melt_single_O() the mole fraction on a single oxygen basis.
Below is how the new parameterisation is added to C_H2S, called “my new amazing parameterisation”:
###################################
# solubility constant for H2S #####
###################################
# C_H2S = wmH2S/fH2S (ppm H2S, fH2S bar)
def C_H2S(PT, melt_wf, models=default_models):
"""
Solubility constant for disolving H2S in the melt: C_H2S = wmH2S/fH2S (ppmw/bar)
Parameters
----------
PT: dict
Dictionary of pressure-temperature conditions
pressure (bars) as "P"
temperature ('C) as "T"
melt_wf: dict
Dictionary of melt composition (SiO2, TiO2, etc.)
Not normally used unless model option requires melt composition.
models: pandas.DataFrame
Minimum requirement is dataframe with index of "hydrogen sulfide" and column
label of "option"
Returns
-------
float
Solubility constant for H2S as <class 'mpfr'>
Model options for hydrogen sulfide
-------------
- 'Basalt_Hughes24' [default] Fig.S6 from Hughes et al. (2024) based on experimental
data Moune et al. (2009) and calculations in Lesne et al. (2011)
- 'BasalticAndesite_Hughes24' Fig.S6 from Hughes et al. (2024) based on experimental
data Moune et al. (2009) and calculations in Lesne et al. (2011)
"""
model = models.loc["hydrogen sulfide", "option"]
# fitted to basalt data from Moune+ 2009 CMP 157:691–707 and Lesne+ 2015 ChemGeol
# 418:104–116
if model == "Basalt_Hughes24":
K = 10.23
# fitted to basaltic andesite data from Moune+ 2009 CMP 157:691–707 and Lesne+ 2015
# ChemGeol 418:104–116
elif model == "BasalticAndesite_Hughes24":
K = 6.82
elif model == "my new amazing parameterisation":
P, T_K = PT['P'], PT['T']+273.15 # bars, K # extracts pressure in bars and temperature in ‘C from PT and convert temperature to K.
tot, Si, Ti, Al, FeT, Fe2, Fe3, Mg, Mn, Ca, Na, K, P_, H, C = mg.melt_cation_proportion(run,melt_wf,setup,species,"no","no") # calculates cation mole fractions from the melt composition without volatiles in the melt and does not consider iron speciation (i.e., all Fe is FeO).
A, B, D = 2.5, 6.8, -9. # sets the constants A, B, and C
C0 = A*Si + B*Mg + D*K # calculates the compositional term from the cation mole fractions
DV = 10. # cm3/mol # set volume change term
DH = -13. # kJ/mol # sets enthalpy change term
P0 = 1. # bar # sets reference pressure
T0 = 1400 # K # sets reference temperature
R = 83.15 # gives value for R the gas constant
B = -((DV/(R*T_K))*(P-P0)) + (DH/R)*((1.0/T0) - (1.0/T_K)) # calculates PT term
K = C0*gp.exp(B) # calculates solubility constant in ppmw/kbar
K = K/1000. # converts fH2S in kbars to fH2S in bars
return K
To use it, the “hydrogen sulfide” option in “models” would be set to “my new amazing parameterisation”.