[CHNOSZ-commits] r654 - in pkg/CHNOSZ: . demo
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Fri Mar 19 09:55:45 CET 2021
Author: jedick
Date: 2021-03-19 09:55:45 +0100 (Fri, 19 Mar 2021)
New Revision: 654
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/demo/Pourbaix.R
pkg/CHNOSZ/demo/contour.R
pkg/CHNOSZ/demo/gold.R
pkg/CHNOSZ/demo/solubility.R
pkg/CHNOSZ/demo/sphalerite.R
Log:
Update more demos for new solubility() arguments
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2021-03-19 08:35:15 UTC (rev 653)
+++ pkg/CHNOSZ/DESCRIPTION 2021-03-19 08:55:45 UTC (rev 654)
@@ -1,6 +1,6 @@
Date: 2021-03-19
Package: CHNOSZ
-Version: 1.4.0-23
+Version: 1.4.0-24
Title: Thermodynamic Calculations and Diagrams for Geochemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/demo/Pourbaix.R
===================================================================
--- pkg/CHNOSZ/demo/Pourbaix.R 2021-03-19 08:35:15 UTC (rev 653)
+++ pkg/CHNOSZ/demo/Pourbaix.R 2021-03-19 08:55:45 UTC (rev 654)
@@ -62,13 +62,13 @@
basis(c(elem_basis, "H2O", "H+", "e-"))
# Find species
-i_cr <- retrieve(element, c("O", "H"), "cr")
-i_aq <- retrieve(element, c("O", "H"), "aq")
+icr <- retrieve(element, c("O", "H"), "cr")
+iaq <- retrieve(element, c("O", "H"), "aq")
# Add solids with unit activity
-species(i_cr, 0)
+species(icr, 0)
# Add aqueous species with activity for first equisolubility line
-species(i_aq, levels[1], add = TRUE)
+species(iaq, levels[1], add = TRUE)
# Calculate affinities of formation of species
# from basis species as a function of Eh and pH
@@ -81,9 +81,9 @@
d_all <- diagram(a_all, names = FALSE, lty = 0, min.area = 0.1, limit.water = limit.water, fill = fill)
# Calculate affinities for minerals
-species(i_cr)
+species(icr)
# Calculate overall solubility (i.e. minimum solubility given all candidate minerals)
-s <- solubility(i_aq, pH = c(pH, res), Eh = c(Eh, res), T = T, P = P, IS = IS, in.terms.of = element)
+s <- solubility(iaq, pH = c(pH, res), Eh = c(Eh, res), T = T, P = P, IS = IS, in.terms.of = element)
# Plot diagram (LAYER 2: equisolubility lines)
diagram(s, levels = levels, contour.method = "flattest", add = TRUE, lwd = 1.5)
@@ -90,7 +90,7 @@
# Calculate affinities for aqueous species
# FIXME: should be able to remove cr species from previous affinity object
-species(i_aq, 0)
+species(iaq, 0)
a_aq <- affinity(pH = c(pH, res), Eh = c(Eh, res), T = T, P = P, IS = IS)
# Plot diagram (LAYER 3: equal-activity lines for aqueous species)
@@ -115,10 +115,10 @@
diagram(a_aq, add = TRUE, lty = 0, col = col, dy = dy, col.names = col.names)
# Add solids with unit activity
-species(i_cr, 0)
+species(icr, 0)
# Add aqueous species with activity for last equisolubility line
# (i.e. unit activity)
-species(i_aq, levels[length(levels)], add = TRUE)
+species(iaq, levels[length(levels)], add = TRUE)
# Calculate affinities of formation of species
# from basis species as a function of Eh and pH
Modified: pkg/CHNOSZ/demo/contour.R
===================================================================
--- pkg/CHNOSZ/demo/contour.R 2021-03-19 08:35:15 UTC (rev 653)
+++ pkg/CHNOSZ/demo/contour.R 2021-03-19 08:55:45 UTC (rev 654)
@@ -1,5 +1,5 @@
# CHNOSZ/demo/contour.R
-# gold solubility contours on logfO2-pH diagram
+# Gold solubility contours on logfO2-pH diagram
# After Williams-Jones et al., 2009, Fig. 3
# doi:10.2113/gselements.5.5.281
# 20181107 initial version
@@ -6,33 +6,36 @@
# 20190415 cleanup for demo
library(CHNOSZ)
-# define temperature (degrees C), pressure (bar), grid resolution
+# Define temperature (degrees C), pressure (bar), grid resolution
T <- 250
P <- 500
res <- 600
-# make smooth (TRUE) or sharp (FALSE) transitions between basis species
+# Make smooth (TRUE) or sharp (FALSE) transitions between basis species
blend <- TRUE
-# set up system
+# Set up system
basis(c("Au", "Cl-", "H2S", "H2O", "oxygen", "H+"))
-species(c("Au(HS)2-", "AuHS", "AuOH", "AuCl2-"))
-# this get us close to total S = 0.01 m
+iaq <- info(c("Au(HS)2-", "AuHS", "AuOH", "AuCl2-"))
+species(iaq)
+# This get us close to total S = 0.01 m
basis("H2S", -2)
-# calculate solution composition for 1 mol/kg NaCl
+# Calculate solution composition for 1 mol/kg NaCl
NaCl <- NaCl(T = T, P = P, m_tot=1)
basis("Cl-", log10(NaCl$m_Cl))
-# calculate affinity with changing basis species
+# Calculate affinity with changing basis species
bases <- c("H2S", "HS-", "HSO4-", "SO4-2")
m <- mosaic(bases, pH = c(2, 10, res), O2 = c(-41, -29, res), T = T, P = P, IS = NaCl$IS, blend = blend)
-# show predominance fields
+# Show predominance fields
diagram(m$A.bases, col = "red", col.names = "red", lty = 2, italic = TRUE)
diagram(m$A.species, add=TRUE, col = "blue", col.names = "blue", lwd = 2, bold = TRUE)
-# calculate and plot solubility
-s <- solubility(m$A.species)
-# convert to ppb
+
+# Calculate and plot solubility of Au (use named 'bases' argument to trigger mosaic calculation)
+species("Au")
+s <- solubility(iaq, bases = bases, pH = c(2, 10, res), O2 = c(-41, -29, res), T = T, P = P, IS = NaCl$IS, blend = blend)
+# Convert to ppb
s <- convert(s, "ppb")
diagram(s, type = "loga.balance", levels = c(1, 10, 100, 1000), add = TRUE)
-# add legend and title
+# Add legend and title
dP <- describe.property(c("T", "P"), c(250, 500))
legend("top", dP, bty = "n", inset = c(0, 0.06))
lx <- lex(lNaCl(1), lS(0.01))
Modified: pkg/CHNOSZ/demo/gold.R
===================================================================
--- pkg/CHNOSZ/demo/gold.R 2021-03-19 08:35:15 UTC (rev 653)
+++ pkg/CHNOSZ/demo/gold.R 2021-03-19 08:55:45 UTC (rev 654)
@@ -56,7 +56,7 @@
basis("H2S", "PPM")
# Calculate solubility of gold
species("Au")
- iaq <- c("Au(HS)2-", "AuHS", "AuOH")
+ iaq <- info(c("Au(HS)2-", "AuHS", "AuOH"))
# (set IS = 0 for diagram to show "log m" instead of "log a")
s <- solubility(iaq, pH = c(3, 8), T = 300, P = 1000, IS = 0)
# Make diagram and show total log molality
@@ -76,14 +76,14 @@
# log(m_Au)-pH diagram similar to Fig. 12b of Stefansson and Seward, 2004
# (doi:10.1016/j.gca.2004.04.006)
Au_pH2 <- function() {
- species(c("Au(HS)2-", "AuHS", "AuOH", "AuCl2-"))
# Apply PPM buffer for fO2 and aH2S
basis("O2", "PPM")
basis("H2S", "PPM")
- # Calculate affinity and solubility
+ # Calculate solubility of gold
# (set IS = 0 for diagram to show "log m" instead of "log a")
- a <- affinity(pH = c(3, 8), T = 450, P = 1000, IS = 0)
- s <- solubility(a)
+ species("Au")
+ iaq <- info(c("Au(HS)2-", "AuHS", "AuOH", "AuCl2-"))
+ s <- solubility(iaq, pH = c(3, 8), T = 450, P = 1000, IS = 0)
# Make diagram and show total log molality
diagram(s, ylim = c(-8, -3), col = col, lwd = 2, lty = 1)
diagram(s, add = TRUE, type = "loga.balance", lwd = 3, lty = 2)
@@ -113,7 +113,6 @@
# log(m_Au)-T diagram like Fig. 2B of Williams-Jones et al., 2009
# (doi:10.2113/gselements.5.5.281)
Au_T1 <- function() {
- species(c("Au(HS)2-", "AuHS", "AuOH", "AuCl2-"))
# Apply PPM buffer for fO2 and aH2S
basis("O2", "PPM")
basis("H2S", "PPM")
@@ -121,9 +120,10 @@
basis("H+", "QMK")
# Estimate solution composition for 1.5 m NaCl and 0.5 m KCl
chl <- chloride(T = seq(150, 550, 10), P = 1000, m_NaCl = 1.5, m_KCl = 0.5)
- # Calculate affinity and solubility
- a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(chl$m_Cl), `K+` = log10(chl$m_K), P = 1000, IS = chl$IS)
- s <- solubility(a)
+ # Calculate solubility of gold
+ species("Au")
+ iaq <- info(c("Au(HS)2-", "AuHS", "AuOH", "AuCl2-"))
+ s <- solubility(iaq, T = seq(150, 550, 10), `Cl-` = log10(chl$m_Cl), `K+` = log10(chl$m_K), P = 1000, IS = chl$IS)
# Make diagram and show total log molality
diagram(s, ylim = c(-10, -3), col = col, lwd = 2, lty = 1)
diagram(s, add = TRUE, type = "loga.balance", lwd = 3, lty = 2)
@@ -143,7 +143,6 @@
# (doi:10.2113/gselements.5.5.281)
# (doi:10.1144/SP402.4)
Au_T2 <- function() {
- species(c("Au(HS)2-", "AuHS", "AuOH", "AuCl2-"))
# Total S = 0.01 m
basis("H2S", -2)
# Apply HM buffer for fO2
@@ -156,8 +155,9 @@
# bases <- c("H2S", "HS-", "SO4-2", "HSO4-")
# m <- mosaic(bases, T = seq(150, 550, 10), `Cl-` = log10(chl$m_Cl), `K+` = log10(chl$m_K), P = 1000, IS = chl$IS)
# s <- solubility(m$A.species)
- a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(chl$m_Cl), `K+` = log10(chl$m_K), P = 1000, IS = chl$IS)
- s <- solubility(a)
+ species("Au")
+ iaq <- info(c("Au(HS)2-", "AuHS", "AuOH", "AuCl2-"))
+ s <- solubility(iaq, T = seq(150, 550, 10), `Cl-` = log10(chl$m_Cl), `K+` = log10(chl$m_K), P = 1000, IS = chl$IS)
# Make diagram and show total log molality
diagram(s, ylim = c(-10, -3), col = col, lwd = 2, lty = 1)
diagram(s, add = TRUE, type = "loga.balance", lwd = 3, lty = 2)
Modified: pkg/CHNOSZ/demo/solubility.R
===================================================================
--- pkg/CHNOSZ/demo/solubility.R 2021-03-19 08:35:15 UTC (rev 653)
+++ pkg/CHNOSZ/demo/solubility.R 2021-03-19 08:55:45 UTC (rev 654)
@@ -22,12 +22,11 @@
IS <- 0
# Start with CO2
-basis(c("carbon dioxide", "H2O", "O2", "H+"))
+basis(c("CO2", "H2O", "O2", "H+"))
# This is ca. atmospheric PCO2
-basis("CO2", -3.5)
-species(c("CO2", "HCO3-", "CO3-2"))
-a <- affinity(pH = c(pH, res), T = T1, IS = IS)
-s <- solubility(a)
+species("carbon dioxide", -3.5)
+iaq <- info(c("CO2", "HCO3-", "CO3-2"))
+s <- solubility(iaq, pH = c(pH, res), T = T1, IS = IS)
# First plot total activity line
diagram(s, ylim = c(-10, 4), type = "loga.balance", lwd = 4, col = "green2")
# Add activities of species
@@ -44,17 +43,16 @@
stopifnot(round(s$loga.balance[c(1, res)])==c(-5, 6))
# CO2 T-pH plot
-a <- affinity(pH = c(pH, res), T = c(T, res), IS = IS)
-s <- solubility(a)
+s <- solubility(iaq, pH = c(pH, res), T = c(T, res), IS = IS)
diagram(s, type = "loga.balance")
title(main = substitute("Solubility of"~what, list(what = expr.species("CO2"))))
# Now do calcite
-basis(c("calcite", "Ca+2", "H2O", "O2", "H+"))
-species(c("CO2", "HCO3-", "CO3-2"))
-a <- affinity(pH = c(pH, res), T = T1, IS = IS)
+basis(c("CO2", "Ca+2", "H2O", "O2", "H+"))
+species("calcite")
+iaq <- info(c("CO2", "HCO3-", "CO3-2"))
# Change this to dissociate = 2 to reproduce straight lines in Fig. 4A of Manning et al., 2013
-s <- solubility(a, dissociate = TRUE)
+s <- solubility(iaq, pH = c(pH, res), T = T1, IS = IS, dissociate = TRUE)
diagram(s, ylim = c(-10, 4), type = "loga.balance", lwd = 4, col = "green2")
diagram(s, add = TRUE, dy = 1)
legend("topright", lty = c(1, 1:3), lwd = c(4, 2, 2, 2),
@@ -66,8 +64,7 @@
stopifnot(round(s$loga.balance[c(1, res)])==c(4, -4))
# Calcite T-pH plot
-a <- affinity(pH = c(pH, res), T = c(T, res), IS = IS)
-s <- solubility(a)
+s <- solubility(iaq, pH = c(pH, res), T = c(T, res), IS = IS, dissociate = TRUE)
diagram(s, type = "loga.balance")
title(main = "Solubility of calcite", font.main = 1)
Modified: pkg/CHNOSZ/demo/sphalerite.R
===================================================================
--- pkg/CHNOSZ/demo/sphalerite.R 2021-03-19 08:35:15 UTC (rev 653)
+++ pkg/CHNOSZ/demo/sphalerite.R 2021-03-19 08:55:45 UTC (rev 654)
@@ -1,35 +1,34 @@
# CHNOSZ/demo/sphalerite.R
-# sphalerite solubility after Akinfiev and Tagirov, 2014, Fig. 13
+# Sphalerite solubility after Akinfiev and Tagirov, 2014, Fig. 13
# 20190526 jmd initial version
library(CHNOSZ)
opar <- par(no.readonly = TRUE)
-# set up chemical system
-basis(c("ZnS", "Cl-", "H2S", "H2O", "O2", "H+"))
-iZn <- retrieve("Zn", c("O", "H", "Cl", "S"), "aq")
-species(iZn)
+# Set up chemical system
+basis(c("Zn+2", "Cl-", "H2S", "H2O", "O2", "H+"))
+species("ZnS")
+iaq <- retrieve("Zn", c("O", "H", "Cl", "S"), "aq")
# a function to make a single plot
plotfun <- function(T = 400, P = 500, m_tot = 0.1, pHmin = 4, logppmmax = 3) {
- # calculate NaCl speciation from simplified model
+ # Calculate NaCl speciation from simplified model
NaCl <- NaCl(T = T, P = P, m_tot = m_tot)
basis("Cl-", log10(NaCl$m_Cl))
basis("H2S", log10(0.05))
- # use mosaic to account for HS- and H2S speciation
- m <- mosaic(c("H2S", "HS-"), pH = c(pHmin, 10), T = T, P = P, IS = NaCl$IS)
- s <- solubility(m$A.species)
+ # Calculate solubility with mosaic (triggered by bases argument) to account for HS- and H2S speciation
+ s <- solubility(iaq, bases = c("H2S", "HS-"), pH = c(pHmin, 10), T = T, P = P, IS = NaCl$IS)
- # convert log activity to log ppm
+ # Convert log activity to log ppm
sp <- convert(s, "logppm")
diagram(sp, ylim = c(-5, logppmmax))
diagram(sp, type = "loga.balance", add = TRUE, lwd = 2, col = "green3")
- # add water neutrality line
+ # Add water neutrality line
pKw <- - subcrt(c("H2O", "OH-", "H+"), c(-1, 1, 1), T = T, P = P)$out$logK
abline(v = pKw / 2, lty = 2, lwd = 2, col = "blue1")
- # add legend
+ # Add legend
l <- lex(lNaCl(m_tot), lTP(T, P))
legend("topright", legend = l, bty = "n")
}
@@ -39,9 +38,9 @@
par(opar)
-### the following code for making multiple plots is not used in the demo ###
+### The following code for making multiple plots is not used in the demo ###
-# a function to make a page of plots
+# A function to make a page of plots
pagefun <- function() {
# set the values of temperature, pressure, and total NaCl
T <- c(400, 400, 250, 250, 100, 100)
@@ -56,7 +55,7 @@
for(i in 1:6) plotfun(T = T[i], P = P[[i]], m_tot = m_tot[i], pHmin = pHmin[i], logppmmax = logppmmax[i])
}
-# a function to make a png file with all the plots
+# A function to make a png file with all the plots
pngfun <- function() {
png("sphalerite.png", width = 1000, height = 1200, pointsize = 24)
pagefun()
@@ -67,6 +66,6 @@
dev.off()
}
-# we don't run these functions in the demo
+# We don't run these functions in the demo
#pagefun()
#pngfun()
More information about the CHNOSZ-commits
mailing list