Creat a solution with dissolved oxygen and dissolved organic carbon

[Zhaoyang]

Hi, everyone,

I want to calculate nitrificaton and denitrification with Phreeqc, but I don't know how to define dissolved oxygen and dissolved organic carbon in a solution. Does anyone know how to define this in Phreeqc? Any help will be appreciated.

Code: [Select]

SOLUTION 0 
        units            mg/l
        temp             25.0
        pH               7.0     charge
        pe               12.5   
        DOC               6
        DO                9
        NO3-              20
        NH4+             0.05
ENDBest regards,
Zhaoyang

[dlparkhurst]

NH3 and dissolved oxygen will not coexist at equilibrium; it will convert to NO3-. Therefore, you must treat Ammonia as a separate "element" from N as is done in the Amm.dat database. To oxidize it to NO3-, you must define a KINETIC reaction.

Similarly, you will need a KINETIC reaction to oxidize Doc to C(4) species.

The file below uses Amm.dat to define Amm (equal to NH3) and defines a new element Doc for dissolved organic carbon. I have defined kinetic reactions that convert Doc to inorganic carbon and Amm to N species, but the rate equations and constants are not based on any real data; you must define rates appropriate for your system.

If you want to calculate denitrification, you will need yet another "element" for N2, say Ntg, and a KINETIC reaction that converts N to Ntg. Ntg (= N2) is defined in Amm.dat, but you also need to disable N2(aq) by using a small log K in a SOLUTION_SPECIES definition or removing it from the database.

Code: [Select]

SOLUTION_MASTER_SPECIES
    Doc           Doc              0     Doc             12
SOLUTION_SPECIES
Doc = Doc
    log_k     0
SOLUTION 0
        units            mg/l
        temp             25.0
        pH               7.0     charge
        pe               12.5   
        Doc               12 mg/L as C
        O(0)              9 O2(g) -0.7 10
        N(5)              20 mg/L as N
        Amm              14 mg/L as N
END
RATES
Amm_oxidation
-start
10 k = 0.1/(24*3600)  # 0.1 per day
20 o2 = TOT("O(0)")
30 rate = k * TOT("Amm") * o2/(1e-10 + o2)
40 moles = rate * TIME
50 SAVE moles
60 END
-end
DOC_oxidation
-start
10 k = 0.01/(24*3600) # 0.01 per day
20 o2 = TOT("O(0)")
30 rate = k * TOT("Doc") * o2/(1e-10 + o2)
40 moles = rate * TIME
50 SAVE moles
60 END
END
USE solution 0
KINETICS
-step 864000 in 10
Amm_oxidation
-formula Amm -1 NH3 1

DOC_oxidation
-formula Doc -1 CH2O 1
EQUILIBRIUM_PHASES
O2(g) -0.7 10

USER_GRAPH 1
    -headings               time O2(aq) Amm Doc
    -axis_titles            "Time, days" "Molality" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / (24*3600)
20 GRAPH_Y TOT("O(0)"), TOT("Amm"), TOT("Doc")
  -end



[Zhaoyang]

Thanks for your code. As mentioned before, I want to account for three reactions, i.e., aerobic respiration, nitrification and denitrification (attached below). Therefore, I modified the code from you to include denitrification. The code runs successfully, but the result show a zero concentration of N2 (Ntg). Could you please provide some help to me? Thanks for your help in advance.

Code: [Select]

SOLUTION_MASTER_SPECIES
    Doc           Doc              0     Doc             12
SOLUTION_SPECIES
Doc = Doc
    log_k     0
Ntg = Ntg # N2
     -log_k 0
     -dw 1.96e-9
     -Vm 7 # Pray et al., 1952, IEC 44. 1146
SOLUTION 0
        units            mg/l
        temp             25.0
        pH               7.0     charge
        pe               12.5   
        Doc               12 mg/L as C
        O(0)              9 O2(g) -0.7 10
        N(5)              20 mg/L as N
        AmmH+              14 mg/L as N
END
RATES
Amm_oxidation
-start
10 k = 0.1/(24*3600)  # 0.1 per day
20 o2 = TOT("O(0)")
30 rate = k * TOT("N(-3)") * o2/(1e-10 + o2)
40 moles = rate * TIME
50 SAVE moles
60 END
-end
DOC_oxidation
-start
10 k = 0.01/(24*3600) # 0.01 per day
20 o2 = TOT("O(0)")
30 rate = k * TOT("Doc") * o2/(1e-10 + o2)
40 moles = rate * TIME
50 SAVE moles
60 END
-end
N_reduction
-start
10 k = 0.005/(24*3600) # 0.005 per day
20 rate = k * TOT("Doc") * TOT("N(-3)")
30 moles = rate * TIME
40 SAVE moles
50 END
-end

USE solution 0
KINETICS
-step 8640000 in 100
Amm_oxidation
-formula AmmH+ -1 O  -2  No3- 1

DOC_oxidation
-formula Doc -1 O  -1   CO2 1

N_reduction
-formula Doc -5 No3- -4  CO2  5  Ntg  2

#EQUILIBRIUM_PHASES
#O2(g) -0.7 10

USER_GRAPH 1
    -headings               time O2(aq) Amm Doc N2 (g)
    -axis_titles            "Time, days" "Molality" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / (24*3600)
20 GRAPH_Y TOT("O(0)"), TOT("N(-3)"), TOT("Doc"),TOT("Ntg")
  -end


[dlparkhurst]

You should pay attention to the warning:

"WARNING: Could not find element in database, AmmH4+."

All of your -formula are incorrect, look at the example that I sent.

Also, I wrote:

"you also need to disable N2(aq) by using a small log K in a SOLUTION_SPECIES definition or removing it from the database. "

[Zhaoyang]

Sorry, I have gone through your example carefully, but I am still confused. Firstly, Amm in the Amm database is NH3, instead of NH4+ in my solution. So why do you define Amm in the solution? Secondly, in your example, Amm is transferred to NH3. This is not the reaction that I want to model, i.e., NH4+ should be transferred to NO3-. I don't understand why all of my -fromula are wrong. Finally, I can find Ntg in the Amm database, but not N2(aq). In this way, do I need to define N2(aq) to disable N2(aq) (attached below)? Thank you very much for your kindly help.

Code: [Select]

SOLUTION_MASTER_SPECIES
    Doc           Doc              0     Doc             12
    N2(aq)        N2(aq)           0     N2(aq)          28
SOLUTION_SPECIES
Doc = Doc
    log_k     0
N2(aq) = N2(aq)  # N2
     -log_k 0.0001

Best regards,
Zhaoyang

[dlparkhurst]

I am close to giving up on you.

This definition of N2(aq) comes from the Amm database (you can also see a definition for N(0) in SOLUTION_MASTER_SPECIES). Use this definition in your input file, which sets the log K to a much smaller number to effectively eliminate the species N2 from the N mole balance:

Code: [Select]

SOLUTION_SPECIES
2 NO3- + 12 H+ + 10 e- = N2 + 6 H2O
# -log_k 207.08
        -log_k       0.0
-delta_h -312.130 kcal
-dw 1.96e-9
-Vm 7 # Pray et al., 1952, IEC 44. 1146

So, there will be mole balance equations for Amm, N (with no N(-3) or N(0)), and Ntg, and there will be no interaction among these "elements" except by KINETIC reaction.

Do you enter "Ca+2", "K+", or "HCO3-" in SOLUTION? No. So do not enter "AmmH+". The concentrations in SOLUTION are the total amounts of elements (or redox states), regardless of the distribution of element among the possible species. The initial solution calculation will distribute Amm among its species (AmmH+ or Amm, depending on pH).

The -formula can be confusing, but you are stating the change in elemental concentrations in the chemical system due to the KINETIC reaction. The formula ignores any charge that is defined, and No is not an element name.

Code: [Select]

-formula Amm -1 NH3 +1

