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") |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
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?