README.Rmd

---
output: github_document
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "man/figures/README-",
  dpi = 600,
  dev = "png",
  out.width = "100%"
)
```

# pbcusol

<!-- badges: start -->
[![R-CMD-check](https://github.com/bentrueman/pbcusol/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/bentrueman/pbcusol/actions/workflows/R-CMD-check.yaml)
[![Codecov test coverage](https://codecov.io/gh/bentrueman/pbcusol/graph/badge.svg)](https://app.codecov.io/gh/bentrueman/pbcusol)
<!-- badges: end -->

`pbcusol` predicts equilibrium lead and copper solubility using the US EPA databases LEADSOL (Schock et al. 1996) and CU2SOL (Schock et al., 1995), along with those available as a part of PHREEQC (Charlton and Parkhurst, 2011; Parkhurst and Appelo, 2013). `pbcusol` uses `tidyphreeqc` (Dunnington, 2019)---a convenient interface for PHREEQC in R---for most solubility computations. It's also available as a [Shiny app](https://bentrueman.shinyapps.io/shinypbcusol/).

## Installation

You can install the development version from [GitHub](https://github.com/) with:

``` r
# install.packages("remotes")
remotes::install_github("bentrueman/pbcusol")
```
## Example

For this example, you will need the `tidyverse` family of packages, along with `plyr::round_any()`.

```{r load, message=FALSE}
library("pbcusol")
library("tidyverse")
library("plyr", include.only = "round_any")
```

```{r theme, echo=FALSE}
theme_set(theme_classic() + theme(legend.position = "right", strip.background = element_blank()))
```

Use `eq_sol()` to calculate the equlibrium solubility of multiple copper phases that occur in drinking water systems, over a wide range of pH values and dissolved inorganic carbon concentrations. (N.B., evaluate on a smaller grid to speed this up!)

```{r cu-sol}
dic_increment_cu <- 1.5
solutions_cu <- list("Cu(OH)2", "Tenorite", "Malachite") %>%
  set_names() %>%
  map_dfr(
    ~eq_sol(
      element = "Cu",
      ph = seq(6, 11, by = .025),
      dic = seq(1, 150, by = dic_increment_cu),
      phase = .x
    ),
    .id = "phase"
  )
```

Plot the data (n.b., plots shown here will differ slightly from those the following code generates). Lytle et al. (2018) have noted that copper solubility predictions are not reliable in the presence of orthophosphate, and they propose an empirical model as an alternative. It is implemented in `eq_sol()` via the argument `empirical = TRUE`.

```{r cu-plot-display, eval=FALSE}
solutions_cu %>% 
  mutate(
    log10_cu_ppb = log10(cu_ppb),
    dic_ppm = plyr::round_any(dic_ppm, dic_increment_cu)
  ) %>% 
  ggplot(aes(x = dic_ppm, y = pH)) + 
  facet_wrap(vars(phase)) + 
  geom_raster(aes(fill = log10_cu_ppb)) +
  geom_contour(aes(z = log10_cu_ppb), col = "white") +
  scale_fill_viridis_c(option = "magma")
```

```{r cu-plot, fig.height=3, echo=FALSE}
solutions_cu %>% 
  mutate(
    phase = fct_recode(phase, "Cu(OH)[2]" = "Cu(OH)2"),
    cu_ppb = log10(cu_ppb),
    dic_ppm = plyr::round_any(dic_ppm, dic_increment_cu)
  ) %>% 
  ggplot(aes(x = dic_ppm, y = pH)) + 
  facet_wrap(vars(phase), labeller = label_parsed) + 
  geom_raster(aes(fill = cu_ppb)) +
  geom_contour(aes(z = cu_ppb), col = "white", linetype = 2, linewidth = .2, breaks = seq(-4, 5, by = .1)) + 
  geom_contour(aes(z = cu_ppb), col = "white", breaks = -3:5) + 
  scale_fill_viridis_c(option = "magma", labels = function(breaks) 1e-3 * (10 ^ breaks)) +
  labs(
    x = expression("[C]"[inorganic]~"(mg L"^"-1"*")"),
    y = "pH",
    fill = expression("[Cu] (mg L"^"-1"*")")
  )
```

Use `eq_sol()` to make the same predictions for lead. The helper function `calculate_dic()` may be useful for converting alkalinity to dissolved inorganic carbon. 

```{r pb-sol}
dic_increment_pb <- .8
solutions_pb <- list("Cerussite", "Hydcerussite", "Hxypyromorphite") %>%
  set_names() %>%
  map_dfr(
    ~ eq_sol(
      element = "Pb",
      ph = seq(6, 11, by = .025),
      dic = seq(1, 80, by = dic_increment_pb),
      phosphate = .16,
      phase = .x
    ),
    .id = "phase"
  )
