-
Notifications
You must be signed in to change notification settings - Fork 2
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Add optional data and update conditions in demo/contour.R
git-svn-id: svn://scm.r-forge.r-project.org/svnroot/chnosz/pkg/CHNOSZ@873 edb9625f-4e0d-4859-8d74-9fd3b1da38cb
- Loading branch information
Showing
4 changed files
with
64 additions
and
40 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,6 +1,6 @@ | ||
Date: 2025-02-15 | ||
Package: CHNOSZ | ||
Version: 2.1.0-44 | ||
Version: 2.1.0-45 | ||
Title: Thermodynamic Calculations and Diagrams for Geochemistry | ||
Authors@R: c( | ||
person("Jeffrey", "Dick", , "[email protected]", role = c("aut", "cre"), | ||
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -1,47 +1,64 @@ | ||
# CHNOSZ/demo/contour.R | ||
# Gold solubility contours on logfO2-pH diagram | ||
# After Williams-Jones et al., 2009, Fig. 3 | ||
# doi:10.2113/gselements.5.5.281 | ||
# After Williams-Jones et al. (2009) doi:10.2113/gselements.5.5.281 | ||
# and Ding et al. (2023) doi:10.1016/j.oregeorev.2022.105231 | ||
# with optional data from Tagirov et al. (2024) doi:10.1016/j.gca.2024.08.022 | ||
# 20181107 initial version | ||
# 20190415 cleanup for demo | ||
# 20250215 add IGEM data and adjust conditions to reproduce plot from Ding et al. | ||
|
||
library(CHNOSZ) | ||
|
||
# Define temperature (degrees C), pressure (bar), grid resolution | ||
T <- 250 | ||
P <- 500 | ||
# Define plot resolution | ||
res <- 600 | ||
# Define temperature (degrees C), pressure (bar) | ||
T <- 250 | ||
P <- 300 | ||
# Define molality of NaCl and total S | ||
m_NaCl <- 1 | ||
sum_S <- 0.01 | ||
# Define ranges of pH and logfO2 | ||
pH <- c(1, 10) | ||
O2 <- c(-40, -30) | ||
# Make smooth (TRUE) or sharp (FALSE) transitions between basis species | ||
blend <- TRUE | ||
|
||
# Set up system | ||
basis(c("Au", "Cl-", "H2S", "H2O", "oxygen", "H+")) | ||
iaq <- info(c("Au(HS)2-", "AuHS", "AuOH", "AuCl2-")) | ||
species(iaq) | ||
# This gets us close to total S = 0.01 m | ||
basis("H2S", -2) | ||
# Calculate solution composition for 1 mol/kg NaCl | ||
NaCl <- NaCl(m_NaCl = 1, T = T, P = P) | ||
basis("Cl-", log10(NaCl$m_Clminus)) | ||
# 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 for S-species | ||
diagram(m$A.bases, col = 8, col.names = 8, lty = 2, italic = TRUE) | ||
# Show predominance fields for Au-species | ||
# Plot multiple times to get deeper color | ||
diagram(m$A.species, add=TRUE, col = 7, col.names = 7, lwd = 2, bold = TRUE) | ||
diagram(m$A.species, add=TRUE, col = 7, col.names = 7, lwd = 2, bold = TRUE) | ||
diagram(m$A.species, add=TRUE, col = 7, col.names = 7, lwd = 2, bold = TRUE) | ||
for(i in 1:2) { | ||
|
||
# Use the default database for the first diagram | ||
if(i == 1) reset() | ||
# Use IGEM data for the second diagram | ||
if(i == 2) add.OBIGT("IGEM") | ||
|
||
# Set up system | ||
basis(c("Au", "Cl-", "H2S", "H2O", "oxygen", "H+")) | ||
iaq <- retrieve("Au", c("H", "O", "S", "Cl"), state = "aq") | ||
species(iaq) | ||
basis("H2S", log10(sum_S)) | ||
# Calculate solution composition for given NaCl molality | ||
NaCl <- NaCl(m_NaCl = m_NaCl, T = T, P = P) | ||
basis("Cl-", log10(NaCl$m_Clminus)) | ||
# Calculate affinity with changing basis species across diagram | ||
bases <- c("H2S", "HS-", "HSO4-", "SO4-2") | ||
m <- mosaic(bases, pH = c(pH, res), O2 = c(O2, res), T = T, P = P, IS = NaCl$IS, blend = blend) | ||
# Show predominance fields for S-species | ||
diagram(m$A.bases, col = 8, col.names = 8, lty = 3, italic = TRUE) | ||
# Show predominance fields for Au-species | ||
d <- diagram(m$A.species, col = 4, col.names = 4, lty = 2, bold = TRUE, add = TRUE) | ||
# Add dot-dash line for water stability limit | ||
water.lines(d, lty = 4) | ||
|
||
# Calculate and plot solubility of Au (use named 'bases' argument to trigger mosaic calculation) | ||
species("Au") | ||
s <- solubility(iaq, bases = bases, pH = c(pH, res), O2 = c(O2, res), T = T, P = P, IS = NaCl$IS, blend = blend) | ||
# Convert to ppb | ||
s <- convert(s, "ppb") | ||
diagram(s, levels = c(1, 10, 100, 1000), col = 1, add = TRUE) | ||
# Add legend and title | ||
dP <- describe.property(c("T", "P"), c(T, P)) | ||
lexpr <- lex(dP, lNaCl(m_NaCl), lS(sum_S)) | ||
legend("topright", lexpr, bty = "n") | ||
if(i == 1) title(main = ("Gold solubility (ppb), after Ding et al., 2023 (default data)"), font.main = 1) | ||
if(i == 2) title(main = ("Gold solubility (ppb), after Ding et al., 2023 (optional data)"), font.main = 1) | ||
|
||
# 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, levels = c(1, 10, 100, 1000), add = TRUE) | ||
# 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)) | ||
legend("topright", lx, bty = "n", inset = c(0.1, 0.05)) | ||
title(main = ("Solubility of gold (ppb), after Williams-Jones et al., 2009, Fig. 3"), font.main = 1) | ||
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters