### 1-D diagrams: logarithms of activities
## Aqueous sulfur species (after Seewald, 1997 and 2001)
basis("CHNOS+")
basis("pH", 5)
species(c("H2S", "S2-2", "S3-2", "S2O3-2", "S2O4-2",
"S3O6-2", "S5O6-2", "S2O6-2", "HSO3-", "SO2", "HSO4-"))
a <- affinity(O2=c(-50, -15), T=325, P=350)
e <- equilibrate(a, loga.balance=-2)
diagram(e, ylim=c(-30, 0), legend.x="topleft")
title(main=paste("Aqueous sulfur speciation, 325 degC, 350 bar\n",
"After Seewald, 1997"))
# try it with and without the loga.balance argument (total activity of
# the balanced quantity, in this case H2S aka sulfur)
## Degrees of formation of ionized forms of glycine
## After Fig. 1 of Aksu and Doyle, 2001
basis("CHNOS+")
species(ispecies <- info(c("glycinium", "glycine", "glycinate")))
a <- affinity(pH=c(0, 14))
e <- equilibrate(a)
diagram(e, alpha=TRUE, lwd=1)
title(main=paste("Degrees of formation of aqueous glycine species\n",
"after Aksu and Doyle, 2001"))
## Degrees of formation of ATP species as a function of
## temperature, after LaRowe and Helgeson, 2007, Fig. 10b
# to make a similar diagram, activity of Mg+2 here is set to
# 10^-4, which is different from LH07, who used 10^-3 total molality
basis(c("CO2", "NH3", "H2O", "H3PO4", "O2", "H+", "Mg+2"),
c(999, 999, 999, 999, 999, -5, -4))
species(c("HATP-3", "H2ATP-2", "MgATP-2", "MgHATP-"))
a <- affinity(T=c(0, 120, 25))
e <- equilibrate(a)
diagram(e, alpha=TRUE)
title(main=paste("Degrees of formation of ATP species,\n",
"pH=5, log(aMg+2)=-3. After LaRowe and Helgeson, 2007"),
cex.main=0.9)
### 2-D diagrams: predominance diagrams
### these use the maximum affinity method
## Fe-S-O at 200 deg C, after Helgeson, 1970
basis(c("Fe", "O2", "S2"))
species(c("iron", "ferrous-oxide", "magnetite",
"hematite", "pyrite", "pyrrhotite"))
# include the high-temperature phase of pyrrhotite
species("pyrrhotite", "cr2")
a <- affinity(S2=c(-50, 0), O2=c(-90, -10), T=200)
diagram(a, fill="heat")
title(main=paste("Fe-S-O, 200 degrees C, 1 bar",
"After Helgeson, 1970", sep="\n"))
## pe-pH diagram for hydrated iron sulfides,
## goethite and pyrite, after Majzlan et al., 2006
# add some of these species to the database
add.obigt()
basis(c("Fe+2", "SO4-2", "H2O", "H+", "e-"),
c(0, log10(3), log10(0.75), 999, 999))
species(c("rhomboclase", "ferricopiapite", "hydronium jarosite",
"goethite", "melanterite", "pyrite"))
a <- affinity(pH=c(-1, 4, 256), pe=c(-5, 23, 256))
d <- diagram(a, main="Fe-S-O-H, after Majzlan et al., 2006")
# the first four species show up along the top of the diagram
stopifnot(all.equal(unique(t(d$predominant)[256,]), 1:4))
water.lines(yaxis="pe")
text(3, 22, describe.basis(thermo$basis[2:3,], digits=2, oneline=TRUE))
text(3, 21, describe.property(c("T", "P"), c(25, 1), oneline=TRUE))
# reset the database
data(thermo)
## Aqueous Aluminum Species F-/OH-, after Tagirov and Schott, 2001
# some of the species are not in the default databse
add.obigt()
# the 999s have no effect on the diagram:
# pH and log_a(F-) are plotting variables
# O2 is not in the formation reactions
# Al+3 is the balanced quantity
basis(c("Al+3", "F-", "H+", "O2", "H2O"), c(rep(999, 4), 0))
species(c("Al+3", "Al(OH)4-", "AlOH+2", "Al(OH)2+", "Al(OH)3",
"AlF+2", "AlF2+", "AlF3", "AlF4-", "Al(OH)2F2-", "Al(OH)2F",
"AlOHF2"), "aq")
a <- affinity(pH=c(0, 10), "F-"=c(-1, -9), T=200)
diagram(a, fill="heat")
title(main=paste("Aqueous aluminium species, T=200 C, P=Psat\n",
"after Tagirov and Schott, 2001"))
# restore thermodynamic database to default
data(thermo)
## Temperature-Pressure: kayanite-sillimanite-andalusite
# cf. Fig. 49 of Helgeson et al., 1978
# this is a system of one component (Al2SiO5), but we need the same
# number of basis species as elements; and avoid using H2O or other
# aqueous species because of T/P limits of the water() calculations;
basis(c("corundum", "quartz", "oxygen"))
species(c("kyanite", "sillimanite", "andalusite"))
# database has transition temperatures of kyanite and andalusite
# at 1 bar only, so we permit calculation at higher temperatures
a <- affinity(T=c(200, 900, 99), P=c(0, 9000, 101), exceed.Ttr=TRUE)
d <- diagram(a, fill=NULL)
bexpr <- sapply(c("Al2O3", "SiO2", "H2O"), expr.species, simplify=FALSE)
btext <- substitute(Al2O3 - SiO2 - H2O, unlist(bexpr))
mtitle(c(as.expression(btext), "after Helgeson et al., 1978"))
# find the approximate position of the triple point
tp <- find.tp(d$predominant)
Ttp <- a$vals[[1]][tp[1, 2]]
Ptp <- rev(a$vals[[2]])[tp[1, 1]]
points(Ttp, Ptp, pch=10, cex=5)
# some testing of the overall geometry
stopifnot(species()$name[d$predominant[1, 1]]=="andalusite")
stopifnot(species()$name[d$predominant[1, 101]]=="kyanite")
stopifnot(species()$name[d$predominant[99, 101]]=="sillimanite")
## calculate the equilibrium logarithm of activity of a
## basis species in different reactions
basis("CHNOS")
species(c("ethanol", "lactic acid", "deoxyribose", "ribose"))
a <- affinity(T=c(0, 150))
diagram(a, what="O2", legend.x="topleft", col=rev(rainbow(4)), lwd=2)
title(main="equilibrium logfO2 for 1e-3 mol/kg of CO2 and ... ")
Run the code above in your browser using DataLab