Modeling nitrogen transport and transformation with PhreeqcRM

[Zhaoyang]

Dear all,

Happy new year!

Recently, I coupled PhreeqcRM with a transport model to investigate nitrogen transport and transformation in coastal unconfined aquifers. However, I encountered a problem.

Attached figure shows the numerical model setup where BCD is the sea boundary subject to tide and OA is the land boundary. I have a total of 38986 cells in the model domain. Three reactions are considered: aerobic respiration, nitrification and denitrification. After running RM_FindComponents in PhreeqcRM, I get 12 species in total. The input file for PhreeqcRM was also attached. I can run the coupling model successfully, but the results seem unreasonable since I cannot see N2 gas production due to denitrification. However, I can see N2 gas production when modeling one cell at the sea boundary with Phreeqc (input file attached). In addition, I am not sure about if the rate expressions I adopted are reasonable. Does anyone has hints on this problem? Thanks in advance.

Best regards,
Zhaoyang

[dlparkhurst]

Here is the way I would change your script.

I am using Amm.dat, which separates ammonia from the rest of the nitrogen system.

I also defined the species N2 to be negligible, so the nitrogen (N) system consists mainly of N(5).

Amm.dat has an element Ntg, which represents unreactive N2(aq).

The formulas for the KINETIC reactions are adjusted to be consistent with the new set of elements: Amm, N, Ntg, and Doc. With these definitions, the dissolved oxygen concentration is calculated without the need to include it explicitly in the formulas.

Code: [Select]

#DATABASE D:\software\PHREEQC\database\phreeqc.dat
DATABASE D:\software\PHREEQC\database\Amm.dat
SOLUTION_MASTER_SPECIES
    Doc           Doc              0       Doc             12
#    Do            Do               0       Do              32
#    Nitrate       Nitrate          0       Nitrate         62
#    Ammh          Ammh             0       Ammh            18
    Salt          Salt             0       Salt             58
SOLUTION_SPECIES
    Doc = Doc
       log_k     0
#    Do = Do
#       log_k     0
#    Nitrate = Nitrate
#       log_k     0
#    Ammh = Ammh
#       log_k     0
    Salt  = Salt
       log_k     0
2 NO3- + 12 H+ + 10 e- = N2 + 6 H2O
-log_k 0

SOLUTION 0
     units            mg/l
     temp             25.0
     density                       1.025
     pH               7.0     charge
     pe               4
     Doc              6 mg/L                # Doc   concentration Measured in the river
#     Do               9 mg/L                # O2    concentration Measured in the river
O(0)            9
#     Nitrate          20 mg/L               # NO3-  concentration Measured in the river
     N(5)             20
#     Ammh             0.05 mg/L             # NH4+  concentration Measured in the river
     Amm              0.05 as NH4
     Salt             35  g/L               # Assumed
END

RATES
DOC_oxidation
-start
10 k = 1.5/(24*3600) # 0.01 per day
#20 o2 = TOT("Do")
20 o2 = MOL("O2")
30 rate = k * TOT("Doc") * o2/(1e-10 + o2)
40 moles = rate * TIME
50 SAVE moles
60 END
-end

Ammh_oxidation
-start
10 k = 0.1/(24*3600)  # 0.1 per day
#20 o2 = TOT("Do")
20 o2 = MOL("O2")
#30 rate = k * TOT("Ammh") * o2/(1e-10 + o2)
30 rate = k * TOT("Amm") * o2/(1e-10 + o2)
40 moles = rate * TIME
50 SAVE moles
60 END
-end

N_reduction
-start
10 k = 5000/(24*3600) # 0.05 per day
#20 o2 = TOT("Do")
20 o2 = MOL("Do")
#30 rate = k * TOT("Doc") * TOT("Nitrate")*1e-10/(o2 + 1e-10)
30 rate = k * TOT("Doc") * TOT("N(5)")*1e-10/(o2 + 1e-10)
40 moles = rate * TIME
50 SAVE moles
60 END
-end

#INCREMENTAL_REACTIONS true
USE solution 0
KINETICS
-step 1036800 in 120
DOC_oxidation
#    -formula Doc -1 Do -1  C 1 H 2  O 3
-formula Doc -1 CH2O +1
Ammh_oxidation
#    -formula Ammh -1 Do -2  N 1 H 4  O 4
-formula Amm -1 NH3 +1
N_reduction
#    -formula Doc -2.5 Nitrate -2 CH2O 2.5   Ntg  1  O  6
-formula Doc -2.5 N -2 CH2O 2.5   Ntg  1 

#EQUILIBRIUM_PHASES
#O2(g) -0.7 10

USER_GRAPH 1
    -headings               time  DO Doc N(5) C(4) N2(g) Amm
    -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("Do"), TOT("Ammh"), TOT("Doc"), TOT("Nitrate"),TOT("C(4)")
20 GRAPH_Y MOL("O2"), TOT("Doc"), TOT("N(5)"), TOT("C(4)"), TOT("Ntg")
30 GRAPH_SY TOT("Amm")
  -end


[Zhaoyang]

Thank you very much for your reply. I have compared these two different scripts carefully. However, I don't understand why there is a great deviation between the results based on these two scripts (attached Figure). As can be seen, DO concentration is not zero for the second script in the end.

In addition, I am confused about how to deal with water concentration. H2O, H and O are included if setting "tf" of RM_SetComponentH2O as 1. My question is that how I can transfer H2O mass calculated from PhreeqcRM to the transport model where the concentration unit is mass/mass. Note that H2O is indicated by pressure in the transport model.

Thanks for your help in advance.

Best regards,
Zhaoyang

[dlparkhurst]

I was careless in the SOLUTION definition and did not create the correct concentrations relative to your definition. However, there are a few other mistakes on both sides that need to be corrected.

Here is the corrected SOLUTION definition for my file:

Code: [Select]

     Doc              6 mg/L                # Doc   concentration Measured in the river
#     Do               9 mg/L                # O2    concentration Measured in the river
O(0)             9 as O
#     Nitrate          20 mg/L               # NO3-  concentration Measured in the river
     N(5)             20 as NO3
#     Ammh             0.05 mg/L             # NH4+  concentration Measured in the river
     Amm              0.05 as NH4
     Salt             35  g/L               # Assumed
END

Note that my SOLUTION needs dissolved oxygen to be converted to mol/kgw of O(0), whereas yours needs to be mol/kgw O2.

Also, there was an error in my file for N_reduction, which should have had the following line:

Code: [Select]

20 o2 = MOL("O2")
On your side, you failed to account correctly for the charge of the ions, which causes the pH calculation to be incorrect. I think these would be the correct definitions:

Code: [Select]

SOLUTION_MASTER_SPECIES
    Doc           Doc              0       Doc             12
    Do            Do               0       Do              32
    Nitrate       Nitrate-          0       Nitrate         62
    Ammh          Ammh+             0       Ammh            18
    Salt          Salt             0       Salt             58
SOLUTION_SPECIES
    Doc = Doc
       log_k     0
    Do = Do
       log_k     0
    Nitrate- = Nitrate-
       log_k     0
    Ammh+ = Ammh+
       log_k     0
    Salt  = Salt
       log_k     0
In addition, your Amm_oxidation reaction incorrectly converts to N, rather than Nitrate. This is the corrected formula:

Code: [Select]

Ammh_oxidation
#    -formula Ammh -1 Do -2  N 1 H 4  O 4
-formula Ammh -1 Do -2  Nitrate 1  H 4  O 1

The two attached files give the same results using the two approaches, but my preference is still phrqc2.pqi.