Expression of simplified organic degradation reaction in PHREEQC

[Daehyun]

Dear PHREEQC Community,

I am an undergraduate student who is not yet very experienced with PHREEQC.
Currently, I am working on developing a geochemical reaction model, and I am still in the early stages of the project.

What I am trying to simulate is the generation of reactive gases through the degradation of organic matter.
At this point, I am focusing on a simple reaction in which formaldehyde (CH₂O) is degraded by oxygen to produce carbon dioxide, as shown below:

  CH₂O + O₂ → CO₂ + H₂O

I would like to know whether the code I have written (included below) reasonably represents this reaction.
Since I was unable to define CH₂O using SOLUTION_MASTER_SPECIES and could not implement the reaction in SOLUTION_SPECIES,
I defined the process using a KINETICS block with the -formula keyword.

Here are the questions I would appreciate help with:

Is it appropriate to define organic matter only through the KINETICS block in this way?

If I want to add more types of organic matter degradation (e.g., with nitrate, sulfate, or methanogenesis),
  what would be the most rational or recommended approach within PHREEQC?

Thank you very much for your help.

Code: [Select]

SOLUTION_MASTER_SPECIES
    Form          Form             0     Form            30   #Formaldehyde
    Ntg           Ntg              0     Ntg             28.013   #N2 gas

SOLUTION_SPECIES
Form = Form   #Formaldehyde
    log_k     0
Ntg = Ntg # N2 gas 
-dw 1.96e-9
-Vm 7 # Pray et al., 1952, IEC 44. 1146

PHASES
Ntg(g)
    Ntg = Ntg
    log_k     -3.1864
    delta_h   -10.4391 kJ
    -analytical_expression -58.453 0.001818 3199 17.909 -27460 0
    -T_c      126.2
    -P_c      33.5
    -Omega    0.039
CO2(g)
    CO2 = CO2
    log_k     -1.468
    delta_h   -4.776 kcal
    -analytical_expression 10.5624 -0.023547 -3972.8 0 587460 1.9194e-05
    -T_c      304.2
    -P_c      72.86
    -Omega    0.225
O2(g)
    O2 = O2
    log_k     -2.8983
    -analytical_expression -7.5001 0.0078981 0 0 200270 0
    -T_c      154.6
    -P_c      49.8
    -Omega    0.021


SOLUTION 1
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    Form      1
    -water    1 # kg
GAS_PHASE 1
    -fixed_pressure
    -pressure 1
    -volume 1
    -temperature 25
    O2(g)     0.7
    CO2(g)    0.3


RATES
Organic_degradation
 -start
10  if (M <= 0) THEN GOTO 200
20  mO2   = MOL("O2")
50  k_O2  = 1.04167E-9    # 1/sec   
80  rate  = k_O2 * mO2/(3e-04 + mO2)
110 moles = rate * M * (M/M0) * TIME
200 SAVE moles
 -end

KINETICS 1

Organic_degradation
    -formula  Form  -1 CH2O  1
    -m        1
    -m0       1
    -tol      1e-008
-steps       100000 in 10 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500

END


[dlparkhurst]

I haven't looked in detail, but kinetic conversion of Form to CH2O is reasonable.

You can simulate the reaction more simply with REACTION, which avoids having to define reaction rates.

If you do want to simulate reaction rates, you can use the same formula but adjust the rate depending on the electron acceptors available. Often the factor you have used is  included to establish a maximum rate for an electron acceptor. Usually, you want to inhibit other electron acceptors until a more favorable electron acceptor is depleted. You can use factors like [MOL("O2") + 1e-5]/1e-5 to inhibit, say S(6) reduction, while dissolved oxygen is present.

You may want a separate kinetic reaction that converts N(5) to Ntg. I think the formula would be

Code: [Select]

-formula Form -1.25 CH2O 1.25 N -1 Ntg 0.5


[Daehyun]

Thank you very much for your helpful response.

I am currently trying to model the degradation of organic matter using kinetic reactions. However, when I define the reaction in the KINETICS block like:
-formula Form -1.25 CH2O 1

it seems to simply add C, H, and O to the aqueous phase as individual elements, rather than treating CH₂O as an actual molecule.

What I want to simulate is the degradation of organic compounds through reactions such as:

