Mn(II)-Ligand oxidation (pH/O2 impact) as a kinetic model

[mymyhihj]

Dear all Phreeqc users,

I'm quite new to Phreeqc, and I'm trying to build a kinetic model for the Mn(II)-ligand oxidation with the effect of pH and O2.
The batch data suggest Mn(II)-A can slowly be oxidized to Mn(III)-A at pH 8. Without ligand, I can successfully build a kinetic model for the oxidation Mn(II) to Mn(III). However, with the presence of ligands A, my graph just shows straight lines and the output just describes the Mn(II)-A and Mn(II) species without any information on Mn(III)-A and Mn(III) species. 
Please find my code below:

SOLUTION_MASTER_SPECIES
    A            A-3              0     C25H46N6O8        560.7   # DFOB ligand
    Mn_di        Mn_di+2          0.0   Mn_di             54.94   # Manganese (II) with atomic weight
    Mn_tri       Mn_tri+3         0.0   Mn_tri            54.94   # Manganese (III) with atomic weight

SOLUTION_SPECIES
A-3 = A-3
    log_k     0
A-3 + H+ = HA-2
    log_k     11
H+ + HA-2 = H2A-
    log_k     9.7
H+ + H2A- = H3A
    log_k     9.2
H+ + H3A = H4A+
    log_k     8.4
Mn_di+2 = Mn_di+2
    log_k     0
Mn_tri+3 = Mn_tri+3
    log_k     0
A-3 + Mn_di+2 = Mn_diA-
    log_k     7.5
HA-2 + Mn_di+2 = Mn_diHA
    log_k     6.8
H2A- + Mn_di+2 = Mn_diH2A+
    log_k     5.2
HA-2 + Mn_tri+3 = Mn_triHA+
    log_k     28.6

# Mn+2 species
H2O + Mn_di+2 = Mn_diOH+ + H+
    log_k     -10.6
    delta_h   55.81 kcal
3H2O + Mn_di+2 = Mn_di(OH)3- + 3H+
    log_k     -34.8

#... and also other Mn+2 species
Cl- + Mn_di+2 = Mn_diCl+
    log_k     0.1
Mn_di+2 + NO3- = Mn_diNO3+
    log_k     0.2
CO3-2 + H+ + Mn_di+2 = Mn_diHCO3+
    log_k     11.629
CO3-2 + Mn_di+2 = Mn_diCO3
    log_k     10.58
Mn_di+2 + SO4-2 = Mn_diSO4
    log_k     -2.58
    delta_h   3.23 kcal
HSO4- + Mn_di+2 = Mn_diHSO4+
    log_k     1.08
HS- + Mn_di+2 = Mn_diS + H+
    log_k     -0.17
3Mn_di+2 + 2PO4-3 = Mn_di3(PO4)2
    log_k     23.8
F- + Mn_di+2 = Mn_diF+
    log_k     1.6

# Mn+3 species
H2O + Mn_tri+3 = Mn_triOH+2 + H+
    log_k     1.0
    #delta_h   10.4 kcal

#... and also other Mn+3 species
2Mn_tri+3 + 3SO4-2 = Mn_tri2(SO4)3
    log_k     5.711
Mn_tri+3 + H2O = Mn_triOH+2 + H+
        log_k   -0.12
Mn_tri+3 + 2 H2O = Mn_tri(OH)2+ + 2 H+
        log_k   0.1
        #delta_h 17.1    kcal
Mn_tri+3 + 3 H2O = Mn_tri(OH)3 + 3 H+
        log_k   -6.3465
        #delta_h 24.8    kcal
Mn_tri+3 + H2A- = Mn_triH2A+2
    log_k     29.6
Mn_tri+3 + A-3 = Mn_triA
    log_k     18.6

CALCULATE_VALUES
  k_homogeneous
  -start
  10 SAVE    1.33e12   
  -end

