Simulation of nitrification and denitrification in porous soil

[Rensburg]

Dear phreeqcuser,
  I am currently working on solute transport modeling for agricultural environments. The model involves four substances: Amm, N, BOD (which is a new element I defined), and dissolved O2. However, I have encountered some issues and would like to seek your advice. For clarity, I have significantly simplified the model.The model is expressed mathematically as follows:

Code: [Select]

dAmm/dt = k1*rate(nitrification) + k2*rate(plant respiration) + ...
dN/dt = k3*rate(nitrification) + k4*rate(denitrification) + k5*rate(plant growth) + ...
dO2/dt = k6*rate(atmosphere) + k7*rate(plant respiration) + k8*rate(nitrification) + k9*rate(BOD degradation) + k10*rate(nitrification) + ...
dBOD/dt = k11*rate(plant death) + k12*rate(denitrification) + ...  Can I define my PHREEQC functions based on the rate of change of a specific substance, similar to the example below, rather than defining KINETICS for a chemical reaction as in the official examples ex9 and ex15?

Code: [Select]

RATES
  Amm_t
  -start
  rate = dAmm/dt
  ...
  -end
  N_t
  -start
  rate = dN/dt
  ...
  -end
  O2_t
  -start
  rate = dO2/dt
  ...
  -end
BOD_t
  -start
  rate = dBOD/dt
  ...
  -end
KINETICS
  Amm_t
  -formula  Amm 1
  N_t
  -formula  N 1
  O2_t
  -formula  O 2
  BOD_t
  -formula  BOD 1
  What's more, since the model involves dissolved oxygen, and the kinetics module treats both O-2 and O2 as O, how should I differentiate between them?For example, dissolved oxygen participates in the nitrification reaction: NH4+ + 2O2 = NO3- + H2O + 2H+. The dissolved oxygen should be reduced to O(-2). However, if I define it in KINETICS as "-formula O 2 O -2", the O cancels out. I don’t know how to represent this change. I think using "Oxg" is not appropriate, because the speciation of O affects the transformation of other elements such as Fe(2), Mn(-7), etc. If I use Oxg, this coupling would not be captured. Additionally, plant photosynthesis produces dissolved oxygen. If I directly use "-formula O 2" in KINETICS, will it add O(-2) to the solution instead of dissolved oxygen?Is there a good solution to distinguish dissolved oxygen while preserving its reactivity with other ions?

   In the SOLUTION module, does defining “O 10” mean the presence of dissolved oxygen in the water?

  Thank you for your help.

[dlparkhurst]

If you want to model all of the redox states of all of the redox elements with KINETICS, you will have to define a separate "element" for each redox state. You can use this database for example.

Code: [Select]

SOLUTION_MASTER_SPECIES
    H             H+              -1     H               1.008
    H(0)          H2               0     H             
    H(1)          H+              -1     0             
    Hzero         Hzero2           0     Hzero           1
    E             e-               0     0               1
    O             H2O              0     O               16
    O(0)          O2               0     O             
    O(-2)         H2O              0     0             
    Ozero         Ozero2           0     Ozero           16
    Ca            Ca+2             0     Ca              40.08
    Mg            Mg+2             0     Mg              24.312
    Na            Na+              0     Na              22.9898
    K             K+               0     K               39.102
    Ferrous       Ferrous+2        0     Ferrous         55.847
    Ferric        Ferric+3        -2     Ferric          55.847
    Manganous     Manganous+2      0     Manganous       54.938
    Manganic      Manganic+3       0     Manganic        54.938
    Al            Al+3             0     Al              26.9815
    Ba            Ba+2             0     Ba              137.34
    Sr            Sr+2             0     Sr              87.62
    Si            H4SiO4           0     SiO2            28.0843
    Cl            Cl-              0     Cl              35.453
    C             CO3-2            2     HCO3            12.0111
    C(4)          CO3-2            2     HCO3            12.0111
    Methane       MethaneH4        0     MethaneH4       12.0111
    Alkalinity    CO3-2            1     Ca0.5(CO3)0.5   50.05
    Sulfate       SulfateO4-2      0     SulfateO4       32.064
    Sulfide       HSulfide-        1     Sulfide         32.064
    Nitrate       NitrateO3-       0     Nitrate         14.0067
    Nitrite       NitriteO2-       0     Nitrite         14.0067
    Nzero         Nzero2           0     Nzero           14.0067
    Amm           AmmH4+           0     Amm             14.0067
    B             H3BO3            0     B               10.81
    P             PO4-3            2     P               30.9738
    F             F-               0     F               18.9984
    Li            Li+              0     Li              6.939
    Br            Br-              0     Br              79.904
    Zn            Zn+2             0     Zn              65.37
    Cd            Cd+2             0     Cd              112.4
    Pb            Pb+2             0     Pb              207.19
    Cupric        Cupric+2         0     Cupric          63.546
    Cuprous       Cuprous+         0     Cuprous         63.546
SOLUTION_SPECIES
H+ = H+
    log_k     0
    -gamma    9 0
    -dw       9.31e-009
e- = e-
    log_k     0
H2O = H2O
    log_k     0
Ca+2 = Ca+2
    log_k     0
    -gamma    5 0.165
    -dw       7.93e-010
    -millero  -19.69 0.1058 -0.001256 1.617 -0.075 0.0008262
Mg+2 = Mg+2
    log_k     0
    -gamma    5.5 0.2
    -dw       7.05e-010
    -millero  -22.32 0.0868 -0.0016 2.017 -0.125 0.001457
Na+ = Na+
    log_k     0
    -gamma    4 0.075
    -dw       1.33e-009
    -millero  -3.46 0.1092 -0.000768 2.698 -0.106 0.001651
K+ = K+
    log_k     0
    -gamma    3.5 0.015
    -dw       1.96e-009
    -millero  7.26 0.0892 -0.000736 2.722 -0.101 0.00151
Ferrous+2 = Ferrous+2
    log_k     0
    -gamma    6 0
    -dw       7.19e-010
Manganous+2 = Manganous+2
    log_k     0
    -gamma    6 0
    -dw       6.88e-010
Al+3 = Al+3
    log_k     0
    -gamma    9 0
    -dw       5.59e-010
Ba+2 = Ba+2
    log_k     0
    -gamma    5 0
    -dw       8.48e-010
Sr+2 = Sr+2
    log_k     0
    -gamma    5.26 0.121
    -dw       7.94e-010
    -millero  -18.44 0.0082 -0.0006 1.727 -0.067 0.00084
H4SiO4 = H4SiO4
    log_k     0
    -dw       1.1e-009
    -millero  56 0 0 0 0 0
Cl- = Cl-
    log_k     0
    -gamma    3.5 0.015
    -dw       2.03e-009
    -millero  16.37 0.0896 -0.001264 -1.494 0.034 -0.000621
CO3-2 = CO3-2
    log_k     0
    -gamma    5.4 0
    -dw       9.55e-010
    -millero  -8.74 0.3 -0.004064 5.65 0 0
SulfateO4-2 = SulfateO4-2
    log_k     0
    -gamma    5 -0.04
    -dw       1.07e-009
    -millero  9.26 0.284 -0.003808 0.4348 -0.0099143 -8.4762e-005
NitrateO3- = NitrateO3-
    log_k     0
    -dw       1.9e-009
    -millero  25.51 0.1888 -0.001984 -0.654 0 0
H3BO3 = H3BO3
    log_k     0
    -dw       1.1e-009
    -millero  36.56 0.13 -0.00081 0 0 0
PO4-3 = PO4-3
    log_k     0
    -gamma    4 0
    -dw       6.12e-010
