Created
August 16, 2018 18:45
-
-
Save paleolimbot/e854c098c6d95330c361014276b5a8db to your computer and use it in GitHub Desktop.
Manganeese 3D pourbaix diagram using PHREEQC and tidyphreeqc
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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") |
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?
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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?