[CHNOSZ-commits] r761 - in pkg/CHNOSZ: . R demo inst man vignettes
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Dec 11 05:21:41 CET 2022
Author: jedick
Date: 2022-12-11 05:21:41 +0100 (Sun, 11 Dec 2022)
New Revision: 761
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/R/NaCl.R
pkg/CHNOSZ/demo/contour.R
pkg/CHNOSZ/demo/gold.R
pkg/CHNOSZ/demo/minsol.R
pkg/CHNOSZ/demo/sphalerite.R
pkg/CHNOSZ/inst/NEWS.Rd
pkg/CHNOSZ/man/NaCl.Rd
pkg/CHNOSZ/vignettes/customizing.Rmd
pkg/CHNOSZ/vignettes/multi-metal.Rmd
Log:
Rewrite NaCl() to include pH dependence
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2022-12-09 06:17:33 UTC (rev 760)
+++ pkg/CHNOSZ/DESCRIPTION 2022-12-11 04:21:41 UTC (rev 761)
@@ -1,6 +1,6 @@
-Date: 2022-12-09
+Date: 2022-12-11
Package: CHNOSZ
-Version: 1.9.9-52
+Version: 1.9.9-53
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/R/NaCl.R
===================================================================
--- pkg/CHNOSZ/R/NaCl.R 2022-12-09 06:17:33 UTC (rev 760)
+++ pkg/CHNOSZ/R/NaCl.R 2022-12-11 04:21:41 UTC (rev 761)
@@ -1,67 +1,77 @@
-# NaCl.R
-# Calculate ion activities and ionic strength
-# for a total molality of NaCl,
-# taking account of ion association: Na+ + Cl- = NaCl(aq)
-# 20181102 jmd first version
-# 20181105 use activity coefficient of Na+
-# 20181106 use activity coefficient of NaCl
-# 20221122 make it work with m_tot = 0
+# CHNOSZ/NaCl.R
+# Calculate ionic strength and molalities of species given a total molality of NaCl
+# 20181102 First version (jmd)
+# 20181106 Use activity coefficients of Na+ and Nacl
+# 20221122 Make it work with m_tot = 0
+# 20221210 Rewritten to include pH dependence;
+# uses affinity() and equilibrate() instead of algebraic equations
-NaCl <- function(T = seq(100, 500, 100), P = 1000, m_tot = 2, ...) {
- # Define a function for the reaction quotient
- logQ <- function(m_Cl, gam_NaCl, gam_Na, gam_Cl) {
- # Starting with Q = a_NaCl / (a_Na+ * a_Cl-),
- # substitute m_tot = m_NaCl + m_Cl and m_Cl = m_Na
- # to write:
- log10( (m_tot - m_Cl) * gam_NaCl / (m_Cl ^ 2 * gam_Na * gam_Cl) )
- }
- # Define a function for affinity = log(K / Q)
- A <- function(m_Cl, gam_NaCl, gam_Na, gam_Cl, logK) logK - logQ(m_Cl, gam_NaCl, gam_Na, gam_Cl)
- # Calculate equilibrium constant at all temperatures (standard conditions: IS = 0)
- logK <- subcrt(c("Na+", "Cl-", "NaCl"), c(-1, -1, 1), T = T, P = P, ...)$out$logK
- # Calculate Debye-Huckel parameters at all temperatures
- wout <- water(c("A_DH", "B_DH"), T = convert(T, "K"), P = P)
- # Initialize output variables
- N <- length(T)
- ISout <- a_Cl <- numeric(N)
- # Initial guess for m_Cl and ionic strength assuming complete dissociation of NaCl
- IS <- m_Cl <- rep(m_tot, N)
- # The corresponding total molality of dissolved species (NaCl + Cl- + Na+)
- m_star <- (m_tot - m_Cl) + 2*m_Cl
- # The species indices for Na+, Cl-, and NaCl(aq)
- ispecies <- info(c("Na+", "Cl-", "NaCl"))
- # Data frame needed for nonideal()
- speciesprops <- rep(list(data.frame(G = numeric(N))), length(ispecies))
- # Only try to find root for non-zero NaCl concentration 20221122
- if(m_tot > 0) {
- # Start by doing calculations for all temperatures
- doit <- !logical(N)
- while(any(doit)) {
- # Calculate activity coefficient at given ionic strength
- gammas <- suppressMessages(nonideal(ispecies, speciesprops, IS = IS, T = convert(T, "K"), P = P, A_DH = wout$A_DH, B_DH = wout$B_DH, m_star = m_star))
- gam_Na <- 10^gammas[[1]]$loggam
- gam_Cl <- 10^gammas[[2]]$loggam
- gam_NaCl <- 10^gammas[[3]]$loggam
- # In case Setchenow calculations are turned off 20191209
- if(length(gam_NaCl) == 0) gam_NaCl <- rep(1, length(T))
- # Solve for m_Cl
- for(i in which(doit)) m_Cl[i] <- uniroot(A, c(0, m_tot), gam_NaCl = gam_NaCl[i], gam_Na = gam_Na[i], gam_Cl = gam_Cl[i], logK = logK[i])$root
- # Calculate new total molality
- m_star <- (m_tot - m_Cl) + 2*m_Cl
- # Calculate new ionic strength and deviation
- ISnew <- m_Cl
- dIS <- ISnew - IS
- # Set new ionic strength
- IS <- ISnew
- # Keep going until the deviation in ionic strength at any temperature is less than 0.01
- doit <- abs(dIS) > 0.01
+NaCl <- function(m_tot = 1, T = 25, P = "Psat", pH = NA, attenuate = FALSE) {
+
+ # Store existing thermo data frame
+ thermo <- get("thermo", CHNOSZ)
+ # Get length of longest variable
+ nTP <- max(length(T), length(P))
+ pH.arg <- pH
+ pH <- rep(pH, length.out = nTP)
+
+ # Start with complete dissociation into Na+ and Cl-,
+ # so ionic strength and molality of Na+ are equal to m_tot
+ m_Na <- IS <- m_tot
+ # Make them same length as T and P
+ IS <- rep(IS, length.out = nTP)
+ m_Na <- rep(m_Na, length.out = nTP)
+ # Set tolerance for convergence to 1/100th of m_tot
+ tolerance <- m_tot / 100
+
+ # If m_tot is 0, return 0 for all variables 20221122
+ zeros <- rep(0, nTP)
+ if(m_tot == 0) return(list(IS = zeros, m_Na = zeros, m_Cl = zeros, m_NaCl = zeros, m_HCl = zeros))
+
+ maxiter <- 100
+ for(i in 1:maxiter) {
+ # Setup chemical system and calculate affinities
+ if(identical(pH.arg, NA)) {
+ basis(c("Na+", "Cl-", "e-"))
+ species(c("Cl-", "NaCl"))
+ a <- suppressMessages(affinity(T = T, P = P, "Na+" = log10(m_Na), IS = IS, transect = TRUE))
+ } else {
+ basis(c("Na+", "Cl-", "H+", "e-"))
+ species(c("Cl-", "NaCl", "HCl"))
+ a <- suppressMessages(affinity(T = T, P = P, pH = pH, "Na+" = log10(m_Na), IS = IS, transect = TRUE))
}
+ # Speciate Cl-
+ e <- suppressMessages(equilibrate(a, loga.balance = log10(m_tot)))
+ # Get molality of each Cl-bearing species
+ m_Cl <- 10^e$loga.equil[[1]]
+ m_NaCl <- 10^e$loga.equil[[2]]
+ if(identical(pH.arg, NA)) m_HCl <- NA else m_HCl <- 10^e$loga.equil[[3]]
+ # Store previous ionic strength and molality of Na+
+ IS_prev <- IS
+ m_Na_prev <- m_Na
+ # Calculate new molality of Na+ and deviation
+ m_Na <- m_tot - m_NaCl
+ # Only go halfway to avoid overshoot
+ if(attenuate) m_Na <- (m_Na_prev + m_Na) / 2
+ dm_Na <- m_Na - m_Na_prev
+ # Calculate ionic strength and deviation
+ IS <- (m_Na + m_Cl) / 2
+ dIS <- IS - IS_prev
+ # Keep going until the deviations in ionic strength and molality of Na+ at all temperatures are less than tolerance
+ converged <- abs(dIS) < tolerance & abs(dm_Na) < tolerance
+ if(all(converged)) {
+ # Add one step without attenuating the deviation of m_Na
+ if(attenuate) attenuate <- FALSE else break
+ }
}
- # Assemble final molality of Cl- and gammas
- gammas <- suppressMessages(nonideal(ispecies, speciesprops, IS = IS, T = convert(T, "K"), P = P, A_DH = wout$A_DH, B_DH = wout$B_DH))
- gam_Na <- 10^gammas[[1]]$loggam
- gam_Cl <- 10^gammas[[2]]$loggam
- gam_NaCl <- 10^gammas[[3]]$loggam
+ if(i == maxiter) {
+ if(attenuate) stop(paste("reached", maxiter, "iterations without converging"))
+ else stop(paste("reached", maxiter, "iterations without converging (try setting attenuate = TRUE)"))
+ }
+
+ # Restore thermo data frame
+ assign("thermo", thermo, CHNOSZ)
# Return the calculated values
- list(IS = IS, m_Cl = m_Cl, gam_Na = gam_Na, gam_Cl = gam_Cl, gam_NaCl = gam_NaCl)
+ list(IS = IS, m_Na = m_Na, m_Cl = m_Cl, m_NaCl = m_NaCl, m_HCl = m_HCl)
+
}
Modified: pkg/CHNOSZ/demo/contour.R
===================================================================
--- pkg/CHNOSZ/demo/contour.R 2022-12-09 06:17:33 UTC (rev 760)
+++ pkg/CHNOSZ/demo/contour.R 2022-12-11 04:21:41 UTC (rev 761)
@@ -20,7 +20,7 @@
# This gets us close to total S = 0.01 m
basis("H2S", -2)
# Calculate solution composition for 1 mol/kg NaCl
-NaCl <- NaCl(T = T, P = P, m_tot=1)
+NaCl <- NaCl(m_tot = 1, T = T, P = P)
basis("Cl-", log10(NaCl$m_Cl))
# Calculate affinity with changing basis species
bases <- c("H2S", "HS-", "HSO4-", "SO4-2")
Modified: pkg/CHNOSZ/demo/gold.R
===================================================================
--- pkg/CHNOSZ/demo/gold.R 2022-12-09 06:17:33 UTC (rev 760)
+++ pkg/CHNOSZ/demo/gold.R 2022-12-11 04:21:41 UTC (rev 761)
@@ -102,7 +102,7 @@
# NaCl solution with total chloride equal to specified NaCl + KCl solution,
# then estimate the molality of K+ in that solution 20181109
chloride <- function(T, P, m_NaCl, m_KCl) {
- NaCl <- NaCl(T = T, P = P, m_tot = m_NaCl + m_KCl)
+ NaCl <- NaCl(m_tot = m_NaCl + m_KCl, T = T, P = P)
# Calculate logK of K+ + Cl- = KCl, adjusted for ionic strength
logKadj <- subcrt(c("K+", "Cl-", "KCl"), c(-1, -1, 1), T = T, P = P, IS = NaCl$IS)$out$logK
# What is the molality of K+ from 0.5 mol KCl in solution with 2 mol total Cl
Modified: pkg/CHNOSZ/demo/minsol.R
===================================================================
--- pkg/CHNOSZ/demo/minsol.R 2022-12-09 06:17:33 UTC (rev 760)
+++ pkg/CHNOSZ/demo/minsol.R 2022-12-11 04:21:41 UTC (rev 761)
@@ -26,7 +26,7 @@
# Molality of NaCl
mNaCl <- 1000 * wNaCl / (mass("NaCl") * (1 - wNaCl))
# Estimate ionic strength and molality of Cl-
-sat <- NaCl(T = T, m_tot = mNaCl)
+sat <- NaCl(m_tot = mNaCl, T = T)
basis("Cl-", log10(sat$m_Cl))
# Add minerals and aqueous species
Modified: pkg/CHNOSZ/demo/sphalerite.R
===================================================================
--- pkg/CHNOSZ/demo/sphalerite.R 2022-12-09 06:17:33 UTC (rev 760)
+++ pkg/CHNOSZ/demo/sphalerite.R 2022-12-11 04:21:41 UTC (rev 761)
@@ -12,7 +12,7 @@
# 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
- NaCl <- NaCl(T = T, P = P, m_tot = m_tot)
+ NaCl <- NaCl(m_tot = m_tot, T = T, P = P)
basis("Cl-", log10(NaCl$m_Cl))
basis("H2S", log10(0.05))
Modified: pkg/CHNOSZ/inst/NEWS.Rd
===================================================================
--- pkg/CHNOSZ/inst/NEWS.Rd 2022-12-09 06:17:33 UTC (rev 760)
+++ pkg/CHNOSZ/inst/NEWS.Rd 2022-12-11 04:21:41 UTC (rev 761)
@@ -12,7 +12,7 @@
% links to vignettes 20220723
\newcommand{\viglink}{\ifelse{html}{\out{<a href="../CHNOSZ/doc/#1.html"><strong>#1.Rmd</strong></a>}}{\bold{#1.Rmd}}}
-\section{Changes in CHNOSZ version 1.9.9-51 (2022-12-08)}{
+\section{Changes in CHNOSZ version 1.9.9-53 (2022-12-11)}{
\subsection{MAJOR USER-VISIBLE CHANGES}{
\itemize{
@@ -205,7 +205,9 @@
\item \code{seq2aa()} now has the \strong{sequence} argument first and a
default of NA for \strong{protein} (the protein name).
- \item \code{NaCl()} now works with \code{mtot = 0}.
+ \item \code{NaCl()} has been rewritten to include pH dependence (i.e.
+ formation of HCl as well as NaCl) by using \code{affinity()} and
+ \code{equilibrate()} instead of algebraic equations.
\item Remove unused \samp{cutoff} value in \code{thermo()$opt}.
Modified: pkg/CHNOSZ/man/NaCl.Rd
===================================================================
--- pkg/CHNOSZ/man/NaCl.Rd 2022-12-09 06:17:33 UTC (rev 760)
+++ pkg/CHNOSZ/man/NaCl.Rd 2022-12-11 04:21:41 UTC (rev 761)
@@ -7,66 +7,66 @@
}
\usage{
- NaCl(T = seq(100, 500, 100), P = 1000, m_tot = 2, ...)
+ NaCl(m_tot = 1, T = 25, P = "Psat", pH = NA, attenuate = FALSE)
}
\arguments{
+ \item{m_tot}{numeric, total molality of NaCl (single value)}
\item{T}{numeric, temperature in \degC}
- \item{P}{numeric, pressure in bar (single value)}
- \item{m_tot}{numeric, total molality of NaCl (single value)}
- \item{...}{additional arguments for \code{\link{subcrt}}}
+ \item{P}{numeric, pressure in bar}
+ \item{pH}{numeric, pH}
+ \item{attenuate}{logical, halve changes of variables in each step?}
}
\details{
-This function calculates speciation (ion activities) and ionic strength in aqueous solutions given a total amount (\code{m_tot}, in mol/kg) of NaCl.
-The function is written for quick calculations along a temperature range (\code{T}) at constant pressure (\code{P}).
-The only reaction considered is Na\S{+} + Cl\S{-} = NaCl(aq).
-The algorithm starts by calculating the equilibrium constant (\emph{K}) of the reaction and assuming complete dissociation of NaCl(aq).
-This also permits calculating the ionic strength from the molalities of Na\S{+} and Cl\S{-}.
-Then, \code{\link{uniroot}} is used to find the equilibrium molality of Cl\S{-}; that is, where the affinity of the reaction (log(\emph{K}/\emph{Q})) becomes zero.
-The activity quotient (\emph{Q}) is evaluated taking account of activity coefficients of Na\S{+}, Cl\S{-}, and NaCl(aq) calculated for the nominal ionic strength (see \code{\link{nonideal}}).
-The calculated molality of Cl\S{-} yields a new estimate of the ionic strength of the system.
-The calculations are iterated until the deviation in ionic strength at all temperatures is less than 0.01.
+Thermodynamic models for metal solubility and speciation involving chloride complexes are commonly specified in terms of amount of NaCl rather than activity (or molality) of Cl\S{-} as an independent variable.
+This function calculates distribution of species and ionic strength in a simple aqueous solution given a total amount (\code{m_tot}, in mol/kg) of NaCl.
+The aqueous Cl-bearing species considered in the system are Cl\S{-}, NaCl, and optionally HCl, and Na\S{+} is present as a basis species; other Na-bearing species such as NaOH are not considered.
+The activity coefficients of charged species are calculated using the Debye-Hückel equation (see \code{\link{nonideal}}) via the \code{IS} argument of \code{\link{affinity}}.
+The function first sets the molality of Na\S{+} and ionic strength equal to \code{m_tot}, then calculates the distribution of Cl-bearing species.
+Based on mass balance of Na atoms, the molality of NaCl is then used to recalculate the molality of Na\S{+}, followed by ionic strength.
+To find a solution, the function iterates until the change of molality of Na\S{+} and ionic strength are both less than \code{m_tot} / 100.
+
+In some cases, the iteration may oscillate around the true values without converging.
+Setting \code{attenuate} to TRUE, which halves the amount of change of these values in each step, may help with convergence.
+If a solution is not found after 100 iterations, the function stops with an error.
+
+If \code{pH} is NA (the default), then HCl is not included in the calculation and its molality in the output is also assigned NA.
+Note that only a single value is accepted for \code{m_tot}, but the other numeric arguments can have length > 1, allowing multiple combinations \code{T}, \code{P}, and \code{pH} in a single function call.
+However, due to limitations in \code{\link{affinity}}, only one of \code{T} and \code{P} can have length > 1.
}
\section{Warning}{
-This function provides only a first-order estimate of the solution composition, and is intended for solubility calculations of relatively insoluble metals in NaCl-dominated solutions.
-The formation of other species such as HCl or NaOH is not accounted for.
+It is important to keep in mind the ionic strength limits of the Debye-Hückel equation, but this function doesn't enforce them.
+Furthermore, metal-ligand complexing is not calculated by this function, so metal solubility and speciation calculations will be accurate only for relatively insoluble metals in NaCl-dominated solutions.
}
\value{
-A list with components \samp{IS} (\dQuote{true} ionic strength from concentrations of unpaired ions), \samp{m_Cl} (molality of Cl\S{-}), \samp{gam_Na}, and \samp{gam_Cl} (activity coefficients of Na\S{+} and Cl\S{-}).
+A list with components \samp{IS} (ionic strength calculated from molalities of Na\S{+} and Cl\S{-}), \samp{m_Cl}, \samp{m_Cl}, \samp{m_NaCl}, and \samp{m_HCl} (molalities of Na\S{+}, Cl\S{-}, NaCl, and HCl).
}
\seealso{
-\code{demo("gold")} for an application of this function.
+This function is used in a few demos (\code{demo("contour")}, \code{demo("gold")}, \code{demo("minsol")}, \code{demo("sphalerite")}).
}
\examples{\dontshow{reset()}
-# ionic strength of solution and activity coefficient of Cl-
-# from HCh version 3.7 (Shvarov and Bastrakov, 1999) at 1000 bar,
-# 100, 200, and 300 degress C, and 1 to 6 molal NaCl
+# Ionic strength calculated with HCh version 3.7 (Shvarov and Bastrakov, 1999)
+# at 1000 bar, 100, 200, and 300 degress C, and 1 to 6 molal NaCl
m.HCh <- 1:6
IS.HCh <- list(`100`=c(0.992, 1.969, 2.926, 3.858, 4.758, 5.619),
`300`=c(0.807, 1.499, 2.136, 2.739, 3.317, 3.875),
`500`=c(0.311, 0.590, 0.861, 1.125, 1.385, 1.642))
-gam_Cl.HCh <- list(`100`=c(0.565, 0.545, 0.551, 0.567, 0.589, 0.615),
- `300`=c(0.366, 0.307, 0.275, 0.254, 0.238, 0.224),
- `500`=c(0.19, 0.137, 0.111, 0.096, 0.085, 0.077))
-# total molality in the calculation with NaCl()
+# Total molality in the calculation with NaCl()
m_tot <- seq(1, 6, 0.5)
N <- length(m_tot)
-# where we'll put the calculated values
+# Where we'll put the calculated values
IS.calc <- data.frame(`100`=numeric(N), `300`=numeric(N), `500`=numeric(N))
-gam_Cl.calc <- data.frame(`100`=numeric(N), `300`=numeric(N), `500`=numeric(N))
# NaCl() is *not* vectorized over m_tot, so we use a loop here
for(i in 1:length(m_tot)) {
- NaCl.out <- NaCl(c(100, 300, 500), P=1000, m_tot=m_tot[i])
+ NaCl.out <- NaCl(m_tot[i], c(100, 300, 500), P = 1000)
IS.calc[i, ] <- NaCl.out$IS
- gam_Cl.calc[i, ] <- NaCl.out$gam_Cl
}
-# plot ionic strength from HCh and NaCl() as points and lines
-opar <- par(mfrow=c(2, 1))
+# Plot ionic strength from HCh and NaCl() as points and lines
col <- c("black", "red", "orange")
plot(c(1,6), c(0,6), xlab="NaCl (mol/kg)", ylab=axis.label("IS"), type="n")
for(i in 1:3) {
@@ -76,22 +76,11 @@
points(m.HCh, IS.HCh[[i]], col=col[i])
lines(m_tot, IS.calc[, i], col=col[i])
}
-# add 1:1 line, legend, and title
-abline(0, 1, lty=3)
+# Add legend and title
dprop <- describe.property(rep("T", 3), c(100, 300, 500))
legend("topleft", dprop, lty=1, pch=1, col=col)
title(main="H2O + NaCl; HCh (points) and 'NaCl()' (lines)")
-# plot activity coefficient (gamma)
-plot(c(1,6), c(0,0.8), xlab="NaCl (mol/kg)", ylab=expression(gamma[Cl^"-"]), type="n")
-for(i in 1:3) {
- points(m.HCh, gam_Cl.HCh[[i]], col=col[i])
- lines(m_tot, gam_Cl.calc[, i], col=col[i])
}
-# we should be fairly close
-maxdiff <- function(x, y) max(abs(y - x))
-stopifnot(maxdiff(unlist(gam_Cl.calc[seq(1,11,2), ]), unlist(gam_Cl.HCh)) < 0.033)
-par(opar)
-}
\references{
Shvarov, Y. and Bastrakov, E. (1999) HCh: A software package for geochemical equilibrium modelling. User's Guide. \emph{Australian Geological Survey Organisation} \bold{1999/25}. \url{http://pid.geoscience.gov.au/dataset/ga/25473}
Modified: pkg/CHNOSZ/vignettes/customizing.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/customizing.Rmd 2022-12-09 06:17:33 UTC (rev 760)
+++ pkg/CHNOSZ/vignettes/customizing.Rmd 2022-12-11 04:21:41 UTC (rev 761)
@@ -467,9 +467,9 @@
### Diagram 1: Constant molality of `r F_`
Now we're ready to make a speciation diagram.
-Our target is to reproduce Fig. 7b of @WWH_21, which is made for 300 °C.
+Our aim is to reproduce Fig. 7b of @WWH_21, which is made for 300 °C.
A constant molality of `r F_` is based on the assumption of complete dissociation of 0.1 m NaF (we'll change this later).
-An ionic strength of 1.2 mol/kg is estimated for a solution with 1.8 m NaCl (use `NaCl(300, m_tot = 1.8)`).
+An ionic strength of 0.9 mol/kg is estimated for a solution with 1.8 m NaCl (use `NaCl(1.8, T = 300)`).
`r NOTE`: because the ionic strength is non-zero, the calculations here refer to molality instead of activity of species (see [An Introduction to CHNOSZ](anintro.html#from-activity-to-molality)).
```{r diagram1, message = FALSE, results = "hide", fig.width = 6, fig.height = 5, out.width = "75%", fig.align = "center", pngquant = pngquant}
@@ -477,7 +477,7 @@
basis("F-", log10(0.1))
iaq <- retrieve("W", c("O", "H", "F"), "aq")
species(iaq)
-a <- affinity(pH = c(2, 7), T = 300, IS = 1.2)
+a <- affinity(pH = c(2, 7), T = 300, IS = 0.9)
e <- equilibrate(a)
col <- c(1, 4, 5, 2)
diagram(e, alpha = TRUE, col = col, lty = 1, lwd = 2, ylab = "Fraction total W")
@@ -495,7 +495,7 @@
```{r a_F}
T <- 300
pH <- seq(2, 7, 0.1)
-logK_HF <- subcrt(c("H+", "F-", "HF"), c(-1, -1, 1), T = T, IS = 1.2)$out$logK
+logK_HF <- subcrt(c("H+", "F-", "HF"), c(-1, -1, 1), T = T, IS = 0.9)$out$logK
F_tot <- 0.1
a_F <- F_tot / (1 + 10^(logK_HF - pH))
```
@@ -505,7 +505,7 @@
basis(c("H+", "WO4-2", "F-", "H2O", "O2"))
iaq <- retrieve("W", c("O", "H", "F"), "aq")
species(iaq)
-a <- affinity(pH = pH, "F-" = log10(a_F), T = T, IS = 1.2)
+a <- affinity(pH = pH, "F-" = log10(a_F), T = T, IS = 0.9)
e <- equilibrate(a)
diagram(e, alpha = TRUE, col = col, lty = 1, lwd = 2, ylab = "Fraction total W")
```
Modified: pkg/CHNOSZ/vignettes/multi-metal.Rmd
===================================================================
--- pkg/CHNOSZ/vignettes/multi-metal.Rmd 2022-12-09 06:17:33 UTC (rev 760)
+++ pkg/CHNOSZ/vignettes/multi-metal.Rmd 2022-12-11 04:21:41 UTC (rev 761)
@@ -555,7 +555,7 @@
# Setup basis species
basis(c("Cu+", "pyrite", "H2S", "oxygen", "H2O", "H+", "Cl-"))
basis("H2S", logmS)
-nacl <- NaCl(T = T, P = "Psat", m_tot = m_NaCl)
+nacl <- NaCl(m_tot = m_NaCl, T = T, P = "Psat")
basis("Cl-", log10(nacl$m_Cl))
# Fe-bearing minerals
species(Fe.cr)
More information about the CHNOSZ-commits
mailing list