F- = F-
    log_k     0
    -gamma    3.5 0
    -dw       1.46e-009
    -millero  -3.05 0.3276 -0.00352 1.271 -0.074 8.857e-005
Li+ = Li+
    log_k     0
    -gamma    6 0
    -dw       1.03e-009
Br- = Br-
    log_k     0
    -gamma    3 0
    -dw       2.01e-009
    -millero  22.98 0.0934 -0.000968 -1.675 0.05 -0.001105
Zn+2 = Zn+2
    log_k     0
    -gamma    5 0
    -dw       7.15e-010
Cd+2 = Cd+2
    log_k     0
    -dw       7.17e-010
Pb+2 = Pb+2
    log_k     0
    -dw       9.45e-010
Cupric+2 = Cupric+2
    log_k     0
    -gamma    6 0
    -dw       7.33e-010
H2O = OH- + H+
    log_k     -14
    delta_h   13.362 kcal
    -analytical_expression -283.971 -0.05069842 13323 102.24447 -1119669 0
    -gamma    3.5 0
    -dw       5.27e-009
2H2O = O2 + 4H+ + 4e-
    log_k     -86.08
    delta_h   134.79 kcal
    -dw       2.35e-009
Ozero2 = Ozero2
    log_k     0
    -dw       2.35e-009
2H+ + 2e- = H2
    log_k     -3.15
    delta_h   -1.759 kcal
    -dw       5.1e-009
Hzero2 = Hzero2
    log_k     0
    -dw       5.1e-009
H+ + CO3-2 = HCO3-
    log_k     10.329
    delta_h   -3.561 kcal
    -analytical_expression 107.8871 0.03252849 -5151.79 -38.92561 563713.9 0
    -gamma    5.4 0
    -dw       1.18e-009
    -millero  21.07 0.185 -0.002248 2.29 -0.006644 -3.667e-006
2H+ + CO3-2 = CO2 + H2O
    log_k     16.681
    delta_h   -5.738 kcal
    -analytical_expression 464.1965 0.09344813 -26986.16 -165.75951 2248628.9 0
    -dw       1.92e-009
MethaneH4 = MethaneH4
    log_k     0
    delta_h   -61.039 kcal
    -dw       1.85e-009
H+ + SulfateO4-2 = HSulfateO4-
    log_k     1.988
    delta_h   3.85 kcal
    -analytical_expression -56.889 0.006473 2307.9 19.8858 0 0
    -dw       1.33e-009
HSulfide- = HSulfide-
    log_k     0
    -gamma    3.5 0
    -dw       1.73e-009
HSulfide- = Sulfide-2 + H+
    log_k     -12.918
    delta_h   12.1 kcal
    -gamma    5 0
    -dw       7.31e-010
H+ + HSulfide- = H2Sulfide
    log_k     6.994
    delta_h   -5.3 kcal
    -analytical_expression -11.17 0.02386 3279 0 0 0
    -dw       2.1e-009
NitriteO2- = NitriteO2-
    log_k     0
    -gamma    3 0
    -dw       1.91e-009
Nzero2 = Nzero2
    log_k     0
    -dw       1.96e-009
AmmH4+ = AmmH3 + H+
    log_k     -9.252
    delta_h   12.48 kcal
    -analytical_expression 0.6322 -0.001225 -2835.76 0 0 0
    -dw       2.28e-009
AmmH4+ = AmmH4+
    log_k     0
    -gamma    2.5 0
    -dw       1.98e-009
    -millero  17.47 -0.0034 0.00076 0 0 0
AmmH4+ + SulfateO4-2 = AmmH4SulfateO4-
    log_k     1.11
H3BO3 = H2BO3- + H+
    log_k     -9.24
    delta_h   3.224 kcal
F- + H3BO3 = BF(OH)3-
    log_k     -0.4
    delta_h   1.85 kcal
2F- + H+ + H3BO3 = BF2(OH)2- + H2O
    log_k     7.63
    delta_h   1.618 kcal
3F- + 2H+ + H3BO3 = BF3OH- + 2H2O
    log_k     13.67
    delta_h   -1.614 kcal
4F- + 3H+ + H3BO3 = BF4- + 3H2O
    log_k     20.28
    delta_h   -1.846 kcal
H+ + PO4-3 = HPO4-2
    log_k     12.346
    delta_h   -3.53 kcal
    -gamma    4 0
    -dw       6.9e-010
2H+ + PO4-3 = H2PO4-
    log_k     19.553
    delta_h   -4.52 kcal
    -gamma    4.5 0
    -dw       8.46e-010
    -millero  33.6 0 0 0 0 0
F- + H+ = HF
    log_k     3.18
    delta_h   3.18 kcal
    -analytical_expression -2.033 0.012645 429.01 0 0 0
2F- + H+ = HF2-
    log_k     3.76
    delta_h   4.55 kcal
Ca+2 + H2O = CaOH+ + H+
    log_k     -12.78
Ca+2 + CO3-2 = CaCO3
    log_k     3.224
    delta_h   3.545 kcal
    -analytical_expression -1228.732 -0.29944 35512.75 485.818 0 0
    -dw       4.46e-010
Ca+2 + H+ + CO3-2 = CaHCO3+
    log_k     11.435
    delta_h   -0.871 kcal
    -analytical_expression 1317.0071 0.34546894 -39916.84 -517.70761 563713.9 0
    -gamma    5.4 0
    -dw       5.06e-010
Ca+2 + SulfateO4-2 = CaSulfateO4
    log_k     2.3
    delta_h   1.65 kcal
    -dw       4.71e-010
Ca+2 + HSulfateO4- = CaHSulfateO4+
    log_k     1.08
Ca+2 + PO4-3 = CaPO4-
    log_k     6.459
    delta_h   3.1 kcal
Ca+2 + HPO4-2 = CaHPO4
    log_k     2.739
    delta_h   3.3 kcal
Ca+2 + H2PO4- = CaH2PO4+
    log_k     1.408
    delta_h   3.4 kcal
Ca+2 + F- = CaF+
    log_k     0.94
    delta_h   4.12 kcal
H2O + Mg+2 = MgOH+ + H+
    log_k     -11.44
    delta_h   15.952 kcal
Mg+2 + CO3-2 = MgCO3
    log_k     2.98
    delta_h   2.713 kcal
    -analytical_expression 0.991 0.00667 0 0 0 0
H+ + Mg+2 + CO3-2 = MgHCO3+
    log_k     11.399
    delta_h   -2.771 kcal
    -analytical_expression 48.6721 0.03252849 -2614.335 -18.00263 563713.9 0
Mg+2 + SulfateO4-2 = MgSulfateO4
    log_k     2.37
    delta_h   4.55 kcal
Mg+2 + PO4-3 = MgPO4-
    log_k     6.589
    delta_h   3.1 kcal
HPO4-2 + Mg+2 = MgHPO4
    log_k     2.87
    delta_h   3.3 kcal
H2PO4- + Mg+2 = MgH2PO4+
    log_k     1.513
    delta_h   3.4 kcal
F- + Mg+2 = MgF+
    log_k     1.82
    delta_h   3.2 kcal
H2O + Na+ = NaOH + H+
    log_k     -14.18
Na+ + CO3-2 = NaCO3-
    log_k     1.27
    delta_h   8.91 kcal
    -dw       5.85e-010
HCO3- + Na+ = NaHCO3
    log_k     -0.25
    -dw       6.73e-010
Na+ + SulfateO4-2 = NaSulfateO4-
    log_k     0.7
    delta_h   1.12 kcal
    -dw       6.18e-010
