Simulate Calcite Solubility with tidyphreeqc

calcite solubility.R

library(tidyverse)
library(tidyphreeqc)
library(ragg)

T_range <- c(5, 10, 15, 20)
CO2_range <- seq(300, 1000, 5)

variables <- expand.grid(T_range, CO2_range)
names(variables) <- c("T [°C]", "CO2 [ppm]")

runconstructor <- function(`T [°C]`, `CO2 [ppm]`){
  phr_input(
    phr_solution(temp = `T [°C]`),
    phr_equilibrium_phases(
      "CO2(g)" = c(log10(`CO2 [ppm]`*10^-6), 10),
      "O2(g)" = c(-0.678, 10),
      "Calcite" = c(0, 100)
    ),
    phr_selected_output(
      pH = TRUE, 
      pe = TRUE,
      equilibrium_phases = c("CO2(g)", "O2(g)", "Calcite")
    )
  )
}

runlist <- pmap(variables, runconstructor)
simres <- map(runlist, phr_run) %>% 
  map_dfr(., as_tibble) %>% 
  filter(state == "react") %>% 
  cbind(variables)

simres$`T [°C]` <- factor(
  simres$`T [°C]`,
  levels = simres$`T [°C]` %>% unique() %>% sort()
)

agg_png(
  filename = "Calcite_solubility.png", 
  width = 16, height = 14, units = "cm",
  res = 600
  )

ggplot(
  data = simres,
  aes(
    x = `CO2 [ppm]`,
    y = d_Calcite *(-1)*100.09*10^3,
    color = `T [°C]`
  )
) +
  geom_line(size = 1.2) +
  geom_vline(xintercept = 410, color = "red", size = 2) +
  annotate(geom = "text", x = 580, y = 24, 
           label="you are here", color= "red", size = 8) +
  geom_segment(x = 500, y = 30, xend = 420, yend = 50,
                    colour= "red", size=1.5,
                    arrow = arrow(length = unit(0.5, "cm"))) +
  viridis::scale_color_viridis(discrete = TRUE, option = "C") +
  coord_cartesian(ylim = c(0, 100)) +
  labs(
    title = "Solubility of Calcite in 1 l pure water",
    y = "Solubilitiy of Calcite [mg/l]"
    ) +
  theme_dark() +
  theme(
    axis.text = element_text(color = "black", size = 12),
    axis.title = element_text(color = "black", size = 14),
    legend.text = element_text(color = "black", size = 12),
    legend.title = element_text(color = "black", size = 14)
    )
invisible(dev.off())