# compare steady state flux through interlayers, standard column and option stagnant.

DATABASE c:\phreeqc\database\phreeqc.dat
SOLUTION_MASTER_SPECIES
#element        species alk     gfw_formula     element_gfw
Hto     Hto             0.0     1       1
Na_tr   Na_tr+          0.0     1       1
Cl_tr   Cl_tr-          0.0     1       1
Sr_tr   Sr_tr+2         0.0     1       1
SOLUTION_SPECIES
Hto = Hto; log_k 0; -gamma 1e6 0; -erm_ddl 1#-dw 2.236e-9
Na_tr+ = Na_tr+; log_k 0; -gamma 1e6 0; #-dw 1.33e-9
Cl_tr- = Cl_tr-; log_k 0; -gamma 1e6 0; #-dw 2.03e-9
Sr_tr+2 = Sr_tr+2; log_k 0; -gamma 1e6 0; #-dw 0.993e-9
EXCHANGE_SPECIES
Na_tr+ + X- = Na_trX; log_k 0.0; -gamma	1e6 0
Sr_tr+2 + 2X- = Sr_trX2; log_k 0.91; -gamma 1e6 0

PHASES
 A_sr
 Sr_trCl2 = Sr_tr+2 + 2Cl-; log_k -14
 A_na
 Na_trCl = Na_tr+ + Cl-; log_k -14
 A_cl
 NaCl_tr = Na+ + Cl_tr-; log_k -14
 A_hto
 Hto = Hto; log_k -15
END

# when doubling the number of cells, decrease the lengths and -water
SOLUTION 0; -water 5.0000e-01
Na 100; Cl 10 charge
Sr_tr 1e-3; Na_tr 1e-3; Hto 1e-3; Cl_tr 1e-3
END

# when doubling the number of cells, decrease the lengths and -water
# the center of the 0-concentration cell should remain at 2.5 m
SOLUTION 1-10; -water 0.5
Na 100; Cl 10 charge
EXCHANGE 1-10; X 0.05; -equil 1
END

EQUILIBRIUM_PHASES 3; A_hto 0 0; A_sr 0 0; A_na 0 0; A_cl 0 0
END

PRINT; -reset false
USER_GRAPH; -head M_HTO M_Na_tr M_Sr_tr/10 analyt_HTO ~~Na_tr ~~Sr_tr
-plot_conc t
-initial_solutions false
-axis_titles "time / years" "mass transfer / mol"
 -start
1 if cell_no = 3 then 19
2 if cell_no < 2 then end
# Mass flow HTO = J * A * DDt = -(Dw * por^n) * (DDc / DDx) * (water / length) * DDt (mol)
#               = -(2e-9 m2/s * 0.6^1) * (1e-3 mol/m3 / -2.5 m) * (0.5e-3 m3 / 1 m) * 5e9 s
#           = 1.2e-6 mol   (== PHRQC)
#      m3                / 1m         mol/m3          / 1m                        / G 
2 A = tot("water") / 1e3 / 1 : DDc = tot("Hto") * 1e3 / 1 : DDt = 5e9 : Dw = 2e-9 / 0.6^-1
3 put(A, 1, 1) : put(DDc, 1, 2) : put(DDt, 1, 3) : put(Dw, 1, 4) : 
#     por_IL/por_w                           mol/kgw * 1e3 dm3/m3 * kgw       / dm3 / 1 m       Dw  / G_IL * por_IL/por_w     
4 Ilw = 0.3 / 0.6 * tot("water") : DDNaX = mol("Na_trX") * 1e3 * tot("water") / Ilw / 1 : DwX = 2e-9 / 3 * 0.3 / 0.6
5 put(DDNaX, 1, 5) : put(DwX, 1, 6)
6 DDSrX = mol("Sr_trX2") * 1e3 * tot("water") / Ilw / 1
7 put(DDSrX, 1, 7)
10 end
19 x = total_time / (24 * 3600 * 365)
20 A_hto = equi("A_hto") : A_sr = equi("A_sr") :  A_na = equi("A_na")
40 J_hto = (A_hto - get(1)) : J_sr = (A_sr - get(2)) : J_na = (A_na - get(3))
60 plot_xy x, J_hto * 1, symbol_size = 0
70 plot_xy x, J_na * 1 , symbol_size = 0
72 plot_xy x, J_sr  *0.1     , symbol_size = 0
80 A = get(1, 1) : DDc = get(1, 2) : DDt = get(1, 3) : Dw = get(1, 4)
90 M_hto = Dw * A * DDc * DDt
100 plot_xy x, M_hto, color = Red, line_width = 0
110 DDNaX = get(1, 5) : DwX = get(1, 6)
120 M_Na = M_hto + DwX * A * DDNaX * DDt
130 plot_xy x, M_Na, color = green, line_width = 0
140 DDSrX = get(1, 7)
150 M_Sr = M_hto + DwX * A * DDSrX * DDt
160 plot_xy x, M_Sr / 10, color = blue, line_width = 0

