Processes > Reactive transport modelling
pH and pe changing in constant boundary solution
(1/1)
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: ---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
--- End code ---
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: ---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
--- End code ---
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.
Navigation
[0] Message Index
Go to full version