Kinetic Modeling Error

[serhatttt]

Hi there,

I have trouble with this code. I want to simulate the mineral change in a reservoir with respect to time, but the code gives the following error;

ERROR: CVode is at maximum calls: 500. Cell: 1. Time: 2.58e+10 s
ERROR: Please increase the maximum calls with -bad_step_max.
Stopping.

Here is the code;


DATABASE c:\phreeqc\database\llnl.dat

SOLUTION_SPECIES
CO3-2 + 2 H+ = CO2 + H2O
     log_k 16.681
     delta_h -5.738 kcal
    -analytical_expression 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
     -analytical_expression 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
Ankerite
CaMg0.3Fe0.7(CO3)2 + 2H+ = Ca++ + 0.3Mg++ + 0.7Fe++ + 2HCO3-
     log_k 1.54
     delta_h 0
     -analytical_expression -1.8649e+03 -2.9583e-01 1.0468e+05 6.7554e+02 -6.0514e+06
END

SOLUTION 1 SK6 Well
      temp 200
      pH   6.25
      pe 4
      redox pe
      units mg/l
      density 1
      C(4)     1468
      Ca       48.29
      Cl       1183
      Mg       4.15
      Na       1073
      S(6)     12.85
      K        58.43
      F        3.01
      B        36.31
      Si       145
      Li       4.59
      Fe       0.488
      Sb       0.01242
      water 1 # kg
      pressure 126.18 atm #Reservoir pressure
EQUILIBRIUM_PHASES 1
      CO2(g) 2.10 0.08345 #Partial pressure of carbondioxide at 50 atm, and injected mol
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.99)*EXP(-87.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_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
Corundum
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_cornd = SI("Corundum")
30 if (M <= 0 and si_cornd < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_cornd > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-12.53)*EXP(-47.61e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.150)
70 k_neut = 10^(-15.70)*EXP(-47.61e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-11.75)*EXP(-47.61e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.550)
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_cornd))
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
Ilmenite
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_ilmen = SI("Ilmenite")
30 if (M <= 0 and si_ilmen < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_ilmen > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-8.35)*EXP(-37.90e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.421)
70 k_neut = 10^(-11.16)*EXP(-37.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_ilmen))
190 moles = r * TIME
200 SAVE moles
   -end
Fluorapatite
   -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_fluorap = SI("Fluorapatite")
30 if (M <= 0 and si_fluorap < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_fluorap > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-3.73)*EXP(-25.00e+04/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.613)
70 k_neut = 10^(-8.00)*EXP(-25.00e+04/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_fluorap))
190 moles = r * TIME
200 SAVE moles
   -end
Kaolinite
   -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_kaol = SI("Kaolinite")
30 if (M <= 0 and si_kaol < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_kaol > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-11.31)*EXP(-65.90e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.777)
70 k_neut = 10^(-13.18)*EXP(-22.20e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-17.05)*EXP(-17.90e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.472)
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_kaol))
190 moles = r * TIME
200 SAVE moles
   -end
Alunite
   -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_alun = SI("Alunite")
30 if (M <= 0 and si_alun < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_alun > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-10.65)*EXP(-32.30e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.133)
70 k_neut = 10^(-12.53)*EXP(-39.00e+03/8.314*(1.0/TK-1.0/298.15))
80 k_base = 10^(-12.53)*EXP(-39.00e+03/8.314*(1.0/TK-1.0/298.15))*ACT("OH-")^(-0.194)
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_alun))
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
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^(-3.19)*EXP(-14.30e+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
Gypsum
   -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_gyps = SI("Gypsum")
30 if (M <= 0 and si_gyps < 0) then goto 200
40 SA = PARM(1)*M
50 if (M = 0 and si_gyps > 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_gyps))
190 moles = r * TIME
200 SAVE moles
   -end
Hematite
    -start
10 REM PARM(1) = MSA (Molar surface area) [m^2/mol]
20 si_hema = SI("Hematite")
30 if (M <= 0 and si_hema < 0) then goto 200
40 SA = PARM(1) * M
50 if (M = 0 and si_hema > 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-9.39)*EXP(-66.20e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(1.000)
70 k_neut = 10^(-14.60)*EXP(-66.20e+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_hema))
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
Stibnite
    -start