CH₂O + O₂ → CO₂ + H₂O
2CH₂O → CH₄ + CO₂

and other similar organic degradation pathways.

To represent this, I believe that CH₂O must be properly defined in the database, but whenever I try to define CH₂O in the SOLUTION_MASTER_SPECIES and SOLUTION_SPECIES blocks, I encounter errors during execution.

As a workaround, I defined CH₂O as a solid phase in the PHASES block, specifying the degradation reaction there, and then used it in the KINETICS block to simulate the reaction.

I would like to ask whether this approach is considered appropriate for representing the degradation of CH₂O, or if there is a better or more recommended way to define and simulate these types of organic matter degradation reactions in PHREEQC.

Thank you again for your time and support.

Code: [Select]

SOLUTION_MASTER_SPECIES
    Ntg           Ntg              0     Ntg             28.013   #N2 gas

SOLUTION_SPECIES
Ntg = Ntg # N2 gas 
-dw 1.96e-9
-Vm 7 # Pray et al., 1952, IEC 44. 1146

PHASES
Ntg(g)
    Ntg = Ntg
    log_k     -3.1864
    delta_h   -10.4391 kJ
    -analytical_expression -58.453 0.001818 3199 17.909 -27460 0
    -T_c      126.2
    -P_c      33.5
    -Omega    0.039
CO2(g)
    CO2 = CO2
    log_k     -1.468
    delta_h   -4.776 kcal
    -analytical_expression 10.5624 -0.023547 -3972.8 0 587460 1.9194e-05
    -T_c      304.2
    -P_c      72.86
    -Omega    0.225
O2(g)
    O2 = O2
    log_k     -2.8983
    -analytical_expression -7.5001 0.0078981 0 0 200270 0
    -T_c      154.6
    -P_c      49.8
    -Omega    0.021
Formalde
    CH2O + O2 = CO2 + H2O
    log_k     -100


SOLUTION 1
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    -water    1 # kg
GAS_PHASE 1
    -fixed_pressure
    -pressure 1
    -volume 1
    -temperature 25
    O2(g)     1
   


RATES
Formalde
 -start
10  if (M <= 0) THEN GOTO 200
20  mO2   = MOL("O2")
50  k_O2  = 1.04167E-9    # 1/sec   
80  rate  = k_O2 * mO2/(3e-04 + mO2)
110 moles = rate * M * (M/M0) * TIME
200 SAVE moles
 -end

KINETICS 1

Formalde
    -formula  Formalde
    -m        1
    -m0       1
    -tol      1e-008
-steps       100000 in 10 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 500

END


[dlparkhurst]

It is correct that C, H 2, and O are added. Redox equilibrium will determine which electron acceptor is used.

[dlparkhurst]

Consider this script.

Code: [Select]


SOLUTION 1
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    -water    1 # kg
O(0) 1  O2(g) -0.7
      S(6) 1
      Na   1    charge

REACTION
CH2O 1
0.004 mole in 100 steps

USER_GRAPH 1
    -headings               rxn O2 SO4 C(4) C(-4)
    -axis_titles            "CH2O added" "log molality" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X RXN
20 GRAPH_Y LOG10(MOL("O2")), LOG10(TOT("S(6)")), LOG10(TOT("C(4)")), LOG10(TOT("C(-4)"))
  -end
    -active                 true

END


[dlparkhurst]

Alternatively, you can define each redox state as a separate "element" and explicitly transfer between each redox state. I think this is a much more difficult undertaking, both in the complexity of the input file, and in trying to create rates for major and minor elements that are thermodynamically consistent.  See the following discussion.

https://phreeqcusers.org/index.php/topic,2781.msg10432.html#msg10432

[Daehyun]

Thank you for the previous discussion and the reference link?it was very helpful.

After testing the inhibitory term for electron acceptor competition in several PHREEQC simulations, I noticed that the two kinetic rate formulations below yield identical results:
in both cases, O₂ is consumed first, and only after it is completely depleted, SO₄?⁻ begins to be consumed.
This behavior persists even when I significantly increase the rate constant for sulfate reduction or adjust the inhibitor expression.

From this, I suspect that the model is only prioritizing based on the rate of CH₂O oxidation (DOC input), and not explicitly recognizing electron acceptor preference based on redox potential.

Code: [Select]

