Calculate logK for Eh reaction and O2 solubility

jedick · Jun 24, 2020 · 22c0c89 · 22c0c89
1 parent 10328cf
commit 22c0c89
Show file tree

Hide file tree

Showing 4 changed files with 52 additions and 26 deletions.
diff --git a/R/readhead.R b/R/readhead.R
@@ -39,6 +39,8 @@ readhead <- function(LINES, quiet = FALSE) {
   ih2o_2 <- grep("* c h2o 2", LINES, fixed = TRUE)
   ih2o_3 <- grep("* c h2o 3", LINES, fixed = TRUE)
   ih2o_4 <- grep("* c h2o 4", LINES, fixed = TRUE)
+  ieh <- grep("* log k for eh reaction", LINES, fixed = TRUE)
+  io2 <- grep("* log k for o2 gas solubility", LINES, fixed = TRUE)
   # find header lines before sections
   jbasis <- grep("basis species$", LINES)
   jredox <- grep("redox couples$", LINES)
@@ -97,6 +99,8 @@ readhead <- function(LINES, quiet = FALSE) {
     ih2o_2 = ih2o_2,
     ih2o_3 = ih2o_3,
     ih2o_4 = ih2o_4,
+    ieh = ieh,
+    io2 = io2,
     ibasis = jbasis,
     iredox = jredox,
     iaqueous = jaqueous,

diff --git a/R/writedat.R b/R/writedat.R
@@ -31,6 +31,9 @@ writedat <- function(outfile, LINES, HEAD, LOGK, ADDS, infile, DH.method, a0_ion
   # calculate coefficients for CO2 and H2O activity 20200623
   CO2 <- co2(LOGK$T)
   H2O <- h2o(LOGK$T)
+  # calculate logK for Eh and O2 solubility reactions 20200624
+  logK_Eh <- suppressMessages(CHNOSZ::subcrt(c("H2O", "oxygen", "H+", "e-"), c(-2, 1, 4, 4), T = LOGK$T, P = LOGK$P)$out$logK)
+  logK_O2 <- suppressMessages(CHNOSZ::subcrt(c("oxygen", "O2"), c(-1, 1), T = LOGK$T, P = LOGK$P)$out$logK)
   # loop over lines of the input file
   for(i in 1:length(LINES)) {
 
@@ -46,35 +49,40 @@ writedat <- function(outfile, LINES, HEAD, LOGK, ADDS, infile, DH.method, a0_ion
         outline <- NA
       }
       # output T and P 20200610
-      if(i == HEAD$iT + 1) outline <- formatline(LOGK$T, 1)
-      if(i == HEAD$iT + 2) outline <- formatline(LOGK$T, 2)
-      if(i == HEAD$iP + 1) outline <- formatline(LOGK$P, 1)
-      if(i == HEAD$iP + 2) outline <- formatline(LOGK$P, 2)
+      if(identical(i, HEAD$iT + 1L)) outline <- formatline(LOGK$T, 1)
+      if(identical(i, HEAD$iT + 2L)) outline <- formatline(LOGK$T, 2)
+      if(identical(i, HEAD$iP + 1L)) outline <- formatline(LOGK$P, 1)
+      if(identical(i, HEAD$iP + 2L)) outline <- formatline(LOGK$P, 2)
       # output Debye-Hückel parameters 20200622
-      if(i == HEAD$iadh + 1) outline <- formatline(DH$adh, 1)
-      if(i == HEAD$iadh + 2) outline <- formatline(DH$adh, 2)
-      if(i == HEAD$ibdh + 1) outline <- formatline(DH$bdh * 1e-8, 1)
-      if(i == HEAD$ibdh + 2) outline <- formatline(DH$bdh * 1e-8, 2)
-      if(i == HEAD$ibdot + 1) outline <- formatline(DH$bdot, 1)
-      if(i == HEAD$ibdot + 2) outline <- formatline(DH$bdot, 2)
+      if(identical(i, HEAD$iadh + 1L)) outline <- formatline(DH$adh, 1)
+      if(identical(i, HEAD$iadh + 2L)) outline <- formatline(DH$adh, 2)
+      if(identical(i, HEAD$ibdh + 1L)) outline <- formatline(DH$bdh * 1e-8, 1)
+      if(identical(i, HEAD$ibdh + 2L)) outline <- formatline(DH$bdh * 1e-8, 2)
+      if(identical(i, HEAD$ibdot + 1L)) outline <- formatline(DH$bdot, 1)
+      if(identical(i, HEAD$ibdot + 2L)) outline <- formatline(DH$bdot, 2)
       # output CO2 parameters 20200623
-      if(i == HEAD$ico2_1 + 1) outline <- formatline(CO2$co2_1, 1, ndec = 6)
-      if(i == HEAD$ico2_1 + 2) outline <- formatline(CO2$co2_1, 2, ndec = 6)
-      if(i == HEAD$ico2_2 + 1) outline <- formatline(CO2$co2_2, 1, ndec = 6)
-      if(i == HEAD$ico2_2 + 2) outline <- formatline(CO2$co2_2, 2, ndec = 6)
-      if(i == HEAD$ico2_3 + 1) outline <- formatline(CO2$co2_3, 1, ndec = 6)
-      if(i == HEAD$ico2_3 + 2) outline <- formatline(CO2$co2_3, 2, ndec = 6)
-      if(i == HEAD$ico2_4 + 1) outline <- formatline(CO2$co2_4, 1, ndec = 6)
-      if(i == HEAD$ico2_4 + 2) outline <- formatline(CO2$co2_4, 2, ndec = 6)
+      if(identical(i, HEAD$ico2_1 + 1L)) outline <- formatline(CO2$co2_1, 1, ndec = 6)
+      if(identical(i, HEAD$ico2_1 + 2L)) outline <- formatline(CO2$co2_1, 2, ndec = 6)
+      if(identical(i, HEAD$ico2_2 + 1L)) outline <- formatline(CO2$co2_2, 1, ndec = 6)
+      if(identical(i, HEAD$ico2_2 + 2L)) outline <- formatline(CO2$co2_2, 2, ndec = 6)
+      if(identical(i, HEAD$ico2_3 + 1L)) outline <- formatline(CO2$co2_3, 1, ndec = 6)
+      if(identical(i, HEAD$ico2_3 + 2L)) outline <- formatline(CO2$co2_3, 2, ndec = 6)
+      if(identical(i, HEAD$ico2_4 + 1L)) outline <- formatline(CO2$co2_4, 1, ndec = 6)
+      if(identical(i, HEAD$ico2_4 + 2L)) outline <- formatline(CO2$co2_4, 2, ndec = 6)
       # output H2O parameters 20200623
-      if(i == HEAD$ih2o_1 + 1) outline <- formatline(H2O$h2o_1, 1, ndec = 6)
-      if(i == HEAD$ih2o_1 + 2) outline <- formatline(H2O$h2o_1, 2, ndec = 6)
-      if(i == HEAD$ih2o_2 + 1) outline <- formatline(H2O$h2o_2, 1, ndec = 6)
-      if(i == HEAD$ih2o_2 + 2) outline <- formatline(H2O$h2o_2, 2, ndec = 6)
-      if(i == HEAD$ih2o_3 + 1) outline <- formatline(H2O$h2o_3, 1, ndec = 6)
-      if(i == HEAD$ih2o_3 + 2) outline <- formatline(H2O$h2o_3, 2, ndec = 6)
-      if(i == HEAD$ih2o_4 + 1) outline <- formatline(H2O$h2o_4, 1, ndec = 6)
-      if(i == HEAD$ih2o_4 + 2) outline <- formatline(H2O$h2o_4, 2, ndec = 6)
+      if(identical(i, HEAD$ih2o_1 + 1L)) outline <- formatline(H2O$h2o_1, 1, ndec = 6)
+      if(identical(i, HEAD$ih2o_1 + 2L)) outline <- formatline(H2O$h2o_1, 2, ndec = 6)
+      if(identical(i, HEAD$ih2o_2 + 1L)) outline <- formatline(H2O$h2o_2, 1, ndec = 6)
+      if(identical(i, HEAD$ih2o_2 + 2L)) outline <- formatline(H2O$h2o_2, 2, ndec = 6)
+      if(identical(i, HEAD$ih2o_3 + 1L)) outline <- formatline(H2O$h2o_3, 1, ndec = 6)
+      if(identical(i, HEAD$ih2o_3 + 2L)) outline <- formatline(H2O$h2o_3, 2, ndec = 6)
+      if(identical(i, HEAD$ih2o_4 + 1L)) outline <- formatline(H2O$h2o_4, 1, ndec = 6)
+      if(identical(i, HEAD$ih2o_4 + 2L)) outline <- formatline(H2O$h2o_4, 2, ndec = 6)
+      # output logKs for Eh reaction and O2 solubility 20200624
+      if(identical(i, HEAD$ieh + 1L)) outline <- formatline(logK_Eh, 1)
+      if(identical(i, HEAD$ieh + 2L)) outline <- formatline(logK_Eh, 2)
+      if(identical(i, HEAD$io2 + 1L)) outline <- formatline(logK_O2, 1)
+      if(identical(i, HEAD$io2 + 2L)) outline <- formatline(logK_O2, 2)
     }
 
     # check if we're in the basis species header

diff --git a/inst/extdata/dates.dat b/inst/extdata/dates.dat
diff --git a/man/logKcalc.Rd b/man/logKcalc.Rd
@@ -105,6 +105,20 @@ Because of a bug with calculating exponents in some versions of K2GWB, the CO\s2
 
 }
 
+\section{Dataset Formats}{
+
+The conversion is compatible with all existing dataset formats of GWB files.
+For dataset formats \samp{feb94} and \samp{oct94}, \code{logKcalc} updates the \logK blocks for the Eh reaction and O\s2 gas solubility.
+(Blocks for H\s2 and N\s2 gas solubility, which are not used by GWB, are not updated.)
+
+These blocks are not present in dataset formats \samp{oct13} and later.
+Instead, the free electron is present as a species; \code{logKcalc} updates its \logK values.
+
+Dataset formats \samp{jul17} and \samp{jan19} specify a fugacity model (tsonopoulos, peng-robinson, or spycher-reed).
+The parameters are located in the data blocks for individual gas species and are not modified by \code{logKcalc}.
+
+}
+
 \seealso{
 \code{\link{modOBIGT}} for preparing the thermodynamic database for the \logK calculation.
 }