HPO4-2 + Na+ = NaHPO4-
    log_k     0.29
F- + Na+ = NaF
    log_k     -0.24
H2O + K+ = KOH + H+
    log_k     -14.46
K+ + SulfateO4-2 = KSulfateO4-
    log_k     0.85
    delta_h   2.25 kcal
    -analytical_expression 3.106 0 -673.6 0 0 0
    -dw       7.46e-010
HPO4-2 + K+ = KHPO4-
    log_k     0.29
Ferrous+2 + H2O = FerrousOH+ + H+
    log_k     -9.5
    delta_h   13.2 kcal
Cl- + Ferrous+2 = FerrousCl+
    log_k     0.14
Ferrous+2 + CO3-2 = FerrousCO3
    log_k     4.38
Ferrous+2 + HCO3- = FerrousHCO3+
    log_k     2
Ferrous+2 + SulfateO4-2 = FerrousSulfateO4
    log_k     2.25
    delta_h   3.23 kcal
Ferrous+2 + HSulfateO4- = FerrousHSulfateO4+
    log_k     1.08
Ferrous+2 + 2HSulfide- = Ferrous(HSulfide)2
    log_k     8.95
Ferrous+2 + 3HSulfide- = Ferrous(HSulfide)3-
    log_k     10.987
Ferrous+2 + HPO4-2 = FerrousHPO4
    log_k     3.6
Ferrous+2 + H2PO4- = FerrousH2PO4+
    log_k     2.7
F- + Ferrous+2 = FerrousF+
    log_k     1
Ferric+3 = Ferric+3
    log_k     0
    -gamma    9 0
Ferric+3 + H2O = FerricOH+2 + H+
    log_k     -2.19
    delta_h   10.4 kcal
Ferric+3 + 2H2O = Ferric(OH)2+ + 2H+
    log_k     -5.67
    delta_h   17.1 kcal
Ferric+3 + 3H2O = Ferric(OH)3 + 3H+
    log_k     -12.56
    delta_h   24.8 kcal
Ferric+3 + 4H2O = Ferric(OH)4- + 4H+
    log_k     -21.6
    delta_h   31.9 kcal
2Ferric+3 + 2H2O = Ferric2(OH)2+4 + 2H+
    log_k     -2.95
    delta_h   13.5 kcal
3Ferric+3 + 4H2O = Ferric3(OH)4+5 + 4H+
    log_k     -6.3
    delta_h   14.3 kcal
Cl- + Ferric+3 = FerricCl+2
    log_k     1.48
    delta_h   5.6 kcal
2Cl- + Ferric+3 = FerricCl2+
    log_k     2.13
3Cl- + Ferric+3 = FerricCl3
    log_k     1.13
Ferric+3 + SulfateO4-2 = FerricSulfateO4+
    log_k     4.04
    delta_h   3.91 kcal
Ferric+3 + HSulfateO4- = FerricHSulfateO4+2
    log_k     2.48
Ferric+3 + 2SulfateO4-2 = Ferric(SulfateO4)2-
    log_k     5.38
    delta_h   4.6 kcal
Ferric+3 + HPO4-2 = FerricHPO4+
    log_k     5.43
    delta_h   5.76 kcal
Ferric+3 + H2PO4- = FerricH2PO4+2
    log_k     5.43
F- + Ferric+3 = FerricF+2
    log_k     6.2
    delta_h   2.7 kcal
2F- + Ferric+3 = FerricF2+
    log_k     10.8
    delta_h   4.8 kcal
3F- + Ferric+3 = FerricF3
    log_k     14
    delta_h   5.4 kcal
H2O + Manganous+2 = ManganousOH+ + H+
    log_k     -10.59
    delta_h   14.4 kcal
Cl- + Manganous+2 = ManganousCl+
    log_k     0.61
2Cl- + Manganous+2 = ManganousCl2
    log_k     0.25
3Cl- + Manganous+2 = ManganousCl3-
    log_k     -0.31
Manganous+2 + CO3-2 = ManganousCO3
    log_k     4.9
HCO3- + Manganous+2 = ManganousHCO3+
    log_k     1.95
Manganous+2 + SulfateO4-2 = ManganousSulfateO4
    log_k     2.25
    delta_h   3.37 kcal
Manganous+2 + 2NitrateO3- = Manganous(NitrateO3)2
    log_k     0.6
    delta_h   -0.396 kcal
F- + Manganous+2 = ManganousF+
    log_k     0.84
Manganic+3 = Manganic+3
    log_k     0
Al+3 + H2O = AlOH+2 + H+
    log_k     -5
    delta_h   11.49 kcal
    -analytical_expression -38.253 0 -656.27 14.327 0 0
Al+3 + 2H2O = Al(OH)2+ + 2H+
    log_k     -10.1
    delta_h   26.9 kcal
    -analytical_expression 88.5 0 -9391.6 -27.121 0 0
Al+3 + 3H2O = Al(OH)3 + 3H+
    log_k     -16.9
    delta_h   39.89 kcal
    -analytical_expression 226.374 0 -18247.8 -73.597 0 0
Al+3 + 4H2O = Al(OH)4- + 4H+
    log_k     -22.7
    delta_h   42.3 kcal
    -analytical_expression 51.578 0 -11168.9 -14.865 0 0
Al+3 + SulfateO4-2 = AlSulfateO4+
    log_k     3.5
    delta_h   2.29 kcal
Al+3 + 2SulfateO4-2 = Al(SulfateO4)2-
    log_k     5
    delta_h   3.11 kcal
Al+3 + HSulfateO4- = AlHSulfateO4+2
    log_k     0.46
Al+3 + F- = AlF+2
    log_k     7
    delta_h   1.06 kcal
Al+3 + 2F- = AlF2+
    log_k     12.7
    delta_h   1.98 kcal
Al+3 + 3F- = AlF3
    log_k     16.8
    delta_h   2.16 kcal
Al+3 + 4F- = AlF4-
    log_k     19.4
    delta_h   2.2 kcal
Al+3 + 5F- = AlF5-2
    log_k     20.6
    delta_h   1.84 kcal
Al+3 + 6F- = AlF6-3
    log_k     20.6
    delta_h   -1.67 kcal
H4SiO4 = H3SiO4- + H+
    log_k     -9.83
    delta_h   6.12 kcal
    -analytical_expression -302.3724 -0.050698 15669.69 108.18466 -1119669 0
H4SiO4 = H2SiO4-2 + 2H+
    log_k     -23
    delta_h   17.6 kcal
    -analytical_expression -294.0184 -0.07265 11204.49 108.18466 -1119669 0
6F- + 4H+ + H4SiO4 = SiF6-2 + 4H2O
    log_k     30.18
    delta_h   -16.26 kcal
Ba+2 + H2O = BaOH+ + H+
    log_k     -13.47
Ba+2 + CO3-2 = BaCO3
    log_k     2.71
    delta_h   3.55 kcal
    -analytical_expression 0.113 0.008721 0 0 0 0
Ba+2 + HCO3- = BaHCO3+
    log_k     0.982
    delta_h   5.56 kcal
    -analytical_expression -3.0938 0.013669 0 0 0 0
Ba+2 + SulfateO4-2 = BaSulfateO4
    log_k     2.7
H2O + Sr+2 = SrOH+ + H+
    log_k     -13.29
    -gamma    5 0
H+ + Sr+2 + CO3-2 = SrHCO3+
    log_k     11.509
    delta_h   2.489 kcal
    -analytical_expression 104.6391 0.04739549 -5151.79 -38.92561 563713.9 0
    -gamma    5.4 0