RATES
Mn_diHA_ox
  -start
  20  Mn_diHA = TOT("Mn_diHA")
  30  pP_O2 = 10^(SI("O2(g)"))
  40  rate = calc_value("k_homogeneous") * (ACT("OH-"))^2 * pP_O2 * Mn_diHA
  50  if rate > Mn_diHA then rate = Mn_diHA
  60  if rate > 0.25*mol("O2") then rate = 0.25*mol("O2")
  70  put(rate,3)
  80  moles = rate*time
  100 SAVE moles
  -end

# Initial condition: 100 ?M Mn(II) and 250 ?M DFOB at pH 8.0, open air
SOLUTION 1
    temp      25
    pH        8.0
    pe        10 O2(g) -0.67
    redox     pe
    units     mmol/kgw
    density   1
    A         0.25
    Cl        10.0 charge
    Mn_di     0.1
    Na        10.0
    #-water    1 # kg

KINETICS 1
Mn_diHA_ox
    -formula  Mn_diHA  -1.0 Mn_triHA+  1.0
    -cvode
    -steps 1e-2 1e-1 1e0 1e1 1e2 1e3 1e4 1e5 1e6 1e7 0.3e8 0.6e8 0.9e8 1e8 1.2e8
INCREMENTAL_REACTIONS true

USER_GRAPH 1 Example 9
    -headings               _time_ Mn_diHA Mn_triHA O2 pH pe
    -axis_titles            "Time, in days" "Micromole per kilogram water" "pH"
    -chart_title            "Oxidation of Mn(II)-DFOB to Mn(III)-DFOB"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600 / 24 
20 GRAPH_Y TOT("Mn_diHA")*1e6, TOT("Mn_triHA+")*1e6, TOT("O2")
30 GRAPH_SY -LA("H+")   
40 GRAPH_SY -LA("e-")   
  -end
    -active                 true

END

I would like to ask for help if anyone knows how to deal with this issue. Any advice is appreciated, and thanks in advance!

Best regards,
My

[dlparkhurst]

I think you were really close. The function TOT only applies to total concentrations of "elements", Mn_di and Mn_tri in your case. The function MOL applies to aqueous (, exchange, and surface) species, like Mn_diHA. So if you replace TOT with MOL in a few places, I think you are good.

I have simplified the graph a little just to emphasize the transition from Mn_di to Mn_tri. Also, I try not to use if statements in RATES to that introduce discontinuities or non-differentiable points. It is better to have smooth, differentiable functions. If you want to have a reaction stop because of lack of oxygen, then a Michalis-Menton factor could be added (O2/(1e-6 + O2), for example).

Code: [Select]

SOLUTION_MASTER_SPECIES
    A            A-3              0     C25H46N6O8        560.7   # DFOB ligand
    Mn_di        Mn_di+2          0.0   Mn_di             54.94   # Manganese (II) with atomic weight
    Mn_tri       Mn_tri+3         0.0   Mn_tri            54.94   # Manganese (III) with atomic weight

SOLUTION_SPECIES
A-3 = A-3
    log_k     0
A-3 + H+ = HA-2
    log_k     11
H+ + HA-2 = H2A-
    log_k     9.7
H+ + H2A- = H3A
    log_k     9.2
H+ + H3A = H4A+
    log_k     8.4
Mn_di+2 = Mn_di+2
    log_k     0
Mn_tri+3 = Mn_tri+3
    log_k     0
A-3 + Mn_di+2 = Mn_diA-
    log_k     7.5
HA-2 + Mn_di+2 = Mn_diHA
    log_k     6.8
H2A- + Mn_di+2 = Mn_diH2A+
    log_k     5.2
HA-2 + Mn_tri+3 = Mn_triHA+
    log_k     28.6

# Mn+2 species
H2O + Mn_di+2 = Mn_diOH+ + H+
    log_k     -10.6
    delta_h   55.81 kcal
3H2O + Mn_di+2 = Mn_di(OH)3- + 3H+
    log_k     -34.8

#... and also other Mn+2 species
Cl- + Mn_di+2 = Mn_diCl+
    log_k     0.1
Mn_di+2 + NO3- = Mn_diNO3+
    log_k     0.2
