DATABASE c:\phreeqc\database\wateq4f.dat
SOLUTION_MASTER_SPECIES; Tetra Tetra 0.0 Tetra 1.0
SOLUTION_SPECIES; Tetra = Tetra; log_k 0.0
SURFACE_MASTER_SPECIES; Sor Sor
SURFACE_SPECIES
 Sor = Sor; log_k 0.0
 Sor + Tetra = SorTetra; log_k -99.52              # log(Kf = 3) - log(m_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)
SOLUTION 1-101                                     # The sand and clay aquifer
SOLUTION 102-202; -water 0.5 # for mcd, specify (relative) amounts of water, used to obtain the interface area for diffusion

SURFACE  102-201                             # Sorption in stagnant cells (clay layer)
 Sor 1e100 1 1
END

SOLUTION 0; Na 1; Cl 1; Tetra 1                          # The polluted water, must be neutral for diffusion of Cl-
# Diffusion is calculated from mixing factors, use with -multi_D true;  -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) : nx = 100
20 for i = 1 to nx
30   mf = 0.953
40   punch eol$ + 'MIX ' + str$(i)  + sc$ + str$(i) + str_f$(mf, 6, 3) + sc$ + str$(i + nx + 1) + str_f$(1 - mf, 6, 3)
50 next i
END

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

TRANSPORT; -cells 100; -shifts 25 1; -timest 0.73e7      # 5 years pollution
 -length 10; -disp 10
 # -stagnant 1 1.55e-9 0.24 0.12      # 1 stagnant layer, alpha, mobile por, immobile por
 -multi_d true 1e-9 1 0.05 1.0 false # por = 1, n = 1, together give tortuosity = 1
 -implicit true 1.5
 -stagnant 1 
END

SOLUTION 0                                         # Clean water enters
TRANSPORT; -shifts 50 1                            # 10 years cleaning
 -punch_fr 50                                      # punch/graph only 50th shift
 -punch  1-100 102-201                             # punch mobile and stagnant cells
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 < 101 then plot_xy x, tot("Tetra")*1e3 : plot_xy x, tot("Cl")*1e3
  30 if cell_no > 101 then plot_xy 975, 0 : plot_xy 975, 0 : plot_xy x, tot("Tetra")*1e3 : plot_xy x, tot("Cl")*1e3
END