Sr+2 + CO3-2 = SrCO3
    log_k     2.81
    delta_h   5.22 kcal
    -analytical_expression -1.019 0.012826 0 0 0 0
Sr+2 + SulfateO4-2 = SrSulfateO4
    log_k     2.29
    delta_h   2.08 kcal
H2O + Li+ = LiOH + H+
    log_k     -13.64
Li+ + SulfateO4-2 = LiSulfateO4-
    log_k     0.64
Cuprous+ = Cuprous+
    log_k     0
    -gamma    2.5 0
Cupric+2 + H2O = CupricOH+ + H+
    log_k     -8
    -gamma    4 0
Cupric+2 + 2H2O = Cupric(OH)2 + 2H+
    log_k     -13.68
Cupric+2 + 3H2O = Cupric(OH)3- + 3H+
    log_k     -26.9
Cupric+2 + 4H2O = Cupric(OH)4-2 + 4H+
    log_k     -39.6
Cupric+2 + SulfateO4-2 = CupricSulfateO4
    log_k     2.31
    delta_h   1.22 kcal
H2O + Zn+2 = ZnOH+ + H+
    log_k     -8.96
    delta_h   13.4 kcal
2H2O + Zn+2 = Zn(OH)2 + 2H+
    log_k     -16.9
3H2O + Zn+2 = Zn(OH)3- + 3H+
    log_k     -28.4
4H2O + Zn+2 = Zn(OH)4-2 + 4H+
    log_k     -41.2
Cl- + Zn+2 = ZnCl+
    log_k     0.43
    delta_h   7.79 kcal
2Cl- + Zn+2 = ZnCl2
    log_k     0.45
    delta_h   8.5 kcal
3Cl- + Zn+2 = ZnCl3-
    log_k     0.5
    delta_h   9.56 kcal
4Cl- + Zn+2 = ZnCl4-2
    log_k     0.2
    delta_h   10.96 kcal
Zn+2 + CO3-2 = ZnCO3
    log_k     5.3
Zn+2 + 2CO3-2 = Zn(CO3)2-2
    log_k     9.63
HCO3- + Zn+2 = ZnHCO3+
    log_k     2.1
SulfateO4-2 + Zn+2 = ZnSulfateO4
    log_k     2.37
    delta_h   1.36 kcal
2SulfateO4-2 + Zn+2 = Zn(SulfateO4)2-2
    log_k     3.28
Cd+2 + H2O = CdOH+ + H+
    log_k     -10.08
    delta_h   13.1 kcal
Cd+2 + 2H2O = Cd(OH)2 + 2H+
    log_k     -20.35
Cd+2 + 3H2O = Cd(OH)3- + 3H+
    log_k     -33.3
Cd+2 + 4H2O = Cd(OH)4-2 + 4H+
    log_k     -47.35
Cd+2 + Cl- = CdCl+
    log_k     1.98
    delta_h   0.59 kcal
Cd+2 + 2Cl- = CdCl2
    log_k     2.6
    delta_h   1.24 kcal
Cd+2 + 3Cl- = CdCl3-
    log_k     2.4
    delta_h   3.9 kcal
Cd+2 + CO3-2 = CdCO3
    log_k     2.9
Cd+2 + 2CO3-2 = Cd(CO3)2-2
    log_k     6.4
Cd+2 + HCO3- = CdHCO3+
    log_k     1.5
Cd+2 + SulfateO4-2 = CdSulfateO4
    log_k     2.46
    delta_h   1.08 kcal
Cd+2 + 2SulfateO4-2 = Cd(SulfateO4)2-2
    log_k     3.5
H2O + Pb+2 = PbOH+ + H+
    log_k     -7.71
2H2O + Pb+2 = Pb(OH)2 + 2H+
    log_k     -17.12
3H2O + Pb+2 = Pb(OH)3- + 3H+
    log_k     -28.06
4H2O + Pb+2 = Pb(OH)4-2 + 4H+
    log_k     -39.7
H2O + 2Pb+2 = Pb2OH+3 + H+
    log_k     -6.36
Cl- + Pb+2 = PbCl+
    log_k     1.6
    delta_h   4.38 kcal
2Cl- + Pb+2 = PbCl2
    log_k     1.8
    delta_h   1.08 kcal
3Cl- + Pb+2 = PbCl3-
    log_k     1.7
    delta_h   2.17 kcal
4Cl- + Pb+2 = PbCl4-2
    log_k     1.38
    delta_h   3.53 kcal
Pb+2 + CO3-2 = PbCO3
    log_k     7.24
Pb+2 + 2CO3-2 = Pb(CO3)2-2
    log_k     10.64
HCO3- + Pb+2 = PbHCO3+
    log_k     2.9
Pb+2 + SulfateO4-2 = PbSulfateO4
    log_k     2.75
Pb+2 + 2SulfateO4-2 = Pb(SulfateO4)2-2
    log_k     3.47
NitrateO3- + Pb+2 = PbNitrateO3+
    log_k     1.17

PHASES
redoxCalcite
    CaCO3 = Ca+2 + CO3-2
    log_k     -8.48
    delta_h   -2.297 kcal
    -analytical_expression -171.9065 -0.077993 2839.319 71.595 0 0
redoxAragonite
    CaCO3 = Ca+2 + CO3-2
    log_k     -8.336
    delta_h   -2.589 kcal
    -analytical_expression -171.9773 -0.077993 2903.293 71.595 0 0
redoxDolomite
    CaMg(CO3)2 = Ca+2 + Mg+2 + 2CO3-2
    log_k     -17.09
    delta_h   -9.436 kcal
redoxSiderite
    FerrousCO3 = Ferrous+2 + CO3-2
    log_k     -10.89
    delta_h   -2.48 kcal
redoxRhodochrosite
    ManganousCO3 = Manganous+2 + CO3-2
    log_k     -11.13
    delta_h   -1.43 kcal
redoxStrontianite
    SrCO3 = Sr+2 + CO3-2
    log_k     -9.271
    delta_h   -0.4 kcal
    -analytical_expression 155.0305 0 -7239.594 -56.58638 0 0
redoxWitherite
    BaCO3 = Ba+2 + CO3-2
    log_k     -8.562
    delta_h   0.703 kcal
    -analytical_expression 607.642 0.121098 -20011.25 -236.4948 0 0
redoxGypsum
    CaSulfateO4:2H2O = Ca+2 + 2H2O + SulfateO4-2
    log_k     -4.58
    delta_h   -0.109 kcal
    -analytical_expression 68.2401 0 -3221.51 -25.0627 0 0
redoxAnhydrite
    CaSulfateO4 = Ca+2 + SulfateO4-2
    log_k     -4.36
    delta_h   -1.71 kcal
    -analytical_expression 197.52 0 -8669.8 -69.835 0 0
redoxCelestite
    SrSulfateO4 = Sr+2 + SulfateO4-2
    log_k     -6.63
    delta_h   -1.037 kcal
    -analytical_expression -14805.9622 -2.4660924 756968.533 5436.3588 -40553604 0
redoxBarite
    BaSulfateO4 = Ba+2 + SulfateO4-2
    log_k     -9.97
    delta_h   6.35 kcal
    -analytical_expression 136.035 0 -7680.41 -48.595 0 0
Hydroxyapatite
    Ca5(PO4)3OH + 4H+ = 5Ca+2 + H2O + 3HPO4-2
    log_k     -3.421
    delta_h   -36.155 kcal
Fluorite
    CaF2 = Ca+2 + 2F-
    log_k     -10.6
    delta_h   4.69 kcal
    -analytical_expression 66.348 0 -4298.2 -25.271 0 0
