[CHNOSZ-commits] r339 - in pkg/CHNOSZ: . R demo inst man
noreply at r-forge.r-project.org
noreply at r-forge.r-project.org
Sun Nov 4 15:44:14 CET 2018
Author: jedick
Date: 2018-11-04 15:44:13 +0100 (Sun, 04 Nov 2018)
New Revision: 339
Added:
pkg/CHNOSZ/R/NaCl.R
pkg/CHNOSZ/man/NaCl.Rd
Modified:
pkg/CHNOSZ/DESCRIPTION
pkg/CHNOSZ/NAMESPACE
pkg/CHNOSZ/demo/gold.R
pkg/CHNOSZ/inst/NEWS
pkg/CHNOSZ/man/solubility.Rd
Log:
add NaCl() for first-order estimate of H2O + NaCl solutions
Modified: pkg/CHNOSZ/DESCRIPTION
===================================================================
--- pkg/CHNOSZ/DESCRIPTION 2018-11-01 07:25:18 UTC (rev 338)
+++ pkg/CHNOSZ/DESCRIPTION 2018-11-04 14:44:13 UTC (rev 339)
@@ -1,6 +1,6 @@
-Date: 2018-11-01
+Date: 2018-11-04
Package: CHNOSZ
-Version: 1.1.3-46
+Version: 1.1.3-47
Title: Thermodynamic Calculations and Diagrams for Geo(bio)chemistry
Authors at R: c(
person("Jeffrey", "Dick", , "j3ffdick at gmail.com", role = c("aut", "cre"),
Modified: pkg/CHNOSZ/NAMESPACE
===================================================================
--- pkg/CHNOSZ/NAMESPACE 2018-11-01 07:25:18 UTC (rev 338)
+++ pkg/CHNOSZ/NAMESPACE 2018-11-04 14:44:13 UTC (rev 339)
@@ -56,7 +56,7 @@
"calculateEpsilon", "calculateQ", "water.DEW", "berman",
"maxdiff", "expect_maxdiff", "Bdot",
# added 20171121 or later
- "dumpdata", "thermo.axis", "solubility"
+ "dumpdata", "thermo.axis", "solubility", "NaCl"
)
# Load shared objects
Added: pkg/CHNOSZ/R/NaCl.R
===================================================================
--- pkg/CHNOSZ/R/NaCl.R (rev 0)
+++ pkg/CHNOSZ/R/NaCl.R 2018-11-04 14:44:13 UTC (rev 339)
@@ -0,0 +1,48 @@
+# NaCl.R
+# calculate ion activities and ionic strength
+# given a total molality of NaCl
+# taking account of ion association: Na+ + Cl- = NaCl(aq)
+# 20181102 jmd first version
+
+NaCl <- function(T=seq(100, 500, 100), P=1000, m_tot=2) {
+ # define a function for the reaction quotient
+ logQ <- function(m_Cl, gamma) {
+ # starting with Q = a_NaCl / (a_Na+ * a_Cl-),
+ # substitute gam_NaCl = 0, m_NaCl + m_Cl = m_tot, m_Cl = m_Na, gam_Cl = gam_Na = gamma
+ # to write:
+ log10( (m_tot - m_Cl) / (m_Cl * gamma) ^ 2 )
+ }
+ # define a function for affinity = log(K / Q)
+ A <- function(m_Cl, gamma, logK) logK - logQ(m_Cl, gamma)
+ # 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 species index for Cl-
+ iCl <- info("Cl-")
+ # we start by doing calculations for all temperatures
+ doit <- !logical(N)
+ while(any(doit)) {
+ # calculate activity coefficient at given ionic strength
+ gamma <- suppressMessages(10^nonideal(iCl, list(data.frame(G=numeric(N))), IS=IS, T=convert(T, "K"), P=P, A_DH=wout$A_DH, B_DH=wout$B_DH)[[1]]$loggam)
+ # solve for m_Cl
+ for(i in which(doit)) m_Cl[i] <- uniroot(A, c(0, m_tot), gamma=gamma[i], logK=logK[i])$root
+ # calculate new ionic strength and deviation
+ ISnew <- m_Cl
+ dIS <- ISnew - IS
+ # set net 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
+ }
+ # assemble final gamma and activity of Cl-
+ gamma <- suppressMessages(10^nonideal(iCl, list(data.frame(G=numeric(N))), IS=IS, T=convert(T, "K"), P=P, A_DH=wout$A_DH, B_DH=wout$B_DH)[[1]]$loggam)
+ a_Cl <- m_Cl * gamma
+ # return the calculated values
+ list(IS=IS, gamma=gamma, a_Cl=a_Cl)
+}
Modified: pkg/CHNOSZ/demo/gold.R
===================================================================
--- pkg/CHNOSZ/demo/gold.R 2018-11-01 07:25:18 UTC (rev 338)
+++ pkg/CHNOSZ/demo/gold.R 2018-11-04 14:44:13 UTC (rev 339)
@@ -31,6 +31,8 @@
# set up system
# use H2S here: it's the predominant species at the pH of the QMK buffer -- see sulfur()
basis(c("Al2O3", "SiO2", "Fe", "Au", "K+", "Cl-", "H2S", "H2O", "oxygen", "H+"))
+# set activity of K+ for 0.5 molal KCl assuming complete dissociation
+basis("K+", log10(0.5))
# create a pH buffer
mod.buffer("QMK", c("quartz", "muscovite", "K-feldspar"), "cr", 0)
@@ -130,23 +132,22 @@
basis("H2S", "PPM")
# apply QMK buffer for pH
basis("H+", "QMK")
- # at 400 degC, 1000 bar, and IS=2, the logarithm of the activity coefficient of Cl- is -0.66:
- loggam <- subcrt("Cl-", T=400, P=1000, IS=2)$out[[1]]$loggam
- # for a total molality of 2 m (1.5 m NaCl and 0.5 m KCl), the activity of Cl- is about 10:
- actCl <- 2/10^loggam
- basis("Cl-", log10(actCl))
+ # calculate solution composition for 2 mol/kg NaCl
+ NaCl <- NaCl(T = seq(150, 550, 10), P = 1000, m_tot=2)
# calculate affinity, equilibrate, solubility
- a <- affinity(T = c(150, 550), P = 1000, IS = 2)
+ a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(NaCl$a_Cl), P = 1000, IS = NaCl$IS)
e <- equilibrate(a)
s <- solubility(e)
# make diagram and show total log molality
diagram(s, ylim = c(-10, -4), col = col, lwd = 2, lty = 1)
diagram(s, add = TRUE, type = "loga.balance", lwd = 3, lty = 2)
# make legend and title
- dprop <- describe.property(c("P", "IS"), c(1000, 2))
- legend("topleft", dprop, bty = "n")
- dbasis <- describe.basis(ibasis = c(6, 9, 7, 10))
- legend("bottomright", dbasis, bty = "n")
+ dP <- describe.property("P", 1000)
+ dNaCl <- expression(NaCl == 2~mol~kg^-1)
+ dK <- describe.basis(ibasis=5)
+ legend("topleft", c(dP, dNaCl, dK), bty = "n")
+ dbasis <- describe.basis(ibasis = c(9, 7, 10))
+ legend("topright", dbasis, bty = "n")
title(main=("After Williams-Jones et al., 2009, Fig. 2B"), font.main = 1)
}
@@ -155,36 +156,29 @@
# (doi:10.1144/SP402.4)
Au_T2 <- function() {
species(c("Au(HS)2-", "Au(HS)", "AuOH", "AuCl2-"))
- # set total activity of H2S
- # TODO: the paper says total S = 0.01 m,
- # but a higher activity makes the diagram here closer to
- # that of Williams-Jones et al., 2009
- basis("H2S", -1)
+ # approximate activity of H2S for total S = 0.01 m
+ basis("H2S", -2)
# apply HM buffer for fO2
basis("O2", "HM")
# apply QMK buffer for pH
basis("H+", "QMK")
- # calculate activity coefficient of Cl- at IS=2
- loggam <- subcrt("Cl-", T = seq(150, 550, 10), P = 1000, IS = 2)$out[[1]]$loggam
- # calculate activity of Cl- given a total molality of 2 m (1.5 m NaCl and 0.5 m KCl)
- actCl <- 2/10^loggam
- # TODO: adjust Cl- activity for increasing logK(T) of Na+ + Cl- = NaCl
+ # calculate solution composition for 2 mol/kg NaCl
+ NaCl <- NaCl(T = seq(150, 550, 10), P = 1000, m_tot=2)
# calculate affinity, equilibrate, solubility
- a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(actCl), P = 1000, IS = 2)
+ a <- affinity(T = seq(150, 550, 10), `Cl-` = log10(NaCl$a_Cl), P = 1000, IS = NaCl$IS)
e <- equilibrate(a)
s <- solubility(e)
# make diagram and show total log molality
diagram(s, ylim = c(-10, -2), col = col, lwd = 2, lty = 1)
diagram(s, add = TRUE, type = "loga.balance", lwd = 3, lty = 2)
# make legend and title
- dprop <- describe.property(c("P", "IS"), c(1000, 2))
- legend("topleft", dprop, bty = "n")
- # show the log molality of Cl-
- basis("Cl-", log10(2))
- dbasis1 <- describe.basis(ibasis = 6, use.molality = TRUE)
- dbasis2 <- describe.basis(ibasis = c(9, 7, 10))
- legend("bottomright", c(dbasis1, dbasis2), bty = "n")
- title(main=("After Williams-Jones et al., 2009, Fig. 2B"), font.main = 1)
+ dP <- describe.property("P", 1000)
+ dNaCl <- expression(NaCl == 2~mol~kg^-1)
+ dK <- describe.basis(ibasis=5)
+ legend("topleft", c(dP, dNaCl, dK), bty = "n")
+ dbasis <- describe.basis(ibasis = c(9, 7, 10))
+ legend("topright", dbasis, bty = "n")
+ title(main=("After Williams-Jones et al., 2009, Fig. 2A"), font.main = 1)
}
# make plots
Modified: pkg/CHNOSZ/inst/NEWS
===================================================================
--- pkg/CHNOSZ/inst/NEWS 2018-11-01 07:25:18 UTC (rev 338)
+++ pkg/CHNOSZ/inst/NEWS 2018-11-04 14:44:13 UTC (rev 339)
@@ -1,4 +1,4 @@
-CHANGES IN CHNOSZ 1.1.3-46 (2018-11-01)
+CHANGES IN CHNOSZ 1.1.3-47 (2018-11-01)
---------------------------------------
NEW FEATURES
@@ -6,10 +6,15 @@
- Add solubility(). Run this after equilibrate() to calculate the
solubility (loga.balance) of the balanced basis species.
-- Add demo/gold.R for calculations of Au solubility (based on diagrams
- from Akinfiev and Zotov, 2001, Stefánsson and Seward, 2004, and
- Williams-Jones et al., 2009).
+- Add NaCl(), implementing a first-order calculation of the speciation
+ of NaCl in water, taking account of activity coefficients and the
+ reaction Na+ + Cl- = NaCl(aq).
+- Add demo/gold.R for calculations of Au solubility in hydrothermal
+ chloride and sulfide solutions (based on diagrams from Akinfiev and
+ Zotov, 2001, Stefánsson and Seward, 2004, and Williams-Jones et al.,
+ 2009).
+
- Revise demo/solubility.R to show solubility calculations for CO2(gas)
and calcite as a function of T and pH.
Added: pkg/CHNOSZ/man/NaCl.Rd
===================================================================
--- pkg/CHNOSZ/man/NaCl.Rd (rev 0)
+++ pkg/CHNOSZ/man/NaCl.Rd 2018-11-04 14:44:13 UTC (rev 339)
@@ -0,0 +1,92 @@
+\encoding{UTF-8}
+\name{NaCl}
+\alias{NaCl}
+\title{Simple NaCl-Water Solution}
+\description{
+Calculate speciation and ionic strength of aqueous solutions with a given molality of NaCl.
+}
+
+\usage{
+ NaCl(T = seq(100, 500, 100), P = 1000, m_tot = 2)
+}
+
+\arguments{
+ \item{T}{numeric, temperature in \degC}
+ \item{P}{numeric, pressure in bar (single value)}
+ \item{m_tot}{numeric, total molality of NaCl (single value)}
+}
+
+\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 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.
+}
+
+\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.
+}
+
+\value{
+A list with components \samp{IS} (\dQuote{true} ionic strength from concentrations of unpaired ions), \samp{gamma} (activity coefficient of Cl\S{-}), \samp{a_Cl} (activity of Cl\S{-}).
+}
+
+\seealso{
+\code{demo("gold")} for an application of this function.
+}
+
+\examples{\dontshow{data(thermo)}
+# ionic strength and activity coefficient of Cl-
+# from HCh (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.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()
+m_tot <- seq(1, 6, 0.5)
+N <- length(m_tot)
+# where we'll put the calculated values
+IS.calc <- data.frame(`100`=numeric(N), `300`=numeric(N), `500`=numeric(N))
+gam.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), m_tot=m_tot[i])
+ IS.calc[i, ] <- NaCl.out$IS
+ gam.calc[i, ] <- NaCl.out$gamma
+}
+# plot ionic strength from HCh and NaCl() as points and lines
+par(mfrow=c(2, 1))
+col <- c("black", "red", "orange")
+plot(c(1,6), c(0,6), xlab="NaCl (mol/kg)", ylab="I (mol/kg)", type="n")
+for(i in 1:3) {
+ 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)
+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,1), xlab="NaCl (mol/kg)", ylab="gamma", type="n")
+for(i in 1:3) {
+ points(m.HCh, gam.HCh[[i]], col=col[i])
+ lines(m_tot, gam.calc[, i], col=col[i])
+}
+}
+
+\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}.
+}
+
+\concept{Extended workflow}
Modified: pkg/CHNOSZ/man/solubility.Rd
===================================================================
--- pkg/CHNOSZ/man/solubility.Rd 2018-11-01 07:25:18 UTC (rev 338)
+++ pkg/CHNOSZ/man/solubility.Rd 2018-11-04 14:44:13 UTC (rev 339)
@@ -38,7 +38,7 @@
}
\seealso{
-\code{demo("solubility")} adds pH-\T diagrams to the CO\s{2} and calcite example here.
+\code{demo("solubility")} adds \T-pH diagrams to the CO\s{2} and calcite example here.
\code{demo("gold")} shows solubility calculations for Au in aqueous solutions with hydroxide, chloride, and hydrosulfide complexes.
}
More information about the CHNOSZ-commits
mailing list