How to determine kinetic time step

[Sidol]

Dear Dr. Parkhurst!

I am trying to modeling soil pore water composition by using distilled water and soil mineral phases in a batch reaction.
In order to avoid the unrealistic dissolution of certain mineral phases (e.g. quartz, feldspars, calcite) I handled them as kinetic reactants.
Here is my code:

Code: [Select]

TITLE A-szint pórusvíz modellje, kinetikus
      "PHASES" blokk:  PHREEQC_ThermoddemV1.10_06Jun2017
SOLUTION 1
    temp      14.1
    pH        7
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    Al        1e-13
    Ca        1e-13
    Cl        1e-13
    F         1e-13
    Fe        1e-13
    K         1e-13
    Mg        1e-13
    N(5)      1e-13
    Na        1e-13
    S         1e-13
    Si        1e-13
    -water    1 # kg


PHASES
Microcline
    K(AlSi3)O8 + 4H+ + 4H2O = Al+3 + 3H4SiO4 + K+
    log_k     0.004
    delta_h   -56.203 kJ
    -analytical_expression -736.77713 -0.12898219 36861.528 270.3714 -1547971 0
VermiculiteSO
    Ca0.445(Si2.778Al1.222)(Al0.216Mg2.475Fe0.254)O10(OH)2 + 10.888H+ = 1.438Al+3 + 0.445Ca+2 + 0.028Fe+2 + 0.226Fe+3 + 0.888H2O + 2.778H4SiO4 + 2.475Mg+2
    log_k     45.888
    delta_h   -441.531 kJ
    -analytical_expression -1922.3485 -0.31254347 118646.07 694.16321 -4816349.5 0
Vermiculite(K)
    K0.86Mg3.00Si3.14Al0.86O10(OH)2 + 9.44H+ + 0.56H2O = 0.86Al+3 + 3.14H4SiO4 + 0.86K+ + 3Mg+2
    log_k     37.445
    delta_h   -335.54 kJ
    -analytical_expression -1693.6279 -0.26466982 102324.37 612.44388 -4326093.9 0
Vermiculite(Mg)
    Mg0.43Mg3.00Si3.14Al0.86O10(OH)2 + 9.44H+ + 0.56H2O = 0.86Al+3 + 3.14H4SiO4 + 3.43Mg+2
    log_k     38.042
    delta_h   -379.809 kJ
    -analytical_expression -1782.4468 -0.27812893 108797 642.42701 -4545784.1 0
Vermiculite(Na)
    Na0.86Mg3.00Si3.14Al0.86O10(OH)2 + 9.44H+ + 0.56H2O = 0.86Al+3 + 3.14H4SiO4 + 3Mg+2 + 0.86Na+
    log_k     38.389
    delta_h   -355.542 kJ
    -analytical_expression -1744.5392 -0.26999731 105980.44 629.66633 -4453583.2 0
EQUILIBRIUM_PHASES 1
    CO2(g)    -2.6 10
    Dolomite  0 0.1966
    Goethite  0 0.6528
    Halloysite 0 0
    Illite    0 0.9118
    Kaolinite 0 2.4165
    O2(g)     -0.686 10
    Vermiculite(K) 0 0.202
    Vermiculite(Mg) 0 0.214
    Vermiculite(Na) 0 0.209
    VermiculiteSO 0 0.2067

KINETICS 1
K-feldspar
    -formula  K-feldspar  1
    -m        0.5661
    -m0       0.5661
    -parms    6.41 0.1
    -tol      1e-08
Albite
    -formula  NaAlSi3O8  1
    -m        0.599
    -m0       0.599
    -parms    6.04 0.1
    -tol      1e-08
Quartz
    -formula  SiO2  1
    -m        20.331
    -m0       20.331
    -parms    0.146 1.5
    -tol      1e-08
Calcite
    -formula  CaCO3  1
    -m        0.3622
    -m0       0.3622
    -parms    167000 0.6
    -tol      1e-08
Pyrite
    -formula  FeS2  1
    -m        0.0604
    -m0       0.0604
    -parms    0.3 0.67 0.5 -0.11
    -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

INCREMENTAL_REACTIONS True
RATES
    K-feldspar
