In bentrueman/pbcusol: Predict Lead and Copper Solubility

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

pbcusol

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.

Installation

You can install the development version from GitHub with:

# 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().

library("pbcusol")
library("tidyverse")
library("plyr", include.only = "round_any")

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!)

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.

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")

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.

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:

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")

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).

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).

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).

# 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.

bentrueman/pbcusol documentation built on Oct. 25, 2024, 1:06 p.m.