Modelling Reactive transport of As(5) Zn Cd Ni Cu with surface complexation

[DM]

Dear PHREEQC users,

I am new to PHREEQC, which I would like to use to model the transport of As(5), Zn, Cd, Ni and Cu under saturated flow conditions in columns (0.12 cm) that contain sand coated with goethite+humic substances (my reactive surface is then goethite-coated with humic acids). SOLUTION 0 is the flushing solution containing the contaminants (As(5), Zn, Cd, Ni, Cu + CaCl2) which was injected for 1436400 s. The initial solution in column entered in SOLUTION 1-2 is 0.005 M CaCl2. After the pulse of contaminants the flushing solution is set back to 0.005 M CaCl2. The number of strong and weak sites for Hfo was adjusted to meet the measured Kd for As(5) from a batch experiment. The input script for sorption of cations to OM I used is based on the example 19B from the PHREEQC Version 3.

Questions:
1. I do not know how to include in my input data that after the injection of 1436400 s of SOLUTION 0, the background solution 0.005 M CaCl2 is injected again for 1184400 s?
2. I am having difficulties with the discretization of the column. If I increase the number of cells to 40 the graph looks very different.

Any advice is highly appreciated.

with kind regards,
Daniela

# DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.4.0-12927\database\wateq4f.dat

INPUT script I used is:
PRINT

-user_print true

SURFACE_MASTER_SPECIES
# Monodentate 60%
H_a H_aH; H_b H_bH; H_c H_cH; H_d H_dH
H_e H_eH; H_f H_fH; H_g H_gH; H_h H_hH

# Bidentate 40%
H_ab H_abH2; H_ad H_adH2; H_af H_afH2; H_ah H_ahH2
H_bc H_bcH2; H_be H_beH2; H_bg H_bgH2; H_cd H_cdH2
H_cf H_cfH2; H_ch H_chH2; H_de H_deH2; H_dg H_dgH2

SURFACE_SPECIES
H_aH = H_aH; log_k 0; H_bH = H_bH; log_k 0; H_cH = H_cH; log_k 0;  \
H_dH = H_dH; log_k 0;
H_eH = H_eH; log_k 0; H_fH = H_fH; log_k 0; H_gH = H_gH; log_k 0;  \
H_hH = H_hH; log_k 0;

H_abH2 = H_abH2; log_k 0; H_adH2 = H_adH2; log_k 0; H_afH2 = H_afH2; log_k 0;
H_ahH2 = H_ahH2; log_k 0; H_bcH2 = H_bcH2; log_k 0; H_beH2 = H_beH2; log_k 0;
H_bgH2 = H_bgH2; log_k 0; H_cdH2 = H_cdH2; log_k 0; H_cfH2 = H_cfH2; log_k 0;
H_chH2 = H_chH2; log_k 0; H_deH2 = H_deH2; log_k 0; H_dgH2 = H_dgH2; log_k 0;

# Protons
H_aH = H_a- + H+; log_k -1.59
H_bH = H_b- + H+; log_k -2.70
H_cH = H_c- + H+; log_k -3.82
H_dH = H_d- + H+; log_k -4.93

H_eH = H_e- + H+; log_k -6.88
H_fH = H_f- + H+; log_k -8.72
H_gH = H_g- + H+; log_k -10.56
H_hH = H_h- + H+; log_k -12.40

H_abH2 = H_abH- + H+; log_k -1.59; H_abH- = H_ab-2 + H+; log_k -2.70
H_adH2 = H_adH- + H+; log_k -1.59; H_adH- = H_ad-2 + H+; log_k -4.93
H_afH2 = H_afH- + H+; log_k -1.59; H_afH- = H_af-2 + H+; log_k -8.72
H_ahH2 = H_ahH- + H+; log_k -1.59; H_ahH- = H_ah-2 + H+; log_k -12.40
H_bcH2 = H_bcH- + H+; log_k -2.70; H_bcH- = H_bc-2 + H+; log_k -3.82
H_beH2 = H_beH- + H+; log_k -2.70; H_beH- = H_be-2 + H+; log_k -6.88
H_bgH2 = H_bgH- + H+; log_k -2.70; H_bgH- = H_bg-2 + H+; log_k -10.56
H_cdH2 = H_cdH- + H+; log_k -3.82; H_cdH- = H_cd-2 + H+; log_k -4.93
H_cfH2 = H_cfH- + H+; log_k -3.82; H_cfH- = H_cf-2 + H+; log_k -8.72
H_chH2 = H_chH- + H+; log_k -3.82; H_chH- = H_ch-2 + H+; log_k -12.40
H_deH2 = H_deH- + H+; log_k -4.93; H_deH- = H_de-2 + H+; log_k -6.88
H_dgH2 = H_dgH- + H+; log_k -4.93; H_dgH- = H_dg-2 + H+; log_k -10.56

