Dithionite redox reaction with Fe(III) minerals

[ktelfeya]

Hi everyone,

I am trying to model a simple batch reaction, but the resulting plot indicates that I am clearly making mistakes in the code. During the experiment, I added 10 mL of a NaS2O4 solution (buffered with 0.1 M NaHCO3) to 5g of sediment with ~1.5% Fe2O3 content with a surface area of ~ 2.2 m2/g. I then sampled the solution over ~40 days. The vessel was kept free from oxygen. I am trying to model the loss of S2O4-- due to reaction with Fe(III)(s). The S2O4-- will also degrade to intermediate S species due to hydrolysis, but the reaction with the sediment should overwhelm this.
The redox reaction of S2O4-- with Fe(III)(s) according to the literature is:
S2O4-- + 2Fe(III)(s) + 2H2O = 2SO3-- + 2Fe(II)(s) + 4H+
And the rate law has been reported as:
-dC/dt = k*Csurf*C, in which k is the rate constant, Csurf is the concentration of available surface sites in the sediment, and C is the concentration of dithionite (from Ammonette et al. (1994): Abiotic reduction of aquifer materials by dithionite: a promising in-situ remediation technology). I’m starting with a rate constant of 1.03e-5 from this report, but this value (and even this rate law) can change.
The experimental data give S2O4-- concentrations of 0.040, 0.019, 0.015, 0.005, and 0.002 M  at 1, 3, 11, 24 and 43 days, respectively (starting concentration 0.1 M).
Here is my code. The resulting S2O4-- concentrations are orders of magnitude too low. It’s definitely possible that I need to change the pe value, but because the redox is in disequilibrium, I’m not sure how to define it.

SOLUTION 1
    temp      25
    pH        8.3
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    Alkalinity 0.1 as HCO3-
    Na        0.15 charge
    S(3)      0.1 #dithionite
    S(4)      1e-06 #sulfite
    -water    0.01 # kg


EQUILIBRIUM_PHASES 1
    Fe(OH)3   0 .00144

KINETICS 1
Dithionite_Fe
    -formula  S2O4-2  1 Fe(OH)3  2
    -m        0.1
    -m0       0.1
    -parms    1.03e-005
    -tol      1e-008
-steps       86400 172800 259200 604800 950400 1382400 2073600 3715200
-step_divide 1
-runge_kutta 3
-bad_step_max 500

RATES
    Dithionite_Fe
-start
10 rem M = current number of moles of dithionite
20 rem M0 = number of moles of dithionite initially present
30 rem k1=PARM(1)
40 Csurf=EQUI("Fe(OH)3")
50 rf = PARM(1)*Csurf*M
60 MOLES = rf * TIME
70 Save MOLES
-end


USER_GRAPH 1
    -axis_titles            "time" "Dithionite (M)" "pH"
    -chart_title            "Dithionite reduction of Fe"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 GRAPH_X TOTAL_TIME/3600/24
20 GRAPH_Y TOT("S(3)")
30 GRAPH_SY -LA("H+")
  -end
    -active                 true


Thank you so much for your help.

Kat

[dlparkhurst]

The polythionates are unstable relative to sulfide and sulfate, so they won't show up in your calculations. You need to set up a system that allows S2O4-- and SO3-- to exist. Here is a potential way to do it using definitions for a "new element" that I have called .  I used the phreeqc.dat database.

You may want to consider the RATES definition. I have taken the approach that Fe(OH)3 is the kinetic reactant rather than an equilibrium phase (or some combination that you defined). Here Fe(OH)3 reacts kinetically and the  reacts accordingly, but there is no connection to Fe(OH)3 equilibrium.  So, my approach may not be the best.

SOLUTION_MASTER_SPECIES
          O3-2         0.0   32 32
(4)       O3-2         0.0   32
(3)       2O4-2        0.0   32

SOLUTION_SPECIES
O3-2 = O3-2
log_k  0.0


 2.0000 H+  + 2.0000 O3--  = 2O4--  + .500 O2  + H2O
#        -llnl_gamma           5.0   
        log_k           -25.2076
   -delta_H   0            # Not possible to calculate enthalpy of reaction   S2O4-2
   -analytic   -2.3172e2   2.0393e-3   -7.1011e0   8.3239e1   9.4155e-1

END
RATES
    Dithionite_Fe
