TITLE Example 11.--Transport and cation exchange.
SOLUTION 0  CaCl2
        units            mmol/kgw
        temp             25.0
        pH               7.0     charge
        pe               12.5    O2(g)   -0.68
        Ca               0.6
        Cl               1.2
SOLUTION 1-40  Initial solution for column
        units            mmol/kgw
        temp             25.0
        pH               7.0     charge
        pe               12.5    O2(g)   -0.68
        Na               1.0
        K                0.2
        N(5)             1.2
END
EXCHANGE 1-40
        -equilibrate 1
        X                0.0011
COPY cell 1 101
END
ADVECTION
        -cells           40
        -shifts          100
        -punch_cells     40
        -punch_frequency 1
        -print_cells     40
        -print_frequency 20
PRINT; -reset false; -status false
SELECTED_OUTPUT
        -file            ex11adv.sel
        -reset           false
        -step
        -totals          Na Cl K Ca
USER_PUNCH
  -heading  Pore_vol
  10 PUNCH (STEP_NO + .5) / 40.
USER_GRAPH 1 Example 11
  -chart_title "Using ADVECTION Data Block"
  -headings Cl Na K Ca
  -axis_titles "Pore volumes" "Millimoles per kilogram water"
  -axis_scale x_axis 0 2.5
  -axis_scale y_axis 0 1.4
  -plot_concentration_vs time
  -start
  10 x = (STEP_NO + 0.5) / cell_no
  20 PLOT_XY x, TOT("Cl")*1000, symbol = None
  30 PLOT_XY x, TOT("Na")*1000, symbol = None
  40 PLOT_XY x, TOT("K")*1000, symbol = None
  50 PLOT_XY x, TOT("Ca")*1000, symbol = None
  60 PUT(1, 1)
  -end
COPY cell 101 1-40
END
USER_GRAPH 1
        -detach
END
TRANSPORT
        -cells           40
        -lengths         40*0.002
        -shifts          100
        -time_step       720.0
        -flow_direction  forward
        -boundary_conditions   flux  flux
        -diffusion_coefficient 0.0e-9
        -dispersivities  40*0.002
        -correct_disp    true
        -punch_cells     40
        -punch_frequency 1
        -print_cells     40
        -print_frequency 20
SELECTED_OUTPUT
        -file            ex11trn.sel
        -reset           false
        -step
        -totals          Na Cl K Ca
        -high_precision true
USER_GRAPH 2 Example 11
  -chart_title "Using TRANSPORT Data Block"
  -headings Cl Na K Ca Cl_analytical
  -axis_titles "Pore volumes" "Millimoles per kilogram water"
  -axis_scale x_axis 0 2.5
  -axis_scale y_axis 0 1.4
#  -batch \temp\11.gif false # After saving, the chart on the monitor is closed.
  -plot_concentration_vs time
  10 x = (STEP_NO + 0.5) / cell_no
  20 PLOT_XY x, TOT("Cl")*1000, symbol = Plus, symbol_size = 2
  30 PLOT_XY x, TOT("Na")*1000, symbol = Plus, symbol_size = 2
  40 PLOT_XY x, TOT("K") *1000, symbol = Plus, symbol_size = 2
  50 PLOT_XY x, TOT("Ca")*1000, symbol = Plus, symbol_size = 2
# calculate Cl_analytical...
  60 DATA 0.254829592, -0.284496736, 1.421413741, -1.453152027, 1.061405429, 0.3275911
  70 READ a1, a2, a3, a4, a5, a6
  80 Peclet = 0.08 / 0.002
  90  z = (1 - x) / SQRT(4 * x / Peclet)
  100 PA = 0
  110 GOSUB 2000 # calculate e_erfc = exp(PA) * erfc(z)
  120 e_erfc1 = e_erfc
  130 z = (1 + x) / SQRT(4 * x / Peclet)
  140 PA = Peclet
  150 GOSUB 2000 # calculate exp(PA) * erfc(z)
  160 y = 0.6 * (e_erfc1 + e_erfc)
  170 PLOT_XY x, y, line_width = 0, symbol = Circle, color = Red
  180 d = (y - TOT("Cl")*1000)^2
  190 IF EXISTS(10) THEN PUT(d + GET(10), 10) ELSE PUT(d, 10)
  200 IF STEP_NO = 2 * CELL_NO THEN print 'SSQD for Cl after 2 Pore Volumes: ', GET(10), '(mmol/L)^2'
  210 END
  2000 REM calculate e_erfc = exp(PA) * erfc(z)...
  2010 sgz = SGN(z)
  2020 z = ABS(z)
  2050 b = 1 / (1 + a6 * z)
  2060 e_erfc = b * (a1 + b * (a2 + b * (a3 + b * (a4 + b * a5)))) * EXP(PA - (z * z))
  2070 IF sgz = -1 THEN e_erfc = 2 * EXP(PA) - e_erfc
  2080 RETURN
END