# Calcium
H_aH + Ca+2 = H_aCa+ + H+; log_k -3.20
H_bH + Ca+2 = H_bCa+ + H+; log_k -3.20
H_cH + Ca+2 = H_cCa+ + H+; log_k -3.20
H_dH + Ca+2 = H_dCa+ + H+; log_k -3.20

H_eH + Ca+2 = H_eCa+ + H+; log_k -6.99
H_fH + Ca+2 = H_fCa+ + H+; log_k -6.99
H_gH + Ca+2 = H_gCa+ + H+; log_k -6.99
H_hH + Ca+2 = H_hCa+ + H+; log_k -6.99

H_abH2 + Ca+2 = H_abCa + 2H+; log_k -6.40
H_adH2 + Ca+2 = H_adCa + 2H+; log_k -6.40
H_afH2 + Ca+2 = H_afCa + 2H+; log_k -7.45
H_ahH2 + Ca+2 = H_ahCa + 2H+; log_k -10.2
H_bcH2 + Ca+2 = H_bcCa + 2H+; log_k -6.40
H_beH2 + Ca+2 = H_beCa + 2H+; log_k -10.2
H_bgH2 + Ca+2 = H_bgCa + 2H+; log_k -10.2
H_cdH2 + Ca+2 = H_cdCa + 2H+; log_k -6.40
H_cfH2 + Ca+2 = H_cfCa + 2H+; log_k -10.2
H_chH2 + Ca+2 = H_chCa + 2H+; log_k -10.2
H_deH2 + Ca+2 = H_deCa + 2H+; log_k -10.2
H_dgH2 + Ca+2 = H_dgCa + 2H+; log_k -10.2

# Cadmium
H_aH + Cd+2 = H_aCd+ + H+; log_k -1.52
H_bH + Cd+2 = H_bCd+ + H+; log_k -1.52
H_cH + Cd+2 = H_cCd+ + H+; log_k -1.52
H_dH + Cd+2 = H_dCd+ + H+; log_k -1.52

H_eH + Cd+2 = H_eCd+ + H+; log_k -5.57
H_fH + Cd+2 = H_fCd+ + H+; log_k -5.57
H_gH + Cd+2 = H_gCd+ + H+; log_k -5.57
H_hH + Cd+2 = H_hCd+ + H+; log_k -5.57

H_abH2 + Cd+2 = H_abCd + 2H+; log_k -3.04
H_adH2 + Cd+2 = H_adCd + 2H+; log_k -3.04
H_afH2 + Cd+2 = H_afCd + 2H+; log_k -7.09
H_ahH2 + Cd+2 = H_ahCd + 2H+; log_k -7.09
H_bcH2 + Cd+2 = H_bcCd + 2H+; log_k -3.04
H_beH2 + Cd+2 = H_beCd + 2H+; log_k -7.09
H_bgH2 + Cd+2 = H_bgCd + 2H+; log_k -7.09
H_cdH2 + Cd+2 = H_cdCd + 2H+; log_k -3.04
H_cfH2 + Cd+2 = H_cfCd + 2H+; log_k -7.09
H_chH2 + Cd+2 = H_chCd + 2H+; log_k -7.09
H_deH2 + Cd+2 = H_deCd + 2H+; log_k -7.09
H_dgH2 + Cd+2 = H_dgCd + 2H+; log_k -7.09

# Nickel
H_aH + Ni+2 = H_aNi+ + H+; log_k -1.4
H_bH + Ni+2 = H_bNi+ + H+; log_k -1.4
H_cH + Ni+2 = H_cNi+ + H+; log_k -1.4
H_dH + Ni+2 = H_dNi+ + H+; log_k -1.4

H_eH + Ni+2 = H_eNi+ + H+; log_k -4.5
H_fH + Ni+2 = H_fNi+ + H+; log_k -4.5
H_gH + Ni+2 = H_gNi+ + H+; log_k -4.5
H_hH + Ni+2 = H_hNi+ + H+; log_k -4.5

