logB.to.OBIGT: Fit Thermodynamic Parameters to Formation Constants (\logB)
In CHNOSZ: Thermodynamic Calculations and Diagrams for Geochemistry
View source: R/logB.to.OBIGT.R
logB.to.OBIGT R Documentation
Fit Thermodynamic Parameters to Formation Constants (\logB)
Description
Fit thermodynamic parameters to experimental formation constants for an aqueous species and add the parameters to OBIGT.
Usage
logB.to.OBIGT(logB, species, coeffs, T, P, npar = 3,
optimize.omega = FALSE, tolerance = 0.05)
Arguments
logB
Values of \logB
species
Species in the formation reaction (the formed species must be last)
coeffs
Coefficients of the formation reaction
T
Temperature (units of T.units)
P
Temperature (‘Psat’ or numerical values with units of P.units)
npar
numeric, number of parameters to use in fitting
optimize.omega
logical, optimize the omega parameter (relevant for charged species)?
tolerance
Tolerance for checking equivalence of input and calculated \logB values
Details
This function updates the OBIGT thermodynamic database with parameters fit to formation constants of aqueous species as a function of temperature.
The logB argument should have decimal logarithm of formation constants for an aqueous complex (\logB).
The formation reaction is defined by species and coeffs.
All species in the formation reaction must be present in OBIGT except for the last species, which is the newly formed species.
The function works as follows.
First, the properties of the known species are combined with \logB to calculate the standard Gibbs energy (G[T]) of the formed species at each value of T and P.
Then, lm is used to fit one or more of the thermodynamic parameters G, S, c1, c2, and omega to the values of G[T].
The first two of these parameters are standard-state values at 25 \degC and 1 bar, and the last three are parameters in the revised HKF equations (e.g. Eq. B25 of Shock et al., 1992).
The fitted parameters for the formed species are then added to OBIGT.
Finally, all.equal is used to test for approximate equivalence of the input values of \logB and calculated equilibrium constants; if the mean absolute difference exceeds tolerance, an error occurs.
To avoid overfitting, only the first three parameters (G, S, and c1) are used by default.
More parameters (up to 5) or fewer (down to 1) can be selected by changing npar.
Volumetric parameters (\a1 to \a4) in the HKF equations currently aren't included, so the resulting fits are valid only at the input pressure values.
Independent of npar, the number of parameters used in the fit is not more than the number of data values (n).
If n is less than 5, then the values of the unused parameters are set to 0 for addition to OBIGT.
For instance, a single value of \logB would be fit only with G, implying that computed values of G[T] have no temperature dependence.
The value of \omega is a constant in the revised HKF equations for uncharged species, but for charged species, it is a function of \T and \P as described by the “g function” (Shock et al., 1992).
An optimization step is available to refine the value of omega at 25 \degC and 1 bar for charged species taking its temperature dependence into account for the fitting.
However, representative datasets are not well represented by variable omega (see first example), so this step is skipped by default.
Furthermore, logB.to.OBIGT by default also sets the z parameter in OBIGT to 0; this enforces a constant \omega for charged species in calculations with subcrt (note that this is a parameter for the HKF equations and does not affect reaction balancing).
Set optimize.omega to TRUE to override the defaults and activate variable \omega for charged species; this takes effect only if npar == 5 and n > 5.
Value
The species index of the new species in OBIGT.
References
Migdisov, Art. A., Zezin, D. and Williams-Jones, A. E. (2011) An experimental study of Cobalt (II) complexation in Cl\S- and H\s2S-bearing hydrothermal solutions. Geochim. Cosmochim. Acta 75, 4065–4079. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.gca.2011.05.003")}
Mei, Y., Sherman, D. M., Liu, W., Etschmann, B., Testemale, D. and Brugger, J. (2015) Zinc complexation in chloride-rich hydrothermal fluids (25–600 \degC): A thermodynamic model derived from ab initio molecular dynamics. Geochim. Cosmochim. Acta 150, 264–284. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.gca.2014.09.023")}
Shock, E. L., Oelkers, E. H., Johnson, J. W., Sverjensky, D. A. and Helgeson, H. C. (1992) Calculation of the thermodynamic properties of aqueous species at high pressures and temperatures: Effective electrostatic radii, dissociation constants and standard partial molal properties to 1000 \degC and 5 kbar. J. Chem. Soc. Faraday Trans. 88, 803–826. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1039/FT9928800803")}
See Also
logB.to.OBIGT calls mod.OBIGT with zap = TRUE to clear the parameters of a formed species if it already exists in the OBIGT database.
If preexisting parameters (e.g. volumetric HKF coefficients) weren't cleared, they would interfere with the fitting done here, which uses only selected parameters.
Examples
## CoHS+ from Migdisov et al. (2011)
logB <- c(6.24, 6.02, 5.84, 5.97, 6.52)
T <- c(120, 150, 200, 250, 300)
P <- "Psat"
species <- c("Co+2", "HS-", "CoHS+")
coeffs <- c(-1, -1, 1)
opar <- par(mfrow = c(2, 1))
for(o.o in c(TRUE, FALSE)) {
# Fit the parameters with or without variable omega
inew <- logB.to.OBIGT(logB, species, coeffs, T, P, npar = 5, optimize.omega = o.o)
# Print the new database entry
info(inew)
# Plot experimental logB
plot(T, logB, "n", c(100, 320), c(5.8, 6.8),
xlab = axis.label("T"), ylab = quote(log~beta))
points(T, logB, pch = 19, cex = 2)
# Plot calculated values
Tfit <- seq(100, 320, 10)
sres <- subcrt(species, coeffs, T = Tfit)
lines(sres$out$T, sres$out$logK, col = 4)
title(describe.reaction(sres$reaction))
legend <- c("Migdisov et al. (2011)",paste0("logB.to.OBIGT(optimize.omega = ",o.o,")"))
legend("top", legend, pch = c(19, NA), lty = c(0, 1), col = c(1, 4),
pt.cex = 2, bg = "#FFFFFFB0")
}
par(opar)
# NB. Optimizing omega leads to unphysical oscillations in the logK (first plot)
## ZnCl+ from Mei et al. (2015)
# Values for 5000 bar calculated with the modified Ryzhenko-Bryzgalin (RB) model
logB <- c(-1.93,-1.16,-0.38,0.45,1.15,1.76,2.30,2.80,3.26,3.70,4.12,4.53,4.92)
species <- c("Zn+2", "Cl-", "ZnCl+")
coeffs <- c(-1, -1, 1)
T <- c(25, 60, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600)
P <- 5000
# Note: ZnCl+ is in the default OBIGT database;
# logB.to.OBIGT() "zaps" the previous parameters before adding the fitted ones
inew <- logB.to.OBIGT(logB, species, coeffs, T, P, npar = 5)
# Plot RB and logB.to.OBIGT values
plot(T, logB, xlab = axis.label("T"), ylab = axis.label("logB"), pch = 19, cex = 2)
Tfit <- seq(25, 600, 25)
sres <- subcrt(species, coeffs, T = Tfit, P = P)
lines(sres$out$T, sres$out$logK, col = 4)
title(describe.reaction(sres$reaction), line = 3)
title("5000 bar", font.main = 1, line = 1)
legend <- c("Mei et al. (2015)", "logB.to.OBIGT()")
legend("topleft", legend, pch = c(19, NA), lty = c(0, 1), col = c(1, 4), pt.cex = 2)
CHNOSZ documentation built on May 29, 2024, 3:30 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.
View source: R/logB.to.OBIGT.R
| logB.to.OBIGT | R Documentation |
Fit Thermodynamic Parameters to Formation Constants (\logB)
Description
Fit thermodynamic parameters to experimental formation constants for an aqueous species and add the parameters to OBIGT.
Usage
logB.to.OBIGT(logB, species, coeffs, T, P, npar = 3,
optimize.omega = FALSE, tolerance = 0.05)
Arguments
logB |
Values of \logB |
species |
Species in the formation reaction (the formed species must be last) |
coeffs |
Coefficients of the formation reaction |
T |
Temperature (units of |
P |
Temperature (‘Psat’ or numerical values with units of |
npar |
numeric, number of parameters to use in fitting |
optimize.omega |
logical, optimize the omega parameter (relevant for charged species)? |
tolerance |
Tolerance for checking equivalence of input and calculated \logB values |
Details
This function updates the OBIGT thermodynamic database with parameters fit to formation constants of aqueous species as a function of temperature.
The logB argument should have decimal logarithm of formation constants for an aqueous complex (\logB).
The formation reaction is defined by species and coeffs.
All species in the formation reaction must be present in OBIGT except for the last species, which is the newly formed species.
The function works as follows.
First, the properties of the known species are combined with \logB to calculate the standard Gibbs energy (G[T]) of the formed species at each value of T and P.
Then, lm is used to fit one or more of the thermodynamic parameters G, S, c1, c2, and omega to the values of G[T].
The first two of these parameters are standard-state values at 25 \degC and 1 bar, and the last three are parameters in the revised HKF equations (e.g. Eq. B25 of Shock et al., 1992).
The fitted parameters for the formed species are then added to OBIGT.
Finally, all.equal is used to test for approximate equivalence of the input values of \logB and calculated equilibrium constants; if the mean absolute difference exceeds tolerance, an error occurs.
To avoid overfitting, only the first three parameters (G, S, and c1) are used by default.
More parameters (up to 5) or fewer (down to 1) can be selected by changing npar.
Volumetric parameters (\a1 to \a4) in the HKF equations currently aren't included, so the resulting fits are valid only at the input pressure values.
Independent of npar, the number of parameters used in the fit is not more than the number of data values (n).
If n is less than 5, then the values of the unused parameters are set to 0 for addition to OBIGT.
For instance, a single value of \logB would be fit only with G, implying that computed values of G[T] have no temperature dependence.
The value of \omega is a constant in the revised HKF equations for uncharged species, but for charged species, it is a function of \T and \P as described by the “g function” (Shock et al., 1992).
An optimization step is available to refine the value of omega at 25 \degC and 1 bar for charged species taking its temperature dependence into account for the fitting.
However, representative datasets are not well represented by variable omega (see first example), so this step is skipped by default.
Furthermore, logB.to.OBIGT by default also sets the z parameter in OBIGT to 0; this enforces a constant \omega for charged species in calculations with subcrt (note that this is a parameter for the HKF equations and does not affect reaction balancing).
Set optimize.omega to TRUE to override the defaults and activate variable \omega for charged species; this takes effect only if npar == 5 and n > 5.
Value
The species index of the new species in OBIGT.
References
Migdisov, Art. A., Zezin, D. and Williams-Jones, A. E. (2011) An experimental study of Cobalt (II) complexation in Cl\S- and H\s2S-bearing hydrothermal solutions. Geochim. Cosmochim. Acta 75, 4065–4079. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.gca.2011.05.003")}
Mei, Y., Sherman, D. M., Liu, W., Etschmann, B., Testemale, D. and Brugger, J. (2015) Zinc complexation in chloride-rich hydrothermal fluids (25–600 \degC): A thermodynamic model derived from ab initio molecular dynamics. Geochim. Cosmochim. Acta 150, 264–284. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1016/j.gca.2014.09.023")}
Shock, E. L., Oelkers, E. H., Johnson, J. W., Sverjensky, D. A. and Helgeson, H. C. (1992) Calculation of the thermodynamic properties of aqueous species at high pressures and temperatures: Effective electrostatic radii, dissociation constants and standard partial molal properties to 1000 \degC and 5 kbar. J. Chem. Soc. Faraday Trans. 88, 803–826. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.1039/FT9928800803")}
See Also
logB.to.OBIGT calls mod.OBIGT with zap = TRUE to clear the parameters of a formed species if it already exists in the OBIGT database.
If preexisting parameters (e.g. volumetric HKF coefficients) weren't cleared, they would interfere with the fitting done here, which uses only selected parameters.
Examples
## CoHS+ from Migdisov et al. (2011)
logB <- c(6.24, 6.02, 5.84, 5.97, 6.52)
T <- c(120, 150, 200, 250, 300)
P <- "Psat"
species <- c("Co+2", "HS-", "CoHS+")
coeffs <- c(-1, -1, 1)
opar <- par(mfrow = c(2, 1))
for(o.o in c(TRUE, FALSE)) {
# Fit the parameters with or without variable omega
inew <- logB.to.OBIGT(logB, species, coeffs, T, P, npar = 5, optimize.omega = o.o)
# Print the new database entry
info(inew)
# Plot experimental logB
plot(T, logB, "n", c(100, 320), c(5.8, 6.8),
xlab = axis.label("T"), ylab = quote(log~beta))
points(T, logB, pch = 19, cex = 2)
# Plot calculated values
Tfit <- seq(100, 320, 10)
sres <- subcrt(species, coeffs, T = Tfit)
lines(sres$out$T, sres$out$logK, col = 4)
title(describe.reaction(sres$reaction))
legend <- c("Migdisov et al. (2011)",paste0("logB.to.OBIGT(optimize.omega = ",o.o,")"))
legend("top", legend, pch = c(19, NA), lty = c(0, 1), col = c(1, 4),
pt.cex = 2, bg = "#FFFFFFB0")
}
par(opar)
# NB. Optimizing omega leads to unphysical oscillations in the logK (first plot)
## ZnCl+ from Mei et al. (2015)
# Values for 5000 bar calculated with the modified Ryzhenko-Bryzgalin (RB) model
logB <- c(-1.93,-1.16,-0.38,0.45,1.15,1.76,2.30,2.80,3.26,3.70,4.12,4.53,4.92)
species <- c("Zn+2", "Cl-", "ZnCl+")
coeffs <- c(-1, -1, 1)
T <- c(25, 60, 100, 150, 200, 250, 300, 350, 400, 450, 500, 550, 600)
P <- 5000
# Note: ZnCl+ is in the default OBIGT database;
# logB.to.OBIGT() "zaps" the previous parameters before adding the fitted ones
inew <- logB.to.OBIGT(logB, species, coeffs, T, P, npar = 5)
# Plot RB and logB.to.OBIGT values
plot(T, logB, xlab = axis.label("T"), ylab = axis.label("logB"), pch = 19, cex = 2)
Tfit <- seq(25, 600, 25)
sres <- subcrt(species, coeffs, T = Tfit, P = P)
lines(sres$out$T, sres$out$logK, col = 4)
title(describe.reaction(sres$reaction), line = 3)
title("5000 bar", font.main = 1, line = 1)
legend <- c("Mei et al. (2015)", "logB.to.OBIGT()")
legend("topleft", legend, pch = c(19, NA), lty = c(0, 1), col = c(1, 4), pt.cex = 2)
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.