SiO2(a)
    SiO2 + 2H2O = H4SiO4
    log_k     -2.71
    delta_h   3.34 kcal
    -analytical_expression -0.26 0 -731 0 0 0
Chalcedony
    SiO2 + 2H2O = H4SiO4
    log_k     -3.55
    delta_h   4.72 kcal
    -analytical_expression -0.09 0 -1032 0 0 0
Quartz
    SiO2 + 2H2O = H4SiO4
    log_k     -3.98
    delta_h   5.99 kcal
    -analytical_expression 0.41 0 -1309 0 0 0
Gibbsite
    Al(OH)3 + 3H+ = Al+3 + 3H2O
    log_k     8.11
    delta_h   -22.8 kcal
Al(OH)3(a)
    Al(OH)3 + 3H+ = Al+3 + 3H2O
    log_k     10.8
    delta_h   -26.5 kcal
Kaolinite
    Al2Si2O5(OH)4 + 6H+ = 2Al+3 + H2O + 2H4SiO4
    log_k     7.435
    delta_h   -35.3 kcal
Albite
    NaAlSi3O8 + 8H2O = Al(OH)4- + 3H4SiO4 + Na+
    log_k     -18.002
    delta_h   25.896 kcal
Anorthite
    CaAl2Si2O8 + 8H2O = 2Al(OH)4- + Ca+2 + 2H4SiO4
    log_k     -19.714
    delta_h   11.58 kcal
K-feldspar
    KAlSi3O8 + 8H2O = Al(OH)4- + 3H4SiO4 + K+
    log_k     -20.573
    delta_h   30.82 kcal
K-mica
    KAl3Si3O10(OH)2 + 10H+ = 3Al+3 + 3H4SiO4 + K+
    log_k     12.703
    delta_h   -59.376 kcal
Chlorite(14A)
    Mg5Al2Si3O10(OH)8 + 16H+ = 2Al+3 + 6H2O + 3H4SiO4 + 5Mg+2
    log_k     68.38
    delta_h   -151.494 kcal
Ca-Montmorillonite
    Ca0.165Al2.33Si3.67O10(OH)2 + 12H2O = 2.33Al(OH)4- + 0.165Ca+2 + 2H+ + 3.67H4SiO4
    log_k     -45.027
    delta_h   58.373 kcal
Talc
    Mg3Si4O10(OH)2 + 6H+ + 4H2O = 4H4SiO4 + 3Mg+2
    log_k     21.399
    delta_h   -46.352 kcal
Illite
    K0.6Mg0.25Al2.3Si3.5O10(OH)2 + 11.2H2O = 2.3Al(OH)4- + 1.2H+ + 3.5H4SiO4 + 0.6K+ + 0.25Mg+2
    log_k     -40.267
    delta_h   54.684 kcal
Chrysotile
    Mg3Si2O5(OH)4 + 6H+ = H2O + 2H4SiO4 + 3Mg+2
    log_k     32.2
    delta_h   -46.8 kcal
    -analytical_expression 13.248 0 10217.1 -6.1894 0 0
Sepiolite
    Mg2Si3O7.5OH:3H2O + 4H+ + 0.5H2O = 3H4SiO4 + 2Mg+2
    log_k     15.76
    delta_h   -10.7 kcal
Sepiolite(d)
    Mg2Si3O7.5OH:3H2O + 4H+ + 0.5H2O = 3H4SiO4 + 2Mg+2
    log_k     18.66
redoxHematite
    Ferric2O3 + 6H+ = 2Ferric+3 + 3H2O
    log_k     -4.008
    delta_h   -30.845 kcal
redoxGoethite
    FerricOOH + 3H+ = Ferric+3 + 2H2O
    log_k     -1
    delta_h   -14.48 kcal
redoxFe(OH)3(a)
    Ferric(OH)3 + 3H+ = Ferric+3 + 3H2O
    log_k     4.891
redoxPyrite
    FerrousSulfide2 + 2H+ + 2e- = Ferrous+2 + 2HSulfide-
    log_k     -18.479
    delta_h   11.3 kcal
redoxFeS(ppt)
    FerrousSulfide + H+ = Ferrous+2 + HSulfide-
    log_k     -3.915
redoxMackinawite
    FerrousSulfide + H+ = Ferrous+2 + HSulfide-
    log_k     -4.648
redoxSulfur
    Sulfide + 2H+ + 2e- = H2Sulfide
    log_k     4.882
    delta_h   -9.5 kcal
redoxVivianite
    Ferrous3(PO4)2:8H2O = 3Ferrous+2 + 8H2O + 2PO4-3
    log_k     -36
redoxHausmannite
    ManganousManganic2O4 + 8H+ = 4H2O + 2Manganic+3 + Manganous+2
    log_k     10.01
redoxManganite
    ManganicOOH + 3H+ = 2H2O + Manganic+3
    log_k     -0.17
redoxPyrochroite
    Manganous(OH)2 + 2H+ = 2H2O + Manganous+2
    log_k     15.2
Halite
    NaCl = Cl- + Na+
    log_k     1.582
    delta_h   0.918 kcal
redoxCO2(g)
    CO2 = CO2
    log_k     -1.468
    delta_h   -4.776 kcal
    -analytical_expression 108.3865 0.01985076 -6919.53 -40.45154 669365 0
redoxO2(g)
    Ozero2 = Ozero2
    log_k     -2.8983
    -analytical_expression -7.5001 0.0078981 0 0 200270 0
H2(g)
    H2 = H2
    log_k     -3.15
    delta_h   -1.759 kcal
redoxH2(g)
    Hzero2 = Hzero2
    log_k     -3.15
    delta_h   -1.759 kcal
H2O(g)
    H2O = H2O
    log_k     1.51
    delta_h   -44.03 kJ
redoxN2(g)
    Nzero2 = Nzero2
    log_k     -3.26
    delta_h   -1.358 kcal
redoxH2S(g)
    H2Sulfide = H2Sulfide
    log_k     -0.997
    delta_h   -4.57 kcal
redoxCH4(g)
    MethaneH4 = MethaneH4
    log_k     -2.86
    delta_h   -3.373 kcal
redoxNH3(g)
    AmmH3 = AmmH3
    log_k     1.77
    delta_h   -8.17 kcal
redoxMelanterite
    FerrousSulfateO4:7H2O = Ferrous+2 + 7H2O + SulfateO4-2
    log_k     -2.209
    delta_h   4.91 kcal
    -analytical_expression 1.447 -0.004153 0 0 -214949 0
redoxAlunite
    KAl3(SulfateO4)2(OH)6 + 6H+ = 3Al+3 + 6H2O + K+ + 2SulfateO4-2
    log_k     -1.4
    delta_h   -50.25 kcal
redoxJarosite-K
    KFerric3(SulfateO4)2(OH)6 + 6H+ = 3Ferric+3 + 6H2O + K+ + 2SulfateO4-2
    log_k     -9.21
    delta_h   -31.28 kcal
Zn(OH)2(e)
    Zn(OH)2 + 2H+ = 2H2O + Zn+2
    log_k     11.5
redoxSmithsonite
    ZnCO3 = Zn+2 + CO3-2
    log_k     -10
    delta_h   -4.36 kcal
redoxSphalerite
    ZnSulfide + H+ = HSulfide- + Zn+2
    log_k     -11.618
    delta_h   8.25 kcal
Willemite
    Zn2SiO4 + 4H+ = H4SiO4 + 2Zn+2
    log_k     15.33
    delta_h   -33.37 kcal
Cd(OH)2
    Cd(OH)2 + 2H+ = Cd+2 + 2H2O
    log_k     13.65
