Kinetics of Fe(2) oxidation and Fe(3) precipitation

[Benkaito]

Please help, I am having trouble generating the graph. Also, the pH doesn't change after simulation. Here is the code:

Code: [Select]

SOLUTION_MASTER_SPECIES
    Fe_di         Fe_di+2          0     Fe_di           55.847
    Fe_tri        Fe_tri+3         0     Fe_tri          55.847

SOLUTION_SPECIES
Fe_di+2 = Fe_di+2
    log_k     0
Fe_di+2 + H2O = Fe_diOH+ + H+
    log_k     -9.5
    delta_h   13.2 kcal
Fe_di+2 + 3H2O = Fe_di(OH)3- + 3H+
    log_k     -31.0
    delta_h   30.3 kcal
Fe_tri+3 = Fe_tri+3
    log_k     0
Fe_tri+3 + H2O = Fe_triOH+2 + H+
    log_k     -2.19
    delta_h   10.4 kcal
Fe_tri+3 + 2H2O = Fe_tri(OH)2+ + 2H+
    log_k   -5.67
    delta_h 17.1   kcal
Fe_tri+3 + 3H2O = Fe_tri(OH)3 +3H+
    log_k     -12.56
    delta_h   24.8 kcal

PHASES
Fe_tri(OH)3(a)
    Fe_tri(OH) 3 + 3H+ = Fe_tri+3 + 3H2O
    log_k     4.891

SOLUTION 1
    temp      25
    pH        3.68
    pe        4 O2(g) -0.67
    redox     pe
    units     mg/l
    density   1
    Fe_di     349
    Fe_tri    5
    Cl        649 charge
    Na        649
    -water    1 # kg

EQUILIBRIUM_PHASES 1
    Fe_tri(OH)3(a) 0 0
    O2(g)     -0.67

RATES
    Fe_di_ox
-start
10 Fe_di = TOT("Fe_di")
20 if(Fe_di<=0) then goto 200
30 mo2 = mol("O2")
40 mNO3 = mol("NO3-")
50 moles = 10^13 * (ACT("OH-"))^2*(mNO3 + mo2)*Fe_di * TIME
200 SAVE moles
-end
    Fe_tri(OH)3(a)
-start
10 Ao = parm(1)
20 V = parm(2)
30 rate = 10^-14*(1-SR("Fe_tri(OH)3(a)")) * A0/V*(m/m0)^0.67
40 moles = rate * time
50 SAVE moles
-end

KINETICS 1
Fe_di_ox
    -formula  Fe_di  -1 Fe_tri  1
Fe_tri(OH)3(a)
    -formula  Fe_tri(OH)  1
    -m        1
    -m0       1
    -parms    22.6 0.162
-steps       1440 in 60 steps # seconds

INCREMENTAL_REACTIONS True

SELECTED_OUTPUT
    -file                 Kinetiks.out.xls
    -reset                false
    -time                 true
    -step                 true
    -ph                   true
USER_PUNCH
-headings
-start
10 PUNCH SIM_TIME/3600/24, TOT("Fe_di")*1000*55.847, TOT("Fe_tri")*1000*55.847, -LA("H+"), SI("Fe_tri(OH)3(a)")
-end
USER_GRAPH 1
    -headings               Days, Fe(2), Fe(3), pH, SI-Fi_tri(OH)3a
    -axis_titles            "Time (days)" "mg/L" "pH"
    -chart_title            "Kinetics of Fe (2) oxidation and Fe(3) precipitation"
    -axis_scale sy_axis     0 7 3.5 3.5
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600 / 24
20 GRAPH_Y TOT("Fe_di")*1000*55.847
30 GRAPH_Y TOT("Fe_tri")*1000*55.847
40 GRAPH_SY -LA("H+")
  -end
END
   

[dlparkhurst]

I think you have a typo with "Ao". Also, you are starting at a low pH, so the ACT("OH-")^2 is a very small number, so your rate is very small. I also am not sure you want to include Fe_tri(OH)3(a) in EQUILIBRIUM_PHASES because you are using KINETICS to simulate that reaction.

