Kinetic modeling of water-CO2-rock at 300 C and 222.67 atm

[gagasutomo]

Good day everyone

I am currently trying to simulate interaction between water-CO2-rock at high temperature and pressure (300 C & 222.67 atm) using KINETICS data block. My PHREEQC code doesn't show result even after 4 hours running. But, when I take out the gas phase (only simulating water-rock interaction), the code show results. Here are my code, Please let me know the problem.

Thank you

DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.3.7-11094\database\llnl.dat
SOLUTION_SPECIES
CO3-2 + 2 H+ = CO2 + H2O
   -log_k   16.681
   -delta_h -5.738   kcal
   -analytic   464.1965   0.09344813   -26986.16   -165.75951   2248628.9
   -dw    1.92e-9
   -Vm   7.29  0.92  2.07  -1.23  -1.60 # ref. 1 + McBride et al. 2015, JCED 60, 171
2CO2 = (CO2)2 # activity correction for CO2 solubility at high P, T
   -log_k -1.8
   -analytical_expression  8.68  -0.0103  -2190
   -Vm   14.58  1.84  4.14  -2.46  -3.20
END
PHASES
CO2(g)
   CO2 = CO2
   -log_k   -1.468
   -delta_h -4.776 kcal
   -analytic   10.5624  -2.3547e-2  -3972.8  0  5.8746e5  1.9194e-5
   -T_c  304.2 # critical T, K
   -P_c   72.86 # critical P, atm
   -Omega 0.225 # acentric factor
END
SOLUTION 1
    temp      300
    pH        5.9
    pe        4
    redox     pe
    units     ppm
    density   1
    Ca        405
    Cl        6741
    K         919
    Mg        67
    Na        3351
    S(6)      288
    Si        70
    C(4)      366
    -water    1 # kg
GAS_PHASE 1
    -fixed_pressure
    -pressure 222.67
    -volume 1
    -temperature 300
    CO2(g)    222.67
RATES
Quartz
   -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_qtz = SI("Quartz")
    30 if (M <= 0 and si_qtz < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_qtz > 0) then SA = 1e-05  #nucleation
    60 k_acid = 0                 
    70 k_neut = 10^(-13.40)*EXP(-90.90e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_qtz))
    190 moles = r * TIME
    200 SAVE moles
     -end
Anorthite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_anort = SI("Anorthite")
   30 if (M <= 0 and si_anort < 0) then goto 200             
   40 SA = PARM(1) * M             
   50 if (M = 0 and si_anort > 0) then SA = 1e-05  #nucleation
   60 k_acid = 10^(-3.5)*EXP(-16.60e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(1.411)                 
   70 k_neut = 10^(-9.82)*EXP(-31.50e+03/8.314*(1.0/TK-1.0/298.15))
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(10^si_anort))
   190 moles = r * TIME
   200 SAVE moles
   -end
Albite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_alb = SI("Albite")
    30 if (M <= 0 and si_alb < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_alb > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-10.16)*EXP(-65.00e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("H+")^(0.457))                 
    70 k_neut = 10^(-12.56)*EXP(-69.80e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 10^(-15.60)*EXP(-71.00e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("OH-")^(-0.572))
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_alb))
    190 moles = r * TIME
    200 SAVE moles
    -end
K-Feldspar
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_kfeld = SI("K-Feldspar")
    30 if (M <= 0 and si_kfeld < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_kfeld > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-10.06)*EXP(-51.70e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)                 
    70 k_neut = 10^(-12.41)*EXP(-38.00e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 10^(-21.20)*EXP(-94.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.823)
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_kfeld))
    190 moles = r * TIME
    200 SAVE moles
    -end
Diopside
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_diopsi = SI("Diopside")
    30 if (M <= 0 and si_diopsi < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_diopsi > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-6.36)*EXP(-96.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.710)                 
    70 k_neut = 10^(-11.11)*EXP(-40.60e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_diopsi))
    190 moles = r * TIME
    200 SAVE moles
    -end