```

Plot the data:

```{r pb-plot-display, eval=FALSE}
solutions_pb %>% 
  mutate(
    log10_pb_ppb = log10(pb_ppb),
    dic_ppm = plyr::round_any(dic_ppm, dic_increment_pb)
  ) %>% 
  ggplot(aes(x = dic_ppm, y = pH)) + 
  facet_wrap(vars(phase)) + 
  geom_raster(aes(fill = log10_pb_ppb)) +
  geom_contour(aes(z = log10_pb_ppb), col = "white") +
  scale_fill_viridis_c(option = "magma")
```

```{r pb-plot, fig.height=3, echo=FALSE}
solutions_pb %>% 
  mutate(
    phase = fct_recode(phase, "Hydroxylpyromorphite" = "Hxypyromorphite", "Hydrocerussite" = "Hydcerussite"),
    pb_ppb = log10(pb_ppb),
    dic_ppm = plyr::round_any(dic_ppm, dic_increment_pb)
  ) %>% 
  ggplot(aes(x = dic_ppm, y = pH)) + 
  facet_wrap(vars(phase)) + 
  geom_raster(aes(fill = pb_ppb)) +
  geom_contour(aes(z = pb_ppb), col = "white", linetype = 2, linewidth = .2, breaks = seq(-4, 5, by = .1)) + 
  geom_contour(aes(z = pb_ppb), col = "white", breaks = -3:5) + 
  scale_fill_viridis_c(option = "magma", labels = function(breaks) 1e-3 * (10 ^ breaks)) +
  labs(
    x = expression("[C]"[inorganic]~"(mg L"^"-1"*")"),
    y = "pH",
    fill = expression("[Pb] (mg L"^"-1"*")")
  )
```

Combining the information in the previous two plots, we can generate prediction surfaces for equilibrium lead and copper solubility that are quite close to those presented in literature, specifically Figure 7 in Schock et al. (1995) and Figure 4-18 in Schock et al. (1996).

```{r combined, message=FALSE, fig.height=3, echo=FALSE}
list(
  "[Pb]" = rename(solutions_pb, conc_ppb = pb_ppb), 
  "[Cu]" = rename(solutions_cu, conc_ppb = cu_ppb)
) %>% 
  bind_rows(.id = "element") %>% 
  mutate(
    dic_ppm = if_else(
      element == "[Pb]", 
      plyr::round_any(dic_ppm, dic_increment_pb), 
      plyr::round_any(dic_ppm, dic_increment_cu)
    )
  ) %>% 
  group_by(element, pH, dic_ppm) %>% 
  summarize(
    phase = phase[which.min(conc_ppb)],
    conc_ppb = min(conc_ppb) %>% log10()
  ) %>% 
  ungroup() %>% 
  ggplot(aes(x = dic_ppm, y = pH)) + 
  facet_wrap(vars(element), scales = "free_x") + 
  geom_point(aes(col = conc_ppb), shape = 15, size = 1) +
  geom_contour(
    aes(z = conc_ppb), col = "white", linetype = 2, 
    linewidth = .2, breaks = seq(-4, 3, by = .1)
  ) + 
  geom_contour(
    data = function(x) x %>% 
      filter(element == "[Cu]"),
    aes(z = conc_ppb), col = "white", breaks = -3:3
  ) + 
  geom_smooth(
    data = function(x) x %>% 
      filter(phase == "Tenorite") %>% 
      group_by(element, dic_ppm) %>% 
      summarize(min = min(pH)),
    aes(x = dic_ppm, y = min),
    linetype = 2, linewidth = 1, col = "white", se = FALSE, span = .4
  ) + 
  geom_smooth(
    data = function(x) x %>% 
      filter(phase == "Hydcerussite") %>% 
      group_by(element, dic_ppm) %>% 
      summarize(min = min(pH)),
    aes(x = dic_ppm, y = min),
    linetype = 2, linewidth = 1, col = "white", se = FALSE, span = .4
  ) +
  geom_smooth(
    data = function(x) x %>% 
      filter(phase == "Hxypyromorphite") %>% 
      group_by(element, dic_ppm) %>% 
      summarize(max = max(pH)) %>% 
      filter(dic_ppm > 35),
    aes(x = dic_ppm, y = max),
    linetype = 2, linewidth = 1, col = "white", se = FALSE, span = .4
  ) +
  scale_color_viridis_c(option = "magma", labels = function(breaks) 1e-3 * (10 ^ breaks)) +
  labs(
    x = expression("[C]"[inorganic]~"(mg L"^"-1"*")"),
    y = "pH",
    col = expression("[element] (mg L"^"-1"*")")
  ) +
  geom_label(
    data = tibble(
      element = c(rep("[Pb]", 3), rep("[Cu]", 2)),
      label = c("PbCO[3]", "Pb[3]*'('*CO[3]*')'[2]*'('*OH*')'[2]", "Pb[5]*'('*PO[4]*')'[3]*OH", 
                "Cu[2]*'('*CO[3]*')'*'('*OH*')'[2]", "CuO"),
      x = c(65, 40, 20, 90, 50),
      y = c(7.6, 10, 7.6, 6.5, 10.5)
    ),
    aes(x = x, y = y, label = label),
    parse = TRUE,
    label.r = unit(0, "cm"), label.size = unit(0, "cm")
  )
