pH and pe changing in constant boundary solution

[Daehyun]

Hello,

I am currently using a PHREEQC transport model in which the initial groundwater is assigned as Solution 0 (constant boundary), and cells 1?9 represent bentonite porewater saturated with groundwater. The goal is to simulate simple diffusion between the boundary groundwater and the bentonite porewater.

When I run the model, the ion concentrations in Solution 0 remain exactly the same as the initial values throughout the simulation. However, pH and pe values in Solution 0 change during the simulation.

Upon further investigation, I suspect this behavior is related to my suppression of CH₄ formation. In the site groundwater data, methane is absent, so I intentionally suppressed methane generation by setting the equilibrium constant (log_k) for the CO₂ → CH₄ reaction to a very low value (e.g., -100). I believe this modification changes the way PHREEQC recalculates pH and pe during each transport step, possibly by altering the dominant redox couples used in the calculations.

I would like to understand:

Is this change in pH and pe an expected consequence of removing or heavily suppressing the CH₄ formation reaction?

What is the recommended approach to keep pH and pe constant at the boundary while preventing CH₄ generation? Would it be better to suppress the reaction using kinetic constraints rather than modifying log_k?

Is there a more appropriate way to represent ?no methane formation? without affecting redox/pH recalculations in the boundary solution?

Any insights or suggestions would be greatly appreciated.

Thank you.

Code: [Select]

SOLUTION_MASTER_SPECIES
    C             CO3-2            2     HCO3            12.0111
    C(4)          CO3-2            2     HCO3           
    Methane       MethaneH4        0     MethaneH4       12.0111
    Doc           Doc              0     C1H2O1          30
    Ntg           Ntg              0     Ntg             28.0134
    Cs            Cs+              0     Cs              132.9054
    U             UO2+2            0     U               238.0289
    I             I-               0     I               126.9045
    Nitrate       NitrateO3-       0     Nitrate         14.0067
    Nitrate(3)    NitrateO2-       0     Nitrate       
    Nitrate(+5)   NitrateO3-       0     Nitrate       

SOLUTION_SPECIES
CO3-2 = CO3-2
    log_k     0
    -gamma    5.4 0
    -dw       9.55e-10
    -millero  -8.74 0.3 -0.004064 5.65 0 0
CO3-2 + H+ = 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-09
    -millero  21.07 0.185 -0.002248 2.29 -0.006644 -3.667e-06
MethaneH4 = MethaneH4
    log_k     0
    delta_h   -61.039 kcal
    -dw       1.85e-09
Doc = Doc
    log_k     0
Ntg = Ntg
    log_k     0
    -dw       1.96e-05
Cs+ = Cs+
    log_k     0
    -gamma    2.5 0
    -dw       2.062e-09
UO2+2 = UO2+2
    log_k     0
    -gamma    4.5 0
I- = I-
    log_k     0
    -gamma    3 0
    -dw       6.6455e-11
Sr+2 = Sr+2
    log_k     0
    -gamma    5.26 0.121
    -dw       5e-09
NitrateO3- = NitrateO3-
-gamma 3 0
-Vm 6.32 6.78 0 -3.06 0.346 0 0.93 0 -0.012 1
-viscosity 8.37e-2 -0.458 1.54e-2 0.34 1.79e-2 5.02e-2 0.7381
-dw 1.9e-9 104 1.11
NitrateO3- + 2H+ + 2e- = NitrateO2- + H2O
-log_k 28.57
-delta_h -43.76 kcal
-gamma 3 0
-Vm 5.5864 5.859 3.4472 -3.0212 1.1847 # supcrt
-dw 1.91e-9
CO3-2 + 10 H+ + 8 e- = CH4 + 3 H2O
-log_k -100
-delta_h -61.039 kcal
-Vm 9.01 -1.11 0 -1.85 -1.5 # Hnedkovsky et al., 1996, JCT 28, 125
-dw 1.85e-9


PHASES
Soc
    Doc = Doc
    log_k     0
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
H2S(g)
    H2S = H+ + HS-
    log_k     -7.93
    delta_h   9.1 kJ
    -analytical_expression -45.07 -0.02418 0 17.9205 0 0
    -T_c      373.2
    -P_c      88.2
    -Omega    0.1