-start
  1   REM Sverdrup and Warfvinge, 1995, mol m^-2 s^-1
  2   REM PARM(1) = Specific area of Kspar m^2/mol Kspar
  3   REM PARM(2) = Adjusts lab rate to field rate
  4   REM temp corr: from A&P, p. 162. E (kJ/mol) / R / 2.303 = H in H*(1/T-1/281)
  5   REM K-Feldspar parameters
 10  DATA 11.7, 0.5, 4e-6, 0.4, 500e-6, 0.15, 14.5, 0.14, 0.15, 13.1, 0.3
 20  RESTORE 10
 30  READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH
 40  DATA 3500, 2000, 2500, 2000
 50  RESTORE 40
 60  READ e_H, e_H2O, e_OH, e_CO2
 70  pk_CO2 = 13
 80  n_CO2 = 0.6
100 REM Generic rate follows
110 dif_temp = 1/TK - 1/281
120 BC       = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2")
130 REM rate by H+
140 pk_H     = pk_H + e_H * dif_temp
150 rate_H   = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC)
160 REM rate by hydrolysis
170 pk_H2O   = pk_H2O + e_H2O * dif_temp
180 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC)
190 REM rate by OH-
200 pk_OH    = pk_OH + e_OH * dif_temp
210 rate_OH  = 10^-pk_OH * ACT("OH-")^o_OH
220 REM rate by CO2
230 pk_CO2   = pk_CO2 + e_CO2 * dif_temp
240 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2
250 rate     = rate_H + rate_H2O + rate_OH + rate_CO2
260 area     = PARM(1) * M0 *(M/M0)^0.67
270 rate     = PARM(2) * area * rate * (1-SR("Microcline"))
280 moles    = rate * TIME
290 SAVE moles
-end
    Albite
-start
  1   REM Sverdrup and Warfvinge, 1995, mol m^-2 s^-1
  2   REM PARM(1) = Specific area of Albite m^2/mol Albite
  3   REM PARM(2) = Adjusts lab rate to field rate
  4   REM temp corr: from A&P, p. 162. E (kJ/mol) / R / 2.303 = H in H*(1/T-1/281)
  5   REM Albite parameters
 10  DATA 11.5, 0.5, 4e-6, 0.4, 500e-6, 0.2, 13.7, 0.14, 0.15, 11.8, 0.3
 20  RESTORE 10
 30  READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH
 40  DATA 3500, 2000, 2500, 2000
 50  RESTORE 40
 60  READ e_H, e_H2O, e_OH, e_CO2
 70  pk_CO2 = 13
 80  n_CO2 = 0.6
100 REM Generic rate follows
110 dif_temp = 1/TK - 1/281
120 BC       = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2")
130 REM rate by H+
140 pk_H     = pk_H + e_H * dif_temp
150 rate_H   = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC)
160 REM rate by hydrolysis
170 pk_H2O   = pk_H2O + e_H2O * dif_temp
180 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC)
190 REM rate by OH-
200 pk_OH    = pk_OH + e_OH * dif_temp
210 rate_OH  = 10^-pk_OH * ACT("OH-")^o_OH
220 REM rate by CO2
230 pk_CO2   = pk_CO2 + e_CO2 * dif_temp
240 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2
250 rate     = rate_H + rate_H2O + rate_OH + rate_CO2
260 area     = PARM(1) * M0 *(M/M0)^0.67
270 rate     = PARM(2) * area * rate * (1-SR("Albite"))
280 moles    = rate * TIME
290 SAVE moles
-end
    Quartz
-start
 1  REM  Specific rate k from Rimstidt and Barnes, 1980, GCA 44,1683
 2  REM  k = 10^-13.7 mol/m2/s (25 C), Ea = 90 kJ/mol
 3  REM  sp. rate * parm(2) due to salts (Dove and Rimstidt, MSA Rev. 29, 259)
 4  REM  PARM(1) = Specific area of Quartz, m^2/mol Quartz
 5  REM  PARM(2) = salt correction: (1 + 1.5 * c_Na (mM)), < 35
10 dif_temp = 1/TK - 1/298
20 pk_w = 13.7 + 4700.4 * dif_temp
40 moles = PARM(1) * M0 * PARM(2) * (M/M0)^0.67 * 10^-pk_w * (1 -  SR("Quartz"))
50 SAVE moles * TIME
-end
    Calcite