```

Use `pb_logk()` and `cu_logk()` to generate tables of the relevant reactions and constants in the LEADSOL and CU2SOL databases. N.B., `pbcusol` uses a modified version of `phreeqc::minteq.dat` that include the data from LEADSOL and CU2SOL. N.B., `phreeqc::minteq.dat` includes two reactions describing complexation of copper and chloride that are not included in Schock et al. (1995).

### Surface complexation (experimental)  

Metal binding to natural organic matter can be modeled using `eq_sol_wham()`, an approximation of the Windermere Humic Aqueous Model (WHAM) (Tipping and Hurley, 1992), as described in Example 19 of Parkhurst and Appelo (2013).

```{r wham}
eq_sol_wham(
  element = "Pb", 
  ph = 7.5, 
  dic = 50, 
  phase = "Cerussite", 
  Na = 4, 
  mass_ha = 3.5e-3
) %>% 
  transmute(
    phase, pH, dic_ppm, 
    solution_pb = pb_ppb, 
    total_pb = mol_Cerussite * 1e6 * 207.21
  )
```

Metal binding to colloidal ferrihydrite can be modeled using `eq_sol_fixed()` (or `eq_sol_wham()`), as described in Example 8 of Parkhurst and Appelo (2013).

```{r ferrihydrite}
# from Example 8:
hfo_surface_area <- 600                                   # m^2 / g
strong_density <- 5e-6 / (hfo_surface_area * .09)         # site density in mol / m^2 
weak_density <- 2e-4 / (hfo_surface_area * .09)           # site density in mol / m^2 
hfo_mass <- 1.7e-4                                        # grams of colloidal ferrihydrite (1e-4 g Fe/L)
hfo_s <- hfo_mass * hfo_surface_area * strong_density     # number of strong binding sites
hfo_w <- hfo_mass * hfo_surface_area * weak_density       # number of weak binding sites

# define surface:
fer_surf <- list(
    Hfo_sOH = c(hfo_s, hfo_surface_area, hfo_mass),
    Hfo_wOH = hfo_w,
    "-equilibrate" = 1,
    "-Donnan"
    )

eq_sol_fixed(
  element = "Pb", 
  ph = 7.5, 
  dic = 5, 
  phosphate = .3, 
  phase = "Hxypyromorphite", 
  surface_components = fer_surf
) %>% 
  transmute(
    phase, pH, dic_ppm, p_ppm, 
    solution_pb = pb_ppb, 
    total_pb = mol_Hxypyromorphite * 1e6 * 207.21 * 5
  )
```


# References 

Charlton, S.R., and D. L. Parkhurst. 2011. Modules based on the geochemical model PHREEQC for use in scripting and programming languages. Computers & Geosciences, v. 37, p. 1653-1663.

Dunnington, D. 2019. tidyphreeqc: Tidy Geochemical Modeling Using PHREEQC. https://github.com/paleolimbot/tidyphreeqc.

Lytle, D. A., M. R. Schock, J. Leo, and  B. Barnes. 2018. A model for estimating the impact of orthophosphate on copper in water. Journal AWWA. 110: E1-E15. https://doi-org.ezproxy.library.dal.ca/10.1002/awwa.1109

Parkhurst, D. L., and C. A. J. Appelo. 2013.  Description of input and examples for PHREEQC version 3--A computer program
for speciation, batch- reaction, one-dimensional transport, and inverse geochemical calculations: U.S. Geological
Survey Techniques and Methods, book 6, chap. A43, 497 p. http://pubs.usgs.gov/tm/06/a43.

Schock, M. R., D. A Lytle, and J. A. Clement. 1995. “Effect of pH, DIC, orthophosphate and sulfate on drinking water cuprosolvency.” National Risk Management Research Lab., Cincinnati, OH (United States).

Schock, M. R., I. Wagner, and R. J. Oliphant. 1996. “Corrosion and solubility of lead in drinking water.” In Internal corrosion of water distribution systems, 2nd ed., p. 131–230. Denver, CO: American Water Works Association Research Foundation.

Tipping, E., and  M. A. Hurley. 1992. A unifying model of cation binding by humic substances. Geochimica et Cosmochimica Acta, v. 56, no. 10, p. 3627-3641.