Manganeese 3D pourbaix diagram using PHREEQC and tidyphreeqc

mn_pourbaix.R

library(tidyverse)
library(gganimate)
library(tidyphreeqc)
phr_use_db_minteq()
theme_set(theme_bw())
library(ggspatial)

result <- phr_run(
  phr_solution_list(
    pH = seq(3, 13, 0.1), 
    pe = seq(-10, 20, 0.1), 
    Mn = 10^2, # mmol / kg water
    temp = seq(0, 70, 1)
  ),
  phr_selected_output(
    activities = c("MnOH+", "Mn+2", "Mn+3", "MnO4-2", "MnO4-"),
    saturation_indices = c(
      "Birnessite", "Bixbyite", "Hausmannite", "Manganite",
      "Nsutite", "Pyrocroite", "Pyrolusite"
    ),
    totals = "Mn",
    temp = TRUE
  )
)

result_df <- result %>%
  as_tibble() %>%
  select(pH, pe, temp = `temp(C)`, Mn = starts_with("Mn"), starts_with("la_"), starts_with("si_")) %>%
  # these values aren't quite equal each time
  mutate(Mn = 10^round(log10(Mn)))

result_activity <- result_df %>%
  select(pH, pe, temp, Mn, starts_with("la_")) %>%
  gather(species, log_activity, starts_with("la_")) %>%
  mutate(species = str_remove(species, "^la_")) %>%
  group_by(pH, pe, temp, Mn) %>%
  summarise(species = species[which.max(log_activity)], log_activity = max(log_activity)) %>%
  ungroup()

result_si <- result_df %>%
  select(pH, pe, temp, Mn, starts_with("si_")) %>%
  gather(species, saturation_index, starts_with("si_")) %>%
  mutate(species = str_remove(species, "^si_")) %>%
  group_by(pH, pe, temp, Mn) %>%
  summarise(species = species[which.max(saturation_index)], saturation_index = max(saturation_index)) %>%
  ungroup()

activity_polygons <- suppressWarnings(
  result_activity %>%
    mutate(species = factor(species), species_int = as.integer(species)) %>%
    group_by(Mn, temp) %>%
    nest() %>%
    mutate(
      rast = map(data, ~raster::rasterFromXYZ(tibble(.$pH, .$pe, .$species_int))),
      poly = map2(
        rast, data, 
        ~sf::st_as_sf(raster::rasterToPolygons(.x, dissolve = TRUE)) %>%
          mutate(species = levels(.y$species)[..species_int])
      )
    ) %>%
    unnest(poly) %>%
    sf::st_as_sf() %>%
    sf::st_simplify(preserveTopology = TRUE)
)

si_polygons <- suppressWarnings(
  result_si %>%
    mutate(species = factor(species), species_int = as.integer(species)) %>%
    group_by(Mn, temp) %>%
    nest() %>%
    mutate(
      rast = map(data, ~raster::rasterFromXYZ(tibble(.$pH, .$pe, .$species_int))),
      poly = map2(
        rast, data, 
        ~sf::st_as_sf(raster::rasterToPolygons(.x, dissolve = TRUE)) %>%
          mutate(species = levels(.y$species)[..species_int])
      )
    ) %>%
    unnest(poly) %>%
    sf::st_as_sf() %>%
    sf::st_simplify(preserveTopology = TRUE)
)

activity_labels <- activity_polygons %>% 
  sf::st_cast("MULTIPOLYGON") %>% 
  sf::st_cast("POLYGON", warn = FALSE) %>%
  mutate(area = sf::st_area(.)) %>%
  filter(area > 0.2) %>%
  sf::st_centroid() %>% 
  df_spatial()

si_labels <- si_polygons %>% 
  sf::st_cast("MULTIPOLYGON") %>% 
  sf::st_cast("POLYGON", warn = FALSE) %>%
  mutate(area = sf::st_area(.)) %>%
  filter(area > 0.2) %>%
  sf::st_centroid() %>% 
  df_spatial()

anim <- ggplot(NULL, aes(x, y)) +
  geom_polygon(
    aes(group = piece_id, fill = species), 
    data = si_polygons %>% df_spatial(), 
    col = NA
  ) +
  geom_polygon(
    aes(group = piece_id), 
    data = activity_polygons %>% df_spatial(), 
    col = "black", fill = NA
  ) +
  # geom_label(
  #   aes(label = species), 
  #   data = si_labels, 
  #   col = "grey50", alpha = NA
  # ) +
  geom_label(
    aes(label = species), 
    data = activity_labels, 
    col = "black", alpha = NA
  ) +
  scale_fill_discrete(
    breaks = c("Birnessite", "Bixbyite", "Hausmannite", "Manganite", "Nsutite", "Pyrocroite", "Pyrolusite"),
    labels = c(
      "Birnessite (MnO2)", "Bixbyite (Mn2O3)", "Hausmannite (Mn3O4)", 
      "Manganite (MnOOH)", "Nsutite (MnO2)", "Pyrocroite (Mn(OH)2)", 
      "Pyrolusite (MnO2)"
    )
  ) +
  labs(
    title = "Dominant species of Mn at 0.1 mol/kg water, {round(frame_time)}°C",
    caption = "As modelled by PHREEQC/minteq database using tidyphreeqc",
    fill = "Most saturated species",
    alpha = "Log Activity",
    x = "pH", y = "pe"
  ) +
  transition_time(temp)

animate(anim, nframes = 300, fps = 10, width = 600, height = 600)
anim_save("mn_animation.gif")

Output:

Hello,
Your code is amazing thank you. It was working fine until a few weeks ago and package update I think, and since then I get this error message, and can't find the answer. Do you have the same problem? Thank you

Error: Can't cast poly$geometry <sfc_MULTIPOLYGON> to poly$geometry <sfc_MULTIPOLYGON>.
Call rlang::last_error() to see a backtrace

rlang::last_error()

message: Can't cast `poly$geometry` to `poly$geometry` . class: `vctrs_error_incompatible_cast` backtrace: 1. base::suppressWarnings(...) 18. vctrs:::vec_cast.default(x = x, to = to, x_arg = x_arg, to_arg = to_arg) 19. vctrs::stop_incompatible_cast(x, to, x_arg = x_arg, to_arg = to_arg) 20. vctrs:::stop_incompatible(...) 21. vctrs:::stop_vctrs(...)

The problem is that sfc objects can't be combined because of the way that unnest() works. Off the top of my head, I don't know what the easiest way around this is, except to continue using an old version of tidyr. You could do so using remotes::install_github("cran/tidyr@61a7c3d"). In the next release of the sf package this should be fixed.

I want to make Pourbaix diagrams that reflect a database I created. Is there a way to load my own database into tidyphreeqc?

Yes! You can do tidyphreeqc::phr_use_db_llnl(readLines("path_to_file.dat"))

Whoops...just tidyphreeqc::phr_use_db(readLines("path_to_file.dat")).

The 3D pourbaix diagram is beautiful and elegant but is a little complex. Will you consider creating a simple 2-D pourbaix diagram (for example, in 25 Celsius degree only), which is more often used. I guess it need convert regions into polygons as well, right?

Please I would like to inquire about how to calculate saturation index (SI) with tidyphreeqc package
I have a set of 50 wells analyzes
What does data entry look like?
How to do it?