SOLUTION_MASTER_SPECIES
    Doc           Doc              0       Doc             12

SOLUTION_SPECIES
Doc = Doc
     log_k     0


SOLUTION 1
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Doc       1000
    O(0)      1
    S(6)      1
    -water    1 # kg




RATES
DOC_oxidation
-start
 10  if (M <= 0) THEN GOTO 200
 20  mO2   = MOL("O2")
 40  mSO4  = MOL("SO4-2")
 50  k_O2  = 1.57e-9
 70  k_SO4 = 1.57e-9
 80  rate  = k_O2 * mO2/(1e-10 + mO2)
100 rate  = rate + k_SO4 * mSO4/(1.e-10 + mSO4) * 1e-6/(1e-6 + mO2)
110 moles = rate * TIME
200 SAVE moles
-end

KINETICS 1
DOC_oxidation
    -formula  CH2O  1 Doc  -1
    -m        1
    -m0       1
    -tol      1e-008
-steps       10 100 1000 10000 100000 1000000  10000000  100000000  1000000000
-step_divide 1
-runge_kutta 3
-bad_step_max 500

INCREMENTAL_REACTIONS
USER_GRAPH 1
    -headings               time O2 SO4
    -axis_scale x_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TIME
20 GRAPH_Y MOL("O2"), MOL("SO4-2")
  -end
    -active                 true


Code: [Select]

SOLUTION_MASTER_SPECIES
    Doc           Doc              0       Doc             12

SOLUTION_SPECIES
Doc = Doc
     log_k     0


SOLUTION 1
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mmol/kgw
    density   1
    Doc       1000
    O(0)      1
    S(6)      1
    -water    1 # kg




RATES
DOC_oxidation
-start
 10  if (M <= 0) THEN GOTO 200
 20  mO2   = MOL("O2")
 40  mSO4  = MOL("SO4-2")
 50  k_O2  = 1.57e-9
 70  k_SO4 = 1.57e-9
 80  rate  = k_O2 * mO2/(1e-10 + mO2)
100 rate  = rate + k_SO4 * mSO4/(1.e-10 + mSO4)
110 moles = rate * TIME
200 SAVE moles
-end

KINETICS 1
DOC_oxidation
    -formula  CH2O  1 Doc  -1
    -m        1
    -m0       1
    -tol      1e-008
-steps       10 100 1000 10000 100000 1000000  10000000  100000000  1000000000
-step_divide 1
-runge_kutta 3
-bad_step_max 500

INCREMENTAL_REACTIONS
USER_GRAPH 1
    -headings               time O2 SO4
    -axis_scale x_axis      auto auto auto auto log
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TIME
20 GRAPH_Y MOL("O2"), MOL("SO4-2")
  -end
    -active                 true

My question is:
Does PHREEQC internally enforce any predefined priority of electron acceptors (such as O₂ > NO₃⁻ > SO₄?⁻), even without inhibition terms?
Or is this outcome solely driven by the kinetic expressions and reactant concentrations I define?

I?d appreciate any clarification on how PHREEQC handles competitive redox pathways when multiple acceptors are present.

[dlparkhurst]

You are reacting CH2O in all cases. You are only changing the rate at which CH2O is added to the solution, so the sequence of reactions will remain the same regardless of the rates or inhibition factors. The approach is sound, in that you could define a faster rate when O2 is present, and a slower rate when O2 is absent and NO3 is present, etc. 

As for your question, there is no explicit definition of the sequence of electron acceptors in the sense that you are asking the question. The sequence of electron acceptors is determined by the log Ks in SOLUTION_SPECIES for the reactions that involve e- (or O2(aq)). You could change the sequence of electron acceptors if you, for instance, increased the log K for the reaction of SO4-2 to S-2, but you would have difficulty justifying such a change.

I'm not sure why this would be, but if you are unhappy with the sequence of reactions (and don't want to change log Ks), you will need to define multiple "elements" to represent the relevant redox states and then write rate equations that explicitly transform one state to another. This would be difficult in that you would have to suppress and invoke reactions with a complicated scheme of inhibition factors and other algorithms to force reactions to occur in a sequence other than the thermodynamic sequence. The forum has a post with a database redox.dat that separates every redox state to a different "element". Kinetic reactions would be used to transfer one redox state to another. There should be further discussion if you search the forum for "redox.dat".