I suggest you systematically start with example 9 and convert it to your file, first by changing the SOLUTION, and then changing the RATES definition. I would eliminate the Fe_tri(OH)3(a), and focus on the oxidation reaction.  Note that you can use PRINT statements in the RATES definition to look at the values that are calculated.

[Benkaito]

Thank you! I will go with Fe(II) oxidation. Now I want to incorporate advective transport into the code but I am not sure how to go about it. I want to simulate the oxidation of Fe(II) into Fe(III) when the solution (written in the code) is passed through a column (that has unlimited amount of oxygen). Also, I added the PRINT but it doesn't add anything that wasn't there before. Thanks again, David. Your suggestions are always insightful and helpful!

Code: [Select]


SOLUTION_MASTER_SPECIES
    Fe_di         Fe_di+2          0     Fe_di           55.847
    Fe_tri        Fe_tri+3         0     Fe_tri          55.847

SOLUTION_SPECIES
Fe_di+2 = Fe_di+2
    log_k     0
Fe_di+2 + H2O = Fe_diOH+ + H+
    log_k     -9.5
    delta_h   13.2 kcal
Fe_di+2 + 3H2O = Fe_di(OH)3- + 3H+
    log_k     -31.0
    delta_h   30.3 kcal
Fe_tri+3 = Fe_tri+3
    log_k     0
Fe_tri+3 + H2O = Fe_triOH+2 + H+
    log_k     -2.19
    delta_h   10.4 kcal
Fe_tri+3 + 2H2O = Fe_tri(OH)2+ + 2H+
    log_k   -5.67
    delta_h 17.1   kcal
Fe_tri+3 + 3H2O = Fe_tri(OH)3 +3H+
    log_k     -12.56
    delta_h   24.8 kcal

PHASES
Goethite
    Fe_triOOH + 3H+ = Fe_tri+3 + 2H2O
    log_k     -1.0
END

SOLUTION 1
    temp      25
    pH        6.2
    pe        4 O2(g) -0.67
    redox     pe
    units     mg/l
    density   1
    Al        43.2
    Cu        1.4
    Fe_di     351
    Mn        21.8
    Zn        5.5
    -water    1 # kg

EQUILIBRIUM_PHASES 1
    O2(g)     -0.67


RATES
    Fe_di_ox
-start
10 Fe_di = TOT("Fe_di")
20 if (Fe_di <= 0) then goto 200
30 p_o2 = 10^(SI("O2(g)"))
40 moles = (2.91e-9 + 1.33e12 * (ACT("OH-"))^2 * p_o2) * Fe_di * TIME
200 SAVE moles
-end
PRINT
    -equilibrium_phases    true
    -kinetics              true
    -saturation_indices    true
    -totals                true
    -warnings              -1

KINETICS 1
Fe_di_ox
    -formula  Fe_di  -1 Fe_tri  1
    -m        1
    -m0       1
    -tol      1e-008
-steps       100 400 10800 21600 50400 86400 172800 172800 172800 172800
-step_divide 1
-runge_kutta 3
-bad_step_max 500

INCREMENTAL_REACTIONS True

SELECTED_OUTPUT
    -file                 selectedbasic.out.xls
USER_PUNCH
-headings Days Fe(2) Fe(3) pH SI_goethite
-start
10 PUNCH SIM_TIME/3600/24 TOT("Fe_di")*1e6, TOT("Fe_tri")*1e6, -LA("H+"), SI("Goethite")
-end
USER_GRAPH 1
    -headings               Time Fe(2) Fe(3) pH
    -axis_titles            "Time (days)" "mmol/kg of water" "pH"
    -chart_title            "Oxidation of ferrous iron"
    -axis_scale sy_axis     0 8 4 4
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600 / 24
20 GRAPH_Y TOT("Fe_di")*1e6, TOT("Fe_tri")*1e6
30 GRAPH_SY -LA("H+")
  -end
    -active                 true