redoxOtavite
    CdCO3 = Cd+2 + CO3-2
    log_k     -12.1
    delta_h   -0.019 kcal
CdSiO3
    CdSiO3 + 2H+ + H2O = Cd+2 + H4SiO4
    log_k     9.06
    delta_h   -16.63 kcal
redoxCdSO4
    CdSulfateO4 = Cd+2 + SulfateO4-2
    log_k     -0.1
    delta_h   -14.74 kcal
redoxCerrusite
    PbCO3 = Pb+2 + CO3-2
    log_k     -13.13
    delta_h   4.86 kcal
redoxAnglesite
    PbSulfateO4 = Pb+2 + SulfateO4-2
    log_k     -7.79
    delta_h   2.15 kcal
Pb(OH)2
    Pb(OH)2 + 2H+ = 2H2O + Pb+2
    log_k     8.15
    delta_h   -13.99 kcal

EXCHANGE_MASTER_SPECIES
    X             X-           
EXCHANGE_SPECIES
X- = X-
    log_k     0
Na+ + X- = NaX
    log_k     0
    -gamma    4 0.075
K+ + X- = KX
    log_k     0.7
    delta_h   -4.3 kJ
    -gamma    3.5 0.015
Li+ + X- = LiX
    log_k     -0.08
    delta_h   1.4 kJ
    -gamma    6 0
AmmH4+ + X- = AmmH4X
    log_k     0.6
    delta_h   -2.4 kJ
    -gamma    2.5 0
Ca+2 + 2X- = CaX2
    log_k     0.8
    delta_h   7.2 kJ
    -gamma    5 0.165
Mg+2 + 2X- = MgX2
    log_k     0.6
    delta_h   7.4 kJ
    -gamma    5.5 0.2
Sr+2 + 2X- = SrX2
    log_k     0.91
    delta_h   5.5 kJ
    -gamma    5.26 0.121
Ba+2 + 2X- = BaX2
    log_k     0.91
    delta_h   4.5 kJ
    -gamma    5 0
Manganous+2 + 2X- = ManganousX2
    log_k     0.52
    -gamma    6 0
Ferrous+2 + 2X- = FerrousX2
    log_k     0.44
    -gamma    6 0
Cupric+2 + 2X- = CupricX2
    log_k     0.6
    -gamma    6 0
2X- + Zn+2 = ZnX2
    log_k     0.8
    -gamma    5 0
Cd+2 + 2X- = CdX2
    log_k     0.8
    -davies 
Pb+2 + 2X- = PbX2
    log_k     1.05
    -davies 
Al+3 + 3X- = AlX3
    log_k     0.41
    -gamma    9 0
AlOH+2 + 2X- = AlOHX2
    log_k     0.89
    -davies 
SURFACE_MASTER_SPECIES
    Hfo_s         Hfo_sOH     
    Hfo_w         Hfo_wOH     
SURFACE_SPECIES
Hfo_sOH = Hfo_sOH
    log_k     0
H+ + Hfo_sOH = Hfo_sOH2+
    log_k     7.29
Hfo_sOH = Hfo_sO- + H+
    log_k     -8.93
Hfo_wOH = Hfo_wOH
    log_k     0
H+ + Hfo_wOH = Hfo_wOH2+
    log_k     7.29
Hfo_wOH = Hfo_wO- + H+
    log_k     -8.93
Ca+2 + Hfo_sOH = Hfo_sOHCa+2
    log_k     4.97
Ca+2 + Hfo_wOH = Hfo_wOCa+ + H+
    log_k     -5.85
Hfo_sOH + Sr+2 = Hfo_sOHSr+2
    log_k     5.01
Hfo_wOH + Sr+2 = Hfo_wOSr+ + H+
    log_k     -6.58
H2O + Hfo_wOH + Sr+2 = Hfo_wOSrOH + 2H+
    log_k     -17.6
Ba+2 + Hfo_sOH = Hfo_sOHBa+2
    log_k     5.46
Ba+2 + Hfo_wOH = Hfo_wOBa+ + H+
    log_k     -7.2
Cd+2 + Hfo_sOH = Hfo_sOCd+ + H+
    log_k     0.47
Cd+2 + Hfo_wOH = Hfo_wOCd+ + H+
    log_k     -2.91
Hfo_sOH + Zn+2 = Hfo_sOZn+ + H+
    log_k     0.99
Hfo_wOH + Zn+2 = Hfo_wOZn+ + H+
    log_k     -1.99
Cupric+2 + Hfo_sOH = Hfo_sOCupric+ + H+
    log_k     2.89
Cupric+2 + Hfo_wOH = Hfo_wOCupric+ + H+
    log_k     0.6
Hfo_sOH + Pb+2 = Hfo_sOPb+ + H+
    log_k     4.65
Hfo_wOH + Pb+2 = Hfo_wOPb+ + H+
    log_k     0.3
Hfo_wOH + Mg+2 = Hfo_wOMg+ + H+
    log_k     -4.6
Hfo_sOH + Manganous+2 = Hfo_sOManganous+ + H+
    log_k     -0.4
Hfo_wOH + Manganous+2 = Hfo_wOManganous+ + H+
    log_k     -3.5
Ferrous+2 + Hfo_sOH = Hfo_sOFerrous+ + H+
    log_k     -0.95
Ferrous+2 + Hfo_wOH = Hfo_wOFerrous+ + H+
    log_k     -2.98
Ferrous+2 + H2O + Hfo_wOH = Hfo_wOFerrousOH + 2H+
    log_k     -11.55
3H+ + Hfo_wOH + PO4-3 = Hfo_wH2PO4 + H2O
    log_k     31.29
2H+ + Hfo_wOH + PO4-3 = Hfo_wHPO4- + H2O
    log_k     25.39
H+ + Hfo_wOH + PO4-3 = Hfo_wPO4-2 + H2O
    log_k     17.72
H3BO3 + Hfo_wOH = Hfo_wH2BO3 + H2O
    log_k     0.62
H+ + Hfo_wOH + SulfateO4-2 = Hfo_wSulfateO4- + H2O
    log_k     7.78
Hfo_wOH + SulfateO4-2 = Hfo_wOHSulfateO4-2
    log_k     0.79
F- + H+ + Hfo_wOH = Hfo_wF + H2O
    log_k     8.7
F- + Hfo_wOH = Hfo_wOHF-
    log_k     1.6
END

You will have to use the appropriate "element" redox states with the approptriate stoichiometry. For example, to define a quarter millimole of dissolved oxygen, you would define

Code: [Select]

SOLUTION
-units mmol/kgw
Ozero 0.5

The nitrification -formula in KINETICS would be

Code: [Select]

-formula Amm -1 Nitrate +1 Ozero -2 O +2

Note that you only need to list the conversions from one redox state to another (with this database), but you need to use the correct stoichiometry. Assuming Doc is CH2O, the denitrification formula would be

Code: [Select]

-formula Nitrate -2 Nzero 2 Doc -2.5 CH2O +2.5

To avoid confusion, you should disallow the formation of H2(aq) and O2(aq) by decreasing the log Ks with the following reactions:

Code: [Select]

SOLUTION_SPECIES
2H2O = O2 + 4H+ + 4e-
   # log_k     -86.08
   log_k -150

2H+ + 2e- = H2
    #log_k     -3.15
    log_k -50
END

Note that it is up to you  to enforce thermodynamic reality. Depending on your reaction network, there are ways to use thermodynamics to implement some of the reactions. For a simple example, if you have Fe(2) in solution, you can oxidize it to Fe(3) simply by adding O2 kinetically.