Montmorillonite-Ca
    Ca0.3Mg0.6Al1.4Si4O10(OH)2:4.45H2O + 6H+ = 1.4Al+3 + 0.3Ca+2 + 0.45H2O + 4H4SiO4 + 0.6Mg+2
    log_k     6.15
    delta_h   -134.134 kJ
    -analytical_expression -17.34927 0 7006.307 0 0 0
    -Vm       220.76 cm3/mol
Albite-high
    NaAlSi3O8 + 4H+ + 4H2O = Al+3 + 3H4SiO4 + Na+
    log_k     4.14
    delta_h   -95.622 kJ
    -analytical_expression -12.61226 0 4994.685 0 0 0
Anorthite
    CaAl2Si2O8 + 8H+ = 2Al+3 + Ca+2 + 2H4SiO4
    log_k     25.31
    delta_h   -314.358 kJ
    -analytical_expression -29.76316 0 16420.06 0 0 0
Quartz
    SiO2 + 2H2O = H4SiO4
    log_k     -16.29
    delta_h   21.166 kJ
    -analytical_expression -0.03187597 0 -1105.577 0 0 0
    -Vm       22.69 cm3/mol
Cristobalite
    SiO2 + 2H2O = H4SiO4
    log_k     -12.31
    delta_h   16.5 kJ
    -analytical_expression -0.2693241 0 -861.855 0 0 0
Calcite
    CaCO3 = CO3-2 + Ca+2
    log_k     -8.48
    delta_h   -10.62 kJ
    -analytical_expression -10.34054 0 554.7212 0 0 0
    -Vm       36.93 cm3/mol
Pyrite
    FeS2 + 2H+ + 2e- = Fe+2 + 2HS-
    log_k     -16.82
    delta_h   50.735 kJ
    -analytical_expression -7.93161 0 -2650.074 0 0 0
    -Vm       23.94 cm3/mol
Kaolinite
    Al2(Si2O5)(OH)4 + 6H+ = 2Al+3 + H2O + 2H4SiO4
    log_k     6.5
    delta_h   -169.718 kJ
    -analytical_expression -23.23332 0 8864.989 0 0 0
Strontianite
    Sr(CO3) = CO3-2 + Sr+2
    log_k     -9.27
    delta_h   -0.366 kJ
    -analytical_expression -9.33412 0 19.11751 0 0 0
Chalcedony
    SiO2 + 2H2O = H4SiO4
    log_k     -3.7281
    delta_h   31.4093 kJ
    -analytical_expression -9.0068 0.0093241 4053.5 -1.083 -750770 0
HSaponite-Ca
    Ca0.17Mg3Al0.34Si3.66O10(OH)2:4.45H2O + 7.36H+ = 0.34Al+3 + 0.17Ca+2 + 1.81H2O + 3.66H4SiO4 + 3Mg+2
    log_k     28.36
    delta_h   -239.662 kJ
    -analytical_expression -13.62698 0 12518.42 0 0 0
    -Vm       223.01 cm3/mol
HSaponite-FeCa
    Ca0.17Mg2FeAl0.34Si3.66O10(OH)2:4.45H2O + 7.36H+ = 0.34Al+3 + 0.17Ca+2 + Fe+2 + 1.81H2O + 3.66H4SiO4 + 2Mg+2
    log_k     27.97
    delta_h   -235.847 kJ
    -analytical_expression -13.34862 0 12319.15 0 0 0
    -Vm       225.59 cm3/mol
HSaponite-FeK
    K0.34Mg2FeAl0.34Si3.66O10(OH)2:1.96H2O + 7.36H+ + 0.68H2O = 0.34Al+3 + Fe+2 + 3.66H4SiO4 + 0.34K+ + 2Mg+2
    log_k     28.11
    delta_h   -242.507 kJ
    -analytical_expression -14.3754 0 12667.02 0 0 0
    -Vm       179.69 cm3/mol
HSaponite-FeMg
    Mg0.17Mg2FeAl0.34Si3.66O10(OH)2:4.61H2O + 7.36H+ = 0.34Al+3 + Fe+2 + 1.97H2O + 3.66H4SiO4 + 2.17Mg+2
    log_k     28.07
    delta_h   -235.257 kJ
    -analytical_expression -13.14526 0 12288.33 0 0 0
    -Vm       223.85 cm3/mol
HSaponite-FeNa
    Na0.34Mg2FeAl0.34Si3.66O10(OH)2:3.84H2O + 7.36H+ = 0.34Al+3 + Fe+2 + 1.2H2O + 3.66H4SiO4 + 2Mg+2 + 0.34Na+
    log_k     27.72
    delta_h   -246.878 kJ
    -analytical_expression -15.53117 0 12895.34 0 0 0
    -Vm       212.99 cm3/mol