Enstatite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_enst = SI("Enstatite")
   30 if (M <= 0 and si_enst < 0) then goto 200             
   40 SA = PARM(1) * M             
   50 if (M = 0 and si_enst > 0) then SA = 1e-05  #nucleation
   60 k_acid = 10^(-9.02)*EXP(-80.00e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.600)                 
   70 k_neut = 10^(-12.72)*EXP(-80.00e+03/8.314*(1.0/TK-1.0/298.15))
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(10^si_enst))
   190 moles = r * TIME
   200 SAVE moles
   -end
Magnetite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_magnet = SI("Magnetite")
   30 if (M <= 0 and si_magnet < 0) then goto 200             
   40 SA = PARM(1) * M             
   50 if (M = 0 and si_magnet > 0) then SA = 1e-05  #nucleation
   60 k_acid = 10^(-8.59)*EXP(-18.60e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.279)                 
   70 k_neut = 10^(-10.78)*EXP(-18.60e+03/8.314*(1.0/TK-1.0/298.15))
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(10^si_magnet))
   190 moles = r * TIME
   200 SAVE moles
   -end
Epidote
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_epidt = SI("Epidote")
    30 if (M <= 0 and si_epidt < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_epidt > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-10.60)*EXP(-71.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.338)                 
    70 k_neut = 10^(-11.99)*EXP(-70.70e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 10^(-17.33)*EXP(-79.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.556)
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_epidt))
    190 moles = r * TIME
    200 SAVE moles
    -end
Calcite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_calc = SI("Calcite")
    30 if (M <= 0 and si_calc < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_calc > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-0.30)*EXP(-14.40e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(1.000)                 
    70 k_neut = 10^(-5.81)*EXP(-23.50e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_carb = 10^(-3.48)*EXP(-35.40e+03/8.314*(1.0/TK-1.0/298.15))*PR_P("CO2(g)")^(1.000)
    90 k_rateconst = k_acid + k_neut + k_carb
    100 r = k_rateconst * SA * (1-(10^si_calc))
    190 moles = r * TIME
    200 SAVE moles
    -end
Illite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_ill = SI("Illite")
    30 if (M <= 0 and si_ill < 0) then goto 200           
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_ill > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-10.12)*EXP(-58.00e+03/8.314*(1.0/TK-1.0/373.15))*ACT("H+")^(0.550)                 
    70 k_neut = 10^(-12.26)*EXP(-54.00e+03/8.314*(1.0/TK-1.0/373.15))
    80 k_base = 10^(-10.60)*EXP(-77.00e+03/8.314*(1.0/TK-1.0/373.15))*ACT("OH-")^(-0.350)
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_ill))
    190 moles = r * TIME
    200 SAVE moles
    -end
Chamosite-7A
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_cham = SI("Chamosite-7A")
    30 if (M <= 0 and si_cham < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_cham > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-11.11)*EXP(-88.00e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)                 
    70 k_neut = 10^(-12.52)*EXP(-88.00e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_cham))
    190 moles = r * TIME
    200 SAVE moles
    -end
Dolomite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_dolo = SI("Dolomite")
    30 if (M <= 0 and si_dolo < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_dolo > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-3.76)*EXP(-56.70e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)                 
    70 k_neut = 10^(-8.60)*EXP(-95.30e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_carb = 10^(-5.37)*EXP(-45.70e+03/8.314*(1.0/TK-1.0/298.15))*PR_P("CO2(g)")^(0.500)
    90 k_rateconst = k_acid + k_neut + k_carb
    100 r = k_rateconst * SA * (1-(10^si_dolo))
    190 moles = r * TIME
    200 SAVE moles
    -end
Anhydrite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_anhyd = SI("Anhydrite")
   30 if (M <= 0 and si_anhyd < 0) then goto 200             
   40 SA = PARM(1)*M             
   50 if (M = 0 and si_anhyd > 0) then SA = 1e-05  #nucleation
   60 k_acid = 0
   70 k_neut = 10^(-2.79)*EXP(-0.00e+03/8.314*(1.0/TK-1.0/298.15))
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(10^si_anhyd))
   190 moles = r * TIME
   200 SAVE moles
   -end
KINETICS 1
Quartz
    -formula  SiO2  1
    -m        32.25
    -m0       32.25
    -parms    48.064
    -tol      1e-12
Anorthite
    -formula  CaAl2(SiO4)2  1
    -m        21.12
    -m0       21.12
    -parms    5.548
    -tol      1e-12