[Rensburg]

  Sorry, I think you misunderstood me. I'm not trying to simulate purely kinetic reactions. I want to distinguish between O(0) and O(-2) because I know that O(0) can oxidize substances like Fe(2+) and As(3+), while O(-2) cannot. My goal is not to simulate nitrification and denitrification per se, but rather to simulate the change in dissolved oxygen under the influence of nitrification and denitrification, in order to study the oxidation capacity of the water. Would the database you mentioned lead to something less intelligent, requiring me to manually define many thermodynamic redox reactions? If you still think that using the database you mentioned is the most suitable approach for me, then what exactly do you mean by "you need to ensure thermodynamic reality yourself"?

  My concern is that when I set up kinetic reactions—for example, for the reaction you mentioned: NH₄⁺ + 2O₂ = NO₃⁻ + H₂O + 2H⁺—if I define it in terms of element gain/loss, using something like "-formula Amm -1 N 1 H 3", such a definition seems unable to reflect the reduction of O₂ in the chemical equation, and therefore cannot represent the weakening of the solution's ability to oxidize Fe(2+) and As(3+). I'm not sure how to handle this. By the way, I'm also wondering whether I should include "H 3" in the definition. Could you give me some advice?

  Similarly, during denitrification, the chemical equation is: NO₃⁻ + 1.25 CH₂O = 0.5 N₂ + 1.25 CO₂ + 0.75 H₂O + OH⁻. Should I define CH₂O as a new element, say "Doc"? In that case, what should the alkalinity of Doc be set to? Also, in this expression, should it be represented as "-formula N -1 Ntg 0.5"?

[dlparkhurst]

>Would the database you mentioned lead to something less intelligent, requiring me to manually define many thermodynamic redox reactions?

Yes, I think so. Thermodynamics would need to be incorporated in the rate equations. Here is a slightly strained example, but, at equilibrium the log partial pressure of O2(g) is about -40 when you have equal amounts of N(5) and N(-3), and N2(aq) is ignored.

Code: [Select]

SOLUTION_SPECIES
2 NO3- + 12 H+ + 10 e- = N2 + 6 H2O
# -log_k 207.08
      -log_k 100
SOLUTION
pe 4 N(5)/N(3)
N(5) 1
N(-3) 1
END

If you start with Amm and have rates to convert to Nitrate and assuming a fixed log P(O2) at -40, then the rate equation needs to know to stop at equilibrium when half the Amm has been converted to Nitrate, which is when the log K for the reaction AmmH4+ + 2O2 = NitrateO3- + 2H+ +  H2O is satisfied.

Amm is the nitrogen part of NH3 or NH4+ in the database that I posted, so AmmH4+ and AmmH3. (However, to add confusion, in the Amm.dat database, Amm represents NH3, so Amm and AmmH+.)

In example 9, Fe_di is converted to Fe_tri and oxygen is consumed automatically; no O is included in the KINETICS formula. Same would happen by converting Amm to Nitrate with the database I posted.

In the redox database I posted Nzero is atoms of N, so the aqueous species is Nzero2. In phreeqc.dat and Amm.dat, Ntg is equivalent to two nitrogens, N2, and you would need two NO3- to make one Ntg. I think the redox database is the better way to go, but that is in retrospect.

Nitrogen is a tough system because it does not react to thermodynamic equilibrium. Nitrification skips N2, and denitrification stops at N2 rather than N(-3). If you can remove N2 (as in the example above), you can simulate rates of oxidation simply by adding O2 to the system kinetically. If you do this, then sulfide, ammonium, iron(2), arsenite, etc will all oxidize in the proper sequence. If you treat these redox active species as separate elements, as in the posted redox.dat, then I think you must add many inhibitor terms and thermodynamic equilibrium constants to the rate equations to simulate plausible thermodynamics. I don't think you want to do this.

If you want to simulate denitrification only, you could use Amm.dat and simply add CH2O kinetically. However, in general, I think you have to develop a kinetic reaction that has a formula that removes N, creates Ntg, and adds CH2O to solution in the proper stoichiometry (using phreeqc.dat or Amm.dat and eliminating N2(aq)).

[Rensburg]

Thank you for your advice. Let me see if I understand you correctly. You mean that when using the database you provided, even if the KINETICS formula does not include O, converting Amm to Nitrate still consumes oxygen. However, to achieve the goal "simulate the change in dissolved oxygen under the influence of nitrification and denitrification, in order to study the oxidation capacity of the water.", using Amm.dat is still the best choice—provided that I define the kinetic reaction with the correct stoichiometry to remove N, produce Ntg, and add CH2O to the solution. Is that correct?

Also, what do you mean by "If you can remove N2 (as in the example above), you can simulate rates of oxidation simply by adding O2 to the system kinetically. If you do this, sulfide, ammonium, iron(2), arsenite, etc will all oxidize in the proper sequence"? Under which database should I perform this operation?

[dlparkhurst]

Here is a script that shows the expected oxidation sequence starting from reducing conditions to equilibrium with atmospheric oxygen.

It also shows the reduction sequence starting from atmospheric oxygen to fully reduced conditions.

It uses phreeqc.dat. In the oxidizing part, N2 is disabled. In the reducing sequence N(-3) is disabled. Alternatively, Amm.dat could be used for the reducing sequence, although not in the same run because only one database is allowed for a phreeqc run.

Now, you will probably say, "But I want to be able to do either oxidizing or reducing without modifying the database". And I will say that this is simply a demonstration of the expected oxidation and reduction sequences. To model denitrification, you can use the definitions in the oxidation sequence, but you will need to have a kinetic reaction that under the appropriate conditions converts NO3- to Ntg (possibly including the conversion of Doc to CH2O) so that all of the N is consumed before conditions are sufficiently reducing that N(-3) becomes the dominant nitrogen species.

Code: [Select]

SOLUTION_SPECIES
2 NO3- + 12 H+ + 10 e- = N2 + 6 H2O
# -log_k 207.08
      -log_k 100
SOLUTION
pH 7
N(-3) 1
Fe(2) 1
S(-2) 1
Cl 1 charge
END
USE solution 1
REACTION
O2 1
0.005 in 20
EQUILIBRIUM_PHASES
O2(g) -0.7 0
USER_GRAPH 1
    -headings               rxn Fe(2) Fe(3) S(-2) S(6) N(-3) N(0) N(5) O2
    -axis_titles            "O2 reacted" "Mol/kgw" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X RXN
20 GRAPH_Y TOT("Fe(2)"), TOT("Fe(3)"), TOT("S(-2)"), TOT("S(6)")
30 GRAPH_Y TOT("N(-3)"), TOT("N(0)"), TOT("N(5)"), MOL("O2")
  -end
    -active                 true
SAVE solution 2
END
SOLUTION_SPECIES
2 NO3- + 12 H+ + 10 e- = N2 + 6 H2O
-log_k 207.08
-delta_h -312.13 kcal
-Vm 7 # Pray et al., 1952, IEC 44 1146
-dw 1.96e-9 -90 # Cadogan et al. 2014, JCED 59, 519
NO3- + 10 H+ + 8 e- = NH4+ + 3 H2O
# -log_k 119.077
      -log_k 0
END
USER_GRAPH 1
-detach
END
USE solution 2
DELETE
    -equilibrium_phases   1
REACTION
CH2O 1
0.005 in 20
USER_GRAPH 2
    -headings               rxn Fe(2) Fe(3) S(-2) S(6) N(-3) N(0) N(5) O2
    -axis_titles            "O2 reacted" "Mol/kgw" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X RXN