HSaponite-K
    K0.34Mg3Al0.34Si3.66O10(OH)2:1.96H2O + 7.36H+ + 0.68H2O = 0.34Al+3 + 3.66H4SiO4 + 0.34K+ + 3Mg+2
    log_k     28.49
    delta_h   -246.322 kJ
    -analytical_expression -14.66376 0 12866.29 0 0 0
    -Vm       177.11 cm3/mol
HSaponite-Mg
    Mg0.17Mg3Al0.34Si3.66O10(OH)2:4.61H2O + 7.36H+ = 0.34Al+3 + 1.97H2O + 3.66H4SiO4 + 3.17Mg+2
    log_k     28.48
    delta_h   -239.062 kJ
    -analytical_expression -13.40186 0 12487.08 0 0 0
    -Vm       221.08 cm3/mol
HSaponite-Na
    Na0.34Mg3Al0.34Si3.66O10(OH)2:3.84H2O + 7.36H+ = 0.34Al+3 + 1.2H2O + 3.66H4SiO4 + 3Mg+2 + 0.34Na+
    log_k     28.03
    delta_h   -250.288 kJ
    -analytical_expression -15.81858 0 13073.45 0 0 0
    -Vm       210.4 cm3/mol
Cu-canister
    Cu = Cu+2 + 2e-
    log_k     -11.49

EXCHANGE_SPECIES
X- = X-
    log_k     0
Na+ + X- = NaX
    log_k     0
    -gamma    4.08 0.082
K+ + X- = KX
    log_k     0.7
    delta_h   -4.3 kJ
    -gamma    3.5 0.015
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
0.5CaX2 + Cs+ = CsX + 0.5Ca+2
    log_k     2.32
CaX2 + Sr+2 = SrX2 + Ca+2
    log_k     -0.05

SURFACE_MASTER_SPECIES
    Hfo_g         Hfo_gOH     
    Hfo_h         Hfo_hOH     

SURFACE_SPECIES
Cs+ + Hfo_gOH = Hfo_gOCs + H+
    log_k     -5.62
Hfo_gOH + Sr+2 = Hfo_gOSr+ + H+
    log_k     -6.85
Hfo_gOH = Hfo_gOH
    log_k     0
Hfo_hOH = Hfo_hOH
    log_k     0
H+ + Hfo_hOH = Hfo_hOH2+
    log_k     8.35
Hfo_hOH = Hfo_hO- + H+
    log_k     -9.59
Hfo_gOH = Hfo_gO- + H+
    log_k     -7.65

RATES
Montmorillonite-Ca
10 rem acid solution parameters
11 a1=1.94984E-13
12 E1=48000
13 n1=0.220
20 a2=3.89045E-15
21 E2=48000
30 rem base solution parameters
31 a3=3.89045E-15
32 E3=48000
33 n2=-0.130
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Montmorillonite-Ca")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1 
85 Rate3=a3*EXP(-E3/R/TK)*ACT("OH-")^n2         
90 Rate=(Rate1+ Rate3)*(1-Sr_mineral)*SA*0.1
100 moles= rate*Time
200 save moles
-end


Cristobalite
-start
20 rem neutral solution parameters
21 a2=5.88844E-13
22 E2=74500
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Cristobalite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
90 Rate=(Rate2)*(1-Sr_mineral)*SA*0.1
100 moles= rate*Time
200 save moles
-end

Albite-high
-start
10 rem acid solution parameters
11 a1=6.91831E-11
12 E1=65000
13 n1=0.457
20 rem neutral solution parameters
21 a2=2.75423E-13
22 E2=69800
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Albite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  #acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
90 Rate=(Rate1+Rate2)*(1-Sr_mineral)*SA*0.1
100 moles= rate*Time
200 save moles
-end

Quartz
-start
20 rem neutral solution parameters
21 a2=3.98107E-14
22 E2=90900
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Quartz")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
90 Rate=(Rate2)*(1-Sr_mineral)*SA*0.1
100 moles= rate*Time
200 save moles
-end


Calcite
-start
10 rem acid solution parameters
11 a1=0.501187
12 E1=14400
13 n1=1
20 rem neutral solution parameters
21 a2=1.54882E-6
22 E2=23500
30 rem basic dependence parameters
31 a3=0.000331
32 E3=35400
33 n2=1
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Calcite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
75 Rate1=a1*EXP(-E1/R/TK)*ACT("H+")^n1  #acid rate expression
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
85 Rate3=a3*EXP(-E3/R/TK)*ACT("OH")^n2    #base rate expression
90 Rate=(Rate1+Rate2+Rate3)*(1-Sr_mineral)*SA*0.1
100 moles= rate*Time
200 save moles
-end

Pyrite
-start
30 rem Neutral dependence parameters
31 a2=2.81838E-05
32 E2=56900
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Pyrite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
80 Rate2=a2*EXP(-E2/R/TK)               #neutral rate expression
90 Rate=Rate2*(1-Sr_mineral)*SA*0.1
100 moles= rate*Time
200 save moles
-end

Kaolinite
-start
10 rem solution parameters
12 k1=1E-14
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Kaolinite")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
90 Rate=k1*(1-Sr_mineral)*SA*0.1
100 moles= rate*Time
200 save moles
-end


Doc_degradation
-start
 2 k_O2 = 45e-6
 3 k_SO4 = 3e-9
10 S_Doc = 1e-8
20 mO2 = MOL("O2")
30 mSO4 = MOL("SO4-2")
40 mDoc = TOT("Doc")
50 rate = (mDoc/(1E-6 + mDoc))*((k_O2*(2.5E-4 + mO2))+(k_SO4*(1e-6+mSO4))*(1e-6/(1e-6 + mO2)))
60 dS = rate * time
70 save dS
-end

Cu-canister
-start
10 rem solution parameters
12 k1=1.1E-8
36 rem rate=0 if no minerals and undersaturated
40 SR_mineral=SR("Cu-canister")
41 if (M<0) then goto 200
42 if (M=0 and SR_mineral<1) then goto 200
43 if (M0<=0) then SA=PARM(1) else SA=PARM(1)*(M/M0)^0.67
50 if (SA<=0) then SA=1
60 R=8.31451
90 Rate=k1*(1-Sr_mineral)*SA*0.1
100 moles= rate*Time
200 save moles
-end

SOLUTION 0
    temp      14.8
    pH        9.05
    pe        -7.4
    redox     pe
    units     mol/kgw
    density   1
    Al        2.728e-06
    C(-4)     0.0012996
    Ca        0.00014222
    Cl        5.049e-05
    Doc       0.01
    F         0.0004263
    Fe        1.2356e-07
    K         8.44e-06
    Mg        1.1932e-05
    Na        0.0016485
    S         6.038e-05
    Si        0.00012483
    Sr        2.659e-06
    U         6.302e-09
    Nitrate   1.0967e-05
    -water    0.2 # kg

SAVE SOLUTION 0
END

#KURT groundwater quility data (DB-3)
SOLUTION 100
    temp      14.8
    pH        9.05
    pe        -7.4
    redox     pe
    units     mol/kgw
    density   1
    Al        2.728e-06
    C(-4)     0.0012996
    Ca        0.00014222
    Cl        5.049e-05
    Doc       0.01
    F         0.0004263
    Fe        1.2356e-07
    K         8.44e-06
    Mg        1.1932e-05
    Na        0.0016485
    S         6.038e-05
    Si        0.00012483
    Sr        2.659e-06
    U         6.302e-09
    Nitrate   1.0967e-05
    -water    0.2 # kg

GAS_PHASE 100
    -fixed_pressure
    -pressure 1
    -volume 0.2
    -temperature 25
    CO2(g)    0.0004
    Ntg(g)    0.79
    O2(g)     0.2
    H2S(g)    0


EQUILIBRIUM_PHASES 100
    Albite-high 0 0.549
    Calcite   0 0.1
    Cristobalite 0 2.147
    Kaolinite 0 0
    Montmorillonite-Ca 0 1.566
    Pyrite    0 0
    Quartz    0 0.333

REACTION_TEMPERATURE 1-10
    60

REACTION_PRESSURE 1-10
    50

EXCHANGE 2-10
    CaX2 0.03085

SURFACE 2-10
    Hfo_gOH    0.00742    70       100
    Hfo_hOH    0.00494

SAVE SOLUTION 1-10

SAVE GAS_PHASE 1-10

END

KINETICS 1
Doc_degradation
    -formula  Doc  -1 CH2O  1
    -m        1e-05
    -m0       1e-05
    -tol      1e-08