This formula states that for every mole of Amm that is removed from the system, 1 mole of NH3 is added to the system; the number of moles of N increases by 1 and H by 3.  Now, there is no species NH3 in the Amm database (and you are minimizing N2(aq) by using a very small log K) , so the N that is added must end up in the N species that are available, namely NO3- and NO2-. Further, to form these species, O2 must be consumed. The reaction will simulate nitrification by conversion from the Amm mole balance to the N mole balance.

You can complain, but the program works the way I think, not necessarily the way that you think.

[Zhaoyang]

Thanks for your reply. I have understood the formula for nitrification according to your explanation, but I still don't know how to define the formula for denitrification based on the formulas for nitrification and aerobic respiration. In addition, why is the formula "Doc -1 CH2O 1" for aerobic respiration? Aerobic respiration will generate CO2. Why is CO2 not included in the formula? In order to figure it out, I also looked at examples given in the manual. However, I failed to work it out (code attached). If possible, could you please further explain the formulas for these the rest two reactions? Thanks in advance.

Best regards,
Zhaoyang

Code: [Select]

SOLUTION_MASTER_SPECIES
    Doc           Doc              0     Doc             12
SOLUTION_SPECIES
    Doc = Doc
       log_k     0
    2 NO3- + 12 H+ + 10 e- = N2 + 6 H2O
       #-log_k 207.08
        -log_k       0
        -delta_h -312.130 kcal
        -dw 1.96e-9
        -Vm 7 # Pray et al., 1952, IEC 44. 1146
SOLUTION 0
     units            mg/l
     temp             25.0
     pH               7.0     charge
     pe               4
     Doc              12 mg/L as C
     O(0)             9 O2(g) -0.68 10
     N(5)             20 mg/L as N
     Amm              12 mg/L as N
END
RATES
Amm_oxidation
-start
10 k = 0.1/(24*3600)  # 0.1 per day
20 o2 = TOT("O(0)")
30 rate = k * TOT("Amm") * o2/(1e-10 + o2)
40 moles = rate * TIME
50 SAVE moles
60 END
-end
DOC_oxidation
-start
10 k = 0.01/(24*3600) # 0.01 per day
20 o2 = TOT("O(0)")
30 rate = k * TOT("Doc") * o2/(1e-10 + o2)
40 moles = rate * TIME
50 SAVE moles
60 END
-end
N_reduction
-start
10 k = 0.005/(24*3600) # 0.005 per day
20 rate = k * TOT("Doc") * TOT("Amm")* o2/(1e-10 + o2)
30 moles = rate * TIME
40 SAVE moles
50 END
-end

USE solution 0
KINETICS
-step 864000 in 20
Amm_oxidation
-formula Amm -1 NH3 1

DOC_oxidation
-formula Doc -1 CH2O 1

N_reduction
-formula Doc -5 Amm -4  CH2O  5  Ntg  2

#EQUILIBRIUM_PHASES
#O2(g) -0.7 10

USER_GRAPH 1
    -headings               time O2(aq) Amm Doc N2 (g)
    -axis_titles            "Time, days" "Molality" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / (24*3600)
20 GRAPH_Y TOT("O(0)"), TOT("Amm"), TOT("Doc"),TOT("Ntg")
  -end


 

[dlparkhurst]

Look at it this way. Given "elements" Amm and N, and that Amm represents NH3 decoupled from the rest of the N system, here is an oxidation reaction for Amm:

Amm + 2O2 = NO3- + H+ + H2O

Now calculate the change in element concentrations for this reaction, products - reactants:

Delta Amm = -1
Delta N       = +1
Delta O      = 0
Delta H      = +3

So the net reaction for the -formula definition is

-formula Amm -1 N +1 H +3 or -formula Amm -1 NH3 +1.

The nitrogen system is a little different because thermodynamic equilibrium is a poor assumption, but for carbon, it produces more reasonable results. So, if you do the same sort of net reaction for Doc (assumed CH2O) reacting to CO2, you will get