180 put(A_hto, 1) : put(A_sr, 2) : put(A_na, 3)
 -end
TRANSPORT
 -shifts 10
 -flow diff; -cells 3; -bcon 1 2;
 -lengths 1
 -time 5e9  1
 -multi_D true 2e-09 0.6 0 1
 -interlayer_D true 0.3 0 3
END

# option stagnant. cell 3 is the constant c cell, 4 stagnant cells.
# mixf = Dp * DDt / (length_tot / internal cells) * (A / V)
# for comparison with regular column:
# mixf = (2e-9 * 0.6) * (5e9) / (2.5  / 4) * (0.5e-3 / 0.001 ) = 4.8
# 2 * mixf for boundary cells 3 ==> 4 and 7 ==> 8
# or mid-way last cell...
# mixf = (2e-9 * 0.6) * (5e9) / (2.5  / 4.5) * (0.5e-3 / 0.001 ) = 5.4
# 2 * mixf for boundary cell 3 ==> 4
PRINT; -reset true
USER_GRAPH 1; -active false
EQUILIBRIUM_PHASES 0-9 # remove from previous simulation
EXCHANGE 0-9
 X 1e-10; -equil 0
SOLUTION 0-2; Na 1
END

SOLUTION 3; -water 5.0000e5
Na 100; Cl 10 charge
Sr_tr 1e-3; Na_tr 1e-3; Hto 1e-3; Cl_tr 1e-3
END

SOLUTION 4-9; -water 0.5
Na 100; Cl 10 charge
EXCHANGE 4-7; -equil 4;  X 0.05
END

# MIX 3; 4 9.6
# MIX 4; 5 4.8
# MIX 5; 6 4.8
# MIX 6; 7 4.8
# MIX 7; 8 9.6
MIX 3; 4 10.8
MIX 4; 5 5.4
MIX 5; 6 5.4
MIX 6; 7 5.4
MIX 7; 8 5.4
USE mix none
END

EQUILIBRIUM_PHASES 8; A_hto 0 0; A_sr 0 0; A_na 0 0; A_cl 0 0
END

USER_GRAPH 1; -head stag
-active true
-connect_simulations false
1 if total_time = 5e9 then put(0, 1) : put(0, 2) : put(0, 3) :
19 x = total_time / (24 * 3600 * 365)
20 A_hto = equi("A_hto") : A_sr = equi("A_sr") :  A_na = equi("A_na")
40 J_hto = (A_hto - get(1)) : J_sr = (A_sr - get(2)) : J_na = (A_na - get(3))
60 plot_xy x, J_hto    
70 plot_xy x, J_na     
72 plot_xy x, J_sr / 10

180 put(A_hto, 1) : put(A_sr, 2) : put(A_na, 3)

TRANSPORT
 -length 10 # for 1 MCDrun for cells 0 and 1
 -war true; -flow diff; -cells 1; -bcon 1 2; -stag 6
 -time 5e9 50 # reduce the mixf to < 0.36
 -initial_time 0
 -punch_c 8
END