CO3-2 + H+ + Mn_di+2 = Mn_diHCO3+
    log_k     11.629
CO3-2 + Mn_di+2 = Mn_diCO3
    log_k     10.58
Mn_di+2 + SO4-2 = Mn_diSO4
    log_k     -2.58
    delta_h   3.23 kcal
HSO4- + Mn_di+2 = Mn_diHSO4+
    log_k     1.08
HS- + Mn_di+2 = Mn_diS + H+
    log_k     -0.17
3Mn_di+2 + 2PO4-3 = Mn_di3(PO4)2
    log_k     23.8
F- + Mn_di+2 = Mn_diF+
    log_k     1.6

# Mn+3 species
H2O + Mn_tri+3 = Mn_triOH+2 + H+
    log_k     1.0
    #delta_h   10.4 kcal

#... and also other Mn+3 species
2Mn_tri+3 + 3SO4-2 = Mn_tri2(SO4)3
    log_k     5.711
Mn_tri+3 + H2O = Mn_triOH+2 + H+
        log_k   -0.12
Mn_tri+3 + 2 H2O = Mn_tri(OH)2+ + 2 H+
        log_k   0.1
        #delta_h 17.1    kcal
Mn_tri+3 + 3 H2O = Mn_tri(OH)3 + 3 H+
        log_k   -6.3465
        #delta_h 24.8    kcal
Mn_tri+3 + H2A- = Mn_triH2A+2
    log_k     29.6
Mn_tri+3 + A-3 = Mn_triA
    log_k     18.6

CALCULATE_VALUES
  k_homogeneous
  -start
  10 SAVE    1.33e12   
  -end

RATES
Mn_diHA_ox
  -start
  #20  Mn_diHA = TOT("Mn_diHA")
  20  Mn_diHA = MOL("Mn_diHA")
  30  pP_O2 = 10^(SI("O2(g)"))
  40  rate = calc_value("k_homogeneous") * (ACT("OH-"))^2 * pP_O2 * Mn_diHA
  #50  if rate > Mn_diHA then rate = Mn_diHA
  #60  if rate > 0.25*mol("O2") then rate = 0.25*mol("O2")
  70  put(rate,3)
  80  moles = rate*time
  100 SAVE moles
110 PRINT moles, rate, pP_O2, Mn_diHA
  -end

# Initial condition: 100 ?M Mn(II) and 250 ?M DFOB at pH 8.0, open air
SOLUTION 1
    temp      25
    pH        8.0
    pe        10 O2(g) -0.67
    redox     pe
    units     mmol/kgw
    density   1
    A         0.25
    Cl        10.0 charge
    Mn_di     0.1
    Na        10.0
    #-water    1 # kg

KINETICS 1
Mn_diHA_ox
    #-formula  Mn_diHA  -1.0 Mn_triHA+  1.0
    -formula  Mn_di  -1.0 Mn_tri  1.0
    -cvode
    -steps 1e-2 1e-1 1e0 1e1 1e2 1e3 1e4 1e5 1e6 1e7 0.3e8 0.6e8 0.9e8 1e8 1.2e8
INCREMENTAL_REACTIONS true

USER_GRAPH 1 Example 9
    #-headings               _time_ Mn_di Mn_tri Mn_diHA Mn_triHA O2 pH pe
    -headings               time Mn_di Mn_tri pH PE
    -axis_titles            "Time, in days" "Micromole per kilogram water" "pH"
    -chart_title            "Oxidation of Mn(II)-DFOB to Mn(III)-DFOB"
    -axis_scale x_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / 3600 / 24
20 GRAPH_Y TOT("Mn_di"), TOT("Mn_tri") #, MOL("Mn_diHA")*1e6, MOL("Mn_triHA+")*1e6, TOT("O2")
30 GRAPH_SY -LA("H+")   
40 GRAPH_SY -LA("e-")   
  -end
    -active                 true

END


[mymyhihj]

Many many thanks Mr. Parkhurst, I finally understand how it works!