0 REM PARM(1) = MSA (Molar surface area) [m^2/mol]   
20 si_stib = SI("Stibnite")
30 if (M <= 0 and si_stib < 0) then goto 200 
40 SA = PARM(1) * M
50 if (M = 0 and si_stib> 0) then SA = 1e-05 #nucleation
60 k_acid = 10^(-8.04)*EXP(-50.8e+03/8.314*(1.0/TK-1.0/298.15))*ACT("H+")^(0.355)
70 k_neut = 0
80 k_base = 0
90 k_rateconst = k_acid + k_neut + k_base
100 r = k_rateconst * SA * (1-(10^si_stib))
190 moles = r * TIME
200 SAVE moles
    -end
KINETICS 1
Quartz
    -formula SiO2 1
    -m 207.05
    -m0 207.05
    -parms 1.36
    -tol 1e-08
   
Anorthite
    -formula CaAl2(SiO4)2 1
    -m 12.20
    -m0 12.20
    -parms 6.03
    -tol 1e-08

Albite
    -formula NaAlSi3O8 1
    -m 0.5650
    -m0 0.5650
    -parms 6.02
    -tol 1e-08

K-Feldspar
    -formula KAlSi3O8 1
    -m 17.11
    -m0 17.11
    -parms 6.523
    -tol 1e-08

Corundum
    -formula Al2O3 1
    -m 23.98
    -m0 23.98
    -parms 1.537
    -tol 1e-08
   
Enstatite
    -formula MgSiO3 1
    -m 10
    -m0 10
    -parms 1.6235
    -tol 1e-08

Ilmenite
    -formula FeTiO3 1
    -m 0.1392
    -m0 0.1392
    -parms 1.916
    -tol 1e-08
   
Hematite
     -formula Fe2O3 1
     -m 10
     -m0 10
     -parms 1.825
     -tol 1e-08   
   
Kaolinite
    -formula Al2Si2O5(OH)4 1
    -m 10
    -m0 10
    -parms 5.98
    -tol 1e-08

Calcite
    -formula CaCO3 1
    -m 10
    -m0 10
    -parms 2.2661
    -tol 1e-08
   
Dolomite
     -formula CaMg(CO3)2 1
     -m 10
     -m0 10
     -parms 3.8685
     -tol 1e-08
     
Anhydrite
    -formula CaSO4 1
    -m 10
    -m0 10
    -parms 2.80
    -tol 1e-08
   
Stibnite
    -m 10
    -m0 10
    -parms 4.40
    -tol 1e-08
    -steps 3 30 300 3000 30000 300000 3000000 30000000 300000000 3000000000 30000000000 300000000000
    -cvode true
   
KNOBS
    -convergence_tolerance 1e-10

INCREMENTAL_REACTIONS True

USER_GRAPH 1
    -headings Time Quartz Anorthite Albite K-Feldspar Corundum Kaolinite Calcite Anhydrite Hematite Dolomite Stibnite
    -axis_titles "Log10 Time" "Mineral (mol)"
    -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("Corundum"), kin("Kaolinite"), kin("Calcite"), kin("Anhydrite"), kin("Hematite"), kin("Dolomite"), kin("Stibnite")
    -end
   
-active true

SELECTED_OUTPUT 1
    -file SK6.sel
    -reset false
    -time true
    -step true
    -ph true
    -pe true
    -totals C(4)
    -molalities K+ Na+ Ca+2 Mg+2 HCO3- SO4-2 Cl- CO3-2 CO2 Sb(OH)3 HSb2S4-
    -kinetic_reactants Quartz Anorthite Albite K-feldspar Corundum Kaolinite Calcite Anhydrite Hematite Dolomite Stibnite
END

[dlparkhurst]

PHREEQC can have problems with multiple kinetic reactants and low solubility minerals.

I have not analyzed your system carefully, but my first guess was that you can complete your calculation if you  should move Hematite from KINETICS to EQUILIBRIUM_PHASES. My cursory look indicated that Hematite was always near equilibrium, so the change should not affect the results. By making this change, the calculation ran to completion.

[serhatttt]

Dear Parkhurst

The code works now. Many thanks for your quick response!