Conditional reaction product in Kinetics

[florisWQ]

Dear all, Hi David,

There's something with the way the KINETICS block works in PHREEQC that I hope you could help me understand. I would like to model the oxidation of DOC by O2 and NO3, while assigning the reacted N to the redox-uncoupled gas Ntg. I'd like to do this for bookkeeping the reactant now, but in the future I may want to add kinetic reactions that include uncoupled N2O and N2.

I have defined a rate for DOC oxidation similar to organic_C in the phreeqc.dat database, see the script snippet below. Both a reaction with O2 and NO3 can occur, but depending on solution composition, both reactions do not need to occur at the same time and I am thinking to inhibit the reaction with NO3 as long as O2 is present.

This leaves me wondering how I could then correctly give the formula of the reaction in the KINETICS block. If both O2 and NO3 reacts, I suppose the reaction mass balance can be given as:
- formula Doc -1 CH2O 1 Ntg +1? But how should I use this block if nitrate is not always reacting at each timestep in the model and in some occasions it would be just O2? Is there a conditional formula I could use in KINETICS with BASIC?

I was also not fully sure one how I should best represent the reduction of just O2 for this reaction? Is - formula Doc -1 CH2O +1 in KINETICS the correct way?

Here are excerpts of my script for use with phreeqc.dat:

Code: [Select]

SOLUTION_MASTER_SPECIES
    Doc           Doc              0       Doc             12

SOLUTION_SPECIES
Doc = Doc
     log_k     0
       
SOLUTION 0
pH ...
temp ...
Doc   1
O(0) 0.2
N(5) 0.1
Na ...
...

RATES
DOC_oxidation
10 k_doc_o2 = 1.57e-9
20 k_doc_NO3 = 1.67e-11 
40 mo2 = MOL("O2")
41 mNO3  = TOT("N(5)")
50 rate = k_doc_o2 * mo2/(1e-10 + mo2)
60 rate = rate + k_doc_NO3 * mNO3/(1.55e-4 + mNO3)
70 moles = rate * TIME
80 SAVE moles

KINETICS
DOC_oxidation; -formula Doc -1 CH2O +1 Ntg +1
#here I do not know how to make Ntg production conditional on NO3 consumption

Many thanks in advance for your help!

Floris

[dlparkhurst]

You have a few choices for your calculation.

(1) You can use Amm.dat, in which case, nitrogen will be reduced only as far as N2. The formula can simply be Doc -1 CH2O +1. I have added an inhibitor factor, so that NO3 will not be reduced until O2 concentration are small.

Code: [Select]

SOLUTION_MASTER_SPECIES
    Doc           Doc              0       Doc             12

SOLUTION_SPECIES
Doc = Doc
     log_k     0
       
SOLUTION 0
Doc   1
O(0) 0.2
N(5) 0.1


RATES
DOC_oxidation
10 k_doc_o2 = 1.57e-9
20 k_doc_NO3 = 1.67e-11
40 mo2 = MOL("O2")
41 mNO3  = TOT("N(5)")
50 rate = k_doc_o2 * mo2/(1e-10 + mo2)
60 rate = rate + k_doc_NO3 * mNO3/(1.55e-4 + mNO3) * 1e-6/(1e-6 + mO2)
70 moles = rate * TIME
80 SAVE moles

KINETICS
DOC_oxidation; -formula Doc -1 CH2O +1
-step 1 1e1 1e2 1e3 1e4 1e5 1e6 1e7
INCREMENTAL_REACTIONS
USER_GRAPH 1
    -headings               time O2 NO3 N2
    -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
20 GRAPH_Y MOL("O2"), TOT("N(5)"), MOL("N2")
  -end
    -active                 true

(2) You can use two RATES definitions and two kinetic reactants using phreeqc.dat, which allows nitrate to be reduced to N(-3). The first rate would be the same as above, but, the second rate would transfer N2 to Ntg (N2 -1 Ntg +1), and could simply be a fast exponential decay for the N2 concentration.

(3) You could use one rate for O2 consumption, and one for NO3 consumption. The NO3 reaction would have a formula equivalent to the following.

Code: [Select]

Doc -1.25 N -1 Ntg +0.5 CH2O +1.25