20 GRAPH_Y TOT("Fe(2)"), TOT("Fe(3)"), TOT("S(-2)"), TOT("S(6)")
30 GRAPH_Y TOT("N(-3)"), TOT("N(0)"), TOT("N(5)"), MOL("O2")
  -end
END


[Rensburg]

You did anticipate my question. I'm not sure if I understand you correctly. Are you saying that I can define nitrification like the oxidation sequence in the phreeqc.dat script you provided, and then use kinetics to simulate denitrification to convert NO₃⁻ to Ntg and Doc to CH₂O? Is that correct? If I use the Amm.dat database, just like modifying SOLUTION_SPECIES in the reduction part of your script, then define kinetics to convert Amm to N? Is that correct?

If I use Amm.dat and simulate both nitrification and denitrification with kinetics, can I do that? Or do you think that only one of nitrification and denitrification can be set up as a kinetic reaction?When you said "under the appropriate conditions convert NO₃⁻ to Ntg (possibly also converting Doc to CH₂O)", do you mean something like Denitrification?

Code: [Select]

SOLUTION_MASTER_SPECIES
Doc Doc 0 CH2O 30
SOLUTION_SPECIES
#NO3- + 10 H+ + 8 e- = NH4+ + 3 H2O (Because Amm.dat does not have this definition, I commented it out)
# -log_k 119.077
#      -log_k 0

2 NO3- + 12 H+ + 10 e- = N2 + 6 H2O
# -log_k 207.08
      -log_k 100

Doc = Doc
      -log_k 0
END
SOLUTION
pH 7
Amm 1
N(-5) 1
Doc 1
Fe(2) 1
S(-2) 1
Cl 1 charge
END
USE solution 1
KINETICS
Nitrification
-formula Amm -1 N 1 H 3
Denitrification
-formula Doc -2.5 N -2 C 2.5 H 5 O 2.5 Ntg 1Would this lead to any unreasonable results?

[dlparkhurst]

Too many questions. Try it.

[Rensburg]

I'm really sorry, I didn't understand what you meant. I tried it this way, and the program runs, but I just don't know whether the results are reasonable.
I've improved my code today after trying hard to understand your ideas. Would you be able to take another look and see if this is reasonable?I have two tasks. One is to correctly simulate the change in dissolved oxygen through the kinetics of nitrification and denitrification. The second is to confirm whether the nitrification kinetics should include "-formula H 3" and, when defining a new species such as Doc, what the alkalinity should be based on, to ensure that the pH value is reasonable during the simulation.
If my current test code can correctly achieve these two goals, I think I will have completed my task. If there is anything unexpected that requires improvement in my code, could you please guide me further? Or is my code completely wrong?

Code: [Select]

SOLUTION_MASTER_SPECIES
    Doc       Doc          0     CH2O         30

SOLUTION_SPECIES
Doc = Doc
    log_k     0
    delta_h   0 kJ
2 NO3- + 12 H+ + 10 e- = N2 + 6 H2O
#-log_k 207.08
       -log_k 100
-delta_h -312.130 kcal
-dw 1.96e-9
-Vm 7 # Pray et al., 1952, IEC 44. 1146
RATES
        Nitrification
        -start
10 k12c = 0.2/(3600*24)
20 k12t = 1.05^(TC-20)
30 Knit = 2
40 f1 = (MOL("O2")*32*1000)/(Knit + MOL("O2")*32*1000)
50 rate = k12c * k12t * f1 * TOT("Amm") * 14*1000
55 rate = rate/14
60 moles = rate * TIME
70 if (M + moles) < 0 then moles = -M
80 SAVE moles
        -end

        Denitrification
        -start
10 k20c = 0.2/(3600*24)
20 k20t = 1.05^(TC-20)
30 Kno3 = 2
40 f1 = Kno3/(Kno3 + MOL("O2")*32*1000)
50 if MOL("Doc") <= 0 then rate = 0 else rate = k20c * k20t * f1 * TOT("N") * 14*1000
60 rate = rate/14
90 moles = rate * TIME
100 if (M + moles) < 0 then moles = -M
110 SAVE moles
        -end

        CBOD_increase
        -start
10 aoc = 2.67
20 K1d = 0.05/(3600*24)
30 plant = 0.005  #unit/mg
40 rate = aoc * K1d * plant
55 rate =  rate/32
50 moles = rate * TIME
60 if (M + moles) < 0 then moles = -M
70 SAVE moles
        -end

        Mineralization_N
10 k71c = 0.1/(3600*24)
20 k71t = 1.05^(TC-20)
30 Kbod = 0.5 # 0~1
40 rate =  k71c * k71t * Kbod * M *1000
45 rate =  rate/14
50 moles = rate * TIME
100 if (M + moles) < 0 then moles = -M
110 SAVE moles
        -end

SOLUTION 1
    temp      3.1
    pH        7.81
    pe        4.1
    redox     pe
    units     mg/l
    density   1
    O(0)      4.47
    Amm       0
    N(5)      0.39
    N(3)      0
    P         0.72
    Doc       0
    Si        6.71
    Fe(2)     0.01
    S(-2)     2#0#3
    Ca        22.165
    Mg        2.739
    Na        0
    K         0
    Cl        1.4481
    S(6)      4.6373
    F         0.09
    -water    1 # kg

SELECTED_OUTPUT
        -file   20260402_2.xlsx
        -high_precision true
USER_PUNCH
        -headings       Nitrification_t Denitrification_t DO_increase_t CBOD_decrease_t CBOD_increase_t Mineralization_N_t Amm N Doc O2 Ntg
        -start
  10 punch KIN_DELTA("Nitrification")  KIN_DELTA("Denitrification") KIN_DELTA("DO_increase") KIN_DELTA("CBOD_decrease") KIN_DELTA("CBOD_increase") KIN_DELTA("Mineralization_N")
  20 punch TOT("Amm")  TOT("N")  MOL("Doc")  MOL("O2") MOL("Ntg")
        -end
USER_GRAPH 1
    -headings               rxn Fe(2) Fe(3) S(-2) S(6) Doc N(-3) N(0) N(5) O2
    -axis_titles            "time"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME
20 GRAPH_Y TOT("Fe(2)"), TOT("Fe(3)"), TOT("S(-2)"), TOT("S(6)"), TOT("Doc")
30 GRAPH_Y TOT("Amm"), TOT("Ntg"), TOT("N(5)"), MOL("O2")/10
  -end
    -active                 true
INCREMENTAL_REACTIONS true

KINETICS 1
Nitrification
    -formula  Amm -1 N 1 H 3
    -tol      1e-08
Denitrification
    -formula  N -1 Doc -1.25 Ntg  0.25 CH2O 1.25
    -tol      1e-08
CBOD_increase
    -formula  Doc  1
    -tol      1e-08
Mineralization_N
    -formula  Amm 1
    -m        0.0004
    -tol      1e-08
-steps       3600 in 60 steps
-step_divide 1
-runge_kutta 3
-bad_step_max 500
END

[dlparkhurst]

You are not helpless. Do the different constituents change in the way you expect? Do you expect production of Amm, N(5), or Ntg?

I will respond to a specific question, preferably on a single part of the calculation; it would make sense to start by considering whether you understand the initial consumption of O2(aq), then proceed to the next expected result. I will not answer whether your script is right or not. You must be the judge of that.

Note that this formula looks wrong, N is not conserved (Ntg*0.25 is equivalent to 0.5 N, whereas the reactant is 1 N):

Code: [Select]

N -1 Doc -1.25 Ntg  0.25 CH2O 1.25


[Rensburg]

I understand. Thank you for reading my code. If the general direction of this setup is correct, I will verify it step by step. If anything unexpected happens, I will consult you again. Thank you very much.