Customizing the thermodynamic database
This vignette was compiled on 2024-11-17 with CHNOSZ version 2.1.0-22.
This vignette will cover some topics about using custom thermodynamic
data in CHNOSZ. The two main functions to remember are
add.OBIGT() to add data from a CSV file
and mod.OBIGT() to add or modify data
through a function interface. A third function,
logB.to.OBIGT(), is provided to fit
thermodynamic parameters to experimental formation constants
(log β).
Before describing the methods to add or modify data, some notes on
the basic structure of the database and data entry conventions are
given. Column names (or parts thereof) are formatted in blue
(e.g. formula), and important notes are
highlighted in yellow (NOTE). Function names in CHNOSZ
are colored red for functions that have side effects (including those
that modify the database;
e.g. add.OBIGT()) and green for
functions that don’t have side effects. Note that
info() is used for querying the OBIGT
thermodynamic database, and subcrt()
is the main function in CHNOSZ for calculating standard thermodynamic
properties as a function of temperature and pressure from the parameters
in the database.
Basic structure of OBIGT
OBIGT is the name of the thermodynamic database in CHNOSZ. The data
are distributed in CSV files in the inst/extdata/OBIGT
directory of the CHNOSZ package. When the package is installed, the
files are copied to the exdata/OBIGT directory of the
installed package location. To find out where this is on your computer,
run the following command.
## [1] "/tmp/RtmpViir5k/Rinst688417e0445/CHNOSZ/extdata/OBIGT"The directory path on your computer will be different. Although possible, it is NOT recommended to edit the data files at that location. This is because they will be overwritten by package updates; moreover, it is good practice to keep all the files needed for your project in a project directory.
This lists the files in the installation directory:
## [1] "AD.csv" "AS04.csv" "Berman_cr.csv"
## [4] "DEW.csv" "GEMSFIT.csv" "H2O_aq.csv"
## [7] "SLOP98.csv" "SUPCRT92.csv" "inorganic_aq.csv"
## [10] "inorganic_cr.csv" "inorganic_gas.csv" "organic_aq.csv"
## [13] "organic_cr.csv" "organic_gas.csv" "organic_liq.csv"
## [16] "refs.csv"Some of these files are used to build the default OBIGT database that is created when CHNOSZ starts up. There are also a number of additional data files that have optional datasets. The OBIGT thermodynamic database vignette summarizes the contents of the default and optional data files. The files can also be opened by a spreadsheet program and used as templates for adding data yourself.
thermo()$OBIGT (hereafter, just OBIGT) is the “live”
version of the database that is assembled from the CSV data files when
CHNOSZ starts up or by using the
reset() or
OBIGT() functions. The OBIGT data frame
is stored in an environment named CHNOSZ that is part of
the namespace of the CHNOSZ package. More specifically, it is part of a
list named thermo, which has the OBIGT database and other
parameters and settings used by CHNOSZ.
reset() restores the entire
thermo object to default values;
OBIGT() restores just the OBIGT data
frame. The latter is useful for seeing the effects of changing the
thermodynamic database on on chemical affinities calculated with
affinity(), without changing the
chemical species.
OBIGT can be modified during an R session; if it couldn’t, some of
the examples in this vignette would not be possible! When you quit R, it
offers the option of saving your workspace so it can be reloaded it when
R is restarted. I always say “no” here; my preference is to load data
into a fresh session every time I start R. This “load saved workspace”
feature means that OBIGT might not be the default database in any given
R session. To ensure that this vignette is run using the default
database, we start by running reset()
to reset OBIGT and the other settings used by CHNOSZ.
## reset: resetting "thermo" object## OBIGT: loading default database with 2018 aqueous, 3570 total speciesthermo() is a convenience function
to access or modify parts of the thermo list object. To see
the first few entries in OBIGT, do this:
## name abbrv formula state ref1 ref2 date model E_units G
## 1 water <NA> H2O liq HGK84 JOH92 2006-10-25 H2O cal NA
## 2 e- <NA> (Z-1) aq electron <NA> 2006-10-28 HKF cal 0
## 3 H+ H+ H+ aq proton <NA> 1997-11-06 HKF cal 0
## 4 Li+ Li+ Li+ aq SH88 <NA> 1997-11-06 HKF cal -69933
## 5 Na+ Na+ Na+ aq SH88 <NA> 1997-11-06 HKF cal -62591
## 6 K+ K+ K+ aq SH88 <NA> 1997-11-06 HKF cal -67510
## H S Cp V a1.a a2.b a3.c a4.d c1.e c2.f
## 1 NA NA NA NA NA NA NA NA NA NA
## 2 0 15.6166 0.00 0.00 0.0000 0.000 0.000 0.0000 0.00 0.000
## 3 0 0.0000 0.00 0.00 0.0000 0.000 0.000 0.0000 0.00 0.000
## 4 -66552 2.7000 14.20 -0.87 -0.0237 -0.069 11.580 -2.7761 19.20 -0.240
## 5 -57433 13.9600 9.06 -1.11 1.8390 -2.285 3.256 -2.7260 18.18 -2.981
## 6 -60270 24.1500 1.98 9.06 3.5590 -1.473 5.435 -2.7120 7.40 -1.791
## omega.lambda z.T
## 1 NA NA
## 2 0.0000 0
## 3 0.0000 0
## 4 0.4862 1
## 5 0.3306 1
## 6 0.1927 1Conventions for data entry in OBIGT
The format of OBIGT is described in the CHNOSZ manual: see
thermo(). Next, we point out some
particular conventions including types of data, required and optional
data, order-of-magnitude scaling. Here are the numbered column names for
reference:
## [1] "1 name" "2 abbrv" "3 formula" "4 state"
## [5] "5 ref1" "6 ref2" "7 date" "8 model"
## [9] "9 E_units" "10 G" "11 H" "12 S"
## [13] "13 Cp" "14 V" "15 a1.a" "16 a2.b"
## [17] "17 a3.c" "18 a4.d" "19 c1.e" "20 c2.f"
## [21] "21 omega.lambda" "22 z.T"Types of data
- Columns 1–9 have character data.
- Column 8 is named model; here are
the two models that are most frequent in the default database:
HKF: Revised Helgeson-Kirkham-Flowers “equation-of-state” parameters for aqueous speciesCGL: Heat capacity coefficients for crystalline, gaseous, and liquid species. The first three terms in the CGL heat capacity equation correspond to the Maier-Kelley equation for heat capacity (Maier and Kelley, 1932); the additional terms are useful for representing heat capacities of minerals (Robie and Hemingway, 1995) and organic gases and liquids (Helgeson et al., 1998).
- Column 8 is named model; here are
the two models that are most frequent in the default database:
- Columns 10–22 have numeric data.
- Columns 10–14 have standard-state thermodynamic properties at 25 °C and 1 bar.
- Columns 15–21 have parameters for calculating thermodynamic
properties at other temperatures and pressures.
- The columns are named by combining the the names of the HKF and CGL coefficients, separated by a dot.
- Column 22 has the charge used in the HKF EOS or the maximum
temperature for CGL species.
- NOTE: The value of charge used in the HKF EOS (in particular, the g function for the temperature derivatives of the ω parameter (Shock et al., 1992)) is taken from this column and not from the chemical formula of the species.
Ranges of HKF and CGL models
To a first approximation, the revised HKF equations of state are applicable within the stability region of liquid water or the supercritical fluid with a density greater than 0.35 g/cm3, and not exceeding the ranges of 0 to 1000 °C and 1 to 5000 bar (see Shock et al., 1992 for details). There are two ways in which these limits are enforced in CHNOSZ:
- The default source of water properties (H2O92D Fortran subroutine modified from SUPCRT92) yields NA values beyonds these T and P ranges. Note that the Deep Earth Water (DEW) model is available to extend the applicable range to pressures of up to 60 kbar (6 GPa) (Sverjensky et al., 2014).
subcrt()generates NA values beyond the density limit; theexceed.rhominargument can be used to enable calculations at lower density in the supercritical region.
The upper temperature limit for validity of the CGL heat-capacity
equation is a species-dependent parameter. This value is stored as a
negative value in the T column of OBIGT
and is used by subcrt() to issue a
warning at temperatures beyond this limit.
Required and optional data
REQUIRED:
- All species need a name and a
state.
- The state can be one of
aq,gas, orcr. - For minerals with higher-temperature polymorphs, they are named
cr2,cr3, etc. - NOTE:
crstands for “crystalline”; this naming convention (which was inherited from SUPCRT92 data files) refers to any solid phases including amorphous SiO2 and other minerals.
- The state can be one of
- A chemical formula is required to do
almost anything useful in CHNOSZ (e.g. check reaction balancing with
subcrt()and add species withbasis()orspecies()).- NOTE: The name of inorganic aqueous species and CH4 in OBIGT is the same as the chemical formula.
- Most minerals, gases, liquids, and organic aqueous species have a
name that is a common name. This permits a
shortcut to identify commonly used species in
subcrt(). - For example,
info("O2")refers to dissolved oxygen, whileinfo("oxygen")orinfo("O2", "gas")refers to the gas.
E_unitsneeds to be defined to perform any calculations of thermodynamic properties.- The value can be
Jfor Joules orcalfor calories.
- The value can be
OPTIONAL: Everything else. Really, it depends on what you need. For
instance, if you just want to use
subcrt() to calculate log K of
a reaction from ΔG° of species at 25 °C, then
G is the only parameter that is
needed.
OPTIONAL but useful:
- abbrv may be an abbreviation
(e.g. Qtz for quartz). It is used by
info()(together with name and formula) to look up species in the database. - date is a timestamp for the data entry (YYYY-MM-DD format in the default OBIGT database).
- ref1 and
ref2 are bibliographic reference keys.
They have matching entries in
extdata/OBIGT/refs.csv, which is used bythermo.refs()to display references, and invignettes/OBIGT.bib, which is used in the OBIGT thermodynamic database vignette to produce a reference list.
NOTE: Other functions in CHNOSZ do not depend on date, ref1, and ref2, so you can put anything there that is convenient for you.
NA or 0?
If a character value (in Columns 1–9) or thermodynamic parameter (in
Columns 10–14) is unknown, use NA. Note that a missing
(blank) value in the file is treated as NA.
- Unknown values for character values (usually abbrv, date, ref1, or ref2) should be NA.
- If you have only two of G,
H, and S,
then the missing one should be NA.
- Do NOT set a missing value of G, H, or S to 0. Zero is a numeric value that is incorrect except for very special cases.
- NOTE:
info()– and, by extension,subcrt()– “know” about the equation ΔG°f = ΔH°f - TΔS°f and the entropies of the elements needed to calculate ΔS°f from values of S in OBIGT. This equation is used to compute a missing value of G, H, or S from the other two, or to cross-check the values if all three are present for any species.
- If you don’t have Cp or
V, then set it to NA.
- If HKF or CGL parameters are present, they will be used to calculate Cp, so thermodynamic properties can be calculated at T > 25 °C.
- If HKF or CGL parameters aren’t present, thermodynamic properties can’t be calculated at T > 25 °C (NAs will propagate to higher T).
If an “equation-of-state” parameter or heat capacity coefficient (Columns 15-21) is unknown, use 0.
- Furthermore, if you would like to assume that
Cp or V is
0, then set it to 0.
- Then, thermodynamic properties will be extrapolated to T > 25 °C and P > 1 bar assuming that Cp and V are 0.
More detail on the inner working of the functions: For both HKF and CGL, if at least one parameter for a species is provided, any NA values of the other parameters are taken to be zero. If all EOS parameters are NA, but values of Cp and/or V are present, they are assumed to be constants for extrapolating thermodynamic properties (e.g. ΔG°) as a function of temperature and pressure.
OOM scaling and info()
HKF parameters in the the CSV files and OBIGT data frame are scaled
by order-of-magnitude (OOM) factors. For these parameters, OOM scaling
is nearly always used in published data tables. See
thermo() for details of the OOM
scaling.
info() provides a simple user
interface to the OBIGT database and is called by other functions in
CHNOSZ to retrieve unscaled values from the database. This is a summary
of its main features:
- Remove OOM scaling. This is used primarily by other functions in CHNOSZ to get a set of unscaled model parameters for calculating thermodynamic properties as a function of T and P.
- Extract the HKF or CGL parts of column names (only if all matching species have the same model).
- Calculate a missing one of G, H, or S if two of them are present.
- Cross-check G,
H, and S if
all of them are present, and print a message if the difference is above
a threshold (see
check.GHS()). - Calculate a missing Cp or V from the model parameters, if possible.
- Cross-check Cp or
V (if present) against the model parameters, if possible, and print a
message if the difference is above a threshold (see
check.EOS()).
NOTE:
info() does NOT
change the units of energy; the values it displays (including possibly
calculated ones) correspond to the E_units
for that species in OBIGT. On the other hand,
subcrt() outputs values in the units
previously selected with the function
E.units().
Case study: NA and 0 in the default database
Use the info() function to look at
the database and subcrt() to
calculate thermodynamic properties.
Let’s look at some minerals first. First use
info() to get the species indices
(i.e. rownumbers) in OBIGT, then pull out the “raw” data (including any
NA values).
## name abbrv formula state ref1 ref2 date model
## 2237 orpiment,amorphous <NA> As2S3 cr NA03 <NA> 2017-10-16 CGL
## 2230 arsenic,alpha <NA> As cr NA03 ZZL+16.1 2017-10-16 CGL
## 2179 tin Sn Sn cr JH85 <NA> 1985-08-00 CGL
## E_units G H S Cp V a1.a a2.b a3.c a4.d c1.e c2.f
## 2237 J -76800 -66900 200.00 NA NA NA NA NA NA NA NA
## 2230 J 0 0 35.63 24.43 12.960 NA NA NA NA NA NA
## 2179 cal 0 0 12.24 NA 16.289 4.42 0.0063 0 0 0 0
## omega.lambda z.T
## 2237 NA NA
## 2230 NA NA
## 2179 0 505.06Based on the values in the Cp column, would you predict that thermodynamic properties at T > 25 °C could be calculated for all of these minerals? Let’s see …
For conciseness we’ll consider a relatively small temperature range
and display only the out part of the
subcrt() output.
## subcrt: 1 species at 3 values of T (ºC) and P (bar) [energy units: J]## T P logK G H S V Cp
## 1 25 1 13.4548 -76800 -66900 200 NA NA
## 2 50 1 NA NA NA NA NA NA
## 3 75 1 NA NA NA NA NA NAThat makes sense; integrating NA Cp to calculate Gibbs energy and other thermodynamic properties would propagate NA, and that is what appears in the output. Now let’s run the calculation for the alpha phase of arsenic.
## subcrt: 1 species at 3 values of T (ºC) and P (bar) [energy units: J]## T P logK G H S V Cp
## 1 25 1 0.000000 0.000 0.00 35.6300 12.96 24.43
## 2 50 1 0.148008 -915.669 610.75 37.5971 12.96 24.43
## 3 75 1 0.281855 -1878.634 1221.50 39.4175 12.96 24.43What happened here? Even though there are no heat capacity coefficients (see above), there is a non-NA value of Cp, and that value is used together with the entropy for calculating Gibbs energy at T > 25 °C. Note that zero for the 25 °C values of G and H in this case is not a placeholder for unknown values (as noted above, unknown values should be represented by NA). Instaed, this is the reference state for the element, for which G and H are by convention equal to zero.
Let’s look at another element in its reference state, tin:
## info.character: found tin(cr) with 1 polymorphic transition## subcrt: 1 species at 3 values of T (ºC) and P (bar) [energy units: J]## subcrt: 2 polymorphs for tin ... polymorph 1 is stable## T P logK G H S V Cp polymorph
## 1 25 1 0.000000 0.00 0.000 51.2122 16.289 26.3523 1
## 2 50 1 0.211327 -1307.40 667.044 53.3602 16.289 27.0113 1
## 3 75 1 0.400148 -2667.09 1350.563 55.3973 16.289 27.6702 1Are you surprised? You might be if you only noticed the NA value for
Cp in OBIGT. However, there are non-NA
values for the heat capacity coefficients, which are used to calculate
Cp° as a function of temperature. When supplied with a
numeric argument (a species index),
info() actually does this to fill in
missing 25 °C values of Cp,
V, and G,
H, or S if
possible, in addition to simplifying column names:
## info.character: found tin(cr) with 1 polymorphic transition## info.numeric: Cp° of tin(cr) is NA; set by EOS parameters to 6.3 cal K-1 mol-1## name abbrv formula state ref1 ref2 date model E_units G H S
## 2179 tin Sn Sn cr JH85 <NA> 1985-08-00 CGL cal 0 0 12.24
## Cp V a b c d e f lambda T
## 2179 6.29834 16.289 4.42 0.0063 0 0 0 0 0 505.06Examples of adding data from a file
Using add.OBIGT() to add data from
optional data files for OBIGT or CSV files you make yourself.
add.OBIGT() with optional data
files
The default database has parameters for many minerals from Berman (1988); a notable exception is sulfide minerals, which are from Helgeson et al. (1978). Besides the different literature sources in ref1, the model column indicates that a different model is used for these minerals (Berman equations or CGL).
## info.numeric: Cp° of pyrite(cr) is NA; set by EOS parameters to 14.84 cal K-1 mol-1## name abbrv formula state ref1 ref2 date model E_units G
## 2073 quartz Qz SiO2 cr Ber88 <NA> 2017-10-01 Berman J -856288
## 2160 pyrite Py FeS2 cr HDNB78 RH95.7 1978-05-05 CGL cal -38293
## H S Cp V a b c d e f lambda T
## 2073 -910700 41.46 44.7423 22.69 NA NA NA NA NA NA NA NA
## 2160 -41000 12.65 14.8425 23.94 17.88 0.00132 -305000 0 0 0 0 1015Sometimes it is useful to load mineral data from the SUPCRT92
database, corresponding largely to the compilation by Helgeson et al. (1978). This can be done with
add.OBIGT(). In this example we load
SUPCRT92 data for just one mineral, quartz.
## add.OBIGT: read 1 rows; made 1 replacements, 0 additions [energy units: cal]## info.numeric: Cp° of quartz(cr) is NA; set by EOS parameters to 10.63 cal K-1 mol-1## name abbrv formula state ref1 ref2 date model E_units G
## 2073 quartz Qtz SiO2 cr HDNB78 <NA> 1978-05-05 CGL cal -204646
## H S Cp V a b c d e f lambda T
## 2073 -217650 9.88 10.6275 22.688 11.22 0.0082 -270000 0 0 0 0 848Here we load all minerals available in the optional SUPCRT92 data
file and then list the names. NOTE:
suppressMessages() is used to suppress messages from
info() about missing parameters, and
unique() is used to list each mineral only once (because
each polymorph has a separate entry).
## add.OBIGT: read 177 rows; made 65 replacements, 112 additions [energy units: cal]## [1] "akermanite" "albite" "albite,high"
## [4] "albite,low" "almandine" "andalusite"
## [7] "andradite" "annite" "anorthite"
## [10] "anthophyllite" "antigorite" "aragonite"
## [13] "boehmite" "brucite" "Ca-Al-pyroxene"
## [16] "calcite" "chrysotile" "clinozoisite"
## [19] "coesite" "corundum" "cristobalite,alpha"
## [22] "cristobalite,beta" "diaspore" "diopside"
## [25] "dolomite" "enstatite" "epidote"
## [28] "fayalite" "ferrosilite" "fluorphlogopite"
## [31] "fluortremolite" "forsterite" "gehlenite"
## [34] "gibbsite" "glaucophane" "grossular"
## [37] "grunerite" "hedenbergite" "hematite"
## [40] "jadeite" "K-feldspar" "kaolinite"
## [43] "kyanite" "lawsonite" "lime"
## [46] "magnesite" "magnetite" "margarite"
## [49] "merwinite" "monticellite" "muscovite"
## [52] "paragonite" "periclase" "phlogopite"
## [55] "prehnite" "pyrope" "pyrophyllite"
## [58] "quartz" "sillimanite" "spinel"
## [61] "talc" "tremolite" "wollastonite"
## [64] "zoisite" "rutile" "aegerine"
## [67] "amesite,14A" "amesite,7A" "amorphous silica"
## [70] "analcime" "analcime,dehydrated" "Ca-phillipsite"
## [73] "celadonite" "chabazite" "chalcedony"
## [76] "chloritoid" "clinochlore,14A" "clinochlore,7A"
## [79] "cordierite,dry" "cordierite,hydrous" "cristobalite"
## [82] "cronstedtite,7A" "cummingtonite" "daphnite,14A"
## [85] "daphnite,7A" "dickite" "dolomite,disordered"
## [88] "dolomite,ordered" "edenite" "epidote,ordered"
## [91] "epistilbite" "ferroedenite" "ferrogedrite"
## [94] "ferropargasite" "ferrotremolite" "fluoredenite"
## [97] "fluorite" "greenalite" "halloysite"
## [100] "hastingsite" "heulandite" "K-phillipsite"
## [103] "kalsilite" "larnite" "laumontite"
## [106] "leonhardite" "magnesiohastingsite" "magnesioriebeckite"
## [109] "microcline,maximum" "minnesotaite" "Na-phillipsite"
## [112] "natrolite" "nepheline" "pargasite"
## [115] "PD-oxyannite" "richterite" "riebeckite"
## [118] "sanidine,high" "sepiolite" "spessartine"
## [121] "staurolite" "stilbite" "wairakite"
## [124] "titanite"add.OBIGT() with other CSV
files
add.OBIGT() can also be used to add
data from a user-specified file to the OBIGT database. The file must be
a CSV (comma separated value) file with column headers that match those
in the default database (i.e., thermo()$OBIGT). As an
example, here are the contents of BZA10.csv, which has
parameters taken from Bazarkina et al. (2010). Missing values are indicated by
NA:
## name abbrv formula state ref1 ref2 date model E_units G H
## 1 CdCl+ NA CdCl+ aq BZA10 NA 2010-07-03 HKF cal -52629 NA
## 2 CdCl2 NA CdCl2 aq BZA10 NA 2010-07-03 HKF cal -84883 NA
## 3 CdCl3- NA CdCl3- aq BZA10 NA 2010-07-03 HKF cal -115399 NA
## 4 CdCl4-2 NA CdCl4-2 aq BZA10 NA 2010-07-03 HKF cal -145583 NA
## S Cp V a1.a a2.b a3.c a4.d c1.e c2.f
## 1 7.06 11.12 2.20 2.2303 -2.3357 6.6681 -2.6824 16.6723 -0.7693
## 2 25.72 116.01 42.21 7.5221 10.5852 1.5895 -3.2166 73.7023 20.5956
## 3 45.15 97.78 63.47 10.8045 18.5994 -1.5605 -3.5479 72.0244 16.8832
## 4 50.61 42.52 81.35 13.8329 25.9938 -4.4669 -3.8536 53.6766 5.6267
## omega.lambda z.T
## 1 0.4372 1
## 2 -0.0495 0
## 3 0.9378 -1
## 4 2.4766 -2Loading the data with add.OBIGT()
produces a message that the new data replace existing species. We can
then use subcrt() to calculate the
equilibrium constant for a reaction involving the new species. Note the
decrease in the stepwise stability constant for the second cadmium
chloride complex with increasing pressure (Bazarkina et al., 2010, Fig.
4).
## add.OBIGT: read 4 rows; made 4 replacements, 0 additions [energy units: cal]## subcrt: 3 species at 2 values of T (ºC) and P (bar) (wet) [energy units: J]## $reaction
## coeff name formula state ispecies model
## 377 -1 CdCl+ CdCl+ aq 377 HKF
## 29 -1 Cl- Cl- aq 29 HKF
## 378 1 CdCl2 CdCl2 aq 378 HKF
##
## $out
## T P rho logK G H S V Cp
## 1 25 1 0.997061 0.641379 -3661.000 2669.61 21.3384 22.6816 561.306
## 2 25 2000 1.071783 0.061853 -353.058 2732.39 10.4541 12.2632 539.264After running reset() we can look up
the source of data in the default OBIGT database (Sverjensky et al.,
1997). Running the reaction with thermodynamic parameters
from the default database, we now see that the equilibrium constant is
not as sensitive to pressure:
## reset: resetting "thermo" object## OBIGT: loading default database with 2018 aqueous, 3570 total species## key author year
## 87 SSH97 D. A. Sverjensky, E. L. Shock and H. C. Helgeson 1997## subcrt: 3 species at 2 values of T (ºC) and P (bar) (wet) [energy units: J]## $reaction
## coeff name formula state ispecies model
## 377 -1 CdCl+ CdCl+ aq 377 HKF
## 29 -1 Cl- Cl- aq 29 HKF
## 378 1 CdCl2 CdCl2 aq 378 HKF
##
## $out
## T P rho logK G H S V Cp
## 1 25 1 0.997061 0.619389 -3535.48 -8522.81 -16.7360 10.04043 214.488
## 2 25 2000 1.071783 0.425896 -2431.02 -9603.90 -24.0664 2.34196 191.416Examples of adding and modifying data with a function
Use mod.OBIGT() to add or modify the
database in the current session. The function requires the name of a
species and one or more properties to change.
mod.OBIGT() for aqueous
species
Let’s add data for CoCl4-2 from Liu et al. (2011). The
values are taken from Table 5 of that paper; note that they are reported
in caloric units, which is rather common for the HKF model. The entry
includes the date in ISO 8601 extended format (e.g. 2020-08-16);
Sys.Date() is used in this example to get the current
date.
mod.OBIGT("CoCl4-2", formula = "CoCl4-2", state = "aq", ref1 = "LBT+11", E_units = "cal",
date = as.character(Sys.Date()), G = -134150, H = -171558, S = 19.55, Cp = 72.09, V = 27.74)## mod.OBIGT: updated CoCl4-2(aq)## [1] 867The function prints a message saying that the species was added and
returns the species index of the new species. Now let’s modify the new
species by adding the HKF coefficients including the OOM multipliers, as
they are usually given in publications. The z at the end
refers to the charge of the species, and is used only for calculating
the “g function” in the revised HKF model, not for balancing
reactions.
mod.OBIGT("CoCl4-2", a1 = 6.5467, a2 = 8.2069, a3 = 2.0130, a4 = -3.1183,
c1 = 76.3357, c2 = 11.6389, omega = 2.9159, z = -2)## mod.OBIGT: no change for CoCl4-2(aq)## [1] 867Let us now calculate the equilibrium constant for the formation of CoCl4-2 from Co+2 and Cl-.
T <- c(25, seq(50, 350, 50))
sres <- subcrt(c("Co+2", "Cl-", "CoCl4-2"), c(-1, -4, 1), T = T)
round(sres$out$logK, 2)## [1] -3.20 -2.96 -2.02 -0.74 0.77 2.50 4.57 7.29The calculated values of logK are identical to those in
Table 9 of Liu et al. (2011), which provides a good indication
that the thermodynamic parameters were entered correctly. Nevertheless,
this isn’t a guarantee that the thermodynamic parameters are consistent
with the provided values of CP° and
V°. We can see this by running
info() to cross-check the parameters
for the new CoCl4-2 species:
## check.EOS: calculated Cp° of CoCl4-2(aq) differs by 1.33 cal K-1 mol-1 from database value## check.EOS: calculated V° of CoCl4-2(aq) differs by -1.19 cm3 mol-1 from database valueThe messages indicate that the given values of CP° and V° differ slightly from those calculated using the HKF parameters.
mod.OBIGT() for minerals
Let’s add data for magnesiochromite from Klemme et al. (2000).
The parameters in this paper are reported in Joules, so we set the
E.units() to J. The value for volume,
in cm3 mol-1, is from Robie
and Hemingway (1995).
H <- -1762000
S <- 119.6
V <- 43.56
mod.OBIGT("magnesiochromite", formula = "MgCr2O4", state = "cr", ref1 = "KOSG00",
date = as.character(Sys.Date()), E_units = "J", H = H, S = S, V = V)## mod.OBIGT: added magnesiochromite(cr) with CGL model and energy units of J## [1] 3571Here are the heat capacity parameters for the Haas-Fisher polynomial equation (Cp = a + bT + cT−2 + dT−0.5 + eT2). As of CHNOSZ 2.0.0, OOM multipliers are not used for these coefficients. 1500 K is a generic value for the high-temperature limit; experimental heat capacities were only reported up to 340 K (Klemme et al., 2000).
a <- 221.4
b <- -0.00102030
c <- -1757210
d <- -1247.9
mod.OBIGT("magnesiochromite", E_units = "J", a = a, b = b, c = c, d = d,
e = 0, f = 0, lambda = 0, T = 1500)## mod.OBIGT: updated magnesiochromite(cr)## [1] 3571NOTE: An additional f term is available, which can have any exponent given in lambda. This offers some flexibility for using heat capacity equations that are different from the Haas-Fisher polynomial.
Now we can use subcrt() to
calculate the heat capacity of magnesiochromite. For this calculation,
we set the temperature units to Kelvin. We also specify a pressure of 1
bar because the default setting of Psat
(liquid-vapor saturation) causes an error below the freezing temperature
of water.
## changed temperature units to K## subcrt: 1 species at 3 values of T (K) and P (bar) [energy units: J]## $species
## name formula state ispecies model
## 3571 magnesiochromite MgCr2O4 cr 3571 CGL
##
## $out
## $out$magnesiochromite
## T P Cp
## 1 250 1 114.105
## 2 300 1 129.522
## 3 340 1 138.175Next we check that the calculated values are within 0.3 J K-1 mol-1 of reference values taken from Fig. 1 of Klemme et al. (2000).
Finally, let’s restore the units setting for later calculations with
subcrt(). (Another way would be to
run reset(), which also resets the
OBIGT database.)
## changed temperature units to CCase study: Formation constants for aqueous tungsten species
Here we use logB.to.OBIGT() to fit
to thermodynamic parameters to experimental formation constants. Some
additional steps are shown to refine a thermodynamic model to generate a
speciation diagram as a function of pH.
Fitting formation constants
logB.to.OBIGT() requires three
things:
- Experimental decimal logarithms of formation constants (log β) as a function of temperature;
- The stoichiometry of the formation reaction in terms of known species (the new species must be last);
- The experimental temperature and pressure.
logB.to.OBIGT() does three
things:
- Combines the formation constants with standard Gibbs energies (ΔG°) of the known species to calculate ΔG° of the new species;
- Fits ΔG° of the new species using 25 °C thermodynamic properties and selected HKF model parameters (i.e., G, S, c1, c2, and omega parameters in OBIGT);
- Adds the parameters to OBIGT for use by other functions in CHNOSZ.
First we set the pressure for all log β data.
Add first species: HWO4- (Wang et al., 2019).
T <- c(250, 300, 350)
logB <- c(5.58, 6.51, 7.99)
species <- c("WO4-2", "H+", "HWO4-")
coeff <- c(-1, -1, 1)
logB.to.OBIGT(logB, species, coeff, T, P)## mod.OBIGT: updated HWO4-(aq)## logB.to.OBIGT: mean difference between logB (experimental) and logK (calculated) is 0## [1] 427Add second species: H3WO4F2- (Wang et al., 2021).
T <- seq(100, 250, 25)
logB <- c(17.00, 17.11, 17.46, 17.75, 18.17, 18.71, 19.23)
# Species and coefficients in the formation reaction
species <- c("H+", "WO4-2", "F-", "H3WO4F2-")
coeff <- c(-3, -1, -2, 1)
logB.to.OBIGT(logB, species, coeff, T, P)## mod.OBIGT: updated H3WO4F2-(aq)## logB.to.OBIGT: mean difference between logB (experimental) and logK (calculated) is 0.0307## [1] 3572Add third species: H2WO4 (Wang et al., 2021). Here we increase the tolerance because there is considerable scatter in the experimental values.
logB <- c(7.12, 7.82, 7.07, 7.76, 7.59, 7.98, 8.28)
species <- c("H+", "WO4-2", "H2WO4")
coeff <- c(-2, -1, 1)
logB.to.OBIGT(logB, species, coeff, T, P, tolerance = 0.3)## mod.OBIGT: updated H2WO4(aq)## logB.to.OBIGT: mean difference between logB (experimental) and logK (calculated) is 0.2035## [1] 906After running, logB.to.OBIGT()
returns the species indices; the low values for
HWO4- (427) and H2WO4 (906)
indicate that the function replaced parameters for these species that
were already present in OBIGT.
Diagram 1: Constant molality of F-
Now we’re ready to make a speciation diagram. Our aim is to reproduce
Fig. 7b of Wang et al. (2021), which is made for 300 °C. A
constant molality of F- is based on the assumption of
complete dissociation of 0.1 m NaF (we’ll change this later). An ionic
strength of 0.9 mol/kg is estimated for a solution with 1.8 m NaCl (use
NaCl(1.8, T = 300)). NOTE: because the ionic
strength is non-zero, the calculations here refer to molality instead of
activity of species (see An Introduction to
CHNOSZ).
basis(c("H+", "WO4-2", "F-", "H2O", "O2"))
basis("F-", log10(0.1))
iaq <- retrieve("W", c("O", "H", "F"), "aq")
species(iaq)
a <- affinity(pH = c(2, 7), T = 300, IS = 0.9)
e <- equilibrate(a)
col <- c(1, 4, 5, 2)
diagram(e, alpha = TRUE, col = col, lty = 1, lwd = 2, ylab = "Fraction total W")
This isn’t quite the diagram we were looking for. The published diagram shows a broad region of coexistence of H3WO4F2- and HWO4- at pH < 5 and increasing abundance of H2WO4 at lower pH.
Diagram 2: Variable molality of F-
In reality, the molality of F- depends strongly on pH
according to the reaction H+ + F- = HF. With a
little algebra, we can calculate the molality of F-
(a_F in the code below) from the equilbrium constant of
this reaction for a given total F concentration (F_tot).
NOTE: It is important to
call subcrt() with a non-zero
IS so that it returns effective equilibrium constants
corrected for ionic strength (try setting IS = 0 yourself
and look at what happens to the diagram).
T <- 300
pH <- seq(2, 7, 0.1)
logK_HF <- subcrt(c("H+", "F-", "HF"), c(-1, -1, 1), T = T, IS = 0.9)$out$logK## info.character: found HF(aq); also available in gas## subcrt: 3 species at 300 ºC and 85.84 bar (wet) [energy units: J]## nonideal: calculations for F- (Bdot equation)## nonideal: calculations for HF (Setchenow equation)Now that we have the molality of F- as a function of pH,
we can provide it in the call to
affinity().
basis(c("H+", "WO4-2", "F-", "H2O", "O2"))
iaq <- retrieve("W", c("O", "H", "F"), "aq")
species(iaq)
a <- affinity(pH = pH, "F-" = log10(a_F), T = T, IS = 0.9)
e <- equilibrate(a)
diagram(e, alpha = TRUE, col = col, lty = 1, lwd = 2, ylab = "Fraction total W")
That’s more like it. We have captured the basic geometry of Fig. 7b in Wang et al. (2021). For instance, in accord with the published diagram, HWO4- plateaus at around 40% of total W, and H2WO4 and H3WO4F2- are nearly equally abundant at pH = 2.
The highest experimental temperature for the formation constants of H2WO4 and H3WO4F2- is 250 °C, but this diagram is drawn for 300 °C. Wang et al. (2021) used the modified Ryzhenko-Bryzgalin (MRB) model to extrapolate to 300 °C. In contrast, we used a different model but obtained quite similar results.
NOTE: The coefficients
in the model used by logB.to.OBIGT()
include 25 °C values of G and
S. These should be conservatively treated
only as fitting parameters and should not be used to compute
thermodynamic properties close to 25 °C unless they were fit to
experimental data in that temperature range.