-start
  1   REM   PARM(1) = specific surface area of calcite, cm^2/mol calcite
  2   REM   PARM(2) = exponent for M/M0
 10  si_cc = SI("Calcite")
 20  IF (M <= 0  and si_cc < 0) THEN GOTO 200
 30  k1 = 10^(0.198 - 444.0 / TK )
 40  k2 = 10^(2.84 - 2177.0 /TK )
 50  IF TC <= 25 THEN k3 = 10^(-5.86 - 317.0 / TK)
 60  IF TC > 25 THEN k3 = 10^(-1.1 - 1737.0 / TK )
 80  IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M
110 rate = area * (k1 * ACT("H+") + k2 * ACT("CO2") + k3 * ACT("H2O"))
120 rate = rate * (1 - 10^(2/3*si_cc))
130 moles = rate * 0.001 * TIME
200 SAVE moles
-end
    Pyrite
-start
  1   REM        Williamson and Rimstidt, 1994
  2   REM        PARM(1) = log10(specific area), log10(m^2 per mole pyrite)
  3   REM        PARM(2) = exp for (M/M0)
  4   REM        PARM(3) = exp for O2
  5   REM        PARM(4) = exp for H+
 10  REM Dissolution in presence of DO
 20  if (M <= 0) THEN GOTO 200
 30  if (SI("Pyrite") >= 0) THEN GOTO 200
 40  log_rate = -8.19 + PARM(3)*LM("O2") + PARM(4)*LM("H+")
 50  log_area = PARM(1) + LOG10(M0) + PARM(2)*LOG10(M/M0)
 60  moles = 10^(log_area + log_rate) * TIME
200 SAVE moles
-end

END
TITLE B-szint pórusvíz modellje, kinetikus
      "PHASES" blokk:  PHREEQC_ThermoddemV1.10_06Jun2017
SOLUTION 2
    temp      14.1
    pH        7
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    Al        1e-13
    Ca        1e-13
    Cl        1e-13
    F         1e-13
    Fe        1e-13
    K         1e-13
    Mg        1e-13
    N(5)      1e-13
    Na        1e-13
    S         1e-13
    Si        1e-13
    -water    1 # kg
PHASES
Microcline
    K(AlSi3)O8 + 4H+ + 4H2O = Al+3 + 3H4SiO4 + K+
    log_k     0.004
    delta_h   -56.203 kJ
    -analytical_expression -736.77713 -0.12898219 36861.528 270.3714 -1547971 0
VermiculiteSO
    Ca0.445(Si2.778Al1.222)(Al0.216Mg2.475Fe0.254)O10(OH)2 + 10.888H+ = 1.438Al+3 + 0.445Ca+2 + 0.028Fe+2 + 0.226Fe+3 + 0.888H2O + 2.778H4SiO4 + 2.475Mg+2
    log_k     45.888
    delta_h   -441.531 kJ
    -analytical_expression -1922.3485 -0.31254347 118646.07 694.16321 -4816349.5 0
Vermiculite(K)
    K0.86Mg3.00Si3.14Al0.86O10(OH)2 + 9.44H+ + 0.56H2O = 0.86Al+3 + 3.14H4SiO4 + 0.86K+ + 3Mg+2
    log_k     37.445
    delta_h   -335.54 kJ
    -analytical_expression -1693.6279 -0.26466982 102324.37 612.44388 -4326093.9 0
Vermiculite(Mg)
    Mg0.43Mg3.00Si3.14Al0.86O10(OH)2 + 9.44H+ + 0.56H2O = 0.86Al+3 + 3.14H4SiO4 + 3.43Mg+2
    log_k     38.042
    delta_h   -379.809 kJ
    -analytical_expression -1782.4468 -0.27812893 108797 642.42701 -4545784.1 0
Vermiculite(Na)
    Na0.86Mg3.00Si3.14Al0.86O10(OH)2 + 9.44H+ + 0.56H2O = 0.86Al+3 + 3.14H4SiO4 + 3Mg+2 + 0.86Na+
    log_k     38.389
    delta_h   -355.542 kJ
    -analytical_expression -1744.5392 -0.26999731 105980.44 629.66633 -4453583.2 0
EQUILIBRIUM_PHASES 2
    CO2(g)    -2.6 10
    Dolomite  0 0.2159
    Goethite  0 1.0755
    Halloysite 0 2.2276
    Illite    0 0.9141
    Kaolinite 0 0
    O2(g)     -0.686 10
    Vermiculite(K) 0 0.2904
    Vermiculite(Mg) 0 0.3077
    Vermiculite(Na) 0 0.3005
    VermiculiteSO 0 0.2972