H_abH2 + Ni+2 = H_abNi + 2H+; log_k -2.8
H_adH2 + Ni+2 = H_adNi + 2H+; log_k -2.8
H_afH2 + Ni+2 = H_afNi + 2H+; log_k -5.9
H_ahH2 + Ni+2 = H_ahNi + 2H+; log_k -5.9
H_bcH2 + Ni+2 = H_bcNi + 2H+; log_k -2.8
H_beH2 + Ni+2 = H_beNi + 2H+; log_k -5.9
H_bgH2 + Ni+2 = H_bgNi + 2H+; log_k -5.9
H_cdH2 + Ni+2 = H_cdNi + 2H+; log_k -2.8
H_cfH2 + Ni+2 = H_cfNi + 2H+; log_k -5.9
H_chH2 + Ni+2 = H_chNi + 2H+; log_k -5.9
H_deH2 + Ni+2 = H_deNi + 2H+; log_k -5.9
H_dgH2 + Ni+2 = H_dgNi + 2H+; log_k -5.9

# Copper
H_aH + Cu+2 = H_aCu+ + H+; log_k -0.63
H_bH + Cu+2 = H_bCu+ + H+; log_k -0.63
H_cH + Cu+2 = H_cCu+ + H+; log_k -0.63
H_dH + Cu+2 = H_dCu+ + H+; log_k -0.63

H_eH + Cu+2 = H_eCu+ + H+; log_k -3.75
H_fH + Cu+2 = H_fCu+ + H+; log_k -3.75
H_gH + Cu+2 = H_gCu+ + H+; log_k -3.75
H_hH + Cu+2 = H_hCu+ + H+; log_k -3.75

H_abH2 + Cu+2 = H_abCu + 2H+; log_k -1.26
H_adH2 + Cu+2 = H_adCu + 2H+; log_k -1.26
H_afH2 + Cu+2 = H_afCu + 2H+; log_k -4.38
H_ahH2 + Cu+2 = H_ahCu + 2H+; log_k -4.38
H_bcH2 + Cu+2 = H_bcCu + 2H+; log_k -1.26
H_beH2 + Cu+2 = H_beCu + 2H+; log_k -4.38
H_bgH2 + Cu+2 = H_bgCu + 2H+; log_k -4.38
H_cdH2 + Cu+2 = H_cdCu + 2H+; log_k -1.26
H_cfH2 + Cu+2 = H_cfCu + 2H+; log_k -4.38
H_chH2 + Cu+2 = H_chCu + 2H+; log_k -4.38
H_deH2 + Cu+2 = H_deCu + 2H+; log_k -4.38
H_dgH2 + Cu+2 = H_dgCu + 2H+; log_k -4.38

# Zinc
H_aH + Zn+2 = H_aZn+ + H+; log_k -1.7
H_bH + Zn+2 = H_bZn+ + H+; log_k -1.7
H_cH + Zn+2 = H_cZn+ + H+; log_k -1.7
H_dH + Zn+2 = H_dZn+ + H+; log_k -1.7

H_eH + Zn+2 = H_eZn+ + H+; log_k -4.9
H_fH + Zn+2 = H_fZn+ + H+; log_k -4.9
H_gH + Zn+2 = H_gZn+ + H+; log_k -4.9
H_hH + Zn+2 = H_hZn+ + H+; log_k -4.9

H_abH2 + Zn+2 = H_abZn + 2H+; log_k -2.4
H_adH2 + Zn+2 = H_adZn + 2H+; log_k -2.4
H_afH2 + Zn+2 = H_afZn + 2H+; log_k -6.6
H_ahH2 + Zn+2 = H_ahZn + 2H+; log_k -6.6
H_bcH2 + Zn+2 = H_bcZn + 2H+; log_k -2.4
H_beH2 + Zn+2 = H_beZn + 2H+; log_k -6.6
H_bgH2 + Zn+2 = H_bgZn + 2H+; log_k -6.6
H_cdH2 + Zn+2 = H_cdZn + 2H+; log_k -2.4
H_cfH2 + Zn+2 = H_cfZn + 2H+; log_k -6.6
H_chH2 + Zn+2 = H_chZn + 2H+; log_k -6.6
H_deH2 + Zn+2 = H_deZn + 2H+; log_k -6.6
H_dgH2 + Zn+2 = H_dgZn + 2H+; log_k -6.6
END

