Simplification of the Oxidation kinetics.

[mloche]

Hi, it’s been some years I’m using PHREEQC and today is my first post on the forum. I don’t think I saw any topics specifically on my issue.

I am currently developing a weathering model. This model will kinetically dissolve a set of primary minerals, as well as kinetically oxidizing the Fe and Mn, and finally precipitating any supersaturated phases at each time step. We are using the BRGM Thermoddem database, which contain a lot of phases.
For the oxidation kinetics, in the PHREEQC example, it’s explained that we need to define again the Master species and solutions species as well as custom phases (here the Goethite), giving new names to these variables.

For example, here is the code following theses rules for the oxidation of Fe.
You can see that in a system containing numerous elements, the number of solution species we need to define again is huge. In addition to this, we will have to define the Mn as well as hundred of minerals with this method.

------------------------------------------------------------------------------------------------------------------------------------------------------------------

TITLE Fe Oxidation
SOLUTION_MASTER_SPECIES
Fe_di              Fe_di+2    0.0     Fe_di              55.847
Fe_tri              Fe_tri+3    0.0     Fe_tri              55.847
SOLUTION_SPECIES
Fe_di+2 = Fe_di+2
        log_k   0.0
Fe_tri+3 = Fe_tri+3
        log_k   0.0
#
# Fe+2 species
#
Fe_di+2 + H2O = Fe_diOH+ + H+
        log_k   -9.5
        delta_h 13.20   kcal
#
#... and also other Fe+2 species
#
Fe_di+2 + Cl- = Fe_diCl+
        log_k   0.14
Fe_di+2 + CO3-2 = Fe_diCO3
        log_k   4.38
Fe_di+2 + HCO3- = Fe_diHCO3+
        log_k   2.0
Fe_di+2 + SO4-2 = Fe_diSO4
        log_k   2.25
        delta_h 3.230   kcal
Fe_di+2 + HSO4- = Fe_diHSO4+
        log_k   1.08
Fe_di+2 + 2HS- = Fe_di(HS)2
        log_k   8.95
Fe_di+2 + 3HS- = Fe_di(HS)3-
        log_k   10.987
Fe_di+2 + HPO4-2 = Fe_diHPO4
        log_k   3.6
Fe_di+2 + H2PO4- = Fe_diH2PO4+
        log_k   2.7
Fe_di+2 + F- = Fe_diF+
        log_k   1.0
#
# Fe+3 species
#
Fe_tri+3 + H2O = Fe_triOH+2 + H+
        log_k   -2.19
        delta_h 10.4    kcal
#
#... and also other Fe+3 species
#
Fe_tri+3 + 2 H2O = Fe_tri(OH)2+ + 2 H+
        log_k   -5.67
        delta_h 17.1    kcal
Fe_tri+3 + 3 H2O = Fe_tri(OH)3 + 3 H+
        log_k   -12.56
        delta_h 24.8    kcal
Fe_tri+3 + 4 H2O = Fe_tri(OH)4- + 4 H+
        log_k   -21.6
        delta_h 31.9    kcal
2 Fe_tri+3 + 2 H2O = Fe_tri2(OH)2+4 + 2 H+
        log_k   -2.95
        delta_h 13.5    kcal
3 Fe_tri+3 + 4 H2O = Fe_tri3(OH)4+5 + 4 H+
        log_k   -6.3
        delta_h 14.3    kcal
Fe_tri+3 + Cl- = Fe_triCl+2
        log_k   1.48
        delta_h 5.6     kcal
Fe_tri+3 + 2 Cl- = Fe_triCl2+
        log_k   2.13
Fe_tri+3 + 3 Cl- = Fe_triCl3
        log_k   1.13
Fe_tri+3 + SO4-2 = Fe_triSO4+
        log_k   4.04
        delta_h 3.91    kcal
Fe_tri+3 + HSO4- = Fe_triHSO4+2
        log_k   2.48
Fe_tri+3 + 2 SO4-2 = Fe_tri(SO4)2-
        log_k   5.38
        delta_h 4.60    kcal
