logKcalc: Calculate Equilibrium Constants for a GWB Data File
In jedick/logKcalc: Equilibrium Constants for GWB Thermodynamic Data Files
logKcalc R Documentation
Calculate Equilibrium Constants for a GWB Data File
Description
Read dissociation reactions from a thermodynamic data file for The Geochemist's Workbench®, calculate equilibrium constants at specified temperature and pressure using the OBIGT database in CHNOSZ, and write a modified GWB file.
Usage
logKcalc(infile = "thermo.tdat", outfile = "thermo_OBIGT.tdat",
T = NULL, P = "Psat", ispecies = NULL, a0_ion = NULL, a0_neutral = 0,
update.formulas = TRUE, DH.method = "bdot", maxprint = Inf)
Arguments
infile
Input GWB data file
outfile
Output GWB data file
T
Temperature (\degC)
P
Pressure (bar)
ispecies
Species index of specie(s) to add to the data file
a0_ion
Ion size parameter for added aqueous ions
a0_neutral
Ion size parameter for added aqueous neutral species
update.formulas
Should formulas for aqueous species and minerals be updated?
DH.method
Method for generating Debye-Hückel parameters, ‘bdot’ or ‘bgamma’
maxprint
Maximum number of missing species to show in print statements
Details
This function reads the input GWB file and assembles lists of species in the file.
Species names are then converted to those used in OBIGT, by automatically matching the species name or by looking in a mapping file (see Names Mapping).
Species that are not available in OBIGT are omitted from the output file.
Species are also omitted if they available in OBIGT but their dissociation reactions in the GWB file are written in terms of basis species or redox couples that are not available.
Two examples are As(OH)\s4\S- and Sn\S++++, which are basis species in the ‘thermo.tdat’ file distributed with GWB but are not present in the default OBIGT database.
addOBIGT can be used to add these species to OBIGT.
Normally, the function obtains values of \T from the input GWB file and calculates values of \P corresponding to \Psat (i.e. pressures corresponding to liquid-vapor saturation).
This is done because some values of \P in existing files are slightly below \Psat, preventing calculation at those conditions.
To use values of \P taken from the input file, set P to NULL.
To use any other values, provide 8 numeric values in T and/or P.
Logarithms (base 10) of equilibrium constants are calculated for each of the following species types: redox, aqueous, electron, mineral, gas.
The corresponding values for each species are inserted into the output GWB file.
Any NA values (e.g. for temperatures above the melting point of a mineral) are represented by ‘500’ in the output file to indicate unavailable data (Bethke and Farrell, 2020).
One or two citations are added for each species based on the references in OBIGT.
If the CHNOSZ version is > 1.3.6, a reference list is added at the end of the output file.
Chemical formulas and reactions
By default, the stoichiometric formulas in OBIGT are added to the ‘formula’ field in the output for aqueous species and minerals.
These formulas may be formatted differently from formulas in existing GWB data files, especially for charged species.
Set update.formulas to FALSE to not add or update formulas, but retain any formulas that may be present in the input file.
This might be required for dataset formats ‘oct13’ or earlier, which do not have formulas for the aqueous species (Bethke and Farrell, 2020).
Note that the names of aqueous species, which may also be formulas, are never changed from the input file.
The mass balance of dissociation reactions is checked (in CHNOSZ::subcrt) using the chemical formulas of species in OBIGT.
To allow for rounding errors in fractional coefficients, reactions are considered unbalanced if the difference for any element exceeds 0.001.
Messages are shown for unbalanced reactions, and these species are omitted from the output file.
This includes the minerals antigorite and pyrrhotite in the GWB ‘thermo.tdat’ file, for which the chemical formulas (and therefore reactions) differ from those in OBIGT.
These minerals can be included by using modOBIGT and the ispecies argument, as shown in the example and thermo.Rmd vignette.
Names Mapping
The names mapping step uses some heuristics to find species in OBIGT.
All names from the GWB file are converted to lower-case, suffixes “(aq)” and “,aq” are removed, and underscores are converted to spaces.
Repeated ‘-’ or ‘+’ signs at the ends of names are changed to numbers (e.g. ‘-2’ or ‘+3’).
Leading ‘1-’ are deleted (e.g. 1-pentanol becomes pentanol).
The names of species in the GWB file, which correspond to chemical formulas for many aqueous species, are used for mapping.
However, any chemical formulas that may be listed in addition to names of minerals or other species are not used for mapping.
Both names and formulas of species in OBIGT are searched, but formulas are searched as-is, so formulas that differ from those used in CHNOSZ (e.g. CH\s3COOH instead of C\s2H\s4O\s2) aren't matched unless they are in the names mapping file.
Species names that are not mapped after these steps are then searched in the names mapping file.
This file was constructed by identifying unmapped species in ‘thermo.tdat’ and GWB files produced by K2GWB and DBCreate.
Some names (for feldspars in particular) have multiple entries in this file, to allow mapping to names used in the Berman or SUPCRT92 datasets, which can be selected using modOBIGT.
Adding Species
By default, only species in the input file that can be mapped to OBIGT are included in the output file.
Other species can be added to the output using the ispecies argument.
The values in this argument are the species index for one or more species, as returned by CHNOSZ::info.
Species can only be added if all of their elements are in the available basis species in the GWB file.
Species with the states ‘aq’, ‘cr’, and ‘gas’ in OBIGT are added to the aqueous, mineral, and gas sections of the output file, respectively; any other states produce an error.
The molecular weight of species is calculated using CHNOSZ::mass, which likely does not correspond precisely to the element masses in the input GWB file.
The default ion size parameter for added aqueous ions (a0_ion) is 4.5 \AA; this setting follows the behavior in UNITHERM for ions without a specified value (Shvarov and Bastrakov, 1999).
The ion size parameter for neutral aqueous species (a0_neutral) has a special meaning in the GWB thermodynamic data file (Bethke and Farrell, 2020).
The default value of 0 corresponds to unit activity coefficient (the default in the UT2K converter used by K2GWB; Cleverley and Bastrakov, 2005).
With a value of -0.5, the activity coefficient is equal to that of CO\s2 in an NaCl solution of the same ionic strength (this is used for B(OH)\s3, O\s2, SiO\s2, CH\s4, and H\s2 in ‘thermo.tdat’).
For values \le -1, the activity coefficient is calculated as the product of the ‘bdot’ parameter and the true ionic strength.
The values of a0_ion and a0_neutral are recycled to the total number of species, including all aqueous ions and neutral species as well as non-aqueous species.
Therefore, to specify particular values of these parameters for different ions or neutral aqueous species, values must be specified for all added species, but will only be used for the respective aqueous species.
Activity and Osmotic Coefficients
Coefficients in the extended Debye-Hückel equation (‘adh’ and ‘bdh’ in the GWB file) are calculated as a function of \T and \P using CHNOSZ::water.
The extended term (‘bdot’ parameter) by default is calculated using a spline fit to values at 25-300 \degC from Helgeson (1969); values at \T > 300 \degC are set to zero.
The DH.method argument can be changed to ‘bgamma’ to use the b-gamma approximation, in which the extended term parameter is a function of temperature and pressure that has been extrapolated to 1000 \degC and 30 kbar (Manning et al., 2013) and provisionally as high as 60 kbar (see Figure S5 of Dick, 2019).
The high-pressure extrapolation is intended to be used in conjunction with the Deep Earth Water (DEW) model for water.
The default for ‘bgamma’ is to set the ion size parameter for all ions in the output file to 3.72 \AA, which is a characteristic value for an NaCl-dominated solution (Helgeson et al., 1981).
This differs from HCh, where ion size parameters for individual ions are calculated from their effective elecrostatic radii (Shvarov and Bastrakov, 1999).
Calling the function with the explicit argument a0_ion = NULL overrides this setting, so that no changes are made to ion size parameters from the input file.
Parameters for calculating the activity coefficient of CO\s2 in a pure NaCl solution and the osmotic coefficient of \H2O (i.e., ‘co2’ and ‘h2o’ blocks of the GWB file) are calculated as described for K2GWB (Cleverley and Bastrakov, 2005).
Because of a bug with calculating exponents in some versions of K2GWB, the CO\s2 parameters (especially in the first block) may not match those in GWB files produced by K2GWB for the same values of \T.
Dataset Formats
The conversion is compatible with all existing dataset formats of GWB files.
For dataset formats ‘feb94’ and ‘oct94’, logKcalc updates the \logK blocks for the Eh reaction and O\s2 gas solubility.
(Blocks for H\s2 and N\s2 gas solubility, which are not used by GWB, are not updated.)
These blocks are not present in dataset formats ‘oct13’ and later.
Instead, the free electron is present as a species; logKcalc updates its \logK values.
Dataset formats ‘jul17’ and ‘jan19’ specify a fugacity model (tsonopoulos, peng-robinson, or spycher-reed).
The parameters are located in the data blocks for individual gas species and are not modified by logKcalc.
References
Bethke, C. M. and Farrell, B. (2020) The Geochemist's Workbench® Release 14 GWB Reference Manual, Aqeuous Solutions, LLC, Champaign, Illinois.
Cleverley, J. S. and Bastrakov, E. N. (2005) K2GWB: Utility for generating thermodynamic data files for The Geochemist's Workbench® at 0-1000 °C and 1-5000 bar from UT2K and the UNITHERM database. Computers & Geosciences 31, 756–767. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.cageo.2005.01.007")}
Dick, J. M. (2019) CHNOSZ: Thermodynamic calculations and diagrams for geochemistry. Front. Earth Sci. 7:180. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.3389/feart.2019.00180")}
Helgeson, H. C., Kirkham, D. H. and Flowers, G. C. (1981) Theoretical prediction of the thermodynamic behavior of aqueous electrolytes at high pressures and temperatures. IV. Calculation of activity coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600\degC and 5 Kb. Am. J. Sci. 281, 1249–1516. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.2475/ajs.281.10.1249")}
Manning, C. E., Shock, E. L. and Sverjensky, D. A. (2013) The chemistry of carbon in aqueous fluids at crustal and upper-mantle conditions: Experimental and theoretical constraints. Rev. Mineral. Geochem. 75, 109–148. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.2138/rmg.2013.75.5")}
Shvarov, Y. and Bastrakov, E. (1999) HCh: A software package for geochemical equilibrium modelling. User's Guide. Australian Geological Survey Organisation 1999/25. http://pid.geoscience.gov.au/dataset/ga/25473
See Also
modOBIGT for preparing the thermodynamic database for the \logK calculation.
Examples
# Show the first few entries in the names mapping file
mapfile <- system.file("extdata/map_names.csv", package = "logKcalc")
head(read.csv(mapfile))
# Process a GWB file
# -- define input and output files
infile <- system.file("extdata/thermo_12elements.tdat", package = "logKcalc")
outfile <- tempfile("thermo", fileext = ".tdat")
# -- make some modifications to OBIGT
# (add minerals from SUPCRT92 database and deprecated H2O(g))
modOBIGT(c("addSUPCRT", "steam"))
# -- identify species to add to the output
ispecies <- info("pyrrhotite")
# -- process the file
logKcalc(infile, outfile, ispecies = ispecies)
jedick/logKcalc documentation built on Jan. 19, 2025, 10:41 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
| logKcalc | R Documentation |
Calculate Equilibrium Constants for a GWB Data File
Description
Read dissociation reactions from a thermodynamic data file for The Geochemist's Workbench®, calculate equilibrium constants at specified temperature and pressure using the OBIGT database in CHNOSZ, and write a modified GWB file.
Usage
logKcalc(infile = "thermo.tdat", outfile = "thermo_OBIGT.tdat",
T = NULL, P = "Psat", ispecies = NULL, a0_ion = NULL, a0_neutral = 0,
update.formulas = TRUE, DH.method = "bdot", maxprint = Inf)
Arguments
infile |
Input GWB data file |
outfile |
Output GWB data file |
T |
Temperature (\degC) |
P |
Pressure (bar) |
ispecies |
Species index of specie(s) to add to the data file |
a0_ion |
Ion size parameter for added aqueous ions |
a0_neutral |
Ion size parameter for added aqueous neutral species |
update.formulas |
Should formulas for aqueous species and minerals be updated? |
DH.method |
Method for generating Debye-Hückel parameters, ‘bdot’ or ‘bgamma’ |
maxprint |
Maximum number of missing species to show in print statements |
Details
This function reads the input GWB file and assembles lists of species in the file. Species names are then converted to those used in OBIGT, by automatically matching the species name or by looking in a mapping file (see Names Mapping). Species that are not available in OBIGT are omitted from the output file.
Species are also omitted if they available in OBIGT but their dissociation reactions in the GWB file are written in terms of basis species or redox couples that are not available.
Two examples are As(OH)\s4\S- and Sn\S++++, which are basis species in the ‘thermo.tdat’ file distributed with GWB but are not present in the default OBIGT database.
addOBIGT can be used to add these species to OBIGT.
Normally, the function obtains values of \T from the input GWB file and calculates values of \P corresponding to \Psat (i.e. pressures corresponding to liquid-vapor saturation).
This is done because some values of \P in existing files are slightly below \Psat, preventing calculation at those conditions.
To use values of \P taken from the input file, set P to NULL.
To use any other values, provide 8 numeric values in T and/or P.
Logarithms (base 10) of equilibrium constants are calculated for each of the following species types: redox, aqueous, electron, mineral, gas. The corresponding values for each species are inserted into the output GWB file. Any NA values (e.g. for temperatures above the melting point of a mineral) are represented by ‘500’ in the output file to indicate unavailable data (Bethke and Farrell, 2020).
One or two citations are added for each species based on the references in OBIGT. If the CHNOSZ version is > 1.3.6, a reference list is added at the end of the output file.
Chemical formulas and reactions
By default, the stoichiometric formulas in OBIGT are added to the ‘formula’ field in the output for aqueous species and minerals.
These formulas may be formatted differently from formulas in existing GWB data files, especially for charged species.
Set update.formulas to FALSE to not add or update formulas, but retain any formulas that may be present in the input file.
This might be required for dataset formats ‘oct13’ or earlier, which do not have formulas for the aqueous species (Bethke and Farrell, 2020).
Note that the names of aqueous species, which may also be formulas, are never changed from the input file.
The mass balance of dissociation reactions is checked (in CHNOSZ::subcrt) using the chemical formulas of species in OBIGT.
To allow for rounding errors in fractional coefficients, reactions are considered unbalanced if the difference for any element exceeds 0.001.
Messages are shown for unbalanced reactions, and these species are omitted from the output file.
This includes the minerals antigorite and pyrrhotite in the GWB ‘thermo.tdat’ file, for which the chemical formulas (and therefore reactions) differ from those in OBIGT.
These minerals can be included by using modOBIGT and the ispecies argument, as shown in the example and thermo.Rmd vignette.
Names Mapping
The names mapping step uses some heuristics to find species in OBIGT. All names from the GWB file are converted to lower-case, suffixes “(aq)” and “,aq” are removed, and underscores are converted to spaces. Repeated ‘-’ or ‘+’ signs at the ends of names are changed to numbers (e.g. ‘-2’ or ‘+3’). Leading ‘1-’ are deleted (e.g. 1-pentanol becomes pentanol).
The names of species in the GWB file, which correspond to chemical formulas for many aqueous species, are used for mapping. However, any chemical formulas that may be listed in addition to names of minerals or other species are not used for mapping. Both names and formulas of species in OBIGT are searched, but formulas are searched as-is, so formulas that differ from those used in CHNOSZ (e.g. CH\s3COOH instead of C\s2H\s4O\s2) aren't matched unless they are in the names mapping file.
Species names that are not mapped after these steps are then searched in the names mapping file.
This file was constructed by identifying unmapped species in ‘thermo.tdat’ and GWB files produced by K2GWB and DBCreate.
Some names (for feldspars in particular) have multiple entries in this file, to allow mapping to names used in the Berman or SUPCRT92 datasets, which can be selected using modOBIGT.
Adding Species
By default, only species in the input file that can be mapped to OBIGT are included in the output file.
Other species can be added to the output using the ispecies argument.
The values in this argument are the species index for one or more species, as returned by CHNOSZ::info.
Species can only be added if all of their elements are in the available basis species in the GWB file.
Species with the states ‘aq’, ‘cr’, and ‘gas’ in OBIGT are added to the aqueous, mineral, and gas sections of the output file, respectively; any other states produce an error.
The molecular weight of species is calculated using CHNOSZ::mass, which likely does not correspond precisely to the element masses in the input GWB file.
The default ion size parameter for added aqueous ions (a0_ion) is 4.5 \AA; this setting follows the behavior in UNITHERM for ions without a specified value (Shvarov and Bastrakov, 1999).
The ion size parameter for neutral aqueous species (a0_neutral) has a special meaning in the GWB thermodynamic data file (Bethke and Farrell, 2020).
The default value of 0 corresponds to unit activity coefficient (the default in the UT2K converter used by K2GWB; Cleverley and Bastrakov, 2005).
With a value of -0.5, the activity coefficient is equal to that of CO\s2 in an NaCl solution of the same ionic strength (this is used for B(OH)\s3, O\s2, SiO\s2, CH\s4, and H\s2 in ‘thermo.tdat’).
For values \le -1, the activity coefficient is calculated as the product of the ‘bdot’ parameter and the true ionic strength.
The values of a0_ion and a0_neutral are recycled to the total number of species, including all aqueous ions and neutral species as well as non-aqueous species.
Therefore, to specify particular values of these parameters for different ions or neutral aqueous species, values must be specified for all added species, but will only be used for the respective aqueous species.
Activity and Osmotic Coefficients
Coefficients in the extended Debye-Hückel equation (‘adh’ and ‘bdh’ in the GWB file) are calculated as a function of \T and \P using CHNOSZ::water.
The extended term (‘bdot’ parameter) by default is calculated using a spline fit to values at 25-300 \degC from Helgeson (1969); values at \T > 300 \degC are set to zero.
The DH.method argument can be changed to ‘bgamma’ to use the b-gamma approximation, in which the extended term parameter is a function of temperature and pressure that has been extrapolated to 1000 \degC and 30 kbar (Manning et al., 2013) and provisionally as high as 60 kbar (see Figure S5 of Dick, 2019).
The high-pressure extrapolation is intended to be used in conjunction with the Deep Earth Water (DEW) model for water.
The default for ‘bgamma’ is to set the ion size parameter for all ions in the output file to 3.72 \AA, which is a characteristic value for an NaCl-dominated solution (Helgeson et al., 1981).
This differs from HCh, where ion size parameters for individual ions are calculated from their effective elecrostatic radii (Shvarov and Bastrakov, 1999).
Calling the function with the explicit argument a0_ion = NULL overrides this setting, so that no changes are made to ion size parameters from the input file.
Parameters for calculating the activity coefficient of CO\s2 in a pure NaCl solution and the osmotic coefficient of \H2O (i.e., ‘co2’ and ‘h2o’ blocks of the GWB file) are calculated as described for K2GWB (Cleverley and Bastrakov, 2005). Because of a bug with calculating exponents in some versions of K2GWB, the CO\s2 parameters (especially in the first block) may not match those in GWB files produced by K2GWB for the same values of \T.
Dataset Formats
The conversion is compatible with all existing dataset formats of GWB files.
For dataset formats ‘feb94’ and ‘oct94’, logKcalc updates the \logK blocks for the Eh reaction and O\s2 gas solubility.
(Blocks for H\s2 and N\s2 gas solubility, which are not used by GWB, are not updated.)
These blocks are not present in dataset formats ‘oct13’ and later.
Instead, the free electron is present as a species; logKcalc updates its \logK values.
Dataset formats ‘jul17’ and ‘jan19’ specify a fugacity model (tsonopoulos, peng-robinson, or spycher-reed).
The parameters are located in the data blocks for individual gas species and are not modified by logKcalc.
References
Bethke, C. M. and Farrell, B. (2020) The Geochemist's Workbench® Release 14 GWB Reference Manual, Aqeuous Solutions, LLC, Champaign, Illinois.
Cleverley, J. S. and Bastrakov, E. N. (2005) K2GWB: Utility for generating thermodynamic data files for The Geochemist's Workbench® at 0-1000 °C and 1-5000 bar from UT2K and the UNITHERM database. Computers & Geosciences 31, 756–767. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.cageo.2005.01.007")}
Dick, J. M. (2019) CHNOSZ: Thermodynamic calculations and diagrams for geochemistry. Front. Earth Sci. 7:180. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.3389/feart.2019.00180")}
Helgeson, H. C., Kirkham, D. H. and Flowers, G. C. (1981) Theoretical prediction of the thermodynamic behavior of aqueous electrolytes at high pressures and temperatures. IV. Calculation of activity coefficients, osmotic coefficients, and apparent molal and standard and relative partial molal properties to 600\degC and 5 Kb. Am. J. Sci. 281, 1249–1516. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.2475/ajs.281.10.1249")}
Manning, C. E., Shock, E. L. and Sverjensky, D. A. (2013) The chemistry of carbon in aqueous fluids at crustal and upper-mantle conditions: Experimental and theoretical constraints. Rev. Mineral. Geochem. 75, 109–148. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.2138/rmg.2013.75.5")}
Shvarov, Y. and Bastrakov, E. (1999) HCh: A software package for geochemical equilibrium modelling. User's Guide. Australian Geological Survey Organisation 1999/25. http://pid.geoscience.gov.au/dataset/ga/25473
See Also
modOBIGT for preparing the thermodynamic database for the \logK calculation.
Examples
# Show the first few entries in the names mapping file
mapfile <- system.file("extdata/map_names.csv", package = "logKcalc")
head(read.csv(mapfile))
# Process a GWB file
# -- define input and output files
infile <- system.file("extdata/thermo_12elements.tdat", package = "logKcalc")
outfile <- tempfile("thermo", fileext = ".tdat")
# -- make some modifications to OBIGT
# (add minerals from SUPCRT92 database and deprecated H2O(g))
modOBIGT(c("addSUPCRT", "steam"))
# -- identify species to add to the output
ispecies <- info("pyrrhotite")
# -- process the file
logKcalc(infile, outfile, ispecies = ispecies)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.