-start
10 rem M = current number of moles of dithionite
20 rem M0 = number of moles of dithionite initially present
30 rem k1=PARM(1)
#40 Csurf=EQUI("Fe(OH)3")
#50 rf = PARM(1)*Csurf*M
50 rf = PARM(1)*M
60 MOLES = rf * TIME
70 Save MOLES
-end
END

SOLUTION 1
    temp      25
    pH        8.3
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    Alkalinity 0.1 as HCO3-
    Na        0.15 charge
    (3)      0.1 #dithionite
    (4)      1e-06 #sulfite
    -water    0.01 # kg
END
EQUILIBRIUM_PHASES 1
    Fe(OH)3(a)   0 .00144
END
USE solution 1
KINETICS 1
Dithionite_Fe
    #-formula  S2O4-2  1 Fe(OH)3  2
    -formula  Fe(OH)3  1
    -m        0.1
    -m0       0.1
    -parms    1.03e-005
    -tol      1e-008
-steps       86400 172800 259200 604800 950400 1382400 2073600 3715200
-step_divide 1
-runge_kutta 3
-bad_step_max 500


USER_GRAPH 1
    -axis_titles            "time" "Dithionite (M)" "pH"
    -chart_title            "Dithionite reduction of Fe"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 GRAPH_X TOTAL_TIME/3600/24
20 GRAPH_Y TOT("(3)")
30 GRAPH_SY -LA("H+")
  -end
    -active                 true

[ktelfeya]

Hi Dr. Parkhurst,
Thank you for your response. It seems that I have to cast the new species in terms of sulfate in order to define them if I use the phreeqc database. Is that correct? Also, the dithionite is still dropping to essentially 0. Is there perhaps something I need to suppress so that loss of dithionite is only due to reaction with Fe?

SOLUTION_MASTER_SPECIES
    S             SO4-2            0     SO4             32
    S(6)          SO4-2            0     SO4           
    S(4)          SO3-2            0     S             
    S(3)          S2O4-2           0     S             

SOLUTION_SPECIES
SO4-2 = SO4-2
    log_k     0
    delta_h   -217.4 kJ
SO4-2 = SO3-2 + 0.5O2
    log_k     -46.6244
    delta_h   267.985 kJ
    -analytical_expression -13.771 0.00065102 -13330 4.7164 -208 0
2H+ + 2SO3-2 = S2O4-2 + H2O + 0.5O2
    log_k     -25.2076
    -analytical_expression -231.72 0.0020393 -7.1011 83.239 0.94155 0
END

RATES
    Dithionite_Fe
-start
10 rem M = current number of moles of dithionite
20 rem M0 = number of moles of dithionite initially present
30 rem k1 = PARM(1)
40 rf=PARM(1)*M
50 MOLES = rf*TIME
60 save MOLES
-end
END

SOLUTION 1
    temp      25
    pH        8.3
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Alkalinity 0.1 as HCO3-
    Na        0.15 charge
    S(3)      0.1
    S(4)      1e-006
    -water    0.01 # kg
END

EQUILIBRIUM_PHASES 1
    Fe(OH)3(a) 0 0.00144
END

USE solution 1

KINETICS 1
Dithionite_Fe
    -formula  Fe(OH)3  1
    -m        0.1
    -m0       0.1
    -parms    1.03e-005
    -tol      1e-008
-steps       86400 172800 259200 604800 950400 1382400 2073600 3715200
-step_divide 1
-runge_kutta 3
-bad_step_max 500

USER_GRAPH 1
    -axis_titles            "time" "Dithionite (M)" "pH"
    -chart_title            "Dithionite reduction of Fe"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  t
  -start
10 GRAPH_X TOTAL_TIME/3600/24
20 GRAPH_Y TOT("S(3)")
30 GRAPH_SY -LA("H+")
  -end
    -active                 true



Thank you!

[dlparkhurst]

I apologize, what got posted was pretty useless. It dropped the brackets from the definitions, so I am attaching a file.  It defines a limited sulfur system with just S(4) and S(3), but defined as the element bracket S bracket. I adjusted the rate a little bit to fit the time series.

I did not finish the definitions in SOLUTION_SPECIES to include any protonation reactions. These may be important for the simulated pH.

I also included EQUILIBRIUM_PHASES to have the reaction stop when Fe(OH)3(a) became saturated.