Albite
    -formula  NaAlSi3O8  1
    -m        27.48
    -m0       27.48
    -parms    5.260
    -tol      1e-12
K-Feldspar
    -formula  KAlSi3O8  1
    -m        10.53
    -m0       10.53
    -parms    86.282
    -tol      1e-12
Diopside
    -formula  CaMgSi2O6  1
    -m        5.66
    -m0       5.66
    -parms    25.122
    -tol      1e-12
Enstatite
    -formula  MgSiO3  1
    -m        17.06
    -m0       17.06
    -parms    933.627
    -tol      1e-12
Magnetite
    -m        2.89
    -m0       2.89
    -parms    23.616
    -tol      1e-12
Epidote
    -m        0.92
    -m0       0.92
    -parms    128.267
    -tol      1e-12
Calcite
    -formula  CaCO3  1
    -m        4.75
    -m0       4.75
    -parms    15.014
    -tol      1e-12
Illite
    -formula  K0.6Mg0.25Al2.3Si3.5O10(OH)2  1
    -m        1.22
    -m0       1.22
    -parms    5956.902
    -tol      1e-12
Chamosite-7A
    -m        7.05
    -m0       7.05
    -parms    107.920
    -tol      1e-12
Dolomite
    -formula  CaMg(CO3)2  1
    -m        0
    -m0       0
    -parms    33.192
    -tol      1e-12
Anhydrite
    -formula  CaSO4  1
    -m        0
    -m0       0
    -parms    16.578
    -tol      1e-12
-steps       315576 3155760 31557600 315576000 3155760000 31557600000 315576000000
-cvode true
-cvode_steps 31557600
-cvode_order 5
INCREMENTAL_REACTIONS True
USER_GRAPH 1
    -headings               Time Quartz Anorthite Albite K-Feldspar Diopside Enstatite Magnetite Epidote Calcite Illite Chamosite-7A Dolomite Anhydrite
    -axis_titles            "Log10 Time" "kin" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X log10(total_time)
20 GRAPH_Y kin("Quartz"), kin("Anorthite"), kin("Albite"), kin("K-Feldspar"), kin("Diopside"), kin("Enstatite"), kin("Magnetite"), kin("Epidote"), kin("Calcite"), kin("Illite"), kin("Chamosite-7A"), kin("Dolomite"), kin("Anhydrite")
  -end
    -active                 true

END

[dlparkhurst]

Try lowering the tolerance and using a sequence of log-increasing time steps. You can check the accuracy by decreasing the tolerance to 1e-9 to see if it makes any difference.

Also, be sure your conceptual model for the gas phase is right. The gas phase is defined as simply 1 L volume at the specified temperature and CO2 partial pressure, and this amount of CO2 partitions between the gas phase and the solution as the kinetic reactions proceed. You may want to consider a specified CO2(g) fugacity with EQUILIBRIUM_PHASES instead.

Code: [Select]

DATABASE C:\Program Files (x86)\USGS\Phreeqc Interactive 3.3.7-11094\database\llnl.dat
SOLUTION_SPECIES
CO3-2 + 2 H+ = CO2 + H2O
   -log_k   16.681
   -delta_h -5.738   kcal
   -analytic   464.1965   0.09344813   -26986.16   -165.75951   2248628.9
   -dw    1.92e-9
   -Vm   7.29  0.92  2.07  -1.23  -1.60 # ref. 1 + McBride et al. 2015, JCED 60, 171
2CO2 = (CO2)2 # activity correction for CO2 solubility at high P, T
   -log_k -1.8
   -analytical_expression  8.68  -0.0103  -2190
   -Vm   14.58  1.84  4.14  -2.46  -3.20
END
PHASES
CO2(g)
   CO2 = CO2
   -log_k   -1.468
   -delta_h -4.776 kcal
   -analytic   10.5624  -2.3547e-2  -3972.8  0  5.8746e5  1.9194e-5
   -T_c  304.2 # critical T, K
   -P_c   72.86 # critical P, atm
   -Omega 0.225 # acentric factor
END
SOLUTION 1
    temp      300
    pH        5.9
    pe        4
    redox     pe
    units     ppm
    density   1
    Ca        405
    Cl        6741
    K         919
    Mg        67
    Na        3351
    S(6)      288
    Si        70
    C(4)      366
    -water    1 # kg