SOLUTION 0 Pulse solution with contaminants
pH 7.25
Ca 0.005  #mol/kg water
Cl 0.01
Ni 2.87e-6
Cu 3.91e-5
Zn 4.77e-5
Cd 3.85e-6
As(5) 6.5e-5
END

SOLUTION 1-2 Background solution without contaminants
pH 7.25
Ca 0.005  #mol/kg water
Cl 0.01

SURFACE 1-2
Hfo_w 0.0001 90 0.0043
Hfo_s 0.0000025

H_a 1.42e-6 43.4e3 2e-3
H_b 1.42e-6; H_c 1.42e-6; H_d 1.42e-6

H_e 7.1e-7; H_f 7.1e-7; H_g 7.1e-7; H_h 7.1e-7

H_ab 4.73e-7; H_ad 4.73e-7; H_af 4.73e-7; H_ah 4.73e-7
H_bc 4.73e-7; H_be 4.73e-7; H_bg 4.73e-7; H_cd 4.73e-7
H_cf 4.73e-7; H_ch 4.73e-7; H_de 4.73e-7; H_dg 4.73e-7
-donnan
-equilibrate surface with solution 1-2
END

PRINT
-reset false

TRANSPORT
-cells 2
-lengths 0.06
-dispersivity 2*0.0006
-shifts 112
-flow_direction forward
-time_step 12876 # the residence time in a cell, which is linked to the velocity; 0.12/4.66e-6 m/s/2 is the pore water velocity of the column 644
-boundary_conditions flux flux # flux boundary conditions at inlet and outlet ends of the column
-correct_disp true
-punch_cells 2
-print_cells 2
-punch_frequency 10
-print_frequency 10


USER_GRAPH
-headings As(5) Zn Ni Cd Cu
-chart_title "Breakthrough curve: As(5) Zn Cd Ni Cu"
-axis_titles "Pore volume" "Concentration in ug/L"
-axis_scale x_axis 0 100
-plot_concentration_vs time
-start
10 plot_XY (step_no + 0.5)/cell_no, tot ("As(5)") * 74.9e6
20 plot_XY (step_no + 0.5)/cell_no, tot ("Zn") * 65.38e6
30 plot_XY (step_no + 0.5)/cell_no, tot ("Ni") * 58.69e6
40 plot_XY (step_no + 0.5)/cell_no, tot ("Cd") * 112.4e6
50 plot_XY (step_no + 0.5)/cell_no, tot ("Cu") * 63.55e6
-end
END

[dlparkhurst]

At the end of your file, append a new SOLUTION 0 definition for the CaCl2 solution, and then append a new TRANSPORT definition; only definitions that change in TRANSPORT (-shifts, for example) need to be defined in the new data block.

If you change the number of cells, you must change the time step to maintain the same velocities. Velocity is cell length divided by time step.

[DM]

Dear David,
I really appreciate your quick response.

The below input I  added at the end of my previous script; however the breakthrough curves are not correctly plotted.  How can I plot the complete breakthrough (adsorption+desorption part)?

SOLUTION 0 Desorption of contaminants
pH 7.25
Ca 0.005  #mol/kg water
Cl 0.01
END

TRANSPORT
-shifts 1839

USER_GRAPH 1
-headings As(5) Zn Ni Cd Cu
-chart_title Breakthrough curve: As(5) Zn Cd Ni Cu
-axis_titles "Pore volume" "Concentration in ug/L"
-axis_scale x_axis 0 100
-plot_concentration_vs time
-start
10 plot_XY (step_no + 0.5)/cell_no, tot ("As(5)") * 74.9e6
20 plot_XY (step_no + 0.5)/cell_no, tot ("Zn") * 65.38e6
30 plot_XY (step_no + 0.5)/cell_no, tot ("Ni") * 58.69e6
40 plot_XY (step_no + 0.5)/cell_no, tot ("Cd") * 112.4e6
50 plot_XY (step_no + 0.5)/cell_no, tot ("Cu") * 63.55e6
-end
END

[dlparkhurst]

If you want to know the amount of an element specifically sorbed on a surface, there is a function SURF. The amount of an element in the diffuse layer (you have defined -donnan) is found with the function EDL or EDL_SPECIES. Concentrations of individual surface species can be found with MOL (mol/kgw).

The total amount of an element in a cell can be found with SYS, which would include equilibrium phases and other reactants.

Look at the section The Basic Interpreter in the manual.