-formula Doc -1 CH2O +1

You could also write the Doc reaction using NO3-, Fe(OH)3, or SO4-2, but you would still get the same reaction Doc -1 CH2O +1.

So, the electron acceptor is not included in the reaction. The log Ks for redox reactions (or indirectly free energy) determine the electron acceptors that are used. Usually, the sequence is O2, NO3-, Fe(3), S(6), but you do not need to impose that; it follows naturally from thermodynamics once you add CH2O to the solution.

[Zhaoyang]

Thanks for your explanation. Now, I see how to define the formula in Kinetics.

However, the numerical model still does not show N2 gas because of denitrification. As you can see from the attached figure, O2 concentration decreases until to zero, but N2 gas concentration is always zero. As expected, denitrification will happen when O2 concentration is zero, but N(5) and Doc concentration are not zero. I have adjusted the constants in Rates block, but N2 gas is always zero. Could you give me some hints (code attached)?

In addition, I have a question about implementing Kinetics in the coupling model. When concentrations are transferred from the transport model to PhreeqcRM, it will do the reaction instantaneously. This seems contradictory with the Kinetics (taking some time).

Thanks in advance.

Code: [Select]

SOLUTION_MASTER_SPECIES
    Doc           Doc              0     Doc             12
SOLUTION_SPECIES
    Doc = Doc
       log_k     0
    2 NO3- + 12 H+ + 10 e- = N2 + 6 H2O
       #-log_k 207.08
        -log_k       0
        -delta_h -312.130 kcal
        -dw 1.96e-9
        -Vm 7 # Pray et al., 1952, IEC 44. 1146
SOLUTION 0
     units            mg/l
     temp             25.0
     pH               7.0     charge
     pe               4
     Doc              6 mg/L               # Doc concentration Measured in the river
     O(0)             9 O2(g) -0.68 10   # O2 concentration Measured in the river
     N(5)             20 mg/L             # No3- concentration Measured in the river
     Amm              0.05 mg/L            # NH4+ concentration Measured in the river
END
RATES
Amm_oxidation
-start
10 k = 0.1/(24*3600)  # 0.1 per day
20 o2 = TOT("O(0)")
30 rate = k * TOT("Amm") * o2/(1e-10 + o2)
40 moles = rate * TIME
50 SAVE moles
60 END
-end
DOC_oxidation
-start
10 k = 0.1/(24*3600) # 0.01 per day
20 o2 = TOT("O(0)")
30 rate = k * TOT("Doc") * o2/(1e-10 + o2)
40 moles = rate * TIME
50 SAVE moles
60 END
-end
N_reduction
-start
10 k = 5/(24*3600) # 0.05 per day
20 rate = k * TOT("Doc") * TOT("N(5)")
30 moles = rate * TIME
40 SAVE moles
50 END
-end

USE solution 0
KINETICS
-step 864000 in 20
Amm_oxidation
    -formula Amm -1 NH3 1

DOC_oxidation
    -formula Doc -1 CH2O 1

N_reduction
    -formula Doc -1 N  -2   CH2O 1   Ntg  1

#EQUILIBRIUM_PHASES
#O2(g) -0.7 10

USER_GRAPH 1
    -headings               time O2(aq) Amm Doc N2(g) N(5)  C(4)
    -axis_titles            "Time, days" "Molality" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / (24*3600)
20 GRAPH_Y TOT("O(0)"), TOT("Amm"), TOT("Doc"),TOT("N(0)"),TOT("N(5)"),TOT("C(4)")
  -end


[dlparkhurst]

No, I don't expect N(0) to accumulate in solution. You will have to think this through for yourself.

And your reaction should be something like

2NO3- + 2.5Doc = Ntg + 2.5CO2 + 2OH- + 1.5H2O

delta N = -2
delta Ntg = 1
delta Doc = -2.5
delta O = 2.5
delta H = 5

or -formula N -2 Ntg 1 Doc -2.5 CH2O 2.5

[Zhaoyang]