GAS_PHASE 1
    -fixed_pressure
    -pressure 222.67
    -volume 1
    -temperature 300
    CO2(g)    222.67
RATES
Quartz
   -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_qtz = SI("Quartz")
    30 if (M <= 0 and si_qtz < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_qtz > 0) then SA = 1e-05  #nucleation
    60 k_acid = 0                 
    70 k_neut = 10^(-13.40)*EXP(-90.90e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_qtz))
    190 moles = r * TIME
    200 SAVE moles
     -end
Anorthite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_anort = SI("Anorthite")
   30 if (M <= 0 and si_anort < 0) then goto 200             
   40 SA = PARM(1) * M             
   50 if (M = 0 and si_anort > 0) then SA = 1e-05  #nucleation
   60 k_acid = 10^(-3.5)*EXP(-16.60e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(1.411)                 
   70 k_neut = 10^(-9.82)*EXP(-31.50e+03/8.314*(1.0/TK-1.0/298.15))
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(10^si_anort))
   190 moles = r * TIME
   200 SAVE moles
   -end
Albite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_alb = SI("Albite")
    30 if (M <= 0 and si_alb < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_alb > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-10.16)*EXP(-65.00e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("H+")^(0.457))                 
    70 k_neut = 10^(-12.56)*EXP(-69.80e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 10^(-15.60)*EXP(-71.00e+03/8.314*(1.0/TK-1.0/298.15))*(ACT("OH-")^(-0.572))
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_alb))
    190 moles = r * TIME
    200 SAVE moles
    -end
K-Feldspar
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_kfeld = SI("K-Feldspar")
    30 if (M <= 0 and si_kfeld < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_kfeld > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-10.06)*EXP(-51.70e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)                 
    70 k_neut = 10^(-12.41)*EXP(-38.00e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 10^(-21.20)*EXP(-94.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.823)
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_kfeld))
    190 moles = r * TIME
    200 SAVE moles
    -end
Diopside
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_diopsi = SI("Diopside")
    30 if (M <= 0 and si_diopsi < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_diopsi > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-6.36)*EXP(-96.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.710)                 
    70 k_neut = 10^(-11.11)*EXP(-40.60e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_diopsi))
    190 moles = r * TIME
    200 SAVE moles
    -end
Enstatite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_enst = SI("Enstatite")
   30 if (M <= 0 and si_enst < 0) then goto 200             
   40 SA = PARM(1) * M             
   50 if (M = 0 and si_enst > 0) then SA = 1e-05  #nucleation
   60 k_acid = 10^(-9.02)*EXP(-80.00e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.600)                 
   70 k_neut = 10^(-12.72)*EXP(-80.00e+03/8.314*(1.0/TK-1.0/298.15))
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(10^si_enst))
   190 moles = r * TIME
   200 SAVE moles
   -end
Magnetite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_magnet = SI("Magnetite")
   30 if (M <= 0 and si_magnet < 0) then goto 200             
   40 SA = PARM(1) * M             
   50 if (M = 0 and si_magnet > 0) then SA = 1e-05  #nucleation
   60 k_acid = 10^(-8.59)*EXP(-18.60e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.279)                 
   70 k_neut = 10^(-10.78)*EXP(-18.60e+03/8.314*(1.0/TK-1.0/298.15))
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(10^si_magnet))
   190 moles = r * TIME
   200 SAVE moles
   -end
Epidote
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_epidt = SI("Epidote")
    30 if (M <= 0 and si_epidt < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_epidt > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-10.60)*EXP(-71.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.338)                 
    70 k_neut = 10^(-11.99)*EXP(-70.70e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 10^(-17.33)*EXP(-79.10e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.556)
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_epidt))
    190 moles = r * TIME
    200 SAVE moles
    -end
