Skip to content

Instantly share code, notes, and snippets.

@paleolimbot
Created March 7, 2019 19:40
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save paleolimbot/aed40d173edc9c1bbc497e6187d6f763 to your computer and use it in GitHub Desktop.
Save paleolimbot/aed40d173edc9c1bbc497e6187d6f763 to your computer and use it in GitHub Desktop.
library(tidyverse)
library(tidyphreeqc)
solution_cacl2 <- phr_solution(
0, "CaCl2",
units = "mmol/kgw",
temp = 25,
pH = "7 charge",
pe = "12.5 O2(g) -0.68",
Ca = 0.6,
Cl = 1.2
)
solution_initial <- phr_solution(
"1-40", "Initial solution for column",
units = "mmol/kgw",
temp = 25,
pH = "7 charge",
pe = "12.5 O2(g) -0.68",
Na = 1,
K = 0.2,
"N(5)" = 1.2
)
# do advection calculation
df_advect <- phr_input(
solution_cacl2,
solution_initial,
phr_input_section("EXCHANGE", "1-40", components = list("-equilibrate" = 1, "X" = 0.0011)),
phr_input_section(
"ADVECTION",
components = list(
"-cells" = 40,
"-shifts" = 100,
"-punch_cells" = 40,
"-punch_frequency" = 1,
"-print_cells" = 40,
"-print_frequency" = 20
)
),
phr_selected_output(totals = c("Na", "Cl", "K", "Ca"), step = TRUE)
) %>%
phr_run() %>%
as_tibble()
# plot advection calculation
df_advect %>%
filter(state == "advect") %>%
ggplot(aes(time, `K(mol/kgw)` * 1000)) +
geom_line()
df_transport <- phr_input(
solution_cacl2,
solution_initial,
phr_input_section("EXCHANGE", "1-40", components = list("-equilibrate" = 1, "X" = 0.0011)),
phr_input_section(
"TRANSPORT",
components = list(
"-cells" = 40,
"-lengths" = 0.002,
"-shifts" = 100,
"-time_step" = 720.0,
"-flow_direction" = "forward",
"-boundary_conditions" = c("flux", "flux"),
"-diffusion_coefficient" = 0,
"-correct_disp" = "true",
"-punch_cells" = 40,
"-punch_frequency" = 1,
"-print_cells" = 40,
"-print_frequency" = 20
)
),
phr_selected_output(totals = c("Na", "Cl", "K", "Ca"))
) %>%
phr_run() %>%
as_tibble()
# plot transport calculation
df_transport %>%
filter(state == "transp") %>%
ggplot(aes(time, `K(mol/kgw)` * 1000)) +
geom_line()
@paleolimbot
Copy link
Author

library(tidyverse)
library(tidyphreeqc)

solution_cacl2 <- phr_solution(
  0, "CaCl2",
  units = "mmol/kgw",
  temp = 25, 
  pH = "7 charge", 
  pe = "12.5 O2(g) -0.68",
  Ca = 0.6,
  Cl = 1.2
)

solution_initial <- phr_solution(
  "1-40", "Initial solution for column",
  units = "mmol/kgw",
  temp = 25, 
  pH = "7 charge", 
  pe = "12.5 O2(g) -0.68",
  Na = 1,
  K = 0.2,
  "N(5)" = 1.2
)

# do advection calculation
df_advect <- phr_input(
  solution_cacl2,
  solution_initial,
  phr_input_section("EXCHANGE", "1-40", components = list("-equilibrate" = 1, "X" = 0.0011)),
  phr_input_section(
    "ADVECTION", 
    components = list(
      "-cells" = 40,
      "-shifts" = 100,
      "-punch_cells" = 40,
      "-punch_frequency" = 1,
      "-print_cells" = 40,
      "-print_frequency" = 20
    )
  ),
  phr_selected_output(totals = c("Na", "Cl", "K", "Ca"), step = TRUE)
) %>%
  phr_run() %>%
  as_tibble()

# plot advection calculation
df_advect %>%
  filter(state == "advect") %>%
  ggplot(aes(time, `K(mol/kgw)` * 1000)) +
  geom_line()

df_transport <- phr_input(
  solution_cacl2,
  solution_initial,
  phr_input_section("EXCHANGE", "1-40", components = list("-equilibrate" = 1, "X" = 0.0011)),
  phr_input_section(
    "TRANSPORT", 
    components = list(
      "-cells" = 40,
      "-lengths" = 0.002,
      "-shifts" = 100,
      "-time_step" = 720.0, 
      "-flow_direction" = "forward",
      "-boundary_conditions" = c("flux", "flux"),
      "-diffusion_coefficient" = 0,
      "-correct_disp" = "true",
      "-punch_cells" = 40,
      "-punch_frequency" = 1,
      "-print_cells" = 40,
      "-print_frequency" = 20
    )
  ),
  phr_selected_output(totals = c("Na", "Cl", "K", "Ca"))
) %>%
  phr_run() %>%
  as_tibble()

# plot transport calculation
df_transport %>%
  filter(state == "transp") %>%
  ggplot(aes(time, `K(mol/kgw)` * 1000)) +
geom_line()

Created on 2019-03-07 by the reprex package (v0.2.1)

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment