DATABASE c:\phreeqc\database\wateq4f.dat # database without diffusion coefficients gives equal diffusivity for all solutes

# Compare transport of tetra and mixing into stagnant zones (Table 10.3 in Appelo and Postma, p. 508)
#   with a multicomponent diffusion calculation.
# Instead of mixing of solutions, mole transfer is calculated with multi_d.
# To obtain concentrations, explicitly define that the stagnant cells have half the water of the mobile cells.

# Define tetra (CCl4) and its sorption in the stagnant zone
SOLUTION_MASTER_SPECIES; Tetra Tetra 0.0 Tetra 1.0
SOLUTION_SPECIES; Tetra = Tetra
SURFACE_MASTER_SPECIES; Sor Sor
SURFACE_SPECIES
 Sor = Sor
# Stagnant water is half of mobile: for retardation = 3, reduce sorption by half.
 Sor + Tetra = SorTetra; log_k -99.82             # log(Kf = 1.5) - log(m_Sor)
END

SOLUTION 1-21                   # mobile water
SOLUTION 22-43; -water 0.5      # stagnant water
SURFACE  22-41                  # sorption calculated with SURFACE
 Sor 1e100 1.0 1e100
END

SOLUTION 0; Na 1; Cl 1; Tetra 1; # For Cl diffusion in multiD, the solutions must be neutral

# Define mixing factors
SELECTED_OUTPUT 1; -reset false; -file st.prn
USER_PUNCH 1
10 sc$ = chr$(59)
20 for i = 1 to 20
30   mf = 0.85
40   punch eol$ + 'MIX ' + str$(i)  + sc$ + str$(i) + str_f$(mf, 6, 3) + sc$ + str$(i + 21) + str_f$(1 - mf, 6, 3)
50 next i
END

SELECTED_OUTPUT 1; -active false
INCLUDE$ st.prn
USE mix none
END

PRINT; -reset false
TRANSPORT; -cells 20; -lengths 50; -shifts 5; -time_step 1 year  # 5 years pollution
-multi_d true 1e-9 1 0.05 1.0 false # por = 1, n = 1, together give tortuosity = 1
-stagnant 1
END

SOLUTION 0              # Cleanup water
TRANSPORT; -shifts 10  # 10 years cleaning
-punch_frequency 10   # punch/graph only 10th shift
-punch  1-20, 22-41     # punch mobile flow and stagnant area
USER_GRAPH 1 multicomponent diffusion...
-headings  CCl4 Cl CCl4_stag Cl_stag
-axis_scale x_axis 0 1000
-axis_titles  "Distance / m" "Concentration, mg/l"
-initial_solutions  false
10 x = dist
20 if cell_no < 21 then plot_xy x, tot("Tetra")*1e3 : plot_xy x, tot("Cl")*1e3
# # to get other line symbols and colors, plot ghost-numbers (975, 0) for the plot-lines defined in line 20
30 if cell_no > 21 then plot_xy 975, 0 : plot_xy 975, 0 : plot_xy x, tot("Tetra")*1e3 : plot_xy x, tot("Cl")*1e3
END

# compare with option stagnant with exchange factor
USER_GRAPH 1; -active false

# Reinitialize the column...
SOLUTION 1-21  
SOLUTION 22-43 # relative amounts of water are defined in option -stagnant
# redefine sorption to give retardation = 3
SURFACE_SPECIES; Sor + Tetra = SorTetra; log_k -99.52 # log ((Kf = 3) / m_Sor) for 1 kg water in the stagnant layer
SURFACE  22-41; Sor 1e100 1.0 1e100
END

SOLUTION 0; Cl 1; Tetra 1; # With -multi_d false, solutions are mixed, the solutions can be charge-imbalanced
# Mixing with mixing factors needs definition of MIX for stagnant cells, use with -stagnant 1 #1.55e-9  0.24  0.12
SELECTED_OUTPUT 1; -reset false; -file st.prn; -active true
USER_PUNCH 1
10 sc$ = chr$(59)
20 for i = 1 to 20
30   mf = 0.85
40   punch eol$ + 'MIX ' + str$(i)  + sc$ + str$(i) + str_f$(mf, 6, 3) + sc$ + str$(i + 21) + str_f$(1 - mf, 6, 3)
42   mf = 0.70 # because the porosity is twice smaller, mix a twice higher fraction into stagnant cells
44   punch eol$ + 'MIX ' + str$(i + 21)  + sc$ + str$(i + 21) + str_f$(mf, 6, 3) + sc$ + str$(i) + str_f$(1 - mf, 6, 3)
50 next i
END

SELECTED_OUTPUT 1; -active false
INCLUDE$ st.prn
USE mix none
END

TRANSPORT; -shifts  5
-multi_d false
-stagnant  1  1.55e-9  0.24  0.12 # 1 stagnant layer, alpha, mobile por, immobile por
END

SOLUTION 0                   # Clean water again
TRANSPORT; -shifts 10       # 10 years cleanup
USER_GRAPH 2 Option -stagnant...
  -headings  CCl4 Cl CCl4_stag Cl_stag
  -axis_scale x_axis 0 1000
  -axis_titles  "Distance / m" "Concentration, mg/l"
  10 x = dist
  20 if cell_no < 21 then plot_xy x, tot("Tetra")*1e3 : plot_xy x, tot("Cl")*1e3
  30 if cell_no > 21 then plot_xy 975, 0 : plot_xy 975, 0 : plot_xy x, tot("Tetra")*1e3 : plot_xy x, tot("Cl")*1e3
END