KINETICS 2-9
Montmorillonite-Ca
    -formula  Ca0.3Mg0.6Al1.4Si4O10(OH)2:4.45H2O  1
    -m        0.1566
    -m0       0.1566
    -parms    5600
    -tol      1e-08
Calcite
    -formula  CaCO3  1
    -m        0.01
    -m0       0.01
    -parms    1
    -tol      1e-08
Albite-high
    -formula  NaAlSi3O8  1
    -m        0.0549
    -m0       0.0549
    -parms    1
    -tol      1e-08
Quartz
    -formula  SiO2  1
    -m        0.0333
    -m0       0.0333
    -parms    1
    -tol      1e-08
Cristobalite
    -formula  SiO2  1
    -m        0.2147
    -m0       0.2147
    -parms    1
    -tol      1e-08
KINETICS 10
Cu-canister
    -formula  Cu  1
    -m        15.74
    -m0       15.74
    -parms    0.00152
    -tol      1e-08
-steps       1 in 1 steps # seconds
-step_divide 1
-runge_kutta 3
-bad_step_max 1000
-cvode true
-cvode_steps 100
-cvode_order 5

TRANSPORT
    -cells                 10
    -shifts                1
    -time_step             63072000 # seconds
    -flow_direction        diffusion_only
    -boundary_conditions   constant closed
    -lengths               10*0.1
    -diffusion_coefficient 1e-10
    -thermal_diffusion     2   1e-10
    -multi_d               false

INCREMENTAL_REACTIONS true
END


[dlparkhurst]

SOLUTION 0 should not change during transport if you don't redefine it and do not have any reactants that are reacting with solution 0.

However, SOLUTION 0 will react one time to produce redox equilibrium with your first call to RUN_CELLS. To see what is happening, you can reproduce that reaction in a PHREEQC run if you include MIX; 0 1.0. With this definition, you can see how SOLUTION 0 changes when it reacts to redox equilibrium.

Removing CH4 from the calculations seems a little extreme. If you post your SOLUTION 0, I could tell why you are getting the changes that you see.

[Daehyun]

hank you for your explanation.
I understand that SOLUTION 0 reacts once to reach redox equilibrium during the first call to RUN_CELLS.

In my case, the field data for SOLUTION 0 comes from ~500 m depth groundwater where methanogenic bacteria are absent, and CH₄ has never been detected in water quality analyses. Methane formation from CO₂ generally requires microbial metabolism, so in my study conditions, the CH₄?HCO₃⁻ interconversion is not realistic.

However, during the simulation, when SOLUTION 0 and SOLUTION 1 mix, PHREEQC still performs a redox adjustment between CH₄ and HCO₃⁻. This results in CH₄ being oxidized to HCO₃⁻ (or vice versa), which alters the dissolved carbon speciation and lowers the pH. I tried setting the log_k for CH₄ formation to a very low value to suppress the reaction, but this only shifted most dissolved carbon to CO₃?⁻ and HCO₃⁻.

My intention is to completely exclude CH₄ from the redox system so that it does not participate in any equilibrium calculations at all, while keeping the pH and carbon speciation realistic.

Could you advise on the best approach?

Code: [Select]

SOLUTION 0
    temp      14.8
    pH        9.05
    pe        -7.4
    redox     pe
    units     mol/kgw
    density   1
    Al        2.728e-06
    C(-4)     0.0012996
    Ca        0.00014222
    Cl        5.049e-05
    Doc       0.01
    F         0.0004263
    Fe        1.2356e-07
    K         8.44e-06
    Mg        1.1932e-05
    Na        0.0016485
    S         6.038e-05
    Si        0.00012483
    Sr        2.659e-06
    U         6.302e-09
    Nitrate   1.0967e-05
    -water    0.2 # kg


[dlparkhurst]

Why did you include C(-4) in SOLUTION 0 if there is no methane? Did you mean C(4)? And your pe is -7, which indicates a strongly reducing environment, and generates about 1 umol/kgw H2(aq).

You have Doc entered, which tells me there is the possibility of redox reactions. Nitrate might be reduced, or U might precipitate.

Your SOLUTION 0 is not consistent with your statements in the post, and I don't understand your redox conditions or the reactions you are trying to simulate. I think you need to do some batch calculations with PHREEQC and possibly some 1D transport modeling to understand clarify your chemical system before you use PhreeqcRM.