Thanks, the reason is that N2 is gas, so I should use SR.

[dlparkhurst]

No.

Also note modification to previous post:

And your reaction should be something like

2NO3- + 2.5Doc = Ntg + 2.5CO2 + 2OH- + 1.5H2O

delta N = -2
delta Ntg = 1
delta Doc = -2.5
delta O = 2.5
delta H = 5

or -formula N -2 Ntg 1 Doc -2.5 CH2O 2.5

[Zhaoyang]

Thanks for reminding about the formula.
N2 does not accumulate in the solution is because of changing a small logk value to block the reaction.
N2 gas is observed to increase when using SI("Ntg(g)") (Attached figure), but this increasing of N2 gas is not caused by denitrification. Attached figure shows that DOC and NO3- concentrations do not decrease when O2 concentration is zero. It doesn't make sense since denitrification is expected to occur.

[dlparkhurst]

You are right that "N(0)" will not accumulate in solution because you have effectively disabled the formation of N2(aq) by the SOLUTION_SPECIES definition.

I not sure that you understand that Ntg does accumulate in solution because it is produced in the KINETIC reaction N_reduction.

You really should look at the output file. Assuming you are using a file similar to the one below, you would see that the N_reduction reaction is proceeding (look at output just before --Solution composition--). However, the rates for  Amm_oxidation and N_reduction are about 100 times smaller than the rate for Doc oxidation (look carefully at the scales of the graphs).

As a demonstration, I have added an inhibitor factor, which inhibits the reduction of nitrate until dissolved oxygen is small.

Code: [Select]

SOLUTION_MASTER_SPECIES
    Doc           Doc              0     Doc             12
SOLUTION_SPECIES
    Doc = Doc
       log_k     0
    2 NO3- + 12 H+ + 10 e- = N2 + 6 H2O
       #-log_k 207.08
        -log_k       0
        -delta_h -312.130 kcal
        -dw 1.96e-9
        -Vm 7 # Pray et al., 1952, IEC 44. 1146
SOLUTION 0
     units            mg/l
     temp             25.0
     pH               7.0     charge
     pe               4
     Doc              6 mg/L               # Doc concentration Measured in the river
     O(0)             9 O2(g) -0.68 10     # O2 concentration Measured in the river
     N(5)             20 mg/L              # No3- concentration Measured in the river
     Amm              0.05 mg/L            # NH4+ concentration Measured in the river
END
RATES
Amm_oxidation
-start
10 k = 0.1/(24*3600)  # 0.1 per day
20 o2 = TOT("O(0)")
30 rate = k * TOT("Amm") * o2/(1e-10 + o2)
40 moles = rate * TIME
50 SAVE moles
60 END
-end
DOC_oxidation
-start
10 k = 0.1/(24*3600) # 0.01 per day
20 o2 = TOT("O(0)")
30 rate = k * TOT("Doc") * o2/(1e-10 + o2)
40 moles = rate * TIME
50 SAVE moles
60 END
-end
N_reduction
-start
10 k = 5/(24*3600) # 0.05 per day
20 rate = k * TOT("Doc") * TOT("N(5)") * 1e-10/(TOT("O(0)") + 1e-10)
30 moles = rate * TIME
40 SAVE moles
50 END
-end
END
INCREMENTAL_REACTIONS true
USE solution 0
KINETICS
-step 864000 in 20
Amm_oxidation
    -formula Amm -1 NH3 1
DOC_oxidation
    -formula Doc -1 CH2O 1
N_reduction
    -formula Doc -2.5 N  -2   CH2O 2.5   Ntg  1
USER_GRAPH 1
    -headings               time O2(aq) Amm Doc N(5)  C(4) Ntg
    -axis_titles            "Time, days" "Molality" "Ntg, molality"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / (24*3600)
20 GRAPH_Y TOT("O(0)"), TOT("Amm"), TOT("Doc"),TOT("N(5)"),TOT("C(4)")
30 GRAPH_SY TOT("Ntg")
  -end
