Nothing
inst/doc/FAQ.R
In CHNOSZ: Thermodynamic Calculations and Diagrams for Geochemistry
## ----setup, include = FALSE---------------------------------------------------
library(CHNOSZ)
options(width = 80)
# Use pngquant to optimize PNG images
library(knitr)
knit_hooks$set(pngquant = hook_pngquant)
pngquant <- "--speed=1 --quality=0-25"
if (!nzchar(Sys.which("pngquant"))) pngquant <- NULL
# To make warnings appear in text box 20230619
# https://selbydavid.com/2017/06/18/rmarkdown-alerts/
knitr::knit_hooks$set(
error = function(x, options) {
paste('\n\n<div class="alert alert-danger">',
gsub('##', '\n', gsub('^##\ Error', '**Error:**', x)),
'</div>', sep = '\n')
},
warning = function(x, options) {
paste('\n\n<div class="alert alert-warning">',
gsub('##', '\n', gsub('^##\ Warning:', '**Warning:**', x)),
'</div>', sep = '\n')
},
message = function(x, options) {
paste('\n\n<div class="alert alert-info">',
gsub('##', '\n', x),
'</div>', sep = '\n')
}
)
# Set dpi 20231129
knitr::opts_chunk$set(
dpi = if(nzchar(Sys.getenv("CHNOSZ_BUILD_LARGE_VIGNETTES"))) 100 else 72
)
## ----echo = F, cache = F------------------------------------------------------
# Merge consecutive messages into a single div 20231114
knitr::knit_hooks$set(document = function(x){
# Not sure why this is needed, but simply computing on 'x' doesn't work
file <- tempfile()
writeLines(x, file)
x <- readLines(file)
# Line numbers of the document with </div>
enddiv <- which(x == "</div>")
# Line numbers with <div class="alert alert-info">
beginalert <- which(x == '<div class="alert alert-info">')
# Find </div> followed <div class="alert alert-info"> (skip empty lines)
removediv <- (enddiv + 2) %in% beginalert
if(any(removediv)) {
# Lines to remove
rmlines1 <- enddiv[removediv]
rmlines2 <- enddiv[removediv] + 1
rmlines3 <- enddiv[removediv] + 2
# Do the removing
x <- x[-c(rmlines1, rmlines2, rmlines3)]
}
x
})
## ----HTML, include = FALSE----------------------------------------------------
NOTE <- '<span style="background-color: yellow;">NOTE</span>'
# CHNOSZ functions
equilibrate_ <- '<code>equilibrate()</code>'
info_ <- '<code>info()</code>'
thermo.refs_ <- '<code>thermo.refs()</code>'
# Math stuff
logK <- "log <i>K</i>"
Hplus <- "H<sup>+</sup>"
HCO2_ <- "HCO<sub>2</sub><sup>−</sup>"
HCO3_ <- "HCO<sub>3</sub><sup>−</sup>"
O2 <- "O<sub>2</sub>"
S2 <- "S<sub>2</sub>"
log <- "log "
aHCO2_ <- "<i>a</i><sub>HCO<sub>2</sub><sup>−</sup></sub>"
aHCO3_ <- "<i>a</i><sub>HCO<sub>3</sub><sup>−</sup></sub>"
logfO2 <- "log <i>f</i><sub>O<sub>2</sub></sub>"
Ctot <- "C<sub>tot</sub>"
C3H5O2_ <- "C<sub>3</sub>H<sub>5</sub>O<sub>2</sub><sup>−</sup>"
a3HCO3_ <- "<i>a</i><sup>3</sup><sub>HCO<sub>3</sub><sup>−</sup></sub>"
aC3H5O2_ <- "<i>a</i><sub>C<sub>3</sub>H<sub>5</sub>O<sub>2</sub><sup>−</sup></sub>"
a2HCO3_ <- "<i>a</i><sup>2</sup><sub>HCO<sub>3</sub><sup>−</sup></sub>"
logCtot <- "log C<sub>tot</sub>"
CO2 <- "CO<sub>2</sub>"
H2O <- "H<sub>2</sub>O"
S3minus <- "S<sub>3</sub><sup>-</sup>"
H2S <- "H<sub>2</sub>S"
SO4__ <- "SO<sub>4</sub><sup>-2</sup>"
Kplus <- "K<sup>+</sup>"
Naplus <- "Na<sup>+</sup>"
Clminus <- "Cl<sup>-</sup>"
H2 <- "H<sub>2</sub>"
## ----alanine_refs, message = FALSE--------------------------------------------
info(info("alanine"))[c("ref1", "ref2")]
thermo.refs(info("alanine"))
## ----PPM_refs, message = FALSE------------------------------------------------
basis(c("pyrite", "pyrrhotite", "oxygen"))
sres <- subcrt("magnetite", 1)
info(sres$reaction$ispecies)[, 1:6]
thermo.refs(sres)
reset()
## ----DEW_Ctot, eval = FALSE---------------------------------------------------
# ### Based on demo/DEW.R in CHNOSZ
#
# # Set up subplots
# par(mfrow = c(1, 2), mar = c(3.0, 3.5, 2.5, 1.0), mgp = c(1.7, 0.3, 0), las = 1, tcl = 0.3, xaxs = "i", yaxs = "i")
#
# # Activate DEW model
# water("DEW")
#
# ###########
# ## logfO2-pH diagram for aqueous inorganic and organic carbon species at high pressure
# ## After Figure 1b of Sverjensky et al., 2014b [SSH14]
# ## (Nature Geoscience, https://doi.org/10.1038/NGEO2291)
# ###########
#
# # Define system
# basis("CHNOS+")
# inorganic <- c("CO2", "HCO3-", "CO3-2", "CH4")
# organic <- c("formic acid", "formate", "acetic acid", "acetate", "propanoic acid", "propanoate")
# species(c(inorganic, organic), 0)
# fill <- c(rep("aliceblue", length(inorganic)), rep("aquamarine", length(organic)))
#
# # A function to make the diagrams
# dfun <- function(T = 600, P = 50000, res = 200, Ctot = 0.03) {
# a <- affinity(pH = c(0, 10, res), O2 = c(-24, -12, res), T = T, P = P)
# # Set total C concentration to 0.03 molal
# # (from EQ3NR model for eclogite [Supporting Information of SSH14])
# e <- equilibrate(a, loga.balance = log10(Ctot))
# diagram(e, fill = fill)
# dp <- describe.property(c(" T", " P"), c(T, P), digits = 0)
# legend("bottomleft", legend = dp, bty = "n")
# }
#
# water("DEW")
# add.OBIGT("DEW")
# dfun(Ctot = 0.03)
# mtitle(c("Inorganic and organic species", "C[total] = 0.03 molal"))
# dfun(Ctot = 3)
# mtitle(c("Inorganic and organic species", "C[total] = 3 molal"))
#
# # Restore default settings for the questions below
# reset()
## ----DEW_Ctot, echo = FALSE, message = FALSE, results = "hide", fig.width = 8, fig.height = 4, out.width = "100%", fig.align = "center", pngquant = pngquant, cache = TRUE----
### Based on demo/DEW.R in CHNOSZ
# Set up subplots
par(mfrow = c(1, 2), mar = c(3.0, 3.5, 2.5, 1.0), mgp = c(1.7, 0.3, 0), las = 1, tcl = 0.3, xaxs = "i", yaxs = "i")
# Activate DEW model
water("DEW")
###########
## logfO2-pH diagram for aqueous inorganic and organic carbon species at high pressure
## After Figure 1b of Sverjensky et al., 2014b [SSH14]
## (Nature Geoscience, https://doi.org/10.1038/NGEO2291)
###########
# Define system
basis("CHNOS+")
inorganic <- c("CO2", "HCO3-", "CO3-2", "CH4")
organic <- c("formic acid", "formate", "acetic acid", "acetate", "propanoic acid", "propanoate")
species(c(inorganic, organic), 0)
fill <- c(rep("aliceblue", length(inorganic)), rep("aquamarine", length(organic)))
# A function to make the diagrams
dfun <- function(T = 600, P = 50000, res = 200, Ctot = 0.03) {
a <- affinity(pH = c(0, 10, res), O2 = c(-24, -12, res), T = T, P = P)
# Set total C concentration to 0.03 molal
# (from EQ3NR model for eclogite [Supporting Information of SSH14])
e <- equilibrate(a, loga.balance = log10(Ctot))
diagram(e, fill = fill)
dp <- describe.property(c(" T", " P"), c(T, P), digits = 0)
legend("bottomleft", legend = dp, bty = "n")
}
water("DEW")
add.OBIGT("DEW")
dfun(Ctot = 0.03)
mtitle(c("Inorganic and organic species", "C[total] = 0.03 molal"))
dfun(Ctot = 3)
mtitle(c("Inorganic and organic species", "C[total] = 3 molal"))
# Restore default settings for the questions below
reset()
## ----pyrrhotite_values, message = FALSE---------------------------------------
# The formula of the new mineral and literature reference
formula <- "Fe0.877S"
ref1 <- "PMW87"
# Because the properties from Pankratz et al. (1987) are listed in calories,
# we need to change the output of subcrt() to calories (the default is Joules)
E.units("cal")
# Use temperature in Kelvin for the calculations below
T.units("K")
# Thermodynamic properties of polymorph 1 at 25 °C (298.15 K)
G1 <- -25543
H1 <- -25200
S1 <- 14.531
Cp1 <- 11.922
# Heat capacity coefficients for polymorph 1
a1 <- 7.510
b1 <- 0.014798
# For volume, use the value from Helgeson et al. (1978)
V1 <- V2 <- 18.2
# Transition temperature
Ttr <- 598
# Transition enthalpy (cal/mol)
DHtr <- 95
# Heat capacity coefficients for polymorph 2
a2 <- -1.709
b2 <- 0.011746
c2 <- 3073400
# Maximum temperature of polymorph 2
T2 <- 1800
## ----pyrrhotite_Ttr, message = FALSE------------------------------------------
DGtr <- 0 # DON'T CHANGE THIS
TDStr <- DHtr - DGtr # TΔS° = ΔH° - ΔG°
DStr <- TDStr / Ttr
## ----pyrrhotite_Cp, results = "hide", message = FALSE-------------------------
mod.OBIGT("pyrrhotite_new", formula = formula, state = "cr", ref1 = ref1,
E_units = "cal", G = 0, H = 0, S = 0, V = V1, Cp = Cp1,
a = a1, b = b1, c = 0, d = 0, e = 0, f = 0, lambda = 0, T = Ttr)
mod.OBIGT("pyrrhotite_new", formula = formula, state = "cr2", ref1 = ref1,
E_units = "cal", G = 0, H = 0, S = 0, V = V2,
a = a2, b = b2, c = c2, d = 0, e = 0, f = 0, lambda = 0, T = T2)
## ----pyrrhotite_info, results = "hide", collapse = TRUE-----------------------
info(info("pyrrhotite_new", "cr"))
info(info("pyrrhotite_new", "cr2"))
## ----pyrrhotite_S, message = FALSE--------------------------------------------
DS1 <- subcrt("pyrrhotite_new", "cr", P = 1, T = Ttr, use.polymorphs = FALSE)$out[[1]]$S
DS2 <- subcrt("pyrrhotite_new", "cr2", P = 1, T = Ttr)$out[[1]]$S
DS298 <- DS1 + DStr - DS2
## ----pyrrhotite_D1_D2, message = FALSE, results = "hide"----------------------
mod.OBIGT("pyrrhotite_new", state = "cr", S = S1)
mod.OBIGT("pyrrhotite_new", state = "cr2", S = S1 + DS298)
D1 <- subcrt("pyrrhotite_new", "cr", P = 1, T = Ttr, use.polymorphs = FALSE)$out[[1]]
D2 <- subcrt("pyrrhotite_new", "cr2", P = 1, T = Ttr)$out[[1]]
## ----pyrrhotite_check_S, results = "hide"-------------------------------------
stopifnot(all.equal(D2$S - D1$S, DStr))
## ----pyrrhotite_DG298_DH298, results = "hide", message = FALSE----------------
DG298 <- D1$G + DGtr - D2$G
DH298 <- D1$H + DHtr - D2$H
mod.OBIGT("pyrrhotite_new", state = "cr", G = G1, H = H1)
mod.OBIGT("pyrrhotite_new", state = "cr2", G = G1 + DG298, H = H1 + DH298)
## ----pyrrhotite_info2, collapse = TRUE----------------------------------------
cr_parameters <- info(info("pyrrhotite_new", "cr"))
stopifnot(abs(check.GHS(cr_parameters)) < 1)
cr2_parameters <- info(info("pyrrhotite_new", "cr2"))
stopifnot(abs(check.GHS(cr2_parameters)) < 1)
## ----pyrrhotite_parameters----------------------------------------------------
cr_parameters
cr2_parameters
## ----pyrrhotite_T, eval = FALSE-----------------------------------------------
# T <- seq(550, 650, 1)
# sout <- subcrt("pyrrhotite_new", T = T, P = 1)$out[[1]]
# par(mfrow = c(1, 4), mar = c(4, 4.2, 1, 1))
# labels <- c(G = "DG0f", H = "DH0f", S = "S0", Cp = "Cp0")
# for(property in c("G", "H", "S", "Cp")) {
# plot(sout$T, sout[, property], col = sout$polymorph,
# xlab = axis.label("T"), ylab = axis.label(labels[property])
# )
# abline(v = Ttr, lty = 3, col = 8)
# if(property == "G")
# legend("bottomleft", c("1", "2"), pch = 1, col = c(1, 2), title = "Polymorph")
# }
#
# # Restore default settings for the questions below
# reset()
## ----pyrrhotite_T, echo = FALSE, message = FALSE, results = "hide", fig.width = 8, fig.height = 2.5, out.width = "100%", fig.align = "center", pngquant = pngquant----
T <- seq(550, 650, 1)
sout <- subcrt("pyrrhotite_new", T = T, P = 1)$out[[1]]
par(mfrow = c(1, 4), mar = c(4, 4.2, 1, 1))
labels <- c(G = "DG0f", H = "DH0f", S = "S0", Cp = "Cp0")
for(property in c("G", "H", "S", "Cp")) {
plot(sout$T, sout[, property], col = sout$polymorph,
xlab = axis.label("T"), ylab = axis.label(labels[property])
)
abline(v = Ttr, lty = 3, col = 8)
if(property == "G")
legend("bottomleft", c("1", "2"), pch = 1, col = c(1, 2), title = "Polymorph")
}
# Restore default settings for the questions below
reset()
## ----trisulfur, eval = FALSE, echo = FALSE------------------------------------
# par(mfrow = c(1, 3))
#
# ## BLOCK 1
# T <- 350
# P <- 5000
# res <- 200
#
# ## BLOCK 2
# wt_percent_S <- 10
# wt_permil_S <- wt_percent_S * 10
# molar_mass_S <- mass("S") # 32.06
# moles_S <- wt_permil_S / molar_mass_S
# grams_H2O <- 1000 - wt_permil_S
# molality_S <- 1000 * moles_S / grams_H2O
# logm_S <- log10(molality_S)
#
# ## BLOCK 3
# basis(c("Ni", "SiO2", "Fe2O3", "H2S", "H2O", "oxygen", "H+"))
# species(c("H2S", "HS-", "SO2", "HSO4-", "SO4-2", "S3-"))
# a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P)
# e <- equilibrate(a, loga.balance = logm_S)
# diagram(e)
#
# ## BLOCK 4
# mod.buffer("NNO", c("nickel", "bunsenite"), state = "cr", logact = 0)
# for(buffer in c("HM", "QFM", "NNO")) {
# basis("O2", buffer)
# logfO2_ <- affinity(return.buffer = TRUE, T = T, P = P)$O2
# abline(h = logfO2_, lty = 3, col = 2)
# text(8.5, logfO2_, buffer, adj = c(0, 0), col = 2, cex = 0.9)
# }
#
# ## BLOCK 5
# pH <- subcrt(c("H2O", "H+", "OH-"), c(-1, 1, 1), T = T, P = P)$out$logK / -2
# abline(v = pH, lty = 2, col = 4)
#
# ## BLOCK 6
# lTP <- describe.property(c("T", "P"), c(T, P))
# lS <- paste(wt_percent_S, "wt% S(aq)")
# ltext <- c(lTP, lS)
# legend("topright", ltext, bty = "n", bg = "transparent")
# title(quote("Parameters for"~S[3]^"-"~"from Pokrovski and Dubessy (2015)"), xpd = NA)
#
# ######## Plot 2: Modify Gibbs energy
#
# oldG <- info(info("S3-"))$G
# newG <- oldG - 12548
# mod.OBIGT("S3-", G = newG)
# a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P)
# e <- equilibrate(a, loga.balance = logm_S)
# diagram(e)
# legend("topright", ltext, bty = "n", bg = "transparent")
# title(quote("Modified"~log*italic(K)~"after Pokrovski and Dubrovinsky (2011)"), xpd = NA)
# OBIGT()
#
# ######## Plot 3: Do it with DEW
#
# T <- 700
# P <- 10000 # 10000 bar = 1 GPa
# oldwat <- water("DEW")
# add.OBIGT("DEW")
# info(species()$ispecies)
# a <- affinity(pH = c(2, 10, res), O2 = c(-18, -10, res), T = T, P = P)
# e <- equilibrate(a, loga.balance = logm_S)
# diagram(e)
# lTP <- describe.property(c("T", "P"), c(T, P))
# ltext <- c(lTP, lS)
# legend("topright", ltext, bty = "n", bg = "transparent")
# title(quote("Deep Earth Water (DEW)"~"model"))
# water(oldwat)
# OBIGT()
## ----trisulfur, eval = FALSE, echo = 3:43-------------------------------------
# par(mfrow = c(1, 3))
#
# ## BLOCK 1
# T <- 350
# P <- 5000
# res <- 200
#
# ## BLOCK 2
# wt_percent_S <- 10
# wt_permil_S <- wt_percent_S * 10
# molar_mass_S <- mass("S") # 32.06
# moles_S <- wt_permil_S / molar_mass_S
# grams_H2O <- 1000 - wt_permil_S
# molality_S <- 1000 * moles_S / grams_H2O
# logm_S <- log10(molality_S)
#
# ## BLOCK 3
# basis(c("Ni", "SiO2", "Fe2O3", "H2S", "H2O", "oxygen", "H+"))
# species(c("H2S", "HS-", "SO2", "HSO4-", "SO4-2", "S3-"))
# a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P)
# e <- equilibrate(a, loga.balance = logm_S)
# diagram(e)
#
# ## BLOCK 4
# mod.buffer("NNO", c("nickel", "bunsenite"), state = "cr", logact = 0)
# for(buffer in c("HM", "QFM", "NNO")) {
# basis("O2", buffer)
# logfO2_ <- affinity(return.buffer = TRUE, T = T, P = P)$O2
# abline(h = logfO2_, lty = 3, col = 2)
# text(8.5, logfO2_, buffer, adj = c(0, 0), col = 2, cex = 0.9)
# }
#
# ## BLOCK 5
# pH <- subcrt(c("H2O", "H+", "OH-"), c(-1, 1, 1), T = T, P = P)$out$logK / -2
# abline(v = pH, lty = 2, col = 4)
#
# ## BLOCK 6
# lTP <- describe.property(c("T", "P"), c(T, P))
# lS <- paste(wt_percent_S, "wt% S(aq)")
# ltext <- c(lTP, lS)
# legend("topright", ltext, bty = "n", bg = "transparent")
# title(quote("Parameters for"~S[3]^"-"~"from Pokrovski and Dubessy (2015)"), xpd = NA)
#
# ######## Plot 2: Modify Gibbs energy
#
# oldG <- info(info("S3-"))$G
# newG <- oldG - 12548
# mod.OBIGT("S3-", G = newG)
# a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P)
# e <- equilibrate(a, loga.balance = logm_S)
# diagram(e)
# legend("topright", ltext, bty = "n", bg = "transparent")
# title(quote("Modified"~log*italic(K)~"after Pokrovski and Dubrovinsky (2011)"), xpd = NA)
# OBIGT()
#
# ######## Plot 3: Do it with DEW
#
# T <- 700
# P <- 10000 # 10000 bar = 1 GPa
# oldwat <- water("DEW")
# add.OBIGT("DEW")
# info(species()$ispecies)
# a <- affinity(pH = c(2, 10, res), O2 = c(-18, -10, res), T = T, P = P)
# e <- equilibrate(a, loga.balance = logm_S)
# diagram(e)
# lTP <- describe.property(c("T", "P"), c(T, P))
# ltext <- c(lTP, lS)
# legend("topright", ltext, bty = "n", bg = "transparent")
# title(quote("Deep Earth Water (DEW)"~"model"))
# water(oldwat)
# OBIGT()
## ----trisulfur_logK, message = FALSE, echo = 1:3------------------------------
species <- c("H2S", "SO4-2", "H+", "S3-", "oxygen", "H2O")
coeffs <- c(-2, -1, -1, 1, 0.75, 2.5)
(calclogK <- subcrt(species, coeffs, T = seq(300, 450, 50), P = 5000)$out$logK)
fcalclogK <- formatC(calclogK, format = "f", digits = 1)
reflogK <- -9.6
dlogK <- calclogK[2] - reflogK
# Put in a test so that we don't get surprised by
# future database updates or changes to this vignette
stopifnot(round(dlogK, 4) == -4.4132)
## ----trisulfur, eval = FALSE, echo = 46:55------------------------------------
# par(mfrow = c(1, 3))
#
# ## BLOCK 1
# T <- 350
# P <- 5000
# res <- 200
#
# ## BLOCK 2
# wt_percent_S <- 10
# wt_permil_S <- wt_percent_S * 10
# molar_mass_S <- mass("S") # 32.06
# moles_S <- wt_permil_S / molar_mass_S
# grams_H2O <- 1000 - wt_permil_S
# molality_S <- 1000 * moles_S / grams_H2O
# logm_S <- log10(molality_S)
#
# ## BLOCK 3
# basis(c("Ni", "SiO2", "Fe2O3", "H2S", "H2O", "oxygen", "H+"))
# species(c("H2S", "HS-", "SO2", "HSO4-", "SO4-2", "S3-"))
# a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P)
# e <- equilibrate(a, loga.balance = logm_S)
# diagram(e)
#
# ## BLOCK 4
# mod.buffer("NNO", c("nickel", "bunsenite"), state = "cr", logact = 0)
# for(buffer in c("HM", "QFM", "NNO")) {
# basis("O2", buffer)
# logfO2_ <- affinity(return.buffer = TRUE, T = T, P = P)$O2
# abline(h = logfO2_, lty = 3, col = 2)
# text(8.5, logfO2_, buffer, adj = c(0, 0), col = 2, cex = 0.9)
# }
#
# ## BLOCK 5
# pH <- subcrt(c("H2O", "H+", "OH-"), c(-1, 1, 1), T = T, P = P)$out$logK / -2
# abline(v = pH, lty = 2, col = 4)
#
# ## BLOCK 6
# lTP <- describe.property(c("T", "P"), c(T, P))
# lS <- paste(wt_percent_S, "wt% S(aq)")
# ltext <- c(lTP, lS)
# legend("topright", ltext, bty = "n", bg = "transparent")
# title(quote("Parameters for"~S[3]^"-"~"from Pokrovski and Dubessy (2015)"), xpd = NA)
#
# ######## Plot 2: Modify Gibbs energy
#
# oldG <- info(info("S3-"))$G
# newG <- oldG - 12548
# mod.OBIGT("S3-", G = newG)
# a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P)
# e <- equilibrate(a, loga.balance = logm_S)
# diagram(e)
# legend("topright", ltext, bty = "n", bg = "transparent")
# title(quote("Modified"~log*italic(K)~"after Pokrovski and Dubrovinsky (2011)"), xpd = NA)
# OBIGT()
#
# ######## Plot 3: Do it with DEW
#
# T <- 700
# P <- 10000 # 10000 bar = 1 GPa
# oldwat <- water("DEW")
# add.OBIGT("DEW")
# info(species()$ispecies)
# a <- affinity(pH = c(2, 10, res), O2 = c(-18, -10, res), T = T, P = P)
# e <- equilibrate(a, loga.balance = logm_S)
# diagram(e)
# lTP <- describe.property(c("T", "P"), c(T, P))
# ltext <- c(lTP, lS)
# legend("topright", ltext, bty = "n", bg = "transparent")
# title(quote("Deep Earth Water (DEW)"~"model"))
# water(oldwat)
# OBIGT()
## ----trisulfur_DEW, message = FALSE, echo = 8:22, collapse = TRUE, fig.keep = "none"----
# The first four lines are stand-ins to make this block runnable and get the right output from info();
# the diagram actually shown in the vignette is made using the 'trisulfur' block above
b <- basis("CHNOS+")
s <- species(c("H2S", "HS-", "SO2", "HSO4-", "SO4-2", "S3-"))
res <- 10
wt_percent_S <- logm_S <- 0
######## End of stand-in code ########
T <- 700
P <- 10000 # 10000 bar = 1 GPa
oldwat <- water("DEW")
add.OBIGT("DEW")
info(species()$ispecies)[, 1:8]
a <- affinity(pH = c(2, 10, res), O2 = c(-18, -10, res), T = T, P = P)
e <- equilibrate(a, loga.balance = logm_S)
diagram(e)
lTP <- describe.property(c("T", "P"), c(T, P))
lS <- paste(wt_percent_S, "wt% S(aq)")
ltext <- c(lTP, lS)
legend("topright", ltext, bty = "n", bg = "transparent")
title(quote("Deep Earth Water (DEW)"~"model"))
water(oldwat)
OBIGT()
## ----trisulfur, echo = FALSE, message = FALSE, results = "hide", fig.width = 10, fig.height = 3.33, out.width = "100%", out.extra='class="full-width"', pngquant = pngquant, cache = TRUE----
par(mfrow = c(1, 3))
## BLOCK 1
T <- 350
P <- 5000
res <- 200
## BLOCK 2
wt_percent_S <- 10
wt_permil_S <- wt_percent_S * 10
molar_mass_S <- mass("S") # 32.06
moles_S <- wt_permil_S / molar_mass_S
grams_H2O <- 1000 - wt_permil_S
molality_S <- 1000 * moles_S / grams_H2O
logm_S <- log10(molality_S)
## BLOCK 3
basis(c("Ni", "SiO2", "Fe2O3", "H2S", "H2O", "oxygen", "H+"))
species(c("H2S", "HS-", "SO2", "HSO4-", "SO4-2", "S3-"))
a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P)
e <- equilibrate(a, loga.balance = logm_S)
diagram(e)
## BLOCK 4
mod.buffer("NNO", c("nickel", "bunsenite"), state = "cr", logact = 0)
for(buffer in c("HM", "QFM", "NNO")) {
basis("O2", buffer)
logfO2_ <- affinity(return.buffer = TRUE, T = T, P = P)$O2
abline(h = logfO2_, lty = 3, col = 2)
text(8.5, logfO2_, buffer, adj = c(0, 0), col = 2, cex = 0.9)
}
## BLOCK 5
pH <- subcrt(c("H2O", "H+", "OH-"), c(-1, 1, 1), T = T, P = P)$out$logK / -2
abline(v = pH, lty = 2, col = 4)
## BLOCK 6
lTP <- describe.property(c("T", "P"), c(T, P))
lS <- paste(wt_percent_S, "wt% S(aq)")
ltext <- c(lTP, lS)
legend("topright", ltext, bty = "n", bg = "transparent")
title(quote("Parameters for"~S[3]^"-"~"from Pokrovski and Dubessy (2015)"), xpd = NA)
######## Plot 2: Modify Gibbs energy
oldG <- info(info("S3-"))$G
newG <- oldG - 12548
mod.OBIGT("S3-", G = newG)
a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P)
e <- equilibrate(a, loga.balance = logm_S)
diagram(e)
legend("topright", ltext, bty = "n", bg = "transparent")
title(quote("Modified"~log*italic(K)~"after Pokrovski and Dubrovinsky (2011)"), xpd = NA)
OBIGT()
######## Plot 3: Do it with DEW
T <- 700
P <- 10000 # 10000 bar = 1 GPa
oldwat <- water("DEW")
add.OBIGT("DEW")
info(species()$ispecies)
a <- affinity(pH = c(2, 10, res), O2 = c(-18, -10, res), T = T, P = P)
e <- equilibrate(a, loga.balance = logm_S)
diagram(e)
lTP <- describe.property(c("T", "P"), c(T, P))
ltext <- c(lTP, lS)
legend("topright", ltext, bty = "n", bg = "transparent")
title(quote("Deep Earth Water (DEW)"~"model"))
water(oldwat)
OBIGT()
## ----pyrrhotite_polymorphs, collapse = TRUE-----------------------------------
subcrt("pyrrhotite", T = c(25, 150, 350), property = "G")$out
## ----pyrite_limit-------------------------------------------------------------
subcrt("pyrite", T = seq(200, 1000, 200), P = 1)
## ----mineral_Ttr, collapse = TRUE---------------------------------------------
file <- system.file("extdata/OBIGT/inorganic_cr.csv", package = "CHNOSZ")
dat <- read.csv(file)
# Reverse rows so highest-T polymorph for each mineral is listed first
dat <- dat[nrow(dat):1, ]
# Remove low-T polymorphs
dat <- dat[!duplicated(dat$name), ]
# Remove minerals with no T limit for phase stability (+ve) or Cp equation (-ve)
dat <- dat[!is.na(dat$z.T), ]
# Keep minerals with phase stability limit
dat <- dat[dat$z.T > 0, ]
# Get names of minerals and put into original order
rev(dat$name)
## ----muscovite_limit, message = FALSE-----------------------------------------
add.OBIGT("SUPCRT92")
subcrt("muscovite", T = 850, P = 4500)
reset()
## ----KMQ_basis_species, message = FALSE---------------------------------------
basis(c("Al+3", "quartz", "K+", "H+", "H2O", "oxygen"))
species(c("kaolinite", "muscovite", "K-feldspar"))
## ----KMQ_m_K, message = FALSE-------------------------------------------------
# Define temperature, pressure, and molality of Cl- (==IS)
T <- 150
P <- 500
IS <- m_Cl <- 1
# Calculate equilibrium constant for Ab-Kfs reaction, corrected for ionic strength
logK_AK <- subcrt(c("albite", "K+", "K-feldspar", "Na+"), c(-1, -1, 1, 1),
T = T, P = P, IS = IS)$out$logK
K_AK <- 10 ^ logK_AK
# Calculate molality of K+
(m_K <- m_Cl / (K_AK + 1))
## ----KMQ_diagram, eval = FALSE, echo = 2:10-----------------------------------
# par(mfrow = c(1, 2))
# basis("K+", log10(m_K))
# a <- affinity(pH = c(2, 10), O2 = c(-55, -38), T = T, P = P, IS = IS)
# diagram(a, srt = 90)
# lTP <- as.expression(c(lT(T), lP(P)))
# legend("topleft", legend = lTP, bty = "n", inset = c(-0.05, 0), cex = 0.9)
# ltxt <- c(quote("Unit molality of Cl"^"-"), "Quartz saturation")
# legend("topright", legend = ltxt, bty = "n", cex = 0.9)
# title("Mineral data from Berman (1988)\nand Sverjensky et al. (1991) (OBIGT default)",
# font.main = 1, cex.main = 0.9)
# add.OBIGT("SUPCRT92")
# a <- affinity(pH = c(2, 10), O2 = c(-55, -38), T = T, P = P, IS = IS)
# diagram(a, srt = 90)
# title("Mineral data from Helgeson et al. (1978)\n(as used in SUPCRT92)",
# font.main = 1, cex.main = 0.9)
# OBIGT()
## ----KMQ_diagram, message = FALSE, fig.width = 8, fig.height = 4, out.width = "100%", results = "hide", echo = FALSE----
par(mfrow = c(1, 2))
basis("K+", log10(m_K))
a <- affinity(pH = c(2, 10), O2 = c(-55, -38), T = T, P = P, IS = IS)
diagram(a, srt = 90)
lTP <- as.expression(c(lT(T), lP(P)))
legend("topleft", legend = lTP, bty = "n", inset = c(-0.05, 0), cex = 0.9)
ltxt <- c(quote("Unit molality of Cl"^"-"), "Quartz saturation")
legend("topright", legend = ltxt, bty = "n", cex = 0.9)
title("Mineral data from Berman (1988)\nand Sverjensky et al. (1991) (OBIGT default)",
font.main = 1, cex.main = 0.9)
add.OBIGT("SUPCRT92")
a <- affinity(pH = c(2, 10), O2 = c(-55, -38), T = T, P = P, IS = IS)
diagram(a, srt = 90)
title("Mineral data from Helgeson et al. (1978)\n(as used in SUPCRT92)",
font.main = 1, cex.main = 0.9)
OBIGT()
## ----KMQ_refs, message = FALSE------------------------------------------------
thermo.refs(species()$ispecies)
## ----KMQ_diagram, eval = FALSE, echo = 11:15----------------------------------
# par(mfrow = c(1, 2))
# basis("K+", log10(m_K))
# a <- affinity(pH = c(2, 10), O2 = c(-55, -38), T = T, P = P, IS = IS)
# diagram(a, srt = 90)
# lTP <- as.expression(c(lT(T), lP(P)))
# legend("topleft", legend = lTP, bty = "n", inset = c(-0.05, 0), cex = 0.9)
# ltxt <- c(quote("Unit molality of Cl"^"-"), "Quartz saturation")
# legend("topright", legend = ltxt, bty = "n", cex = 0.9)
# title("Mineral data from Berman (1988)\nand Sverjensky et al. (1991) (OBIGT default)",
# font.main = 1, cex.main = 0.9)
# add.OBIGT("SUPCRT92")
# a <- affinity(pH = c(2, 10), O2 = c(-55, -38), T = T, P = P, IS = IS)
# diagram(a, srt = 90)
# title("Mineral data from Helgeson et al. (1978)\n(as used in SUPCRT92)",
# font.main = 1, cex.main = 0.9)
# OBIGT()
Try the CHNOSZ package in your browser
Any scripts or data that you put into this service are public.
CHNOSZ documentation built on Feb. 12, 2024, 3 p.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.
## ----setup, include = FALSE---------------------------------------------------
library(CHNOSZ)
options(width = 80)
# Use pngquant to optimize PNG images
library(knitr)
knit_hooks$set(pngquant = hook_pngquant)
pngquant <- "--speed=1 --quality=0-25"
if (!nzchar(Sys.which("pngquant"))) pngquant <- NULL
# To make warnings appear in text box 20230619
# https://selbydavid.com/2017/06/18/rmarkdown-alerts/
knitr::knit_hooks$set(
error = function(x, options) {
paste('\n\n<div class="alert alert-danger">',
gsub('##', '\n', gsub('^##\ Error', '**Error:**', x)),
'</div>', sep = '\n')
},
warning = function(x, options) {
paste('\n\n<div class="alert alert-warning">',
gsub('##', '\n', gsub('^##\ Warning:', '**Warning:**', x)),
'</div>', sep = '\n')
},
message = function(x, options) {
paste('\n\n<div class="alert alert-info">',
gsub('##', '\n', x),
'</div>', sep = '\n')
}
)
# Set dpi 20231129
knitr::opts_chunk$set(
dpi = if(nzchar(Sys.getenv("CHNOSZ_BUILD_LARGE_VIGNETTES"))) 100 else 72
)
## ----echo = F, cache = F------------------------------------------------------
# Merge consecutive messages into a single div 20231114
knitr::knit_hooks$set(document = function(x){
# Not sure why this is needed, but simply computing on 'x' doesn't work
file <- tempfile()
writeLines(x, file)
x <- readLines(file)
# Line numbers of the document with </div>
enddiv <- which(x == "</div>")
# Line numbers with <div class="alert alert-info">
beginalert <- which(x == '<div class="alert alert-info">')
# Find </div> followed <div class="alert alert-info"> (skip empty lines)
removediv <- (enddiv + 2) %in% beginalert
if(any(removediv)) {
# Lines to remove
rmlines1 <- enddiv[removediv]
rmlines2 <- enddiv[removediv] + 1
rmlines3 <- enddiv[removediv] + 2
# Do the removing
x <- x[-c(rmlines1, rmlines2, rmlines3)]
}
x
})
## ----HTML, include = FALSE----------------------------------------------------
NOTE <- '<span style="background-color: yellow;">NOTE</span>'
# CHNOSZ functions
equilibrate_ <- '<code>equilibrate()</code>'
info_ <- '<code>info()</code>'
thermo.refs_ <- '<code>thermo.refs()</code>'
# Math stuff
logK <- "log <i>K</i>"
Hplus <- "H<sup>+</sup>"
HCO2_ <- "HCO<sub>2</sub><sup>−</sup>"
HCO3_ <- "HCO<sub>3</sub><sup>−</sup>"
O2 <- "O<sub>2</sub>"
S2 <- "S<sub>2</sub>"
log <- "log "
aHCO2_ <- "<i>a</i><sub>HCO<sub>2</sub><sup>−</sup></sub>"
aHCO3_ <- "<i>a</i><sub>HCO<sub>3</sub><sup>−</sup></sub>"
logfO2 <- "log <i>f</i><sub>O<sub>2</sub></sub>"
Ctot <- "C<sub>tot</sub>"
C3H5O2_ <- "C<sub>3</sub>H<sub>5</sub>O<sub>2</sub><sup>−</sup>"
a3HCO3_ <- "<i>a</i><sup>3</sup><sub>HCO<sub>3</sub><sup>−</sup></sub>"
aC3H5O2_ <- "<i>a</i><sub>C<sub>3</sub>H<sub>5</sub>O<sub>2</sub><sup>−</sup></sub>"
a2HCO3_ <- "<i>a</i><sup>2</sup><sub>HCO<sub>3</sub><sup>−</sup></sub>"
logCtot <- "log C<sub>tot</sub>"
CO2 <- "CO<sub>2</sub>"
H2O <- "H<sub>2</sub>O"
S3minus <- "S<sub>3</sub><sup>-</sup>"
H2S <- "H<sub>2</sub>S"
SO4__ <- "SO<sub>4</sub><sup>-2</sup>"
Kplus <- "K<sup>+</sup>"
Naplus <- "Na<sup>+</sup>"
Clminus <- "Cl<sup>-</sup>"
H2 <- "H<sub>2</sub>"
## ----alanine_refs, message = FALSE--------------------------------------------
info(info("alanine"))[c("ref1", "ref2")]
thermo.refs(info("alanine"))
## ----PPM_refs, message = FALSE------------------------------------------------
basis(c("pyrite", "pyrrhotite", "oxygen"))
sres <- subcrt("magnetite", 1)
info(sres$reaction$ispecies)[, 1:6]
thermo.refs(sres)
reset()
## ----DEW_Ctot, eval = FALSE---------------------------------------------------
# ### Based on demo/DEW.R in CHNOSZ
#
# # Set up subplots
# par(mfrow = c(1, 2), mar = c(3.0, 3.5, 2.5, 1.0), mgp = c(1.7, 0.3, 0), las = 1, tcl = 0.3, xaxs = "i", yaxs = "i")
#
# # Activate DEW model
# water("DEW")
#
# ###########
# ## logfO2-pH diagram for aqueous inorganic and organic carbon species at high pressure
# ## After Figure 1b of Sverjensky et al., 2014b [SSH14]
# ## (Nature Geoscience, https://doi.org/10.1038/NGEO2291)
# ###########
#
# # Define system
# basis("CHNOS+")
# inorganic <- c("CO2", "HCO3-", "CO3-2", "CH4")
# organic <- c("formic acid", "formate", "acetic acid", "acetate", "propanoic acid", "propanoate")
# species(c(inorganic, organic), 0)
# fill <- c(rep("aliceblue", length(inorganic)), rep("aquamarine", length(organic)))
#
# # A function to make the diagrams
# dfun <- function(T = 600, P = 50000, res = 200, Ctot = 0.03) {
# a <- affinity(pH = c(0, 10, res), O2 = c(-24, -12, res), T = T, P = P)
# # Set total C concentration to 0.03 molal
# # (from EQ3NR model for eclogite [Supporting Information of SSH14])
# e <- equilibrate(a, loga.balance = log10(Ctot))
# diagram(e, fill = fill)
# dp <- describe.property(c(" T", " P"), c(T, P), digits = 0)
# legend("bottomleft", legend = dp, bty = "n")
# }
#
# water("DEW")
# add.OBIGT("DEW")
# dfun(Ctot = 0.03)
# mtitle(c("Inorganic and organic species", "C[total] = 0.03 molal"))
# dfun(Ctot = 3)
# mtitle(c("Inorganic and organic species", "C[total] = 3 molal"))
#
# # Restore default settings for the questions below
# reset()
## ----DEW_Ctot, echo = FALSE, message = FALSE, results = "hide", fig.width = 8, fig.height = 4, out.width = "100%", fig.align = "center", pngquant = pngquant, cache = TRUE----
### Based on demo/DEW.R in CHNOSZ
# Set up subplots
par(mfrow = c(1, 2), mar = c(3.0, 3.5, 2.5, 1.0), mgp = c(1.7, 0.3, 0), las = 1, tcl = 0.3, xaxs = "i", yaxs = "i")
# Activate DEW model
water("DEW")
###########
## logfO2-pH diagram for aqueous inorganic and organic carbon species at high pressure
## After Figure 1b of Sverjensky et al., 2014b [SSH14]
## (Nature Geoscience, https://doi.org/10.1038/NGEO2291)
###########
# Define system
basis("CHNOS+")
inorganic <- c("CO2", "HCO3-", "CO3-2", "CH4")
organic <- c("formic acid", "formate", "acetic acid", "acetate", "propanoic acid", "propanoate")
species(c(inorganic, organic), 0)
fill <- c(rep("aliceblue", length(inorganic)), rep("aquamarine", length(organic)))
# A function to make the diagrams
dfun <- function(T = 600, P = 50000, res = 200, Ctot = 0.03) {
a <- affinity(pH = c(0, 10, res), O2 = c(-24, -12, res), T = T, P = P)
# Set total C concentration to 0.03 molal
# (from EQ3NR model for eclogite [Supporting Information of SSH14])
e <- equilibrate(a, loga.balance = log10(Ctot))
diagram(e, fill = fill)
dp <- describe.property(c(" T", " P"), c(T, P), digits = 0)
legend("bottomleft", legend = dp, bty = "n")
}
water("DEW")
add.OBIGT("DEW")
dfun(Ctot = 0.03)
mtitle(c("Inorganic and organic species", "C[total] = 0.03 molal"))
dfun(Ctot = 3)
mtitle(c("Inorganic and organic species", "C[total] = 3 molal"))
# Restore default settings for the questions below
reset()
## ----pyrrhotite_values, message = FALSE---------------------------------------
# The formula of the new mineral and literature reference
formula <- "Fe0.877S"
ref1 <- "PMW87"
# Because the properties from Pankratz et al. (1987) are listed in calories,
# we need to change the output of subcrt() to calories (the default is Joules)
E.units("cal")
# Use temperature in Kelvin for the calculations below
T.units("K")
# Thermodynamic properties of polymorph 1 at 25 °C (298.15 K)
G1 <- -25543
H1 <- -25200
S1 <- 14.531
Cp1 <- 11.922
# Heat capacity coefficients for polymorph 1
a1 <- 7.510
b1 <- 0.014798
# For volume, use the value from Helgeson et al. (1978)
V1 <- V2 <- 18.2
# Transition temperature
Ttr <- 598
# Transition enthalpy (cal/mol)
DHtr <- 95
# Heat capacity coefficients for polymorph 2
a2 <- -1.709
b2 <- 0.011746
c2 <- 3073400
# Maximum temperature of polymorph 2
T2 <- 1800
## ----pyrrhotite_Ttr, message = FALSE------------------------------------------
DGtr <- 0 # DON'T CHANGE THIS
TDStr <- DHtr - DGtr # TΔS° = ΔH° - ΔG°
DStr <- TDStr / Ttr
## ----pyrrhotite_Cp, results = "hide", message = FALSE-------------------------
mod.OBIGT("pyrrhotite_new", formula = formula, state = "cr", ref1 = ref1,
E_units = "cal", G = 0, H = 0, S = 0, V = V1, Cp = Cp1,
a = a1, b = b1, c = 0, d = 0, e = 0, f = 0, lambda = 0, T = Ttr)
mod.OBIGT("pyrrhotite_new", formula = formula, state = "cr2", ref1 = ref1,
E_units = "cal", G = 0, H = 0, S = 0, V = V2,
a = a2, b = b2, c = c2, d = 0, e = 0, f = 0, lambda = 0, T = T2)
## ----pyrrhotite_info, results = "hide", collapse = TRUE-----------------------
info(info("pyrrhotite_new", "cr"))
info(info("pyrrhotite_new", "cr2"))
## ----pyrrhotite_S, message = FALSE--------------------------------------------
DS1 <- subcrt("pyrrhotite_new", "cr", P = 1, T = Ttr, use.polymorphs = FALSE)$out[[1]]$S
DS2 <- subcrt("pyrrhotite_new", "cr2", P = 1, T = Ttr)$out[[1]]$S
DS298 <- DS1 + DStr - DS2
## ----pyrrhotite_D1_D2, message = FALSE, results = "hide"----------------------
mod.OBIGT("pyrrhotite_new", state = "cr", S = S1)
mod.OBIGT("pyrrhotite_new", state = "cr2", S = S1 + DS298)
D1 <- subcrt("pyrrhotite_new", "cr", P = 1, T = Ttr, use.polymorphs = FALSE)$out[[1]]
D2 <- subcrt("pyrrhotite_new", "cr2", P = 1, T = Ttr)$out[[1]]
## ----pyrrhotite_check_S, results = "hide"-------------------------------------
stopifnot(all.equal(D2$S - D1$S, DStr))
## ----pyrrhotite_DG298_DH298, results = "hide", message = FALSE----------------
DG298 <- D1$G + DGtr - D2$G
DH298 <- D1$H + DHtr - D2$H
mod.OBIGT("pyrrhotite_new", state = "cr", G = G1, H = H1)
mod.OBIGT("pyrrhotite_new", state = "cr2", G = G1 + DG298, H = H1 + DH298)
## ----pyrrhotite_info2, collapse = TRUE----------------------------------------
cr_parameters <- info(info("pyrrhotite_new", "cr"))
stopifnot(abs(check.GHS(cr_parameters)) < 1)
cr2_parameters <- info(info("pyrrhotite_new", "cr2"))
stopifnot(abs(check.GHS(cr2_parameters)) < 1)
## ----pyrrhotite_parameters----------------------------------------------------
cr_parameters
cr2_parameters
## ----pyrrhotite_T, eval = FALSE-----------------------------------------------
# T <- seq(550, 650, 1)
# sout <- subcrt("pyrrhotite_new", T = T, P = 1)$out[[1]]
# par(mfrow = c(1, 4), mar = c(4, 4.2, 1, 1))
# labels <- c(G = "DG0f", H = "DH0f", S = "S0", Cp = "Cp0")
# for(property in c("G", "H", "S", "Cp")) {
# plot(sout$T, sout[, property], col = sout$polymorph,
# xlab = axis.label("T"), ylab = axis.label(labels[property])
# )
# abline(v = Ttr, lty = 3, col = 8)
# if(property == "G")
# legend("bottomleft", c("1", "2"), pch = 1, col = c(1, 2), title = "Polymorph")
# }
#
# # Restore default settings for the questions below
# reset()
## ----pyrrhotite_T, echo = FALSE, message = FALSE, results = "hide", fig.width = 8, fig.height = 2.5, out.width = "100%", fig.align = "center", pngquant = pngquant----
T <- seq(550, 650, 1)
sout <- subcrt("pyrrhotite_new", T = T, P = 1)$out[[1]]
par(mfrow = c(1, 4), mar = c(4, 4.2, 1, 1))
labels <- c(G = "DG0f", H = "DH0f", S = "S0", Cp = "Cp0")
for(property in c("G", "H", "S", "Cp")) {
plot(sout$T, sout[, property], col = sout$polymorph,
xlab = axis.label("T"), ylab = axis.label(labels[property])
)
abline(v = Ttr, lty = 3, col = 8)
if(property == "G")
legend("bottomleft", c("1", "2"), pch = 1, col = c(1, 2), title = "Polymorph")
}
# Restore default settings for the questions below
reset()
## ----trisulfur, eval = FALSE, echo = FALSE------------------------------------
# par(mfrow = c(1, 3))
#
# ## BLOCK 1
# T <- 350
# P <- 5000
# res <- 200
#
# ## BLOCK 2
# wt_percent_S <- 10
# wt_permil_S <- wt_percent_S * 10
# molar_mass_S <- mass("S") # 32.06
# moles_S <- wt_permil_S / molar_mass_S
# grams_H2O <- 1000 - wt_permil_S
# molality_S <- 1000 * moles_S / grams_H2O
# logm_S <- log10(molality_S)
#
# ## BLOCK 3
# basis(c("Ni", "SiO2", "Fe2O3", "H2S", "H2O", "oxygen", "H+"))
# species(c("H2S", "HS-", "SO2", "HSO4-", "SO4-2", "S3-"))
# a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P)
# e <- equilibrate(a, loga.balance = logm_S)
# diagram(e)
#
# ## BLOCK 4
# mod.buffer("NNO", c("nickel", "bunsenite"), state = "cr", logact = 0)
# for(buffer in c("HM", "QFM", "NNO")) {
# basis("O2", buffer)
# logfO2_ <- affinity(return.buffer = TRUE, T = T, P = P)$O2
# abline(h = logfO2_, lty = 3, col = 2)
# text(8.5, logfO2_, buffer, adj = c(0, 0), col = 2, cex = 0.9)
# }
#
# ## BLOCK 5
# pH <- subcrt(c("H2O", "H+", "OH-"), c(-1, 1, 1), T = T, P = P)$out$logK / -2
# abline(v = pH, lty = 2, col = 4)
#
# ## BLOCK 6
# lTP <- describe.property(c("T", "P"), c(T, P))
# lS <- paste(wt_percent_S, "wt% S(aq)")
# ltext <- c(lTP, lS)
# legend("topright", ltext, bty = "n", bg = "transparent")
# title(quote("Parameters for"~S[3]^"-"~"from Pokrovski and Dubessy (2015)"), xpd = NA)
#
# ######## Plot 2: Modify Gibbs energy
#
# oldG <- info(info("S3-"))$G
# newG <- oldG - 12548
# mod.OBIGT("S3-", G = newG)
# a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P)
# e <- equilibrate(a, loga.balance = logm_S)
# diagram(e)
# legend("topright", ltext, bty = "n", bg = "transparent")
# title(quote("Modified"~log*italic(K)~"after Pokrovski and Dubrovinsky (2011)"), xpd = NA)
# OBIGT()
#
# ######## Plot 3: Do it with DEW
#
# T <- 700
# P <- 10000 # 10000 bar = 1 GPa
# oldwat <- water("DEW")
# add.OBIGT("DEW")
# info(species()$ispecies)
# a <- affinity(pH = c(2, 10, res), O2 = c(-18, -10, res), T = T, P = P)
# e <- equilibrate(a, loga.balance = logm_S)
# diagram(e)
# lTP <- describe.property(c("T", "P"), c(T, P))
# ltext <- c(lTP, lS)
# legend("topright", ltext, bty = "n", bg = "transparent")
# title(quote("Deep Earth Water (DEW)"~"model"))
# water(oldwat)
# OBIGT()
## ----trisulfur, eval = FALSE, echo = 3:43-------------------------------------
# par(mfrow = c(1, 3))
#
# ## BLOCK 1
# T <- 350
# P <- 5000
# res <- 200
#
# ## BLOCK 2
# wt_percent_S <- 10
# wt_permil_S <- wt_percent_S * 10
# molar_mass_S <- mass("S") # 32.06
# moles_S <- wt_permil_S / molar_mass_S
# grams_H2O <- 1000 - wt_permil_S
# molality_S <- 1000 * moles_S / grams_H2O
# logm_S <- log10(molality_S)
#
# ## BLOCK 3
# basis(c("Ni", "SiO2", "Fe2O3", "H2S", "H2O", "oxygen", "H+"))
# species(c("H2S", "HS-", "SO2", "HSO4-", "SO4-2", "S3-"))
# a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P)
# e <- equilibrate(a, loga.balance = logm_S)
# diagram(e)
#
# ## BLOCK 4
# mod.buffer("NNO", c("nickel", "bunsenite"), state = "cr", logact = 0)
# for(buffer in c("HM", "QFM", "NNO")) {
# basis("O2", buffer)
# logfO2_ <- affinity(return.buffer = TRUE, T = T, P = P)$O2
# abline(h = logfO2_, lty = 3, col = 2)
# text(8.5, logfO2_, buffer, adj = c(0, 0), col = 2, cex = 0.9)
# }
#
# ## BLOCK 5
# pH <- subcrt(c("H2O", "H+", "OH-"), c(-1, 1, 1), T = T, P = P)$out$logK / -2
# abline(v = pH, lty = 2, col = 4)
#
# ## BLOCK 6
# lTP <- describe.property(c("T", "P"), c(T, P))
# lS <- paste(wt_percent_S, "wt% S(aq)")
# ltext <- c(lTP, lS)
# legend("topright", ltext, bty = "n", bg = "transparent")
# title(quote("Parameters for"~S[3]^"-"~"from Pokrovski and Dubessy (2015)"), xpd = NA)
#
# ######## Plot 2: Modify Gibbs energy
#
# oldG <- info(info("S3-"))$G
# newG <- oldG - 12548
# mod.OBIGT("S3-", G = newG)
# a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P)
# e <- equilibrate(a, loga.balance = logm_S)
# diagram(e)
# legend("topright", ltext, bty = "n", bg = "transparent")
# title(quote("Modified"~log*italic(K)~"after Pokrovski and Dubrovinsky (2011)"), xpd = NA)
# OBIGT()
#
# ######## Plot 3: Do it with DEW
#
# T <- 700
# P <- 10000 # 10000 bar = 1 GPa
# oldwat <- water("DEW")
# add.OBIGT("DEW")
# info(species()$ispecies)
# a <- affinity(pH = c(2, 10, res), O2 = c(-18, -10, res), T = T, P = P)
# e <- equilibrate(a, loga.balance = logm_S)
# diagram(e)
# lTP <- describe.property(c("T", "P"), c(T, P))
# ltext <- c(lTP, lS)
# legend("topright", ltext, bty = "n", bg = "transparent")
# title(quote("Deep Earth Water (DEW)"~"model"))
# water(oldwat)
# OBIGT()
## ----trisulfur_logK, message = FALSE, echo = 1:3------------------------------
species <- c("H2S", "SO4-2", "H+", "S3-", "oxygen", "H2O")
coeffs <- c(-2, -1, -1, 1, 0.75, 2.5)
(calclogK <- subcrt(species, coeffs, T = seq(300, 450, 50), P = 5000)$out$logK)
fcalclogK <- formatC(calclogK, format = "f", digits = 1)
reflogK <- -9.6
dlogK <- calclogK[2] - reflogK
# Put in a test so that we don't get surprised by
# future database updates or changes to this vignette
stopifnot(round(dlogK, 4) == -4.4132)
## ----trisulfur, eval = FALSE, echo = 46:55------------------------------------
# par(mfrow = c(1, 3))
#
# ## BLOCK 1
# T <- 350
# P <- 5000
# res <- 200
#
# ## BLOCK 2
# wt_percent_S <- 10
# wt_permil_S <- wt_percent_S * 10
# molar_mass_S <- mass("S") # 32.06
# moles_S <- wt_permil_S / molar_mass_S
# grams_H2O <- 1000 - wt_permil_S
# molality_S <- 1000 * moles_S / grams_H2O
# logm_S <- log10(molality_S)
#
# ## BLOCK 3
# basis(c("Ni", "SiO2", "Fe2O3", "H2S", "H2O", "oxygen", "H+"))
# species(c("H2S", "HS-", "SO2", "HSO4-", "SO4-2", "S3-"))
# a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P)
# e <- equilibrate(a, loga.balance = logm_S)
# diagram(e)
#
# ## BLOCK 4
# mod.buffer("NNO", c("nickel", "bunsenite"), state = "cr", logact = 0)
# for(buffer in c("HM", "QFM", "NNO")) {
# basis("O2", buffer)
# logfO2_ <- affinity(return.buffer = TRUE, T = T, P = P)$O2
# abline(h = logfO2_, lty = 3, col = 2)
# text(8.5, logfO2_, buffer, adj = c(0, 0), col = 2, cex = 0.9)
# }
#
# ## BLOCK 5
# pH <- subcrt(c("H2O", "H+", "OH-"), c(-1, 1, 1), T = T, P = P)$out$logK / -2
# abline(v = pH, lty = 2, col = 4)
#
# ## BLOCK 6
# lTP <- describe.property(c("T", "P"), c(T, P))
# lS <- paste(wt_percent_S, "wt% S(aq)")
# ltext <- c(lTP, lS)
# legend("topright", ltext, bty = "n", bg = "transparent")
# title(quote("Parameters for"~S[3]^"-"~"from Pokrovski and Dubessy (2015)"), xpd = NA)
#
# ######## Plot 2: Modify Gibbs energy
#
# oldG <- info(info("S3-"))$G
# newG <- oldG - 12548
# mod.OBIGT("S3-", G = newG)
# a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P)
# e <- equilibrate(a, loga.balance = logm_S)
# diagram(e)
# legend("topright", ltext, bty = "n", bg = "transparent")
# title(quote("Modified"~log*italic(K)~"after Pokrovski and Dubrovinsky (2011)"), xpd = NA)
# OBIGT()
#
# ######## Plot 3: Do it with DEW
#
# T <- 700
# P <- 10000 # 10000 bar = 1 GPa
# oldwat <- water("DEW")
# add.OBIGT("DEW")
# info(species()$ispecies)
# a <- affinity(pH = c(2, 10, res), O2 = c(-18, -10, res), T = T, P = P)
# e <- equilibrate(a, loga.balance = logm_S)
# diagram(e)
# lTP <- describe.property(c("T", "P"), c(T, P))
# ltext <- c(lTP, lS)
# legend("topright", ltext, bty = "n", bg = "transparent")
# title(quote("Deep Earth Water (DEW)"~"model"))
# water(oldwat)
# OBIGT()
## ----trisulfur_DEW, message = FALSE, echo = 8:22, collapse = TRUE, fig.keep = "none"----
# The first four lines are stand-ins to make this block runnable and get the right output from info();
# the diagram actually shown in the vignette is made using the 'trisulfur' block above
b <- basis("CHNOS+")
s <- species(c("H2S", "HS-", "SO2", "HSO4-", "SO4-2", "S3-"))
res <- 10
wt_percent_S <- logm_S <- 0
######## End of stand-in code ########
T <- 700
P <- 10000 # 10000 bar = 1 GPa
oldwat <- water("DEW")
add.OBIGT("DEW")
info(species()$ispecies)[, 1:8]
a <- affinity(pH = c(2, 10, res), O2 = c(-18, -10, res), T = T, P = P)
e <- equilibrate(a, loga.balance = logm_S)
diagram(e)
lTP <- describe.property(c("T", "P"), c(T, P))
lS <- paste(wt_percent_S, "wt% S(aq)")
ltext <- c(lTP, lS)
legend("topright", ltext, bty = "n", bg = "transparent")
title(quote("Deep Earth Water (DEW)"~"model"))
water(oldwat)
OBIGT()
## ----trisulfur, echo = FALSE, message = FALSE, results = "hide", fig.width = 10, fig.height = 3.33, out.width = "100%", out.extra='class="full-width"', pngquant = pngquant, cache = TRUE----
par(mfrow = c(1, 3))
## BLOCK 1
T <- 350
P <- 5000
res <- 200
## BLOCK 2
wt_percent_S <- 10
wt_permil_S <- wt_percent_S * 10
molar_mass_S <- mass("S") # 32.06
moles_S <- wt_permil_S / molar_mass_S
grams_H2O <- 1000 - wt_permil_S
molality_S <- 1000 * moles_S / grams_H2O
logm_S <- log10(molality_S)
## BLOCK 3
basis(c("Ni", "SiO2", "Fe2O3", "H2S", "H2O", "oxygen", "H+"))
species(c("H2S", "HS-", "SO2", "HSO4-", "SO4-2", "S3-"))
a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P)
e <- equilibrate(a, loga.balance = logm_S)
diagram(e)
## BLOCK 4
mod.buffer("NNO", c("nickel", "bunsenite"), state = "cr", logact = 0)
for(buffer in c("HM", "QFM", "NNO")) {
basis("O2", buffer)
logfO2_ <- affinity(return.buffer = TRUE, T = T, P = P)$O2
abline(h = logfO2_, lty = 3, col = 2)
text(8.5, logfO2_, buffer, adj = c(0, 0), col = 2, cex = 0.9)
}
## BLOCK 5
pH <- subcrt(c("H2O", "H+", "OH-"), c(-1, 1, 1), T = T, P = P)$out$logK / -2
abline(v = pH, lty = 2, col = 4)
## BLOCK 6
lTP <- describe.property(c("T", "P"), c(T, P))
lS <- paste(wt_percent_S, "wt% S(aq)")
ltext <- c(lTP, lS)
legend("topright", ltext, bty = "n", bg = "transparent")
title(quote("Parameters for"~S[3]^"-"~"from Pokrovski and Dubessy (2015)"), xpd = NA)
######## Plot 2: Modify Gibbs energy
oldG <- info(info("S3-"))$G
newG <- oldG - 12548
mod.OBIGT("S3-", G = newG)
a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P)
e <- equilibrate(a, loga.balance = logm_S)
diagram(e)
legend("topright", ltext, bty = "n", bg = "transparent")
title(quote("Modified"~log*italic(K)~"after Pokrovski and Dubrovinsky (2011)"), xpd = NA)
OBIGT()
######## Plot 3: Do it with DEW
T <- 700
P <- 10000 # 10000 bar = 1 GPa
oldwat <- water("DEW")
add.OBIGT("DEW")
info(species()$ispecies)
a <- affinity(pH = c(2, 10, res), O2 = c(-18, -10, res), T = T, P = P)
e <- equilibrate(a, loga.balance = logm_S)
diagram(e)
lTP <- describe.property(c("T", "P"), c(T, P))
ltext <- c(lTP, lS)
legend("topright", ltext, bty = "n", bg = "transparent")
title(quote("Deep Earth Water (DEW)"~"model"))
water(oldwat)
OBIGT()
## ----pyrrhotite_polymorphs, collapse = TRUE-----------------------------------
subcrt("pyrrhotite", T = c(25, 150, 350), property = "G")$out
## ----pyrite_limit-------------------------------------------------------------
subcrt("pyrite", T = seq(200, 1000, 200), P = 1)
## ----mineral_Ttr, collapse = TRUE---------------------------------------------
file <- system.file("extdata/OBIGT/inorganic_cr.csv", package = "CHNOSZ")
dat <- read.csv(file)
# Reverse rows so highest-T polymorph for each mineral is listed first
dat <- dat[nrow(dat):1, ]
# Remove low-T polymorphs
dat <- dat[!duplicated(dat$name), ]
# Remove minerals with no T limit for phase stability (+ve) or Cp equation (-ve)
dat <- dat[!is.na(dat$z.T), ]
# Keep minerals with phase stability limit
dat <- dat[dat$z.T > 0, ]
# Get names of minerals and put into original order
rev(dat$name)
## ----muscovite_limit, message = FALSE-----------------------------------------
add.OBIGT("SUPCRT92")
subcrt("muscovite", T = 850, P = 4500)
reset()
## ----KMQ_basis_species, message = FALSE---------------------------------------
basis(c("Al+3", "quartz", "K+", "H+", "H2O", "oxygen"))
species(c("kaolinite", "muscovite", "K-feldspar"))
## ----KMQ_m_K, message = FALSE-------------------------------------------------
# Define temperature, pressure, and molality of Cl- (==IS)
T <- 150
P <- 500
IS <- m_Cl <- 1
# Calculate equilibrium constant for Ab-Kfs reaction, corrected for ionic strength
logK_AK <- subcrt(c("albite", "K+", "K-feldspar", "Na+"), c(-1, -1, 1, 1),
T = T, P = P, IS = IS)$out$logK
K_AK <- 10 ^ logK_AK
# Calculate molality of K+
(m_K <- m_Cl / (K_AK + 1))
## ----KMQ_diagram, eval = FALSE, echo = 2:10-----------------------------------
# par(mfrow = c(1, 2))
# basis("K+", log10(m_K))
# a <- affinity(pH = c(2, 10), O2 = c(-55, -38), T = T, P = P, IS = IS)
# diagram(a, srt = 90)
# lTP <- as.expression(c(lT(T), lP(P)))
# legend("topleft", legend = lTP, bty = "n", inset = c(-0.05, 0), cex = 0.9)
# ltxt <- c(quote("Unit molality of Cl"^"-"), "Quartz saturation")
# legend("topright", legend = ltxt, bty = "n", cex = 0.9)
# title("Mineral data from Berman (1988)\nand Sverjensky et al. (1991) (OBIGT default)",
# font.main = 1, cex.main = 0.9)
# add.OBIGT("SUPCRT92")
# a <- affinity(pH = c(2, 10), O2 = c(-55, -38), T = T, P = P, IS = IS)
# diagram(a, srt = 90)
# title("Mineral data from Helgeson et al. (1978)\n(as used in SUPCRT92)",
# font.main = 1, cex.main = 0.9)
# OBIGT()
## ----KMQ_diagram, message = FALSE, fig.width = 8, fig.height = 4, out.width = "100%", results = "hide", echo = FALSE----
par(mfrow = c(1, 2))
basis("K+", log10(m_K))
a <- affinity(pH = c(2, 10), O2 = c(-55, -38), T = T, P = P, IS = IS)
diagram(a, srt = 90)
lTP <- as.expression(c(lT(T), lP(P)))
legend("topleft", legend = lTP, bty = "n", inset = c(-0.05, 0), cex = 0.9)
ltxt <- c(quote("Unit molality of Cl"^"-"), "Quartz saturation")
legend("topright", legend = ltxt, bty = "n", cex = 0.9)
title("Mineral data from Berman (1988)\nand Sverjensky et al. (1991) (OBIGT default)",
font.main = 1, cex.main = 0.9)
add.OBIGT("SUPCRT92")
a <- affinity(pH = c(2, 10), O2 = c(-55, -38), T = T, P = P, IS = IS)
diagram(a, srt = 90)
title("Mineral data from Helgeson et al. (1978)\n(as used in SUPCRT92)",
font.main = 1, cex.main = 0.9)
OBIGT()
## ----KMQ_refs, message = FALSE------------------------------------------------
thermo.refs(species()$ispecies)
## ----KMQ_diagram, eval = FALSE, echo = 11:15----------------------------------
# par(mfrow = c(1, 2))
# basis("K+", log10(m_K))
# a <- affinity(pH = c(2, 10), O2 = c(-55, -38), T = T, P = P, IS = IS)
# diagram(a, srt = 90)
# lTP <- as.expression(c(lT(T), lP(P)))
# legend("topleft", legend = lTP, bty = "n", inset = c(-0.05, 0), cex = 0.9)
# ltxt <- c(quote("Unit molality of Cl"^"-"), "Quartz saturation")
# legend("topright", legend = ltxt, bty = "n", cex = 0.9)
# title("Mineral data from Berman (1988)\nand Sverjensky et al. (1991) (OBIGT default)",
# font.main = 1, cex.main = 0.9)
# add.OBIGT("SUPCRT92")
# a <- affinity(pH = c(2, 10), O2 = c(-55, -38), T = T, P = P, IS = IS)
# diagram(a, srt = 90)
# title("Mineral data from Helgeson et al. (1978)\n(as used in SUPCRT92)",
# font.main = 1, cex.main = 0.9)
# OBIGT()
Try the CHNOSZ package in your browser
Any scripts or data that you put into this service are public.
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.