[dlparkhurst]

No, I meant you could add a Basic line in RATES to look at the details and help debug the rate calculation, like

100 PRINT moles, ACT("OH-")^2, ...

As for transport, you need a TRANSPORT data block. I would start with few cells, say four, and few shifts, and look closely at the results.

[Benkaito]

Thanks! I expect to see Fe(2) increase and Fe(3) decrease with pore volumes but that is not what I am getting when I run this code. I think I am missing something. I have attached the graphs
t

Code: [Select]

SOLUTION_MASTER_SPECIES
    Fe_di         Fe_di+2          0     Fe_di           55.847
    Fe_tri        Fe_tri+3         0     Fe_tri          55.847

SOLUTION_SPECIES
Fe_di+2 = Fe_di+2
    log_k     0
Fe_di+2 + H2O = Fe_diOH+ + H+
    log_k     -9.5
    delta_h   13.2 kcal
Fe_di+2 + 3H2O = Fe_di(OH)3- + 3H+
    log_k     -31.0
    delta_h   30.3 kcal
Fe_tri+3 = Fe_tri+3
    log_k     0
Fe_tri+3 + H2O = Fe_triOH+2 + H+
    log_k     -2.19
    delta_h   10.4 kcal
Fe_tri+3 + 2H2O = Fe_tri(OH)2+ + 2H+
    log_k   -5.67
    delta_h 17.1   kcal
Fe_tri+3 + 3H2O = Fe_tri(OH)3 +3H+
    log_k     -12.56
    delta_h   24.8 kcal

PHASES
Goethite
    Fe_triOOH + 3H+ = Fe_tri+3 + 2H2O
    log_k     -1.0
END

SOLUTION 0
    temp      25
    pH        6.5
    pe        4 O2(g) -0.67
    redox     pe
    units     mg/l
    density   1
    Al        43.2
    Cu        1.4
    Fe_di     351
    Mn        21.8
    Zn        5.5
    -water    1 # kg 

EQUILIBRIUM_PHASES 1
    O2(g)     -0.67

RATES
    Fe_di_ox
-start
10 Fe_di = TOT("Fe_di")
20 if (Fe_di <= 0) then goto 200
30 p_o2 = SR("O2(g)")
40 moles = (2.91e-9 + 1.33e12 * (ACT("OH-"))^2* p_o2)* Fe_di*TIME
100 PRINT moles,(2.91e-9 + 1.33e12 * (ACT("OH-"))^2* p_o2)* Fe_di*TIME
200 SAVE moles
-end

SELECTED_OUTPUT
    -file                 KPedited.out.xls

INCREMENTAL_REACTIONS True

KINETICS 1
Fe_di_ox
    -formula  Fe_di  -1 Fe_tri  1
    -m        1
    -m0       1
    -tol      1e-008
-steps       640800 in 200 steps # seconds
-step_divide 0.0001
-runge_kutta 3
-bad_step_max 500

 
SOLUTION 1-4
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mg/l
    density   1
    O(0)      2.1
    -water    1 # kg

TRANSPORT
    -cells                 4
    -shifts                10
    -time_step             300 # seconds
    -dispersivities        4*0.005
    -correct_disp          true
    -diffusion_coefficient 0

USER_PUNCH
-headings Days Fe(2) Fe(3) pH SI_goethite
-start
10 PUNCH SIM_TIME/3600/24 TOT("Fe_di")*1e6, TOT("Fe_tri")*1e6, -LA("H+"), SI("Goethite")
-end
USER_GRAPH 1
    -headings               Time Fe(2) Fe(3) pH
    -axis_titles            "Time (days)" "mmol/kg of water" "pH"
    -chart_title            "Oxidation of ferrous iron"
    -axis_scale sy_axis     0 8 4 4
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600 / 24
20 GRAPH_Y TOT("Fe_di")*1e6, TOT("Fe_tri")*1e6
30 GRAPH_SY -LA("H+")
  -end
    -active                 true