INCREMENTAL_REACTIONS True
KINETICS 2
K-feldspar
    -formula  K-feldspar  1
    -m        0.5661
    -m0       0.5661
    -parms    6.41 0.1
    -tol      1e-08
Albite
    -formula  NaAlSi3O8  1
    -m        0.599
    -m0       0.599
    -parms    6.04 0.1
    -tol      1e-08
Quartz
    -formula  SiO2  1
    -m        20.331
    -m0       20.331
    -parms    0.146 1.5
    -tol      1e-08
Calcite
    -formula  CaCO3  1
    -m        0.3622
    -m0       0.3622
    -parms    167000 0.6
    -tol      1e-08
Pyrite
    -formula  FeS2  1
    -m        0.0604
    -m0       0.0604
    -parms    0.3 0.67 0.5 -0.11
    -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
RATES
    K-feldspar
-start
  1   REM Sverdrup and Warfvinge, 1995, mol m^-2 s^-1
  2   REM PARM(1) = Specific area of Kspar m^2/mol Kspar
  3   REM PARM(2) = Adjusts lab rate to field rate
  4   REM temp corr: from A&P, p. 162. E (kJ/mol) / R / 2.303 = H in H*(1/T-1/281)
  5   REM K-Feldspar parameters
 10  DATA 11.7, 0.5, 4e-6, 0.4, 500e-6, 0.15, 14.5, 0.14, 0.15, 13.1, 0.3
 20  RESTORE 10
 30  READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH
 40  DATA 3500, 2000, 2500, 2000
 50  RESTORE 40
 60  READ e_H, e_H2O, e_OH, e_CO2
 70  pk_CO2 = 13
 80  n_CO2 = 0.6
100 REM Generic rate follows
110 dif_temp = 1/TK - 1/281
120 BC       = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2")
130 REM rate by H+
140 pk_H     = pk_H + e_H * dif_temp
150 rate_H   = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC)
160 REM rate by hydrolysis
170 pk_H2O   = pk_H2O + e_H2O * dif_temp
180 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC)
190 REM rate by OH-
200 pk_OH    = pk_OH + e_OH * dif_temp
210 rate_OH  = 10^-pk_OH * ACT("OH-")^o_OH
220 REM rate by CO2
230 pk_CO2   = pk_CO2 + e_CO2 * dif_temp
240 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2
250 rate     = rate_H + rate_H2O + rate_OH + rate_CO2
260 area     = PARM(1) * M0 *(M/M0)^0.67
270 rate     = PARM(2) * area * rate * (1-SR("Microcline"))
280 moles    = rate * TIME
290 SAVE moles
-end
    Albite
-start
  1   REM Sverdrup and Warfvinge, 1995, mol m^-2 s^-1
  2   REM PARM(1) = Specific area of Albite m^2/mol Albite
  3   REM PARM(2) = Adjusts lab rate to field rate
  4   REM temp corr: from A&P, p. 162. E (kJ/mol) / R / 2.303 = H in H*(1/T-1/281)
  5   REM Albite parameters
 10  DATA 11.5, 0.5, 4e-6, 0.4, 500e-6, 0.2, 13.7, 0.14, 0.15, 11.8, 0.3
 20  RESTORE 10
 30  READ pK_H, n_H, lim_Al, x_Al, lim_BC, x_BC, pK_H2O, z_Al, z_BC, pK_OH, o_OH
 40  DATA 3500, 2000, 2500, 2000
 50  RESTORE 40
 60  READ e_H, e_H2O, e_OH, e_CO2
 70  pk_CO2 = 13
 80  n_CO2 = 0.6
