---
title: "Regressing thermodynamic data"
subtitle: "EOSregress in [CHNOSZ](http://www.chnosz.net)"
author: "Jeffrey M. Dick"
date: "`r Sys.Date()`"
output:
tufte::tufte_html:
tufte_features: ["background"]
css: "vig.css"
toc: true
mathjax: null
tufte::tufte_handout:
citation_package: natbib
latex_engine: xelatex
tufte::tufte_book:
citation_package: natbib
latex_engine: xelatex
vignette: >
%\VignetteIndexEntry{Regressing thermodynamic data}
%\VignetteEngine{knitr::rmarkdown}
\usepackage[utf8]{inputenc}
bibliography: vig.bib
link-citations: yes
csl: elementa.csl
---
```{r width80, include=FALSE}
options(width = 80)
```
```{r digits6, include=FALSE}
options(digits = 6)
```
```{r HTML, include=FALSE}
## some frequently used HTML expressions
V0 <- "*V*°"
Cp0 <- "*C*_{P}°"
c1 <- "*c*_{1}"
c2 <- "*c*_{2}"
a1 <- "*a*_{1}"
a2 <- "*a*_{2}"
a3 <- "*a*_{3}"
a4 <- "*a*_{4}"
h4sio4 <- "H_{4}SiO_{4}"
sio2 <- "SiO_{2}"
h2o <- "H_{2}O"
ch4 <- "CH_{4}"
wPrTr <- "ω_{Pr,Tr}"
```
```{r setup, include=FALSE}
library(knitr)
# invalidate cache when the tufte version changes
opts_chunk$set(tidy = FALSE, cache.extra = packageVersion('tufte'))
options(htmltools.dir.version = FALSE)
# adjust plot margins
knit_hooks$set(small.mar = function(before, options, envir) {
if (before) par(mar = c(4.2, 4.2, .1, .1)) # smaller margin on top and right
})
# use pngquant to optimize PNG images
knit_hooks$set(pngquant = hook_pngquant)
# pngquant isn't available on R-Forge ...
if (!nzchar(Sys.which("pngquant"))) {
pngquant <- NULL
# save space by using a lower resolution
dpi <- 50
} else {
pngquant <- "--speed=1 --quality=0-50"
dpi <- 72
}
## colorize messages 20171031
## adapted from https://gist.github.com/yihui/2629886#file-knitr-color-msg-rnw
color_block = function(color) {
function(x, options) sprintf('%s

', color, x)
}
knit_hooks$set(warning = color_block('magenta'), error = color_block('red'), message = color_block('blue'))
```
This [vignette](index.html) demonstrates `EOSregress()` and related functions for the regression of heat capacity and volumetric data to obtain "equations of state" coefficients.
# A note on the equations
The CHNOSZ thermodynamic database uses the revised Helgeson-Kirkham-Flowers equations of state (EOS) for aqueous species.
Different terms in these equations give the "non-solvation" and "solvation" contributions to the standard molal heat capacity (`r Cp0`) and volume (`r V0`) as a function of temperature (*T*) and pressure (*P*).
The equations were originally described by @HKF81; for `r Cp0` the equation is
> `r Cp0` = `r c1` + `r c2` / (*T* - Θ)^{2} + ω*TX*
Here, `r c1` and `r c2` are the non-solvation parameters, and ω is the solvation parameter.
Θ is equal to 228 K, and *X* is one of the "Born functions" that relates the solvation process to the dielectric properties of water.
For neutral species, all of the parameters, including the "effective" value of ω, are constant.
However, **for charged species, ω has non-zero derivatives at *T* > 100 °C, increasing with temperature**, that depend on the charge and ionic radius.
Revisions by @TH88 and @SOJSH92 refined the equations to improve their high-*T* and *P* behavior.
# A note on the algorithms
The regression functions are essentially a wrapper around the [`water()`](http://www.chnosz.net/manual/water.html)-related functions in CHNOSZ and [`lm()` in R](http://stat.ethz.ch/R-manual/R-patched/library/stats/html/lm.html).
The former provide values of the Born functions.
Accordingly, numerical values for all of the *variables* in the equations (e.g. 1/(*T* - Θ)^{2} and *TX*) can be calculated at each experimental *T* and *P*.
Applying a linear model (`lm`), the *coefficients* in the model can be obtained.
In this case, they correspond to the HKF parameters, e.g. `c1` (the intercept), `c2` and ω.
The coefficients **in this model** are constants, restricing application of the model to neutral (uncharged) species, for which the high-temperature derivatives of ω are 0.
Further below, a procedure is described to iteratively solve the equations for charged species.
# An example for neutral species
This is from the first example from `EOSregress()`.
```{marginfigure}
The `?` indicates a documentation topic in R. To view it, type `?EOSregress` or `help(EOSregress)` at the R prompt.
```
Here, we regress experimental heat capacities of aqueous methane (`r ch4`) using the revised HKF equations.
First, load CHNOSZ and its database:
```{r library_CHNOSZ}
library(CHNOSZ)
data(thermo)
```
Now, read a data file with experimental measurements from @HW97.
Here we also convert J to cal and MPa to bar:
```{r Cpdat}
file <- system.file("extdata/cpetc/Cp.CH4.HW97.csv", package = "CHNOSZ")
Cpdat <- read.csv(file)
Cpdat$Cp <- convert(Cpdat$Cp, "cal")
Cpdat$P <- convert(Cpdat$P, "bar")
```
Next, we specify the terms in the HKF equations and perform the regression using `EOSregress()`.
In the second call to `EOSregress()` below, only for *T* < 600 K are included in the regression.
The coefficients here correspond to `r c1`, `r c2` and ω in the HKF equations.
```{r EOSregress}
var <- c("invTTheta2", "TXBorn")
Cplm_high <- EOSregress(Cpdat, var)
Cplm_low <- EOSregress(Cpdat, var, T.max = 600)
Cplm_low$coefficients
```
```{r EOSplot, fig.margin = TRUE, fig.cap = "Heat capacity of aqueous methane.", fig.width=3.5, fig.height=3.5, cache=TRUE, results="hide", message=FALSE, echo=FALSE, dpi=dpi, out.width=672, out.height=336, pngquant=pngquant}
EOSplot(Cpdat, coefficients = round(Cplm_low$coefficients, 1))
EOSplot(Cpdat, coeficients = Cplm_high, add = TRUE, lty = 3)
PS01_data <- EOScoeffs("CH4", "Cp")
EOSplot(Cpdat, coefficients = PS01_data, add = TRUE, lty=2, col="blue")
```
We can use `EOSplot()` to plot the data and fitted lines and show the coefficients in the legend.
The solid line shows the fit to the lower-temperature data.
The fit to all data, represented by the dotted line, doesn't capture the low-temperature trends in the data.
```{r EOSplot, eval=FALSE}
```
```{marginfigure}
Be aware that the lines shown by `EOSplot()` are calculated for a single pressure only, despite the temperature- and pressure-dependence of the data and regressions.
```
`EOScoeffs()` is a small function that is used to retrieve the HKF parameters in the database in CHNOSZ.
The dashed blue line shows calculated values for methane using these parameters, which are from @PS01.
Compare the database values with the regressed values shown in the legend of figure above.
Some differences are expected as the values in the database are derived from different regression techniques applied to different sets of data:
```{r Cpcoeffs, message=FALSE}
EOScoeffs("CH4", "Cp")
```
## Setting the value of omega
Given both high-temperature volumetric and calorimetric data for neutral species, the effective value of ω is most reliably regressed from the latter [@SSW01].
Let's regress volumetric data using a value of omega taken from the heat capacity regression.
First, read the data from @HWM96.
```{r Vdat}
file <- system.file("extdata/cpetc/V.CH4.HWM96.csv", package = "CHNOSZ")
Vdat <- read.csv(file)
Vdat$P <- convert(Vdat$P, "bar")
```
Compressibilities of species (measured or estimated) are implied by the full set of HKF volumetric parameters (`r a1`, `r a2`, `r a3`, `r a4`).
In this example we model volumes at nearly constant *P*.
Therefore, we can use a simpler equation for `r V0` written in terms of the "isobaric fit parameters" (Tanger and Helgeson 1988, p. 35) σ and ξ, together with the solvation contribution that depends on the *Q* Born function:
> `r V0` = σ + ξ / (*T* - Θ) - ω*Q*
```{marginfigure}
`EOSvar()` actually returns the negative of *Q*, so the omega symbol here carries no negative sign.
```
Now we calculate the *Q* Born function using `EOSvar()` and multiply by ω (from the heat capacity regression) to get the solvation volume at each experimental temperature and pressure.
Subtract the solvation volume from the experimental volumes and create a new data frame holding the calculated "non-solvation" volume.
```{marginfigure}
Because we are dealing with volumes, the units of ω are converted according to 1 cal = 41.84 cm^{3}∙bar.
```
```{r Vdat_non}
QBorn <- EOSvar("QBorn", T = Vdat$T, P = Vdat$P)
Vomega <- convert(Cplm_low$coefficients[["TXBorn"]], "cm3bar")
V_sol <- Vomega * QBorn
V_non <- Vdat$V - V_sol
Vdat_non <- data.frame(T = Vdat$T, P = Vdat$P, V = V_non)
```
Next, regress the non-solvation volume using the non-solvation terms in the HKF model.
As with `r Cp0`, also get the values of the parameters from the database for comparison with the regression results.
```{r Vdat_non_regress, message=FALSE}
var <- "invTTheta"
Vnonlm <- EOSregress(Vdat_non, var, T.max = 450)
Vcoeffs <- round(c(Vnonlm$coefficients, QBorn = Vomega), 1)
Vcoeffs_database <- EOScoeffs("CH4", "V")
```
```{r Vplot, fig.margin=TRUE, results="hide", message=FALSE, echo=FALSE, fig.width=3.5, fig.height=7, fig.cap="Volume of aqueous methane.", dpi=dpi, out.width=672, out.height=672, pngquant=pngquant}
par(mfrow = c(2, 1))
# plot 1
EOSplot(Vdat, coefficients = Vcoeffs)
EOSplot(Vdat, coefficients = Vcoeffs_database, add = TRUE, lty = 2)
# plot 2
EOSplot(Vdat, coefficients = Vcoeffs_database, T.plot = 600, lty = 2)
EOSplot(Vdat, coefficients = Vcoeffs, add = TRUE)
```
Finally, plot the data and regressions.
The first plot shows all the data, and the second the low-temperature subset (*T* < 600 K).
The solid line is the two-term fit for σ and ξ (using ω from the heat capacity regression), and the dashed line shows the volumes calculated using the parameters in the database.
The plot legends give the parameters from the two-term fit (first plot), or from the database (second plot).
```{r Vplot, eval=FALSE}
```
The equation for `r V0` provides a reasonable approximation of the trend of lowest-temperature data (*T* < 450 K).
However, the equation does not closely reproduce the trend of higher-temperature `r V0` data (*T* < 600 K), nor behavior in the critical region.
Because of these issues, some researchers are exploring alternatives to the HKF model for aqueous nonelectrolytes.
(See also an example in `?EOSregress`.)
# An example for charged species
For this example, let's generate synthetic data for Na^{+} using its parameters in the database.
In the call to `subcrt()` below, `convert = FALSE` means to take *T* in units of K.
```{r Nadat}
T <- convert(seq(0, 600, 50), "K")
P <- 1000
prop.PT <- subcrt("Na+", T = T, P = P, grid = "T", convert = FALSE)$out[[1]]
Nadat <- prop.PT[, c("T", "P", "Cp")]
```
As noted above, ω for electrolytes is not a constant.
What happens if we apply the constant-ω model anyway, knowing it's not applicable (especially at high temperature)?
```{r Nalm, fig.margin=TRUE, fig.width=3.5, fig.height=3.5, fig.cap="Heat capacity of Na^{+} (inapplicable: constant ω).", dpi=dpi, out.width=672, out.height=336, pngquant=pngquant}
var <- c("invTTheta2", "TXBorn")
Nalm <- EOSregress(Nadat, var, T.max = 600)
EOSplot(Nadat, coefficients = Nalm$coefficients, fun.legend = NULL)
EOSplot(Nadat, add = TRUE, lty = 3)
```
As before, the solid line is a fit to relatively low-temperature (*T* < 600 K) data, and the dotted line a fit to the entire temperature range of the data.
The fits using constant ω are clearly not acceptable.
There is, however, a way out.
A different variable, `Cp_s_var`, can be used to specify the calculation of the "solvation" heat capacity in the HKF equations using the temperature- and pressure-dependent corrections for charged species.
To use this variable, the values of `r wPrTr` (omega at the reference temperature and pressure) and *Z* (charge) must be given, in addition to *T* and *P*.
Of course, right now we *don't know* the value of `r wPrTr`---it is the purpose of the regression to find it!
But we can make a first guess using the value of ω found above.
```{r Navars1}
var1 <- c("invTTheta2", "Cp_s_var")
omega.guess <- coef(Nalm)[3]
```
Then, we can use an iterative procedure that refines successive guesses of `r wPrTr`.
The convergence criterion is measured by the difference in sequential regressed values of ω.
```{r Nawhile, fig.margin=TRUE, fig.width=3.5, fig.height=3.5, fig.cap="Heat capacity of Na^{+} (variable ω).", dpi=dpi, out.width=672, out.height=336, pngquant=pngquant}
diff.omega <- 999
while(abs(diff.omega) > 1) {
Nalm1 <- EOSregress(Nadat, var1, omega.PrTr = tail(omega.guess, 1), Z = 1)
omega.guess <- c(omega.guess, coef(Nalm1)[3])
diff.omega <- tail(diff(omega.guess), 1)
}
EOSplot(Nadat, coefficients = signif(coef(Nalm1), 6),
omega.PrTr = tail(omega.guess, 1), Z = 1)
EOScoeffs("Na+", "Cp")
```
Alrighty! We managed to obtain HKF coefficients from synthetic data for Na^{+}.
The regressed values of the HKF coefficients (shown in the plot legend) are very close to the database values (printed by the call to `EOScoeffs()`) used to generate the synthetic data.
## Doing it for volume
Just like above, but using synthetic `r V0` data.
Note that the regressed value of ω has volumetric units (cm^{3}∙bar/mol), while `omega.PrTr` is in caloric units (cal/mol).
Compared to `r Cp0`, the regression of `r V0` is very finicky.
Given a starting guess of `r wPrTr` of 1400000 cm^{3}∙bar/mol, the iteration converges on 1394890 instead of the "true" database value of 1383230 (represented by dashed line in the plot).
```{r NaVolume, fig.margin=TRUE, fig.width=3.5, fig.height=3.5, fig.cap="Volume of Na^{+} (variable ω).", results="hide", message=FALSE, echo=FALSE, dpi=dpi, out.width=672, out.height=336, pngquant=pngquant}
T <- convert(seq(0, 600, 25), "K")
P <- 1000
prop.PT <- subcrt("Na+", T = T, P = P, grid = "T", convert = FALSE)$out[[1]]
NaVdat <- prop.PT[, c("T", "P", "V")]
var1 <- c("invTTheta", "V_s_var")
omega.guess <- 1400000
diff.omega <- 999
while(abs(diff.omega) > 1) {
NaVlm1 <- EOSregress(NaVdat, var1,
omega.PrTr = tail(convert(omega.guess, "calories"), 1), Z = 1)
omega.guess <- c(omega.guess, coef(NaVlm1)[3])
diff.omega <- tail(diff(omega.guess), 1)
}
EOSplot(NaVdat, coefficients = signif(coef(NaVlm1), 6),
omega.PrTr = tail(convert(omega.guess, "calories"), 1), Z = 1,
fun.legend = "bottomleft")
coefficients <- EOScoeffs("Na+", "V", P = 1000)
names(coefficients)[3] <- "V_s_var"
EOSplot(NaVdat, coefficients = coefficients, Z = 1, add = TRUE, lty = 2,
omega.PrTr = convert(coefficients["V_s_var"], "calories"))
```
```{r NaVolume, eval=FALSE}
```
# Making a pseudospecies: `r h4sio4`
Some mineral stability diagrams use the activity of `r h4sio4` as a variable.
However, the primary species for dissolved silica in CHNOSZ is `r sio2`(aq) [@SHS89].
The pseudo-reaction with zero properties, `r h4sio4` = `r sio2` + 2 `r h2o`, defines the properties of the pseudospecies `r h4sio4`.
First we go about calculating the properties of `r sio2` + 2 `r h2o`.
We do this over a range of *T* and *P*, but include many points near 25 °C to improve the fit of the regression in that region:
`r op <- options(warn = -1)`
```{r SiO2_2H2O, message=FALSE}
s_25C <- subcrt(c("SiO2", "H2O"), c(1, 2), T = 25)$out
s_near25 <- subcrt(c("SiO2", "H2O"), c(1, 2), T = seq(20, 30, length.out=50))$out
s_lowT <- subcrt(c("SiO2", "H2O"), c(1, 2), T = seq(0, 100, 10))$out
s_Psat <- subcrt(c("SiO2", "H2O"), c(1, 2))$out
s_P500 <- subcrt(c("SiO2", "H2O"), c(1, 2), T = seq(0, 1000, 100), P = 500)$out
s_P1000 <- subcrt(c("SiO2", "H2O"), c(1, 2), T = seq(0, 1000, 100), P = 1000)$out
```
`r options(op)`
Now we can start making the new species, with thermodynamic properties calculated at 25 °C:
```{r new_H4SiO4}
mod.obigt("new-H4SiO4", formula = "H4SiO4", ref1 = "this_vignette",
date = today(), G = s_25C$G, H = s_25C$H, S = s_25C$S,
Cp = s_25C$Cp, V = s_25C$V, z = 0)
```
To prepare for the regression, combine the calculated data and convert °C to K:
```{r substuff}
substuff <- rbind(s_near25, s_lowT, s_Psat, s_P500, s_P1000)
substuff$T <- convert(substuff$T, "K")
```
Now let's run a `r Cp0` regression and update the new species with the regressed HKF coefficients.
Note that we apply order-of-magnitude scaling to the coefficients (see `?thermo`):
```{r Cp_H4SiO4, results="hide"}
Cpdat <- substuff[, c("T", "P", "Cp")]
var <- c("invTTheta2", "TXBorn")
Cplm <- EOSregress(Cpdat, var)
Cpcoeffs <- Cplm$coefficients
mod.obigt("new-H4SiO4", c1 = Cpcoeffs[1],
c2 = Cpcoeffs[2]/10000, omega = Cpcoeffs[3]/100000)
```
Let's get ready to regress `r V0` data.
We use the strategy shown above to calculate non-solvation volume using ω from the `r Cp0` regression:
```{r V_H4SiO4_nonsolvation}
Vdat <- substuff[, c("T", "P", "V")]
QBorn <- EOSvar("QBorn", T = Vdat$T, P = Vdat$P)
Vomega <- convert(Cplm$coefficients[["TXBorn"]], "cm3bar")
V_sol <- Vomega * QBorn
V_non <- Vdat$V - V_sol
Vdat$V <- V_non
```
Here's the `r V0` regression for the pseudospecies.
We specify the variables for the `r a1`, `r a2`, `r a3`, and `r a4` terms in the HKF equations.
```{r V_H4SiO4, results="hide"}
var <- c("invPPsi", "invTTheta", "invPPsiTTheta")
Vlm <- EOSregress(Vdat, var)
Vcoeffs <- convert(Vlm$coefficients, "calories")
mod.obigt("new-H4SiO4", a1 = Vcoeffs[1]*10, a2 = Vcoeffs[2]/100,
a3 = Vcoeffs[3], a4 = Vcoeffs[4]/10000)
```
We just made a new `r h4sio4` pseudospecies.
In February 2017, the database in CHNOSZ was updated with a `pseudo-H4SiO4` that was made using this procedure.
For comparison, the parameters of `r h4sio4` given by @Ste01, which are contained in an optional data file, can be read using `add.obigt()`:
```{r width180, include=FALSE}
options(width = 180)
```
```{r info_H4SiO4, message=FALSE}
add.obigt("SUPCRTBL", "H4SiO4")
info(info(c("new-H4SiO4", "pseudo-H4SiO4", "H4SiO4")))
```
```{r width80, include=FALSE}
```
```{r subcrt_H4SiO4, fig.margin=TRUE, fig.width=4, fig.height=4, small.mar=TRUE, echo=FALSE, results="hide", message=FALSE, dpi=dpi, out.width="100%", cache=TRUE, fig.cap="Comparison of H_{4}SiO_{4} pseudospecies.", pngquant=pngquant}
s1 <- subcrt(c("pseudo-H4SiO4", "SiO2", "H2O"), c(-1, 1, 2))
plot(s1$out$T, s1$out$G, type = "l", ylim = c(-100, 600),
xlab = axis.label("T"), ylab = axis.label("DG0"))
s2 <- subcrt(c("H4SiO4", "SiO2", "H2O"), c(-1, 1, 2))
lines(s2$out$T, s2$out$G, lty = 2)
abline(h = 0, lty = 3)
legend("topright", legend = c("pseudo-H4SiO4 (CHNOSZ)",
"H4SiO4 (Stefánsson, 2001)"), lty = c(1, 2, NA), bty = "n")
text(225, 250, describe.reaction(s1$reaction))
```
Let's compare CHNOSZ's `pseudo-H4SiO4` (with the same properties as the `new-H4SiO4` we just made) and `H4SiO4` from Stefánsson (2001) with the calculated properties of `r sio2` + 2 `r h2o` as a function of temperature:
```{r subcrt_H4SiO4, eval=FALSE}
```
The large differences for `H4SiO4` at low temperature agree with the comparison to Shock et al. (1989) shown in Figure 3 of Stefánsson (2001).
Therefore, Stefánsson's `r h4sio4` should not be used with the _default_ database in CHNOSZ for making mineral activity diagrams.
Instead, the pseudospecies with properties calculated here, `pseudo-H4SiO4`, is preferable for use with the _default_ database in CHNOSZ---see `?transfer` and [`demo(activity_ratios)`](../demo).
On the other hand, Stefánsson's `r h4sio4` is compatible with `r sio2` in the _optional_ SUPCRTBL data file [@ZZL_16].
# Other possibilities
These functions are limited to treatment of calorimetric and volumetric data.
Other software tools have been described recently for deriving HKF parameters and other thermodynamic properties from different types of experimental data [@MKDW15; @Shv15].