USER_GRAPH 2
    -headings               Fe(2) Fe(3)
    -axis_titles            "Pore volumes" "Millimoles per kilogram water"
    -chart_title            "Using TRANSPORT Data Block"
    -axis_scale x_axis      0 3
    -axis_scale y_axis      0 1.5
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 x = (STEP_NO + 0.5) / cell_no
20 PLOT_XY x, TOT("Fe_di")*1000, symbol = plus, symbol_size = 2
30 PLOT_XY x, TOT("Fe_di")*1000, symbol = plus, symbol_size = 2
  -end
    -active                 true

SELECTED_OUTPUT
    -file                 Transport.out.xls
    -reset                false
    -step
    -totals               Fe_di  Fe_tri
    -high_precision       true
END


[dlparkhurst]

Fe(2) should decrease and Fe(3) increase. The column should be more or less steady state after 4 shifts.

Code: [Select]

SOLUTION_MASTER_SPECIES
    Fe_di         Fe_di+2          0     Fe_di           55.847
    Fe_tri        Fe_tri+3         0     Fe_tri          55.847
SOLUTION_SPECIES
Fe_di+2 = Fe_di+2
    log_k     0
Fe_di+2 + H2O = Fe_diOH+ + H+
    log_k     -9.5
    delta_h   13.2 kcal
Fe_di+2 + 3H2O = Fe_di(OH)3- + 3H+
    log_k     -31.0
    delta_h   30.3 kcal
Fe_tri+3 = Fe_tri+3
    log_k     0
Fe_tri+3 + H2O = Fe_triOH+2 + H+
    log_k     -2.19
    delta_h   10.4 kcal
Fe_tri+3 + 2H2O = Fe_tri(OH)2+ + 2H+
    log_k   -5.67
    delta_h 17.1   kcal
Fe_tri+3 + 3H2O = Fe_tri(OH)3 +3H+
    log_k     -12.56
    delta_h   24.8 kcal
PHASES
Goethite
    Fe_triOOH + 3H+ = Fe_tri+3 + 2H2O
    log_k     -1.0
RATES
    Fe_di_ox
-start
10 Fe_di = TOT("Fe_di")
20 if (Fe_di <= 0) then goto 200
30 p_o2 = SR("O2(g)")
40 moles = (2.91e-9 + 1.33e12 * (ACT("OH-"))^2* p_o2)* Fe_di*TIME
100 PRINT moles,(2.91e-9 + 1.33e12 * (ACT("OH-"))^2* p_o2)* Fe_di*TIME
200 SAVE moles
-end
END
SOLUTION 0
    temp      25
    pH        6.5
    pe        4 O2(g) -0.67
    redox     pe
    units     mg/l
    density   1
    Al        43.2
    Cu        1.4
    Fe_di     351
    Mn        21.8
    Zn        5.5
    -water    1 # kg
END
EQUILIBRIUM_PHASES 1
    O2(g)     -0.67
END
KINETICS 1-4
Fe_di_ox
    -formula  Fe_di  -1 Fe_tri  1
    -m        1
    -m0       1
    -tol      1e-008
END
SOLUTION 1-4
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mg/l
    density   1
    O(0)      2.1
    -water    1 # kg
END
TRANSPORT
    -cells                 4
    -shifts                10
    -time_step             300 # seconds
    -dispersivities        4*0.005
    -correct_disp          true
    -diffusion_coefficient 0
    -punch_frequency 10
USER_GRAPH 1
    -headings               Time Fe(2) Fe(3) pH
    -axis_titles            "Time (days)" "umol/kg of water" "pH"
    -chart_title            "Oxidation of ferrous iron"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X dist
20 GRAPH_Y TOT("Fe_di")*1e6, TOT("Fe_tri")*1e6
30 GRAPH_SY -LA("H+")
  -end
    -active                 true
END


[Benkaito]

THANK YOU SO MUCH!!!