In paleolimbot/easyphreeqc: Model Reactions Using USGS Phreeqc

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "README-"
)

easyphreeqc

The goal of easyphreeqc is to provide a more useful interface to the existing phreeqc package.

Installation

You can install easyphreeqc from github with:

# install.packages("devtools")
devtools::install_github("paleolimbot/easyphreeqc")

Example

Running PHREEQC is accomplished using the phreeqc() function, which calls the program and generates the output. The function accepts character vectors of input, which can be generated using intput helper functions such as solution(), selected_output(), equilibrium_phases(), and reaction_temperature() (or roll your own input using phreeqc_input()). To get the raw output of PHREEQC, you will have to pass quiet = FALSE.

library(easyphreeqc)
phreeqc(
  solution(pH = 7, temp = 25),
  quiet = FALSE
)

To get the results as a data frame, we need to supply a selected_output() to the input file.

phreeqc(
  solution(pH = 7, temp = 25),
  selected_output(pH = TRUE, temp = TRUE, activities = c("OH-", "H+", "O2"))
)

To find the distribution of a few solutions, you can add solutions to the input (you will have to number them).

phreeqc(
  solution(1, pH = 6, temp = 25),
  solution(2, pH = 7, temp = 25),
  solution(3, pH = 8, temp = 25),
  selected_output(pH = TRUE, temp = TRUE, activities = c("OH-", "H+", "O2"))
)

Databases

Some elements (for example, mercury) aren't included in the base database. There are a number of databases included in the PHREEQC package, that you can choose by specifying the db argument of phreeqc(). One that includes mercury is the "minteq" database.

phreeqc(
  solution(1, pH = 7, temp = 25, Hg = 0.1),
  solution(2, pH = 7, temp = 25, Hg = 0.2),
  solution(3, pH = 7, temp = 25, Hg = 0.3),
  selected_output(activities = c("Hg", "Hg2+2", "Hg(OH)2", "Hg(OH)2",
                                 "HgOH+", "Hg(OH)3-")),
  db = "minteq"
)

(Note that I called phreeqc() with quiet = FALSE to first find the species included in selected_output())

Multiple solutions

The input for phreeqc() can accept list() objects, so you can use something like plyr::mlply() to create a list of solutions like the one above to pass to phreeqc(). You could use expand.grid() to generate a large number of solutions in this way.

# create data.frame with arguments to solution()
solution_info <- data.frame(
  .number = 1:3,
  pH = 7,
  temp = 25,
  Hg = c(0.1, 0.2, 0.3)
)

# use mlply to turn this into a list()
solution_list <- plyr::mlply(solution_info, solution)

# pass to phreeqc()
phreeqc(
  solution_list,
  selected_output(activities = c("Hg", "Hg2+2", "Hg(OH)2", "Hg(OH)2",
                                 "HgOH+", "Hg(OH)3-")),
  db = "minteq"
)