USER_GRAPH 2
    -headings               time  DOC_oxidation Amm_oxidation N_reduction
    -axis_titles            "Days" "Rate DOC_ox, mol/d" "Rate Amm_ox and N_red, mol/d"
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X TOTAL_TIME / (24*3600)
20 delta = KIN_TIME / (24*3600)
30 GRAPH_Y KIN_DELTA("DOC_oxidation") / delta
40 GRAPH_SY KIN_DELTA("Amm_oxidation") / delta
50 GRAPH_SY KIN_DELTA("N_reduction") / delta 
  -end
    -active                 true
END


[Zhaoyang]

Thanks for your reply. I have got reasonable results with your help.

Two quick questions. The manual states "The -steps identifier is used only during batch-reaction calculations; it is not needed for transport calculations." This means I don't need to indicate the steps in KINETICS when implementing it in transport calculations. Do I understand it correctly? In other words, do I need to set this time step for the steps in KINETICS?

In addition, do I need to define a rate for O2 consumption? I see in some papers where Multiple-Monod kinetic is used to model these three reactions including a rate for O2 consumption (attached figure).

Thanks for your help in advance.

Best regards,
Zhaoyang

[dlparkhurst]

I'm glad you read the manual. It is correct.

It is easy to accidentally define a batch reaction by including KINETICS and SOLUTION between END statements. Without an explicit -steps definition in KINETICS, the time step defaults to 1 sec. If you are careful and do not define a batch reaction with KINETICS, no time step is needed. The time step for kinetics during TRANSPORT is defined by the time step definition in the TRANSPORT data block.

The way the ammonia oxidation reaction is defined, oxygen is automatically consumed to generate NO3-. Similarly, the addition of CH2O automatically consumes the favored electron acceptor, O2 if it is available. In denitrification CH2O will consume the oxygen remaining from NO3- when you remove N. If no other electron acceptor is available, addition of CH2O will generate CO2 and CH4.

If you were to separate O2(aq) into a separate redox state, say Ox, then you would need to include Ox explicitly in rate equations. However, it will get more and more difficult to keep the reactions consistent with thermodynamics. Adding CH2O to a solution of O2, Fe(OH)3, MnO2, H3AsO4, UO2+2, NO3-, and SO4-2, with no alteration of the databases, will mostly reduce O2, but thermodynamics says you will start to generate reduced forms of the other elements. As O2(aq) is completely removed, it will be difficult to know which electron acceptors are used next if a KINETIC approach is used for all valence states of all redox elements.

[Zhaoyang]

Thanks. Actually, I am new in the reaction model. I uesd to work on modeling groundwater flow and salt transport in coastal unconfined aquifers. To extend my research, I want to inculde the reaction model.

For the first point, I assume it is also suitable when coupling a trasport model with PhreeqcRM. In other words, the -steps in KINETICS is still the same as that in the transport model when KINETICS is ture.

For the second point, Phreeqc adopts a different strategy to consider aerobic respiration, nitrification and denitrification, instead of Multiple-Monod kinetic. In this way, I was wondering if the rate expressions in Phreeqc are the same as those defined by Multiple-Monod kinetic. If not, could you recommend some references where I can find the rate expressions for these three reactions?

Thanks a lot and have a nice day.

[dlparkhurst]

No, I think you are on your own now.

You cited the Monod reactions in a previous post. I haven't worked seriously with the nitrogen system in years, so you will have to search the literature.

You can either implement the Monod oxygen reactions in PHREEQC by adding another new element, or you can find another program to do the Monod calculations for comparison.

[MichaelZ20]

http://dx.doi.org/10.1016/j.chemgeo.2013.12.003
Integrated modeling of biogeochemical reactions and associated isotope fractionations at batch scale: A tool to monitor enhanced biodenitrification applications
P Rodr?guez-Escales, BM van Breukelen, G Vidal-Gavilan, A Soler, ...
Chemical Geology 365, 20-29   

[Zhaoyang]

Thank you so much. I have figured it out.