Calcite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_calc = SI("Calcite")
    30 if (M <= 0 and si_calc < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_calc > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-0.30)*EXP(-14.40e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(1.000)                 
    70 k_neut = 10^(-5.81)*EXP(-23.50e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_carb = 10^(-3.48)*EXP(-35.40e+03/8.314*(1.0/TK-1.0/298.15))*PR_P("CO2(g)")^(1.000)
    90 k_rateconst = k_acid + k_neut + k_carb
    100 r = k_rateconst * SA * (1-(10^si_calc))
    190 moles = r * TIME
    200 SAVE moles
    -end
Illite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_ill = SI("Illite")
    30 if (M <= 0 and si_ill < 0) then goto 200           
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_ill > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-10.12)*EXP(-58.00e+03/8.314*(1.0/TK-1.0/373.15))*ACT("H+")^(0.550)                 
    70 k_neut = 10^(-12.26)*EXP(-54.00e+03/8.314*(1.0/TK-1.0/373.15))
    80 k_base = 10^(-10.60)*EXP(-77.00e+03/8.314*(1.0/TK-1.0/373.15))*ACT("OH-")^(-0.350)
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_ill))
    190 moles = r * TIME
    200 SAVE moles
    -end
Chamosite-7A
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_cham = SI("Chamosite-7A")
    30 if (M <= 0 and si_cham < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_cham > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-11.11)*EXP(-88.00e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)                 
    70 k_neut = 10^(-12.52)*EXP(-88.00e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_base = 0
    90 k_rateconst = k_acid + k_neut + k_base
    100 r = k_rateconst * SA * (1-(10^si_cham))
    190 moles = r * TIME
    200 SAVE moles
    -end
Dolomite
    -start
    10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
    20 si_dolo = SI("Dolomite")
    30 if (M <= 0 and si_dolo < 0) then goto 200             
    40 SA = PARM(1) * M             
    50 if (M = 0 and si_dolo > 0) then SA = 1e-05  #nucleation
    60 k_acid = 10^(-3.76)*EXP(-56.70e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.500)                 
    70 k_neut = 10^(-8.60)*EXP(-95.30e+03/8.314*(1.0/TK-1.0/298.15))
    80 k_carb = 10^(-5.37)*EXP(-45.70e+03/8.314*(1.0/TK-1.0/298.15))*PR_P("CO2(g)")^(0.500)
    90 k_rateconst = k_acid + k_neut + k_carb
    100 r = k_rateconst * SA * (1-(10^si_dolo))
    190 moles = r * TIME
    200 SAVE moles
    -end
Anhydrite
   -start
   10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
   20 si_anhyd = SI("Anhydrite")
   30 if (M <= 0 and si_anhyd < 0) then goto 200             
   40 SA = PARM(1)*M             
   50 if (M = 0 and si_anhyd > 0) then SA = 1e-05  #nucleation
   60 k_acid = 0
   70 k_neut = 10^(-2.79)*EXP(-0.00e+03/8.314*(1.0/TK-1.0/298.15))
   80 k_base = 0
   90 k_rateconst = k_acid + k_neut + k_base
   100 r = k_rateconst * SA * (1-(10^si_anhyd))
   190 moles = r * TIME
   200 SAVE moles
   -end
KINETICS 1
Quartz
    -formula  SiO2  1
    -m        32.25
    -m0       32.25
    -parms    48.064
    #-tol      1e-12
Anorthite
    -formula  CaAl2(SiO4)2  1
    -m        21.12
    -m0       21.12
    -parms    5.548
    #-tol      1e-12
Albite
    -formula  NaAlSi3O8  1
    -m        27.48
    -m0       27.48
    -parms    5.260
    #-tol      1e-12
K-Feldspar
    -formula  KAlSi3O8  1
    -m        10.53
    -m0       10.53
    -parms    86.282
    #-tol      1e-12
Diopside
    -formula  CaMgSi2O6  1
    -m        5.66
    -m0       5.66
    -parms    25.122
    #-tol      1e-12
Enstatite
    -formula  MgSiO3  1
    -m        17.06
    -m0       17.06
    -parms    933.627
    #-tol      1e-12
Magnetite
    -m        2.89
    -m0       2.89
    -parms    23.616
    #-tol      1e-12
Epidote
    -m        0.92
    -m0       0.92
    -parms    128.267
    #-tol      1e-12
Calcite
    -formula  CaCO3  1
    -m        4.75
    -m0       4.75
    -parms    15.014
    #-tol      1e-12
Illite
    -formula  K0.6Mg0.25Al2.3Si3.5O10(OH)2  1
    -m        1.22
    -m0       1.22
    -parms    5956.902
    #-tol      1e-12
Chamosite-7A
    -m        7.05
    -m0       7.05
    -parms    107.920
    #-tol      1e-12
Dolomite
    -formula  CaMg(CO3)2  1
    -m        0
    -m0       0
    -parms    33.192
    #-tol      1e-12
Anhydrite
    -formula  CaSO4  1
    -m        0
    -m0       0
    -parms    16.578
    #-tol      1e-12
-steps 0.001 0.003 0.01 0.03 0.1 0.3 1 3 10 30 100 300 1000 3000 10000 30000 100000 300000 1e6 3e6 1e7 3e7 1e8 3e8 1e9 3e9 1e10 3e10 1e11 3e11
#-steps 315576 3155760 31557600 315576000 3155760000 31557600000 315576000000
-cvode true
#-cvode_steps 31557600
#-cvode_order 5
KNOBS
-conv 1e-10
INCREMENTAL_REACTIONS True
USER_GRAPH 1
    -headings               Time Quartz Anorthite Albite K-Feldspar Diopside Enstatite Magnetite Epidote Calcite Illite Chamosite-7A Dolomite Anhydrite
    -axis_titles            "Log10 Time" "kin" ""
    -initial_solutions      false
    -connect_simulations    true
    -plot_concentration_vs  x
  -start
10 GRAPH_X log10(total_time)
20 GRAPH_Y kin("Quartz"), kin("Anorthite"), kin("Albite"), kin("K-Feldspar"), kin("Diopside"), kin("Enstatite"), kin("Magnetite"), kin("Epidote"), kin("Calcite"), kin("Illite"), kin("Chamosite-7A"), kin("Dolomite"), kin("Anhydrite")
  -end
    -active                 true

END


[gagasutomo]

Thank you for your advice Mr. Parkhurst.

I tried lowering the tolerance level to 1e-7 and it worked. Although it is higher than the tolerance level you recommended (1e-9).

Also, I didn't understand the difference between specifying CO2 gas phase in GAS_PHASE and EQUILIBRIUM_PHASE, could you please en light me. Additionally, how can I specify CO2(g) fugacity in The EQUILIBRIUM_PHASE?

[dlparkhurst]

When you define a GAS_PHASE, the initial number of moles of each component is determined by the specified partial pressures, and the temperature, pressure, and volume assigned to the GAS_PHASE. If CO2(g) is assigned a pressure of 10 atm. Let's assume an ideal gas, so the fugacity is equal to the partial pressure. When the GAS_PHASE reacts with a solution, CO2(g) will partition between the solution and the solution, and in general, the partial pressure of CO2(g) in the gas phase will change and be equal to the calculated partial pressure in the solution.

If you specify a partial pressure in EQUILIBRIUM_PHASES (target saturation index is equal to the log fugacity that is desired), then when the equilibrium phase reacts with a solution, the solution will attain the specified partial pressure.

So, with a gas phase, the partial pressure will vary, whereas with an equilibrium phase the partial pressure will remain constant (as long as moles of the component remain in the equilibrium phase).

[gagasutomo]

Thank you again for your enlightenment Mr. Parkhurst

I also tried equilibrium modeling between water-CO2-rock at 300 C and 222.67 atm by specifying:
 1. SOLUTION (for water, with specified T)
 2. GAS_PHASE (for CO2, with specified P, T and V)
 3. EQUILIBRIUM_PHASE (for rock minerals)

What would be the difference between specifying the CO2 in the GAS_PHASE and in the EQUILIBRIUM_PHASE in this case? In kinetic modeling, minerals are assigned separately in the KINETICS command. While in equilibrium modeling, it is specified in EQUILIBRIUM_PHASE. How can I assigned the desired log fugacity when other minerals is also in EQUILIBRIUM_PHASE? What is your suggestion?

[dlparkhurst]

All of my previous comments apply. With GAS_PHASE, the partial pressure will vary with reaction (and, with a fixed pressure gas phase, the gas phase could disappear.)  With EQUILIBRIUM_PHASES, the partial pressure (fugacity) will be fixed.

You can have CO2(g) in EQUILIBRIUM_PHASES in addition to other minerals.