DATABASE ../database/phreeqc.dat
#PRINT
#   -reset   false

TITLE Example 12b.--Compare traditional and multicomponent diffusive transport of heat and solutes.
      Multicomponent diffusion coefficients decrease with the viscosity of the solution, changing more as temperature decreases.

SOLUTION 0   Fixed temp 24C, and NaCl conc (first type boundary cond) at inlet
    temp 24; Na 24; Cl 24
SOLUTION 1-31   24 mM KBr, initial temp 0C
    temp  0; K  24; Br 24
EXCHANGE_SPECIES
    K+ + X- = KX
EXCHANGE 1-31
    KX    0.048
END
TRANSPORT
   -cells 10
   -lengths 0.3333333
   -shifts 1
   -flow_direction         diffusion
   -boundary_conditions    constant  closed
   -thermal_diffusion      3.0    # heat is retarded equal to Na
   -diffusion_coefficient  0.3e-9  # m^2/s
# -multi_d true 0.3e-9 1 0.05 1.0 false
# -implicit true 
   -time_step              1.0e+9  1 # max_mixf = 2/9 = Dt_max * De / Dx^2. Dt_max = 8.2140e+07 seconds, Number of mixes = 1e10 / 8.214e7 = 122
USER_GRAPH 1 Example 12b
   -headings Tradit:Na Cl TC Analyt
   # -headings TC Analyt
   -chart_title "Compare traditional and multicomponent diffusive transport"
   -axis_titles "DISTANCE, IN METERS" "MILLIMOLES PER KILOGRAM WATER", "DEGREES CELSIUS"
   -initial_solutions false
   -plot_concentration_vs x
   -axis_scale sy_axis 0
   -start
  10 x = DIST
  20 PLOT_XY x, TOT("Na")*1000, symbol = Plus
  30 PLOT_XY x, TOT("Cl")*1000, symbol = Plus
  40 PLOT_XY x, TC, symbol = XCross, y-axis 2 : print tc, diff_c("Cl-")
  50 if (x > 10 OR SIM_TIME <= 0) THEN END
  60 DATA 0.254829592, -0.284496736, 1.421413741, -1.453152027, 1.061405429, 0.3275911
  70 READ a1, a2, a3, a4, a5, a6
# Calculate and plot Cl analytical...
  80 z = x / (2 * SQRT(3e-10 * SIM_TIME / 1.0))
  90 GOSUB 2000
  100 PLOT_XY x, 24 * erfc, color = Green, symbol = Square, symbol_size = 8,\
                            line_width = 0
# Calculate and plot 3 times retarded Na and temperature analytical...
  110 z = z * SQRT(3.0)
  120 GOSUB 2000
  130 PLOT_XY x, 24 * erfc, color = Blue, symbol = Square, symbol_size = 8,\
                            line_width = 0
  140 END
  2000 REM calculate erfc...
  2050 b = 1 / (1 + a6 * z)
  2060 erfc = b * (a1 + b * (a2 + b * (a3 + b * (a4 + b * a5)))) * EXP(-(z * z))
  2080 RETURN
   -end
END

# Reinitialize the column...
copy cell 31 1-30
END
TRANSPORT
-shifts 1
   -multi_d true 2.33e-9 1 0.05 0 false # will give the traditional results when tc = 25 throughout
   -thermal_diffusion  3.0 2.33e-9  # define the diffusion coefficient for heat equal to Na
USER_GRAPH 1 Example 12b
   -headings MultiD&Visc:Na Cl TC
   -start
  10 x = DIST
  20 PLOT_XY x, TOT("Na")*1000, symbol = Circle, line_width = 0, symbol_size = 5, color = Red
  30 PLOT_XY x, TOT("Cl")*1000, symbol = Circle, line_width = 0, symbol_size = 5, color = Green
  40 PLOT_XY x, TC,             symbol = Circle, line_width = 0, symbol_size = 5, y-axis 2, color = Blue
END
# Reinitialize the column...
copy cell 31 1-30
USER_GRAPH 1; -connect_simulations false
END
TRANSPORT
-shifts 1
   -multi_d true 2.33e-9 1 0.05 0 false # will give the traditional results when tc = 25 throughout
   -thermal_diffusion  3.0 2.33e-9  # define the diffusion coefficient for heat equal to Na
   -implicit true 3 -12             # max_mixf = 3, min_dif_LM = -12
USER_GRAPH 1 Example 12b
   -headings MultiD&Visc&Implicit:Na Cl TC
   -start
  10 x = DIST
  20 PLOT_XY x, TOT("Na")*1000, symbol = XCross, line_width = 0, symbol_size = 9, color = Red
  30 PLOT_XY x, TOT("Cl")*1000, symbol = XCross, line_width = 0, symbol_size = 9, color = Green
  40 PLOT_XY x, TC,             symbol = XCross, line_width = 0, symbol_size = 9, y-axis 2, color = Blue
END