Fe_tri+3 + HPO4-2 = Fe_triHPO4+
        log_k   5.43
        delta_h 5.76    kcal
Fe_tri+3 + H2PO4- = Fe_triH2PO4+2
        log_k   5.43
Fe_tri+3 + F- = Fe_triF+2
        log_k   6.2
        delta_h 2.7     kcal
Fe_tri+3 + 2 F- = Fe_triF2+
        log_k   10.8
        delta_h 4.8     kcal
Fe_tri+3 + 3 F- = Fe_triF3
        log_k   14.0
        delta_h 5.4     kcal
PHASES
Goethite
        Fe_triOOH + 3 H+ = Fe_tri+3 + 2 H2O
        log_k   -1.0
# Fake phase used to fix pH
FIX_H+
      H+ = H+
      logK 0

RATES
# rate (Singer & Stumm 1970)
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
KINETICS 1
Fe_di_ox
        -Formula  Fe_di  -1.0   Fe_tri  1.0
        -steps 3.15e7 in 100 steps
        -step_divide 1e-4

SOLUTION 1 
temp 20
Fe_di 1.01e-3 #Same as Tabata 2021

equilibrium_phases
O2(g) -3.12 # pH 6.5 500mb 1500 ppm
CO2(g) -0.75
# Add NaOH to pH 10
FIX_H+  -6.5 NaOH 10

INCREMENTAL_REACTIONS true
------------------------------------------------------------------------------------------------------------------------------------------------------------------
The idea here is to simplify this system to use the already built-in variable for the different master species, solutions species, and phases. Doing this way, we would be able to inject directly the Fe(3) at each time step for speciation and precipitation calculation. The following code is an example of this “simplified code”, which should give the same results as the first one except we don’t need to declare all the solutions species, master species and phases.

TITLE Iron oxidation

PHASES
# Fake phase used to fix pH
FIX_H+
      H+ = H+
      logK 0


RATES
# rate (Singer & Stumm 1970)
Fe_(2)_ox
-start
10  Fe_di = TOT("Fe(2)")
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

KINETICS 1
Fe_(2)_ox
        -Formula Fe(2)  -1.0   Fe(3)  1.0
        -steps 3.15e7 in 100 steps
        -step_divide 1e-4

SOLUTION 1 
temp 20
Fe(2) 1.01e-3 #Same as Tabata 2021


equilibrium_phases
O2(g) -3.12 # pH 6.5 500mb 1500 ppm
CO2(g) -0.75
# Add NaOH to pH 10
FIX_H+  -6.5 NaOH 10


INCREMENTAL_REACTIONS true
END

------------------------------------------------------------------------------------------------------------------------------------------------------------------

Unfortunately, this code is not working properly, returning the following error:

ERROR: Parsing error in get_elts_in_species, unexpected character, 2.
ERROR: Error in -formula: Fe(2)
ERROR: Parsing error in get_elts_in_species, unexpected character, 3.
ERROR: Error in -formula: Fe(3)
ERROR: Program stopping due to input errors.

My understanding is that without declaring again new master species, solution species and phases, PHREEQC seems unable to perform Kinetic calculation and then using the thermodynamic database, not recognizing the built in variable such as  Fe(2) and Fe(3). 
Does anyone know what I am doing wrong (or what I don’t understand!), and how to make it work, if possible!

Thanks you !

[dlparkhurst]

PHREEQC does not distinguish valence state in the -formula for a kinetic reaction, so you get the error message about the unexpected character.

Unfortunately, you will have to use the first method that you present, where iron is separated into two valence states. You will need to define all iron phases with reactions involving Fe_di and Fe_tri, and use these two "elements" in your -formula definitions.

You probably do not need to do as much for Mn. Mn+2 is the predominant aqueous species under most conditions, and, in an oxidizing environment, Mn can precipitate as Mn(IV) minerals. However, you can simply use KINETICS for the precipitation of the minerals, without having to define separate aqueous valence states.