100 REM Generic rate follows
110 dif_temp = 1/TK - 1/281
120 BC       = ACT("Na+") + ACT("K+") + ACT("Mg+2") + ACT("Ca+2")
130 REM rate by H+
140 pk_H     = pk_H + e_H * dif_temp
150 rate_H   = 10^-pk_H * ACT("H+")^n_H / ((1 + ACT("Al+3") / lim_Al)^x_Al * (1 + BC / lim_BC)^x_BC)
160 REM rate by hydrolysis
170 pk_H2O   = pk_H2O + e_H2O * dif_temp
180 rate_H2O = 10^-pk_H2O / ((1 + ACT("Al+3") / lim_Al)^z_Al * (1 + BC / lim_BC)^z_BC)
190 REM rate by OH-
200 pk_OH    = pk_OH + e_OH * dif_temp
210 rate_OH  = 10^-pk_OH * ACT("OH-")^o_OH
220 REM rate by CO2
230 pk_CO2   = pk_CO2 + e_CO2 * dif_temp
240 rate_CO2 = 10^-pk_CO2 * (SR("CO2(g)"))^n_CO2
250 rate     = rate_H + rate_H2O + rate_OH + rate_CO2
260 area     = PARM(1) * M0 *(M/M0)^0.67
270 rate     = PARM(2) * area * rate * (1-SR("Albite"))
280 moles    = rate * TIME
290 SAVE moles
-end
    Quartz
-start
 1  REM  Specific rate k from Rimstidt and Barnes, 1980, GCA 44,1683
 2  REM  k = 10^-13.7 mol/m2/s (25 C), Ea = 90 kJ/mol
 3  REM  sp. rate * parm(2) due to salts (Dove and Rimstidt, MSA Rev. 29, 259)
 4  REM  PARM(1) = Specific area of Quartz, m^2/mol Quartz
 5  REM  PARM(2) = salt correction: (1 + 1.5 * c_Na (mM)), < 35
10 dif_temp = 1/TK - 1/298
20 pk_w = 13.7 + 4700.4 * dif_temp
40 moles = PARM(1) * M0 * PARM(2) * (M/M0)^0.67 * 10^-pk_w * (1 -  SR("Quartz"))
50 SAVE moles * TIME
-end
    Calcite
-start
  1   REM   PARM(1) = specific surface area of calcite, cm^2/mol calcite
  2   REM   PARM(2) = exponent for M/M0
 10  si_cc = SI("Calcite")
 20  IF (M <= 0  and si_cc < 0) THEN GOTO 200
 30  k1 = 10^(0.198 - 444.0 / TK )
 40  k2 = 10^(2.84 - 2177.0 /TK )
 50  IF TC <= 25 THEN k3 = 10^(-5.86 - 317.0 / TK)
 60  IF TC > 25 THEN k3 = 10^(-1.1 - 1737.0 / TK )
 80  IF M0 > 0 THEN area = PARM(1)*M0*(M/M0)^PARM(2) ELSE area = PARM(1)*M
110 rate = area * (k1 * ACT("H+") + k2 * ACT("CO2") + k3 * ACT("H2O"))
120 rate = rate * (1 - 10^(2/3*si_cc))
130 moles = rate * 0.001 * TIME
200 SAVE moles
-end
    Pyrite
-start
  1   REM        Williamson and Rimstidt, 1994
  2   REM        PARM(1) = log10(specific area), log10(m^2 per mole pyrite)
  3   REM        PARM(2) = exp for (M/M0)
  4   REM        PARM(3) = exp for O2
  5   REM        PARM(4) = exp for H+
 10  REM Dissolution in presence of DO
 20  if (M <= 0) THEN GOTO 200
 30  if (SI("Pyrite") >= 0) THEN GOTO 200
 40  log_rate = -8.19 + PARM(3)*LM("O2") + PARM(4)*LM("H+")
 50  log_area = PARM(1) + LOG10(M0) + PARM(2)*LOG10(M/M0)
 60  moles = 10^(log_area + log_rate) * TIME
200 SAVE moles
-end
END

I wold like to ask that what kinetic time step should I choose int his case?
If I use linear steps at KINETIC block (attached KINETIC.jpg) what it is advisable to enter as Total time (seconds) and Steps?

Thank You for Your help in advance.

Krisztian

[dlparkhurst]

I'm not sure I agree with your general approach. I don't know your system, but I think it is difficult to select the minerals that react, much less the rates of reaction. I'm also skeptical about your selection of EQUILIBRIUM_PHASES versus KINETICS. But if you are going this way, you can select whatever time step you want. I'd add USER_GRAPH to plot saturation indices of the KINETIC definitions and simply try longer and longer time steps until you reach equilibrium.

If you are trying to match a water composition, you can try INVERSE_MODELING to look for sets of minerals and gases that account for the evolution of rain, pure water, to the observed water composition.