Impurity of CO2 is not affecting final result

[nurlan]

Dear Mr Parkhurst,
I am checking if the impurity of CO2 is affecting the mineral dissolution and precipitation,
For that I chose the H2S, CH4, SO2 and NO2 (the last two didn't work), each having portion of 5% and 95 % CO2. Simply took the 5% of Partial pressure for each mixing gas. Below I attached the model code. I think CH4 or H2S is not reacting or at very low amount because delta is too small. I checked many other papers, even 1% of impurity is affecting their results. COuldn't get where I should change, maybe I need to specify some reactions in the database.
Also, I couldn't figure out why SOx, NOx gases can't react.
Thank you a lot.

Code: [Select]

DATABASE c:\phreeqc\database\llnl.dat
SOLUTION 1 #solution into which CO2 is dissolved

temp 40
pH 6.55
pe -3.1
redox pe
units mol/kgw
density 1.1596
pressure 88.8

Al          1.153e-08 
C                 9.941e-03 
Ca                3.857e-02 
Cl                5.237e+00 
Fe                6.842e-04 
K                 4.024e-02 
Mg                6.319e-02 
Na                4.903e+00 
S                 1.170e-06 
Si                1.843e-04 
Sr                7.824e-04 
water 1
#Equilibrium_phases 1
# dawsonite 0 0
#CO2(g) 1.727
#ch4(g) 0.647
GAS_PHASE 1
-fixed_pressure
-pressure 88.8
-volume 10
CO2(g) 84.36
ch4(g) 4.44

 KINETICS 1
-cvode
Calcite
 -parms 1.549e-6 5.012e-1
 -m0 21.90
 -m 21.9
Quartz
-m0    158.2
-parms  60.08  0.111 #https://escholarship.org/content/qt3h77d5bj/qt3h77d5bj_noSplash_925369b5813359fe46019e017e739492.pdf
Kaolinite
-m0    3.1
-parms  258.16   13.2 #https://escholarship.org/content/qt3h77d5bj/qt3h77d5bj_noSplash_925369b5813359fe46019e017e739492.pdf

Albite
-m0 9
-parms 263.02  0.164 #https://escholarship.org/content/qt3h77d5bj/qt3h77d5bj_noSplash_925369b5813359fe46019e017e739492.pdf
#Dawsonite
#-m0     0.0
       # -parms 144.0  9.8 #hallevang
Magnesite
  -m0     0.177
  -parms   84.31  0.8 # https://www.sciencedirect.com/science/article/pii/S0883292714002315#b0230

Muscovite
-m0    0
-parms   398.71 0.68
-steps 10*315569520
INCREMENTAL_REACTIONS true

RATES
calcite
-start
 40 SR_ca = SR("Calcite")
 45 if (M <= 0 and SR_ca < 1) then goto 200
 50 ek1 = exp(-23.5e3/8.314472*(1/TK - 1/298.15))
 52 ek2 = exp(-14.4e3/8.314472*(1/TK - 1/298.15))
 56 SA = 0.00004
 60 r1=ek1*PARM(1)
 65 r2=ek2*PARM(2)*act("H+")
 70 rate = SA*(M/M0)^0.67 *(1 - SR_ca)*(r1+r2)
 80 moles=rate*TIME
 90 PUT(rate,3)
 200 SAVE moles
 -end
Quartz
-start
10 reacS = parm(1)*parm(2)*m
13 n = 0.0
15 logkH = -13.99
    16 Ea = 87700
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = (act("H+")^n)*kHT/1
23 if SR("Quartz") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Quartz")))^2)
    24 k_n = 1
25 gamma_n = 5e10
26  lamda = 1
    30  if SR("Quartz") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Quartz"))*time
    40  if SR("Quartz") > 1 then moles = (-reacS*lamda*k_p*(SR("Quartz") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end
Kaolinite
-start
10 reacS = parm(1)*parm(2)*m
13 n = 0.777
15 logkH = -11.31
    16 Ea = 65900
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = (act("H+")^n)*kHT/10
23 if SR("Kaolinite") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Kaolinite")))^2)
    24 k_n = 1
25 gamma_n = 2e10
26  lamda = 1
    30  if SR("Kaolinite") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Kaolinite"))*time
    40  if SR("Kaolinite") > 1 then moles = (-reacS*lamda*k_p*(SR("Kaolinite") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end


Albite
-start
10 reacS = parm(1)*parm(2)*m
13 n = 0.457
15 logkH = -10.16
    16 Ea = 65000
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = (act("H+")^n)*kHT/1
23 if SR("Albite") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Albite")))^2)
    24 k_n = 1
25 gamma_n = 2e10
26  lamda = 1
    30  if SR("Albite") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Albite"))*time
    40  if SR("Albite") > 1 then moles = (-reacS*lamda*k_p*(SR("Albite") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end

Dawsonite
-start
10 reacS = parm(1)*parm(2)*m
13 n = 0.98
15 logkH = -4.5
    16 Ea = 63820
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = (act("H+")^n)*kHT
23 if SR("Dawsonite") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Dawsonite")))^2)
    24 k_n = 1
25 gamma_n = 2e10
26  lamda = 1
    30  if SR("Dawsonite") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Dawsonite"))*time
    40  if SR("Dawsonite") > 1 then moles = (-reacS*lamda*k_p*(SR("Dawsonite") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end
Magnesite
-start
10 reacS = parm(1)*parm(2)*m
13 n = 1.0
15 logkH = -6.38
    16 Ea = 14400
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = 7.36E-14
23 if SR("Magnesite") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Magnesite")))^2)
    24 k_n = 1
25 gamma_n = 4e10
26  lamda = 1
    30  if SR("Magnesite") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Magnesite"))*time
    40  if SR("Magnesite") > 1 then moles = (-reacS*lamda*k_p*(SR("Magnesite") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end
Muscovite  #Rate parameters from muscovite
-start
10 reacS = parm(1)*parm(2)*m
13 n = 0.37
15 logkH = -11.85
    16 Ea = 22000
    17 R = 8.3145
    20 logkHT = logkH - (Ea/(2.303*R))*(1/Tk - 1/298.15)
    21 kHT = 10^logkH
22 k_p = (act("H+")^n)*kHT/1
23 if SR("Muscovite") > 1 then expTerm = 1/((Tk^3)*(LOG(SR("Muscovite")))^2)
    24 k_n = 1
25 gamma_n = 2e10
26  lamda = 1
    30  if SR("Muscovite") <= 1 then moles = reacS*kHT*(act("H+")^n)*(1 - SR("Muscovite"))*time
    40  if SR("Muscovite") > 1 then moles = (-reacS*lamda*k_p*(SR("Muscovite") - 1)^2 - k_n*exp(-gamma_n*expTerm))*time
    100 save moles
-end

KNOBS
       -iterations             150
       -convergence_tolerance  1e-8
       -tolerance              1e-10
      -step_size              10
     -pe_step_size           5
    -diagonal_scale         TRUE
    -debug_diffuse_layer    TRUE
      -debug_inverse          TRUE
      -debug_model            TRUE
      -debug_prep             TRUE
      -debug_set              TRUE
      -logfile                TRUE
SELECTED_OUTPUT
-file ch4-impure-100y.sel
-pe true
-totals Br Ca Mg Na K Si C(4) Al Fe O(0) Fe(3) Fe(2) S(6)
  -si O2(g) CO2(g) Quartz Kaolinite  Albite  Dawsonite Magnesite  calcite
-kinetic_reactants Quartz Kaolinite  Albite  Dawsonite Magnesite  calcite
#-equilibrium_phases dawsonite co2(g) ch4(g)
-gases co2(g) no2(g)
END     
     
     
     

[dlparkhurst]

Looks like the CH4(g) reacted to me. After one step, some methane has dissolved and sulfate has been reduced to sulfide.

Most of your reactions are not redox sensitive, so there is not a large effect on the calculation.

If you add SOx or NOx, then either CH4(g) will be oxidized or SOx and NOx will be reduced, depending on the relative amounts of the gases. I don't think they will coexist at equilibrium.

Note you are using llnl.dat, which only allows for ideal gases. Phreeqc.dat, Amm.dat, and pitzer.dat has the Peng-Robinson equation of state for gases, which may start to be important at 88 atm.