## ----setup, include = FALSE--------------------------------------------------- library(CHNOSZ) options(width = 80) # Use pngquant to optimize PNG images library(knitr) knit_hooks$set(pngquant = hook_pngquant) pngquant <- if(nzchar(Sys.getenv("CHNOSZ_BUILD_LARGE_VIGNETTES"))) "--speed=1 --quality=0-25" else "--speed=1 --quality=0-10" if (!nzchar(Sys.which("pngquant"))) pngquant <- NULL # To make warnings appear in text box 20230619 # https://selbydavid.com/2017/06/18/rmarkdown-alerts/ knitr::knit_hooks$set( error = function(x, options) { paste('\n\n
', gsub('##', '\n', gsub('^##\ Error', '**Error:**', x)), '
', sep = '\n') }, warning = function(x, options) { paste('\n\n
', gsub('##', '\n', gsub('^##\ Warning:', '**Warning:**', x)), '
', sep = '\n') }, message = function(x, options) { paste('\n\n
', gsub('##', '\n', x), '
', sep = '\n') } ) # Set dpi 20231129 knitr::opts_chunk$set( dpi = if(nzchar(Sys.getenv("CHNOSZ_BUILD_LARGE_VIGNETTES"))) 100 else 72 ) ## ----echo = F, cache = F------------------------------------------------------ # Merge consecutive messages into a single div 20231114 knitr::knit_hooks$set(document = function(x){ # Not sure why this is needed, but simply computing on 'x' doesn't work file <- tempfile() writeLines(x, file) x <- readLines(file) # Line numbers of the document with enddiv <- which(x == "") # Line numbers with
beginalert <- which(x == '
') # Find
followed
(skip empty lines) removediv <- (enddiv + 2) %in% beginalert if(any(removediv)) { # Lines to remove rmlines1 <- enddiv[removediv] rmlines2 <- enddiv[removediv] + 1 rmlines3 <- enddiv[removediv] + 2 # Do the removing x <- x[-c(rmlines1, rmlines2, rmlines3)] } x }) ## ----HTML, include = FALSE---------------------------------------------------- NOTE <- 'NOTE' # CHNOSZ functions equilibrate_ <- 'equilibrate()' info_ <- 'info()' thermo.refs_ <- 'thermo.refs()' # Math stuff logK <- "log K" Hplus <- "H+" HCO2_ <- "HCO2−" HCO3_ <- "HCO3−" O2 <- "O2" S2 <- "S2" log <- "log " aHCO2_ <- "aHCO2−" aHCO3_ <- "aHCO3−" logfO2 <- "log fO2" Ctot <- "Ctot" C3H5O2_ <- "C3H5O2−" a3HCO3_ <- "a3HCO3−" aC3H5O2_ <- "aC3H5O2−" a2HCO3_ <- "a2HCO3−" logCtot <- "log Ctot" CO2 <- "CO2" H2O <- "H2O" S3minus <- "S3-" H2S <- "H2S" SO4__ <- "SO4-2" Kplus <- "K+" Naplus <- "Na+" Clminus <- "Cl-" H2 <- "H2" Tmax <- "Tmax" ## ----alanine_refs, message = FALSE-------------------------------------------- info(info("alanine"))[c("ref1", "ref2")] thermo.refs(info("alanine")) ## ----PPM_refs, message = FALSE------------------------------------------------ basis(c("pyrite", "pyrrhotite", "oxygen")) sres <- subcrt("magnetite", 1) info(sres$reaction$ispecies)[, 1:6] thermo.refs(sres) reset() ## ----DEW_Ctot, eval = FALSE--------------------------------------------------- # ### Based on demo/DEW.R in CHNOSZ # # # Set up subplots # par(mfrow = c(1, 2), mar = c(3.0, 3.5, 2.5, 1.0), mgp = c(1.7, 0.3, 0), las = 1, tcl = 0.3, xaxs = "i", yaxs = "i") # # # Activate DEW model # water("DEW") # # ########### # ## logfO2-pH diagram for aqueous inorganic and organic carbon species at high pressure # ## After Figure 1b of Sverjensky et al., 2014b [SSH14] # ## (Nature Geoscience, https://doi.org/10.1038/NGEO2291) # ########### # # # Define system # basis("CHNOS+") # inorganic <- c("CO2", "HCO3-", "CO3-2", "CH4") # organic <- c("formic acid", "formate", "acetic acid", "acetate", "propanoic acid", "propanoate") # species(c(inorganic, organic), 0) # fill <- c(rep("aliceblue", length(inorganic)), rep("aquamarine", length(organic))) # # # A function to make the diagrams # dfun <- function(T = 600, P = 50000, res = 200, Ctot = 0.03) { # a <- affinity(pH = c(0, 10, res), O2 = c(-24, -12, res), T = T, P = P) # # Set total C concentration to 0.03 molal # # (from EQ3NR model for eclogite [Supporting Information of SSH14]) # e <- equilibrate(a, loga.balance = log10(Ctot)) # diagram(e, fill = fill) # dp <- describe.property(c(" T", " P"), c(T, P), digits = 0) # legend("bottomleft", legend = dp, bty = "n") # } # # water("DEW") # add.OBIGT("DEW") # dfun(Ctot = 0.03) # mtitle(c("Inorganic and organic species", "C[total] = 0.03 molal")) # dfun(Ctot = 3) # mtitle(c("Inorganic and organic species", "C[total] = 3 molal")) # # # Restore default settings for the questions below # reset() ## ----DEW_Ctot, echo = FALSE, message = FALSE, results = "hide", fig.width = 8, fig.height = 4, out.width = "100%", fig.align = "center", pngquant = pngquant, cache = TRUE---- ### Based on demo/DEW.R in CHNOSZ # Set up subplots par(mfrow = c(1, 2), mar = c(3.0, 3.5, 2.5, 1.0), mgp = c(1.7, 0.3, 0), las = 1, tcl = 0.3, xaxs = "i", yaxs = "i") # Activate DEW model water("DEW") ########### ## logfO2-pH diagram for aqueous inorganic and organic carbon species at high pressure ## After Figure 1b of Sverjensky et al., 2014b [SSH14] ## (Nature Geoscience, https://doi.org/10.1038/NGEO2291) ########### # Define system basis("CHNOS+") inorganic <- c("CO2", "HCO3-", "CO3-2", "CH4") organic <- c("formic acid", "formate", "acetic acid", "acetate", "propanoic acid", "propanoate") species(c(inorganic, organic), 0) fill <- c(rep("aliceblue", length(inorganic)), rep("aquamarine", length(organic))) # A function to make the diagrams dfun <- function(T = 600, P = 50000, res = 200, Ctot = 0.03) { a <- affinity(pH = c(0, 10, res), O2 = c(-24, -12, res), T = T, P = P) # Set total C concentration to 0.03 molal # (from EQ3NR model for eclogite [Supporting Information of SSH14]) e <- equilibrate(a, loga.balance = log10(Ctot)) diagram(e, fill = fill) dp <- describe.property(c(" T", " P"), c(T, P), digits = 0) legend("bottomleft", legend = dp, bty = "n") } water("DEW") add.OBIGT("DEW") dfun(Ctot = 0.03) mtitle(c("Inorganic and organic species", "C[total] = 0.03 molal")) dfun(Ctot = 3) mtitle(c("Inorganic and organic species", "C[total] = 3 molal")) # Restore default settings for the questions below reset() ## ----pyrrhotite_values, message = FALSE--------------------------------------- # The formula of the new mineral and literature reference formula <- "Fe0.877S" ref1 <- "PMW87" # Because the properties from Pankratz et al. (1987) are listed in calories, # we need to change the output of subcrt() to calories (the default is Joules) E.units("cal") # Use temperature in Kelvin for the calculations below T.units("K") # Thermodynamic properties of polymorph 1 at 25 °C (298.15 K) G1 <- -25543 H1 <- -25200 S1 <- 14.531 Cp1 <- 11.922 # Heat capacity coefficients for polymorph 1 a1 <- 7.510 b1 <- 0.014798 # For volume, use the value from Helgeson et al. (1978) V1 <- V2 <- 18.2 # Transition temperature Ttr <- 598 # Transition enthalpy (cal/mol) DHtr <- 95 # Heat capacity coefficients for polymorph 2 a2 <- -1.709 b2 <- 0.011746 c2 <- 3073400 # Maximum temperature of polymorph 2 T2 <- 1800 ## ----pyrrhotite_Ttr, message = FALSE------------------------------------------ DGtr <- 0 # DON'T CHANGE THIS TDStr <- DHtr - DGtr # TΔS° = ΔH° - ΔG° DStr <- TDStr / Ttr ## ----pyrrhotite_Cp, results = "hide", message = FALSE------------------------- mod.OBIGT("pyrrhotite_new", formula = formula, state = "cr", ref1 = ref1, E_units = "cal", G = 0, H = 0, S = 0, V = V1, Cp = Cp1, a = a1, b = b1, c = 0, d = 0, e = 0, f = 0, lambda = 0, T = Ttr) mod.OBIGT("pyrrhotite_new", formula = formula, state = "cr2", ref1 = ref1, E_units = "cal", G = 0, H = 0, S = 0, V = V2, a = a2, b = b2, c = c2, d = 0, e = 0, f = 0, lambda = 0, T = T2) ## ----pyrrhotite_info, results = "hide", collapse = TRUE----------------------- info(info("pyrrhotite_new", "cr")) info(info("pyrrhotite_new", "cr2")) ## ----pyrrhotite_S, message = FALSE-------------------------------------------- DS1 <- subcrt("pyrrhotite_new", "cr", P = 1, T = Ttr, use.polymorphs = FALSE)$out[[1]]$S DS2 <- subcrt("pyrrhotite_new", "cr2", P = 1, T = Ttr)$out[[1]]$S DS298 <- DS1 + DStr - DS2 ## ----pyrrhotite_D1_D2, message = FALSE, results = "hide"---------------------- mod.OBIGT("pyrrhotite_new", state = "cr", S = S1) mod.OBIGT("pyrrhotite_new", state = "cr2", S = S1 + DS298) D1 <- subcrt("pyrrhotite_new", "cr", P = 1, T = Ttr, use.polymorphs = FALSE)$out[[1]] D2 <- subcrt("pyrrhotite_new", "cr2", P = 1, T = Ttr)$out[[1]] ## ----pyrrhotite_check_S, results = "hide"------------------------------------- stopifnot(all.equal(D2$S - D1$S, DStr)) ## ----pyrrhotite_DG298_DH298, results = "hide", message = FALSE---------------- DG298 <- D1$G + DGtr - D2$G DH298 <- D1$H + DHtr - D2$H mod.OBIGT("pyrrhotite_new", state = "cr", G = G1, H = H1) mod.OBIGT("pyrrhotite_new", state = "cr2", G = G1 + DG298, H = H1 + DH298) ## ----pyrrhotite_info2, collapse = TRUE---------------------------------------- cr_parameters <- info(info("pyrrhotite_new", "cr")) stopifnot(abs(check.GHS(cr_parameters)) < 1) cr2_parameters <- info(info("pyrrhotite_new", "cr2")) stopifnot(abs(check.GHS(cr2_parameters)) < 1) ## ----pyrrhotite_parameters---------------------------------------------------- cr_parameters cr2_parameters ## ----pyrrhotite_T, eval = FALSE----------------------------------------------- # T <- seq(550, 650, 1) # sout <- subcrt("pyrrhotite_new", T = T, P = 1)$out[[1]] # par(mfrow = c(1, 4), mar = c(4, 4.2, 1, 1)) # labels <- c(G = "DG0f", H = "DH0f", S = "S0", Cp = "Cp0") # for(property in c("G", "H", "S", "Cp")) { # plot(sout$T, sout[, property], col = sout$polymorph, # xlab = axis.label("T"), ylab = axis.label(labels[property]) # ) # abline(v = Ttr, lty = 3, col = 8) # if(property == "G") # legend("bottomleft", c("1", "2"), pch = 1, col = c(1, 2), title = "Polymorph") # } # # # Restore default settings for the questions below # reset() ## ----pyrrhotite_T, echo = FALSE, message = FALSE, results = "hide", fig.width = 8, fig.height = 2.5, out.width = "100%", fig.align = "center", pngquant = pngquant---- T <- seq(550, 650, 1) sout <- subcrt("pyrrhotite_new", T = T, P = 1)$out[[1]] par(mfrow = c(1, 4), mar = c(4, 4.2, 1, 1)) labels <- c(G = "DG0f", H = "DH0f", S = "S0", Cp = "Cp0") for(property in c("G", "H", "S", "Cp")) { plot(sout$T, sout[, property], col = sout$polymorph, xlab = axis.label("T"), ylab = axis.label(labels[property]) ) abline(v = Ttr, lty = 3, col = 8) if(property == "G") legend("bottomleft", c("1", "2"), pch = 1, col = c(1, 2), title = "Polymorph") } # Restore default settings for the questions below reset() ## ----trisulfur, eval = FALSE, echo = FALSE------------------------------------ # par(mfrow = c(1, 3)) # # ## BLOCK 1 # T <- 350 # P <- 5000 # res <- 200 # # ## BLOCK 2 # wt_percent_S <- 10 # wt_permil_S <- wt_percent_S * 10 # molar_mass_S <- mass("S") # 32.06 # moles_S <- wt_permil_S / molar_mass_S # grams_H2O <- 1000 - wt_permil_S # molality_S <- 1000 * moles_S / grams_H2O # logm_S <- log10(molality_S) # # ## BLOCK 3 # basis(c("Ni", "SiO2", "Fe2O3", "H2S", "H2O", "oxygen", "H+")) # species(c("H2S", "HS-", "SO2", "HSO4-", "SO4-2", "S3-")) # a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P) # e <- equilibrate(a, loga.balance = logm_S) # diagram(e) # # ## BLOCK 4 # mod.buffer("NNO", c("nickel", "bunsenite"), state = "cr", logact = 0) # for(buffer in c("HM", "QFM", "NNO")) { # basis("O2", buffer) # logfO2_ <- affinity(return.buffer = TRUE, T = T, P = P)$O2 # abline(h = logfO2_, lty = 3, col = 2) # text(8.5, logfO2_, buffer, adj = c(0, 0), col = 2, cex = 0.9) # } # # ## BLOCK 5 # pH <- subcrt(c("H2O", "H+", "OH-"), c(-1, 1, 1), T = T, P = P)$out$logK / -2 # abline(v = pH, lty = 2, col = 4) # # ## BLOCK 6 # lTP <- describe.property(c("T", "P"), c(T, P)) # lS <- paste(wt_percent_S, "wt% S(aq)") # ltext <- c(lTP, lS) # legend("topright", ltext, bty = "n", bg = "transparent") # title(quote("Parameters for"~S[3]^"-"~"from Pokrovski and Dubessy (2015)"), xpd = NA) # # ######## Plot 2: Modify Gibbs energy # # oldG <- info(info("S3-"))$G # newG <- oldG - 12548 # mod.OBIGT("S3-", G = newG) # a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P) # e <- equilibrate(a, loga.balance = logm_S) # diagram(e) # legend("topright", ltext, bty = "n", bg = "transparent") # title(quote("Modified"~log*italic(K)~"after Pokrovski and Dubrovinsky (2011)"), xpd = NA) # OBIGT() # # ######## Plot 3: Do it with DEW # # T <- 700 # P <- 10000 # 10000 bar = 1 GPa # oldwat <- water("DEW") # add.OBIGT("DEW") # info(species()$ispecies) # a <- affinity(pH = c(2, 10, res), O2 = c(-18, -10, res), T = T, P = P) # e <- equilibrate(a, loga.balance = logm_S) # diagram(e) # lTP <- describe.property(c("T", "P"), c(T, P)) # ltext <- c(lTP, lS) # legend("topright", ltext, bty = "n", bg = "transparent") # title(quote("Deep Earth Water (DEW)"~"model")) # water(oldwat) # OBIGT() ## ----trisulfur, eval = FALSE, echo = 3:43------------------------------------- # par(mfrow = c(1, 3)) # # ## BLOCK 1 # T <- 350 # P <- 5000 # res <- 200 # # ## BLOCK 2 # wt_percent_S <- 10 # wt_permil_S <- wt_percent_S * 10 # molar_mass_S <- mass("S") # 32.06 # moles_S <- wt_permil_S / molar_mass_S # grams_H2O <- 1000 - wt_permil_S # molality_S <- 1000 * moles_S / grams_H2O # logm_S <- log10(molality_S) # # ## BLOCK 3 # basis(c("Ni", "SiO2", "Fe2O3", "H2S", "H2O", "oxygen", "H+")) # species(c("H2S", "HS-", "SO2", "HSO4-", "SO4-2", "S3-")) # a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P) # e <- equilibrate(a, loga.balance = logm_S) # diagram(e) # # ## BLOCK 4 # mod.buffer("NNO", c("nickel", "bunsenite"), state = "cr", logact = 0) # for(buffer in c("HM", "QFM", "NNO")) { # basis("O2", buffer) # logfO2_ <- affinity(return.buffer = TRUE, T = T, P = P)$O2 # abline(h = logfO2_, lty = 3, col = 2) # text(8.5, logfO2_, buffer, adj = c(0, 0), col = 2, cex = 0.9) # } # # ## BLOCK 5 # pH <- subcrt(c("H2O", "H+", "OH-"), c(-1, 1, 1), T = T, P = P)$out$logK / -2 # abline(v = pH, lty = 2, col = 4) # # ## BLOCK 6 # lTP <- describe.property(c("T", "P"), c(T, P)) # lS <- paste(wt_percent_S, "wt% S(aq)") # ltext <- c(lTP, lS) # legend("topright", ltext, bty = "n", bg = "transparent") # title(quote("Parameters for"~S[3]^"-"~"from Pokrovski and Dubessy (2015)"), xpd = NA) # # ######## Plot 2: Modify Gibbs energy # # oldG <- info(info("S3-"))$G # newG <- oldG - 12548 # mod.OBIGT("S3-", G = newG) # a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P) # e <- equilibrate(a, loga.balance = logm_S) # diagram(e) # legend("topright", ltext, bty = "n", bg = "transparent") # title(quote("Modified"~log*italic(K)~"after Pokrovski and Dubrovinsky (2011)"), xpd = NA) # OBIGT() # # ######## Plot 3: Do it with DEW # # T <- 700 # P <- 10000 # 10000 bar = 1 GPa # oldwat <- water("DEW") # add.OBIGT("DEW") # info(species()$ispecies) # a <- affinity(pH = c(2, 10, res), O2 = c(-18, -10, res), T = T, P = P) # e <- equilibrate(a, loga.balance = logm_S) # diagram(e) # lTP <- describe.property(c("T", "P"), c(T, P)) # ltext <- c(lTP, lS) # legend("topright", ltext, bty = "n", bg = "transparent") # title(quote("Deep Earth Water (DEW)"~"model")) # water(oldwat) # OBIGT() ## ----trisulfur_logK, message = FALSE, echo = 1:3------------------------------ species <- c("H2S", "SO4-2", "H+", "S3-", "oxygen", "H2O") coeffs <- c(-2, -1, -1, 1, 0.75, 2.5) (calclogK <- subcrt(species, coeffs, T = seq(300, 450, 50), P = 5000)$out$logK) fcalclogK <- formatC(calclogK, format = "f", digits = 1) reflogK <- -9.6 dlogK <- calclogK[2] - reflogK # Put in a test so that we don't get surprised by # future database updates or changes to this vignette stopifnot(round(dlogK, 4) == -4.4132) ## ----trisulfur, eval = FALSE, echo = 46:55------------------------------------ # par(mfrow = c(1, 3)) # # ## BLOCK 1 # T <- 350 # P <- 5000 # res <- 200 # # ## BLOCK 2 # wt_percent_S <- 10 # wt_permil_S <- wt_percent_S * 10 # molar_mass_S <- mass("S") # 32.06 # moles_S <- wt_permil_S / molar_mass_S # grams_H2O <- 1000 - wt_permil_S # molality_S <- 1000 * moles_S / grams_H2O # logm_S <- log10(molality_S) # # ## BLOCK 3 # basis(c("Ni", "SiO2", "Fe2O3", "H2S", "H2O", "oxygen", "H+")) # species(c("H2S", "HS-", "SO2", "HSO4-", "SO4-2", "S3-")) # a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P) # e <- equilibrate(a, loga.balance = logm_S) # diagram(e) # # ## BLOCK 4 # mod.buffer("NNO", c("nickel", "bunsenite"), state = "cr", logact = 0) # for(buffer in c("HM", "QFM", "NNO")) { # basis("O2", buffer) # logfO2_ <- affinity(return.buffer = TRUE, T = T, P = P)$O2 # abline(h = logfO2_, lty = 3, col = 2) # text(8.5, logfO2_, buffer, adj = c(0, 0), col = 2, cex = 0.9) # } # # ## BLOCK 5 # pH <- subcrt(c("H2O", "H+", "OH-"), c(-1, 1, 1), T = T, P = P)$out$logK / -2 # abline(v = pH, lty = 2, col = 4) # # ## BLOCK 6 # lTP <- describe.property(c("T", "P"), c(T, P)) # lS <- paste(wt_percent_S, "wt% S(aq)") # ltext <- c(lTP, lS) # legend("topright", ltext, bty = "n", bg = "transparent") # title(quote("Parameters for"~S[3]^"-"~"from Pokrovski and Dubessy (2015)"), xpd = NA) # # ######## Plot 2: Modify Gibbs energy # # oldG <- info(info("S3-"))$G # newG <- oldG - 12548 # mod.OBIGT("S3-", G = newG) # a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P) # e <- equilibrate(a, loga.balance = logm_S) # diagram(e) # legend("topright", ltext, bty = "n", bg = "transparent") # title(quote("Modified"~log*italic(K)~"after Pokrovski and Dubrovinsky (2011)"), xpd = NA) # OBIGT() # # ######## Plot 3: Do it with DEW # # T <- 700 # P <- 10000 # 10000 bar = 1 GPa # oldwat <- water("DEW") # add.OBIGT("DEW") # info(species()$ispecies) # a <- affinity(pH = c(2, 10, res), O2 = c(-18, -10, res), T = T, P = P) # e <- equilibrate(a, loga.balance = logm_S) # diagram(e) # lTP <- describe.property(c("T", "P"), c(T, P)) # ltext <- c(lTP, lS) # legend("topright", ltext, bty = "n", bg = "transparent") # title(quote("Deep Earth Water (DEW)"~"model")) # water(oldwat) # OBIGT() ## ----trisulfur_DEW, message = FALSE, echo = 8:22, collapse = TRUE, fig.keep = "none"---- # The first four lines are stand-ins to make this block runnable and get the right output from info(); # the diagram actually shown in the vignette is made using the 'trisulfur' block above b <- basis("CHNOS+") s <- species(c("H2S", "HS-", "SO2", "HSO4-", "SO4-2", "S3-")) res <- 10 wt_percent_S <- logm_S <- 0 ######## End of stand-in code ######## T <- 700 P <- 10000 # 10000 bar = 1 GPa oldwat <- water("DEW") add.OBIGT("DEW") info(species()$ispecies)[, 1:8] a <- affinity(pH = c(2, 10, res), O2 = c(-18, -10, res), T = T, P = P) e <- equilibrate(a, loga.balance = logm_S) diagram(e) lTP <- describe.property(c("T", "P"), c(T, P)) lS <- paste(wt_percent_S, "wt% S(aq)") ltext <- c(lTP, lS) legend("topright", ltext, bty = "n", bg = "transparent") title(quote("Deep Earth Water (DEW)"~"model")) water(oldwat) OBIGT() ## ----trisulfur, echo = FALSE, message = FALSE, results = "hide", fig.width = 10, fig.height = 3.33, out.width = "100%", out.extra='class="full-width"', pngquant = pngquant, cache = TRUE---- par(mfrow = c(1, 3)) ## BLOCK 1 T <- 350 P <- 5000 res <- 200 ## BLOCK 2 wt_percent_S <- 10 wt_permil_S <- wt_percent_S * 10 molar_mass_S <- mass("S") # 32.06 moles_S <- wt_permil_S / molar_mass_S grams_H2O <- 1000 - wt_permil_S molality_S <- 1000 * moles_S / grams_H2O logm_S <- log10(molality_S) ## BLOCK 3 basis(c("Ni", "SiO2", "Fe2O3", "H2S", "H2O", "oxygen", "H+")) species(c("H2S", "HS-", "SO2", "HSO4-", "SO4-2", "S3-")) a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P) e <- equilibrate(a, loga.balance = logm_S) diagram(e) ## BLOCK 4 mod.buffer("NNO", c("nickel", "bunsenite"), state = "cr", logact = 0) for(buffer in c("HM", "QFM", "NNO")) { basis("O2", buffer) logfO2_ <- affinity(return.buffer = TRUE, T = T, P = P)$O2 abline(h = logfO2_, lty = 3, col = 2) text(8.5, logfO2_, buffer, adj = c(0, 0), col = 2, cex = 0.9) } ## BLOCK 5 pH <- subcrt(c("H2O", "H+", "OH-"), c(-1, 1, 1), T = T, P = P)$out$logK / -2 abline(v = pH, lty = 2, col = 4) ## BLOCK 6 lTP <- describe.property(c("T", "P"), c(T, P)) lS <- paste(wt_percent_S, "wt% S(aq)") ltext <- c(lTP, lS) legend("topright", ltext, bty = "n", bg = "transparent") title(quote("Parameters for"~S[3]^"-"~"from Pokrovski and Dubessy (2015)"), xpd = NA) ######## Plot 2: Modify Gibbs energy oldG <- info(info("S3-"))$G newG <- oldG - 12548 mod.OBIGT("S3-", G = newG) a <- affinity(pH = c(2, 10, res), O2 = c(-34, -22, res), T = T, P = P) e <- equilibrate(a, loga.balance = logm_S) diagram(e) legend("topright", ltext, bty = "n", bg = "transparent") title(quote("Modified"~log*italic(K)~"after Pokrovski and Dubrovinsky (2011)"), xpd = NA) OBIGT() ######## Plot 3: Do it with DEW T <- 700 P <- 10000 # 10000 bar = 1 GPa oldwat <- water("DEW") add.OBIGT("DEW") info(species()$ispecies) a <- affinity(pH = c(2, 10, res), O2 = c(-18, -10, res), T = T, P = P) e <- equilibrate(a, loga.balance = logm_S) diagram(e) lTP <- describe.property(c("T", "P"), c(T, P)) ltext <- c(lTP, lS) legend("topright", ltext, bty = "n", bg = "transparent") title(quote("Deep Earth Water (DEW)"~"model")) water(oldwat) OBIGT() ## ----pyrrhotite_polymorphs, collapse = TRUE----------------------------------- subcrt("pyrrhotite", T = c(25, 150, 350), property = "G")$out ## ----pyrite_limit------------------------------------------------------------- subcrt("pyrite", T = seq(200, 1000, 200), P = 1) ## ----mineral_Ttr, collapse = TRUE--------------------------------------------- file <- system.file("extdata/OBIGT/inorganic_cr.csv", package = "CHNOSZ") dat <- read.csv(file) dat$name[dat$model == "CGL_Ttr"] ## ----muscovite_limit, message = FALSE----------------------------------------- add.OBIGT("SUPCRT92") subcrt("muscovite", T = 850, P = 4500) reset() ## ----KMQ_basis_species, message = FALSE--------------------------------------- basis(c("Al+3", "quartz", "K+", "H+", "H2O", "oxygen")) species(c("kaolinite", "muscovite", "K-feldspar")) ## ----KMQ_m_K, message = FALSE------------------------------------------------- # Define temperature, pressure, and molality of Cl- (==IS) T <- 150 P <- 500 IS <- m_Cl <- 1 # Calculate equilibrium constant for Ab-Kfs reaction, corrected for ionic strength logK_AK <- subcrt(c("albite", "K+", "K-feldspar", "Na+"), c(-1, -1, 1, 1), T = T, P = P, IS = IS)$out$logK K_AK <- 10 ^ logK_AK # Calculate molality of K+ (m_K <- m_Cl / (K_AK + 1)) ## ----KMQ_diagram, eval = FALSE, echo = 2:10----------------------------------- # par(mfrow = c(1, 2)) # basis("K+", log10(m_K)) # a <- affinity(pH = c(2, 10), O2 = c(-55, -38), T = T, P = P, IS = IS) # diagram(a, srt = 90) # lTP <- as.expression(c(lT(T), lP(P))) # legend("topleft", legend = lTP, bty = "n", inset = c(-0.05, 0), cex = 0.9) # ltxt <- c(quote("Unit molality of Cl"^"-"), "Quartz saturation") # legend("topright", legend = ltxt, bty = "n", cex = 0.9) # title("Mineral data from Berman (1988)\nand Sverjensky et al. (1991) (OBIGT default)", # font.main = 1, cex.main = 0.9) # add.OBIGT("SUPCRT92") # a <- affinity(pH = c(2, 10), O2 = c(-55, -38), T = T, P = P, IS = IS) # diagram(a, srt = 90) # title("Mineral data from Helgeson et al. (1978)\n(as used in SUPCRT92)", # font.main = 1, cex.main = 0.9) # OBIGT() ## ----KMQ_diagram, message = FALSE, fig.width = 8, fig.height = 4, out.width = "100%", results = "hide", echo = FALSE---- par(mfrow = c(1, 2)) basis("K+", log10(m_K)) a <- affinity(pH = c(2, 10), O2 = c(-55, -38), T = T, P = P, IS = IS) diagram(a, srt = 90) lTP <- as.expression(c(lT(T), lP(P))) legend("topleft", legend = lTP, bty = "n", inset = c(-0.05, 0), cex = 0.9) ltxt <- c(quote("Unit molality of Cl"^"-"), "Quartz saturation") legend("topright", legend = ltxt, bty = "n", cex = 0.9) title("Mineral data from Berman (1988)\nand Sverjensky et al. (1991) (OBIGT default)", font.main = 1, cex.main = 0.9) add.OBIGT("SUPCRT92") a <- affinity(pH = c(2, 10), O2 = c(-55, -38), T = T, P = P, IS = IS) diagram(a, srt = 90) title("Mineral data from Helgeson et al. (1978)\n(as used in SUPCRT92)", font.main = 1, cex.main = 0.9) OBIGT() ## ----KMQ_refs, message = FALSE------------------------------------------------ thermo.refs(species()$ispecies) ## ----KMQ_diagram, eval = FALSE, echo = 11:15---------------------------------- # par(mfrow = c(1, 2)) # basis("K+", log10(m_K)) # a <- affinity(pH = c(2, 10), O2 = c(-55, -38), T = T, P = P, IS = IS) # diagram(a, srt = 90) # lTP <- as.expression(c(lT(T), lP(P))) # legend("topleft", legend = lTP, bty = "n", inset = c(-0.05, 0), cex = 0.9) # ltxt <- c(quote("Unit molality of Cl"^"-"), "Quartz saturation") # legend("topright", legend = ltxt, bty = "n", cex = 0.9) # title("Mineral data from Berman (1988)\nand Sverjensky et al. (1991) (OBIGT default)", # font.main = 1, cex.main = 0.9) # add.OBIGT("SUPCRT92") # a <- affinity(pH = c(2, 10), O2 = c(-55, -38), T = T, P = P, IS = IS) # diagram(a, srt = 90) # title("Mineral data from Helgeson et al. (1978)\n(as used in SUPCRT92)", # font.main = 1, cex.main = 0.9) # OBIGT() ## ----PyPo_logfO2, message = FALSE, results = "hide"--------------------------- # Constants a_S <- 0.001 T <- 150 # Get the pKa of H2S (note the minus sign!) pKa_H2S <- - subcrt(c("H2S", "HS-", "H+"), c(-1, 1, 1), T = T)$out$logK # Reaction 1 between pyrite and pyrrhotite with H2S basis(c("pyrite", "H2S", "oxygen", "H2O", "H+")) PyPo_H2S <- subcrt("pyrrhotite", 1, T = T) # Reaction 1 law of mass action (LMA) # logf_O2 = 2 logK_1 - 2 loga_H2S logf_O2_1 <- function(loga_H2S) 2 * PyPo_H2S$out$logK - 2 * loga_H2S # Reaction 2 between pyrite and pyrrhotite with HS- basis(c("pyrite", "HS-", "oxygen", "H2O", "H+")) PyPo_HS <- subcrt("pyrrhotite", 1, T = T) # Reaction 2 LMA # logf_O2 = 2 pH + 2 logK_2 - 2 loga_HS- logf_O2_2 <- function(pH, loga_HS) 2 * pH + 2 * PyPo_HS$out$logK - 2 * loga_HS # The logf_O2 for each reaction is the same at the pKa of H2S stopifnot(all.equal(logf_O2_1(log10(a_S)), logf_O2_2(pKa_H2S, log10(a_S)))) stopifnot(all.equal(logf_O2_1(log10(a_S / 2)), logf_O2_2(pKa_H2S, log10(a_S / 2)))) ## ----PyPo_plot, eval = FALSE-------------------------------------------------- # par(mfrow = c(1, 2)) # # # Let's draw this out # thermo.plot.new(c(2, 10), c(-56, -46), xlab = "pH", ylab = axis.label("O2")) # lines(c(2, pKa_H2S), c(logf_O2_1(log10(a_S)), logf_O2_1(log10(a_S))), col = 2) # lines(c(pKa_H2S, 10), c(logf_O2_2(pKa_H2S, log10(a_S)), logf_O2_2(10, log10(a_S))), col = 4) # abline(v = pKa_H2S, lty = 2, col = 8) # text(6, -46.5, expr.species("H2S"), adj = c(0.5, 1)) # text(7, -46.5, expr.species("HS-"), adj = c(0.5, 1)) # text(4, -54, "pyrite") # text(4, -55.5, "pyrrhotite") # points(pKa_H2S, logf_O2_1(log10(a_S)), pch = 0) # points(pKa_H2S, logf_O2_1(log10(a_S / 2)), pch = 19) # legend("bottomright", c("Yes", "No"), pch = c(19, 0), title = as.expression(quote(sum(S) == "constant")), bty = "n") # legend("topleft", legend = lT(T), bty = "n") # title(quote("Straight lines for" ~ italic(a)[H[2]*S] ~ "=" ~ italic(a)[HS^"-"] == "0.001"), font.main = 1) # # # Do it with mosaic() # basis(c("pyrite", "H2S", "oxygen", "H2O", "H+")) # # This sets the activity of H2S, which is used for total S when mosaic is called with blend = TRUE # basis("H2S", -3) # # This defines the basis species to swap through # bases <- c("H2S", "HS-") # species(c("pyrite", "pyrrhotite")) # m <- mosaic(bases, pH = c(2, 10, 500), O2 = c(-56, -46, 500), T = T) # diagram(m$A.species, dx = c(-1, -3), dy = c(-3, -1.5), col = 6) # abline(v = pKa_H2S, lty = 2, col = 8) # text(6, -46.5, expr.species("H2S"), adj = c(0.5, 1)) # text(7, -46.5, expr.species("HS-"), adj = c(0.5, 1)) # points(pKa_H2S, logf_O2_1(log10(a_S)), pch = 0) # points(pKa_H2S, logf_O2_1(log10(a_S / 2)), pch = 19) # legend("bottomright", c("Yes", "No"), pch = c(19, 0), title = as.expression(quote(sum(S) == "constant")), bty = "n") # legend("topleft", legend = lT(T), bty = "n") # title(quote("Curved line for" ~ italic(a)[H[2]*S] + italic(a)[HS^"-"] == "0.001"), font.main = 1) ## ----PyPo_plot, echo = FALSE, message = FALSE, results = "hide", fig.width = 9, fig.height = 4.5, out.width = "100%", fig.align = "center", pngquant = pngquant, cache = TRUE---- par(mfrow = c(1, 2)) # Let's draw this out thermo.plot.new(c(2, 10), c(-56, -46), xlab = "pH", ylab = axis.label("O2")) lines(c(2, pKa_H2S), c(logf_O2_1(log10(a_S)), logf_O2_1(log10(a_S))), col = 2) lines(c(pKa_H2S, 10), c(logf_O2_2(pKa_H2S, log10(a_S)), logf_O2_2(10, log10(a_S))), col = 4) abline(v = pKa_H2S, lty = 2, col = 8) text(6, -46.5, expr.species("H2S"), adj = c(0.5, 1)) text(7, -46.5, expr.species("HS-"), adj = c(0.5, 1)) text(4, -54, "pyrite") text(4, -55.5, "pyrrhotite") points(pKa_H2S, logf_O2_1(log10(a_S)), pch = 0) points(pKa_H2S, logf_O2_1(log10(a_S / 2)), pch = 19) legend("bottomright", c("Yes", "No"), pch = c(19, 0), title = as.expression(quote(sum(S) == "constant")), bty = "n") legend("topleft", legend = lT(T), bty = "n") title(quote("Straight lines for" ~ italic(a)[H[2]*S] ~ "=" ~ italic(a)[HS^"-"] == "0.001"), font.main = 1) # Do it with mosaic() basis(c("pyrite", "H2S", "oxygen", "H2O", "H+")) # This sets the activity of H2S, which is used for total S when mosaic is called with blend = TRUE basis("H2S", -3) # This defines the basis species to swap through bases <- c("H2S", "HS-") species(c("pyrite", "pyrrhotite")) m <- mosaic(bases, pH = c(2, 10, 500), O2 = c(-56, -46, 500), T = T) diagram(m$A.species, dx = c(-1, -3), dy = c(-3, -1.5), col = 6) abline(v = pKa_H2S, lty = 2, col = 8) text(6, -46.5, expr.species("H2S"), adj = c(0.5, 1)) text(7, -46.5, expr.species("HS-"), adj = c(0.5, 1)) points(pKa_H2S, logf_O2_1(log10(a_S)), pch = 0) points(pKa_H2S, logf_O2_1(log10(a_S / 2)), pch = 19) legend("bottomright", c("Yes", "No"), pch = c(19, 0), title = as.expression(quote(sum(S) == "constant")), bty = "n") legend("topleft", legend = lT(T), bty = "n") title(quote("Curved line for" ~ italic(a)[H[2]*S] + italic(a)[HS^"-"] == "0.001"), font.main = 1) ## ----Fe-S-O-C, message = FALSE, results = "hide", fig.width = 5, fig.height = 5, out.width = "60%", fig.align = "center", pngquant = pngquant, cache = TRUE---- basis(c("FeO", "SO4-2", "CO3-2", "H2O", "H+", "oxygen")) basis("SO4-2", -3) basis("CO3-2", -0.6) species(c("hematite", "pyrite", "pyrrhotite", "magnetite", "siderite")) bases <- list( c("SO4-2", "HSO4-", "HS-", "H2S"), c("CO3-2", "HCO3-", "CO2") ) m <- mosaic(bases, pH = c(0, 14, 500), O2 = c(-57, -35, 500), T = 150, IS = 0.77) d <- diagram(m$A.species, fill = "terrain", dx = c(0, 0, 0, 0, 0.3)) water.lines(d)