reset()
## each of the available methods
nonideal("Alberty")
nonideal("bgamma0")
nonideal("bgamma")
nonideal("Bdot0")
nonideal("Bdot") # the default
## what's the activity coefficient of Na+ at
## 25 degC and 1 bar and an ionic strength of 0.7?
sres <- subcrt("Na+", T = 25, IS = 0.7)
# exponentiate to convert log10(gamma) to gamma
10^sres$out[[1]]$loggam
# now use a different method
nonideal("bgamma")
sres <- subcrt("Na+", T = 25, IS = 0.7)
10^sres$out[[1]]$loggam
## what are activity coefficients of -3, -2, -1, 0, +1, +2, +3 charged species
## as a function of ionic strength and temperature?
# first choose the method
nonideal("Bdot")
# define the ionic strength and temperature increments
IS <- c(0.001, 0.01, 0.1, 0.7)
T <- seq(0, 100, 25)
# use species charged -3, -2, -1, 0, +1, +2, +3
species <- c("PO4-3", "HPO4-2", "H2PO4-", "H3PO4", "Na+", "Ca+2", "Al+3")
# initialize empty output table for T (rows) and charge (columns)
gamtab <- matrix(nrow = length(T), ncol = length(species))
rownames(gamtab) <- T
colnames(gamtab) <- -3:3
# make a list of tables to hold the activity coefficients, one for each IS
gamma <- rep(list(gamtab), length(IS))
names(gamma) <- IS
# loop over the values of ionic strength
for(i in seq_along(IS)) {
# calculate properties of species, including logarithm of activity coefficient
sres <- subcrt(species, T = T, IS = IS[i])
# exponentiate to convert log10(gamma) to gamma, and put the values into the tables
for(j in seq_along(species)) gamma[[i]][, j] <- 10^sres$out[[j]]$loggam
}
# print the output and make a plot
print(gamma)
matplot(T, gamma$`0.001`, type = "l")
title(main = "activity coefficients of -3, -2, -1, 0, +1, +2, +3 charged species")
## Alberty, 2003 p. 16 Table 1.3: adjusted pKa of acetic acid
## we use the 'IS' argument in subcrt() to calculate adjusted thermodynamic properties
# set ideal.H to FALSE to calculate activity coefficients for the proton
# (needed for replication of the values in Alberty's book)
nonideal("Alberty")
thermo("opt$ideal.H" = FALSE)
sres <- subcrt(c("acetate", "H+", "acetic acid"), c(-1, -1, 1),
IS=c(0, 0.1, 0.25), T=25, property="logK")
# we're within 0.01 of Alberty's pK values
Alberty_logK <- c(4.75, 4.54, 4.47)
# The maximum (absolute) pairwise difference between x and y
max(abs(Alberty_logK - sres$out$logK)) # 0.0072
# reset option to default
thermo("opt$ideal.H" = TRUE)
## An example using IS with affinity():
## speciation of phosphate as a function of ionic strength
opar <- par(mfrow = c(2, 1))
basis("CHNOPS+")
Ts <- c(25, 100)
species(c("PO4-3", "HPO4-2", "H2PO4-"))
for(T in Ts) {
a <- affinity(IS = c(0, 0.14), T = T)
e <- equilibrate(a)
if(T==25) diagram(e, ylim = c(-3.0, -2.6), legend.x = NULL)
else diagram(e, add = TRUE, names = FALSE, col = "red")
}
title(main="Non-ideality model for phosphate species")
dp <- describe.property(c("pH", "T", "T"), c(7, Ts))
legend("topright", lty = c(NA, 1, 1), col = c(NA, "black", "red"), legend = dp)
text(0.07, -2.76, expr.species("HPO4-2"))
text(0.07, -2.90, expr.species("H2PO4-"))
## phosphate predominance f(IS,pH)
a <- affinity(IS = c(0, 0.14), pH = c(6, 13), T = Ts[1])
d <- diagram(a, fill = NULL)
a <- affinity(IS = c(0, 0.14), pH = c(6, 13), T = Ts[2])
d <- diagram(a, add = TRUE, names = FALSE, col = "red")
par(opar)
## activity coefficients for monovalent ions at 700 degC, 10 kbar
# after Manning, 2013, Fig. 7
# here we use the b_gamma equation
nonideal("bgamma")
IS <- c(0.001, 0.01, 0.1, 1, 2, 2.79)
# we're above 5000 bar, so need to use IAPWS-95 or DEW
oldwat <- water("DEW")
sres <- subcrt("Na+", T = 700, P = 10000, IS = IS)
water(oldwat)
# compare the calculated activity coefficient to values from Manning's figure
gamma <- 10^sres$out[[1]]$loggam
Manning_gamma <- c(0.93, 0.82, 0.65, 0.76, 1.28, 2)
gamma - Manning_gamma
## Plot the data and splines used for calculating b_gamma
## (extended term parameter)
bgamma(showsplines = "T")
bgamma(showsplines = "P")
## a longer example, using nonideal() directly
# Alberty, 2003 p. 273-276: activity coefficient (gamma)
# as a function of ionic strength and temperature
nonideal("Alberty")
IS <- seq(0, 0.25, 0.005)
T <- c(0, 25, 40)
lty <- 1:3
species <- c("H2PO4-", "HADP-2", "HATP-3", "ATP-4")
col <- rainbow(4)
thermo.plot.new(xlim = range(IS), ylim = c(0, 1),
xlab = axis.label("IS"), ylab = "gamma")
for(j in 1:3) {
# use subcrt to generate speciesprops
speciesprops <- subcrt(species, T = rep(T[j], length(IS)))$out
# use nonideal to calculate loggamma; this also adjusts G, H, S, Cp,
# but we don't use them here
nonidealprops <- nonideal(species, speciesprops, IS = IS, T = convert(T[j], "K"))
for(i in 1:4) lines(IS, 10^(nonidealprops[[i]]$loggam), lty=lty[j], col=col[i])
}
t1 <- "Activity coefficient (gamma) of -1,-2,-3,-4 charged species"
t2 <- quote("at 0, 25, and 40 "*degree*"C, after Alberty, 2003")
mtitle(as.expression(c(t1, t2)))
legend("topright", lty=c(NA, 1:3), bty="n",
legend=c(as.expression(axis.label("T")), 0, 25, 40))
legend("top", lty=1, col=col, bty="n",
legend = as.expression(lapply(species, expr.species)))
## reset method to default
nonideal("Bdot") # or reset()
Run the code above in your browser using DataLab