KINETICS RATE from PALANDRI AND KHARAKA (2004) and others

[dkonan01]

Hello everyone,

I am pretty new in using PHREEQC and I would like some explanations from someone for improving my understanding. I started learning PHREEQC 3 months ago as part of a CCUS project and I reached a point where I am a bit confused regarding the kinetic rate constant for mineral dissolution and precipitation.
My concern might be basic but I still don't get it until now. In Palandri's paper the rate is expressed in term of log k (25 C), while in other papers the rate constant is just written in term of k 25. I would like to know the difference and how can I used both values to obtain the same result?

 My 2nd concern is that i tried to reproduce a simulation of the kinetic dissolution and precipitation from a published paper but after 3 days it was still running without any result. I concluded that i made some mistakes that i would also to get some explanations if possible. Thank you for your help

Code: [Select]

TITLE Kinetic dissolution and precipitation of Minerals
SOLUTION 1
    temp      25
    pH        7
    pe        4
    redox     pe
    units     mol/kgw
    density   1
    -water    1 # kg
EQUILIBRIUM_PHASES 1
    Albite    0 0.025
    Anhydrite 0 0.007
    Ca-Montmorillonite 0 0.005
    Calcite   0 0.133
    Chlorite(14A) 0 0.001
    Dolomite  0 0.083
    Illite    0 0.021
    K-feldspar 0 0.025
    Kaolinite 0 0.036
    Pyrite    0 0.011
    Quartz    0 0.588
PHASES
Albite
    NaAlSi3O8 + 8H2O = Al(OH)4- + 3H4SiO4 + Na+
    -analytical_expression -11.694 0.014429 13784 -7.2866 -1613600 0
    -Vm       0.00010131 m3/mol
Anhydrite
    CaSO4 = Ca+2 + SO4-2
    -analytical_expression -209.86 -0.078823 5096.9 85.642 79.594 0
    -Vm       4.61e-05 m3/mol
Ca-Montmorillonite
    Ca0.165Al2.33Si3.67O10(OH)2 + 12H2O = 2.33Al(OH)4- + 0.165Ca+2 + 2H+ + 3.67H4SiO4
    -analytical_expression 6.0725 0.010644 16024 -16.334 -1798200 0
    -Vm       0.00015616 m3/mol
Calcite
    CaCO3 = CO3-2 + Ca+2
    -analytical_expression -149.78 -0.04837 4897.4 60.458 76.464 0
    -Vm       3.69e-05 m3/mol
Chlorite(14A)
    Mg5Al2Si3O10(OH)8 + 16H+ = 2Al+3 + 6H2O + 3H4SiO4 + 5Mg+2
    log_k     68.38
    delta_h   -151.494 kcal
    -Vm       0.00019713 m3/mol
Dolomite
    CaMg(CO3)2 = 2CO3-2 + Ca+2 + Mg+2
    -analytical_expression -317.82 -0.098179 10845 126.57 169.32 0
    -Vm       6.45e-05 m3/mol
Illite
    K0.6Mg0.25Al2.3Si3.5O10(OH)2 + 11.2H2O = 2.3Al(OH)4- + 1.2H+ + 3.5H4SiO4 + 0.6K+ + 0.25Mg+2
    -analytical_expression 26.069 -0.0012553 13670 -20.232 -1120400 0
    -Vm       0.00014148 m3/mol
K-feldspar
    KAlSi3O8 + 8H2O = Al(OH)4- + 3H4SiO4 + K+
    -analytical_expression -1.0684 0.013111 11671 -9.9129 -1585500 0
    -Vm       0.00010815 m3/mol
Kaolinite
    Al2Si2O5(OH)4 + 6H+ = 2Al+3 + H2O + 2H4SiO4
    -analytical_expression 16.835 -0.0078939 7763.6 -12.19 -323540 0
    -Vm       9.935e-05 m3/mol
Pyrite
    FeS2 + 2H+ + 2e- = Fe+2 + 2HS-
    -analytical_expression -241.95 -0.087948 -629.11 99.248 -9.7454 0
    -Vm       2.348e-05 m3/mol
Quartz
    SiO2 + 2H2O = H4SiO4
    -analytical_expression 0.077698 0.010612 3465.1 -4.3551 -721380 0
    -Vm       2.267e-05 m3/mol
CO2(g)
    CO2 = CO2
    -analytical_expression 10.5624 -0.023547 -3972.8 0 587460 1.9194e-05
    -T_c      304.2
    -P_c      72.86
    -Omega    0.225
REACTION 1
    CO2(g)     1
    10 moles
REACTION_TEMPERATURE 1
    50
REACTION_PRESSURE 1
    79
KINETICS 1
Albite
    -formula  NaAlSi3O8  1
    -m        0.025
    -m0       0.025
    -parms    980 1e-12 0.00010131 244.6189
    -tol      1e-08
Anhydrite
    -formula  CaSO4  1
    -m        0.007
    -m0       0.007
    -parms    50000 1.5e-06 4.61e-05 162.6931
    -tol      1e-08
Ca-Montmorillonite
    -formula  Ca0.165Al2.33Si3.67O10(OH)2  1
    -m        0.005
    -m0       0.005
    -parms    10000000 4e-13 0.00015616 10.4996
    -tol      1e-08
Calcite
    -formula  CaCO3  1
    -m        0.133
    -m0       0.133
    -parms    50000 1.5e-06 3.69e-05 162.6931
    -tol      1e-08
Chlorite_(14A)
    -formula  Mg5Al2Si3O10(OH)8  1
    -m        0.001
    -m0       0.001
    -parms    980 2.51e-12 0.00019713 162.6931
    -tol      1e-08
Dolomite
    -formula  CaMg(CO3)2  1
    -m        0.083
    -m0       0.083
    -parms    127000 1e-09 6.45e-05 162.6931
    -tol      1e-08
Illite
    -formula  K0.6Mg0.25Al2.3Si3.5O10(OH)2  1
    -m        0.021
    -m0       0.021
    -parms    461000 1e-13 0.00014148 5.9523
    -tol      1e-08
K-feldspar
    -formula  KAlSi3O8  1
    -m        0.025
    -m0       0.025
    -parms    111300 1.778e-10 0.00010815 66.1464
    -tol      1e-08
Kaolinite
    -formula  Al2Si2O5(OH)4  1
    -m        0.036
    -m0       0.036
    -parms    10000000 4e-13 9.935e-05 10.4996
    -tol      1e-08
Pyrite
    -formula  FeS2  1
    -m        0.011
    -m0       0.011
    -parms    1290 4e-11 2.348e-05 162.6931
    -tol      1e-08
Quartz
    -formula  SiO2  1
    -m        0.588
    -m0       0.588
    -parms    57410 1.035e-14 2.267e-05 1225.096
    -tol      1e-08
-steps       432000 in 5 steps # seconds
-step_divide 100
-runge_kutta 3
-bad_step_max 500
RATES
    Albite
-start
 10 R=8.3143
 20 K=PARM(2)* PARM(4)
 30 if m>m0 then m=m0
 40 rate=K*PARM(1)*(1+1.41)
 60 moles=rate*TIME
200 SAVE moles
-end
    Anhydrite
-start
 10 r=8.3143
 20 k1=PARM(2)*PARM (4)
 40 if m>m0 then m=m0
 60 rate=k1*PARM(1)*(1+3.45)
 70 moles=rate*TIME
200 SAVE moles
-end
    Ca-Montmorillonite
-start
 10 r1=8.3143
 20 k2=PARM(2)*PARM(4)
 40 if m>m0 then m=m0
 60 rate=k2*PARM(1)*(1+2.29)
 70 moles=rate*TIME
200 SAVE moles
-end
    Calcite
-start
 10 R1=8.3143
 20 k3=PARM(2)*PARM (4)
 40 if m>m0 then m=m0
 60 rate=k3*PARM(1)*(1-0.00)
 70 moles=rate*TIME
200 SAVE moles
-end
    Chlorite_(14A)
-start
 10 R2=8.3143
 20 k4=PARM(2)*PARM(4)
 40 if m>m0 then m=m0
 60 rate=k4*PARM(1)*(1-0.00)
 70 moles=rate*TIME
200 SAVE moles
-end
    Dolomite
-start
 10 R3=8.3143
 20 K5=PARM(2)*PARM(4)
 40 if m>m0 then m=m0
 60 rate=K5*PARM(1)*(1-0.00)
 70 moles=rate*TIME
200 SAVE moles
-end
    Illite
-start
 10 R4=8.3143
 20 K6=PARM(2)*PARM(4)
 40 if m>m0 then m=m0
 60 rate=K6*PARM(1)*(1+1.04)
 70 moles=rate*TIME
200 SAVE moles
-end
    K-feldspar
-start
 10 R5=8.3143
 20 K7=PARM(2)*PARM(4)
 40 if m>m0 then m=m0
 60 rate= K7*PARM(1)*(1-0.00)
 70 moles=rate*TIME
200 SAVE moles
-end
    Pyrite
-start
 10 R6=8.3143
 20 K8=PARM(2)*PARM(4)
 40 if m>m0 then m=m0
 60 rate=K8*PARM(1)*(1-0.00)
 70 moles=rate*TIME
200 SAVE moles
-end
    Quartz
-start
 10 R7=8.3143
 20 K9=PARM(2)*PARM(4)
 40 if m>m0 then m=m0
 60 rate=K9*PARM(1)*(1-0.00)
 70 moles=rate*TIME
200 SAVE moles
-end
    Kaolinite
-start
 10 R8=8.3143
 20 K10=PARM(2)*PARM(4)
 40 if m>m0 then m=m0
 60 rate=K10*PARM(1)*(1-0.00)
 70 moles=rate*TIME
200 SAVE moles
-end
INCREMENTAL_REACTIONS True
END

[dlparkhurst]

Patience is a virtue, but you have too much. You should make sure things are working the way you expect before you invest three days in a calculation.

The main problem here is that you have defined both EQUILIBRIUM_PHASES and KINETICS for the same minerals. You should define one or the other for each mineral; either the mineral is always in equilibrium, or it is reacting toward equilibrium with KINETICS. I would probably use equilibrium for Calcite, dolomite, and anhydrite, but follow the description in the paper you are trying to reproduce.

To test rates calculation, I suggest you start by running for 1, 10, 100, ... seconds to get a feel for the relative rates of reactions.

Your rate definitions should have a factor that causes the rates to approach zero at equilibrium. Often that factor is [1 - SR("xxx")].

You cannot set M in the RATES definitions. Remove all instances of "if m>m0 then m=m0". It is not necessary to check for removal of all of a mineral in the RATES definition.

Installation files for PHREEQC version 3.8.0 programs are https://github.com/usgs-coupled/. Look in the Releases section on the right side of the page for each program.

In version 3.8.0, there is a new rates database from Oelkers and coworkers named kinec.v2.dat. It has RATES definitions for a large set of minerals.

There is also a database phreeqc_rates.dat that has Palandri and Kharaka rate parameters in the data block RATE_PARAMETERS_PK. There is a Basic function RATE_PK that calculates the Palandri and Kharakha rate for a mineral based on the parameters (mol/m^2/s). However, you would need to provide the RATES definitions that include the surface area and affinity factors. See comments within the database for example RATES definitions.

Alternatively, you can get files for Palandri and Kharaka (with some modifications) from Chen Zhu's compilation at the University of Indiana (https://hydrogeochem.earth.indiana.edu/software/index.html).

[dkonan01]

Thank you very much sir for your help.

[dlparkhurst]

Here is a script with some rate examples using phreeqc_rates.dat.

Code: [Select]

DATABASE ../database/phreeqc_rates.dat

RATES
Albite_PK # Palandri and Kharaka, 2004
5  REM PARMS: 1 affinity, 2 m^2/mol, 3 roughness, 4 exponent
10 if parm(1) = 1 then affinity = 1 else affinity = 1 - SR("Albite") : if affinity < parm(1) then SAVE 0 : END
20 rate = RATE_PK("Albite")
30 IF M > 0 THEN area = M * parm(2) * parm(3) * (M/M0)^parm(4) ELSE area = 0
40 SAVE area * rate * affinity * TIME
-end

KINETICS 1
Albite_PK
-formula NaAlSi3O8; -parms  0  1  1  0.67
-m0 1; -time 1 # default
END

SOLUTION 1
PHASES
  Fix_pH; H+ = H+
  LiBr; LiBr = Li+ + Br-; -log_k -20  # (very) unsoluble phase with base cation and acid anion, permits to use HBr or LiOH as reactant
SELECTED_OUTPUT 1
        -file kinetic_rates_pH.inc
        -reset false
USER_PUNCH 1    # write out the pH's to equilibrate...
        10 FOR i = 0 to 14 STEP 0.5
        20      punch EOL$ + 'USE solution 1'
        30      punch EOL$ + 'EQUILIBRIUM_PHASES 1'
        40      punch EOL$ + ' LiBr'
        50      punch EOL$ + ' Fix_pH ' + TRIM(STR$(-i)) + ' LiOH 10' # ...or HBr as reactant
        60      punch EOL$ + 'USE kinetics 1'
        70      punch EOL$ + 'END'
        80 NEXT i
END

PRINT; -reset false
SELECTED_OUTPUT 1; -active false
USER_GRAPH 1; -headings pH Palandri
-axis_titles  pH "log10(initial rate / (mol / m2 / s))"
-axis_scale x_axis 0 14
10 graph_x -la("H+")
20 graph_sy log10(tot("Al"))
INCLUDE$ kinetic_rates_pH.inc
END

RATES
Albite_Svd # Sverdrup, 2019
5  REM PARMS: 1 affinity, 2 m^2/mol, 3 roughness, 4 exponent
10 if parm(1) = 1 then affinity = 1 else affinity = 1 - SR("Albite") : if affinity < parm(1) then SAVE 0 : END
20 rate = RATE_SVD("Albite")
30 IF M > 0 THEN area = M * parm(2) * parm(3) * (M/M0)^parm(4) ELSE area = 0
40 SAVE area * rate * affinity * TIME
-end

KINETICS 1
Albite_Svd
-formula NaAlSi3O8; -parms  0  1  20  0.67 # roughness = 20
USER_GRAPH 1; -headings pH Sverdup*20
INCLUDE$ kinetic_rates_pH.inc
END

KINETICS 1
Albite # from Sverdrup and Warfvinge, 1995
-formula NaAlSi3O8; -parms  1  20 # roughness = 20
USER_GRAPH 1; -headings pH Sverdup`95*20
INCLUDE$ kinetic_rates_pH.inc
END

RATES
Albite_Hermanska # 2022
5  REM PARMS: 1 affinity, 2 m^2/mol, 3 roughness, 4 exponent
10 if parm(1) = 1 then affinity = 1 else affinity = 1 - SR("Albite") : if affinity < parm(1) then SAVE 0 : END
20 rate = RATE_HERMANSKA("Albite")
30 IF M > 0 THEN area = M * parm(2) * parm(3) * (M/M0)^parm(4) ELSE area = 0
40 SAVE area * rate * affinity * TIME
-end

KINETICS 1
Albite_Hermanska
-formula NaAlSi3O8; -parms  0  1  1  0.67
USER_GRAPH 1; -headings pH Hermanska
INCLUDE$ kinetic_rates_pH.inc
END

USE solution 1
REACTION_TEMPERATURE 1; 25 25 in 21
USER_GRAPH 1; -headings Albite_data
10 data 1.1, 2.05, 2.45, 2.9, 3, 3.5, 4.1, 5.1, 5.35, 5.47, 5.63, 5.63, 5.73, 7.73, 9.95, 9.95, 9.95, 10.6, 11.2, 11.55, 12.3
20 data -10.25, -10.55, -10.82, -11.25, -11.1, -11.4, -11.47, -11.82, -11.75, -11.65, -11.83, -11.92, -11.92, -11.83, -10.97, -11.05, -11.13, -10.95, -10.55, -10.6, -10.38 # Chou, L., Wollast, R., 1985. Steady-state kinetics and dissolution mechanisms of albite. Am. J. Sci. 285, 963?993.
30 restore 10 : dim ph(21) : for i = 1 to step_no : read ph(i) : next i
40 restore 20 : dim lk(21) : for i = 1 to step_no : read lk(i) : next i
50 i = step_no : plot_xy ph(i), lk(i), line_width = 0, color = Black, y_axis = 2, symbol_size = 10, symbol = Circle
END

# compare rates for calcite dissolution
#   of Palandri and Kharaka, 2004 and Plummer, Wigley and Parkhurst, 1978
#   at different initial CO2 concentrations.
# =====================================

USER_GRAPH 1; -active false

RATES
Calcite_PK # Palandri and Kharaka, 2004
5  REM PARMS: 1 affinity, 2 m^2/mol, 3 roughness, 4 exponent
10 if parm(1) = 1 then affinity = 1 else affinity = 1 - SR("calcite") : if affinity < parm(1) then SAVE 0 : END
20 rate = RATE_PK("calcite")
30 IF M > 0 THEN area = M * parm(2) * parm(3) * (M/M0)^parm(4) ELSE area = 0
40 SAVE area * rate * affinity * TIME
-end

SOLUTION 1
pH 7 charge; C(4) 1 CO2(g) -2.5
KINETICS 1
calcite_PK
-formula CaCO3; -parms  0  1e-2  1  0.67
-time 0.1 10*1 hour
INCREMENTAL_REACTIONS true
USER_GRAPH 2; -headings h Palandri_SI(CO2_g).=.-2.5
-axis_titles "time / hours" "Calcite dissolved / (mmol/kgw)"
10 graph_x total_time / 3600 : graph_sy tot("Ca") * 1e3
END

USE solution 1
KINETICS 1
Calcite
-parms  1e2  0.67  # cm^2/mol calcite, exp factor
-time  0.1 10*1 hour
USER_GRAPH 2; -headings h Plummer.Wigley.Parkhurst
END

SOLUTION 1
pH 7 charge; C(4) 1 CO2(g) -1.5
KINETICS 1
calcite_PK
-formula CaCO3
-parms  0  1e-2  1  0.67
-time 0.1 10*1 hour
USER_GRAPH 2; -headings h Palandri_SI(CO2_g).=.-1.5
END

USE solution 1
KINETICS 1
Calcite
-parms  1e2  0.67
-time  0.1 10*1 hour
USER_GRAPH 2; -headings h Plummer.Wigley.Parkhurst
END

# compare rates for quartz dissolution
#   and the effect of NaCl
# =====================================

USER_GRAPH 2; -active false

RATES
Quartz_PK # Palandri and Kharaka, 2004
5  REM PARMS: 1 affinity, 2 m^2/mol, 3 roughness, 4 exponent
10 if parm(1) = 1 then affinity = 1 else affinity = 1 - SR("Quartz") : if affinity < parm(1) then SAVE 0 : END
20 rate = RATE_PK("Quartz")
30 IF M > 0 THEN area = M * parm(2) * parm(3) * (M/M0)^parm(4) ELSE area = 0
40 SAVE area * rate * affinity * TIME
-end

SOLUTION 1
pH 7 charge
KINETICS 1
Quartz_PK
-formula SiO2
-parms  0  6  1  0.67
-time 0.1 10*1 year
INCREMENTAL_REACTIONS true
USER_GRAPH 3; -headings h Palandri
-axis_titles "time / years" "Quartz dissolved / (mmol/kgw)"
10 graph_x total_time / 3.15e7 : graph_sy tot("Si") * 1e3
END

RATES
Quartz_Hermanska #
5  REM PARMS: 1 affinity, 2 m^2/mol, 3 roughness, 4 exponent
10 if parm(1) = 1 then affinity = 1 else affinity = 1 - SR("Quartz") : if affinity < parm(1) then SAVE 0 : END
20 rate = RATE_HERMANSKA("Quartz")
30 IF M > 0 THEN area = M * parm(2) * parm(3) * (M/M0)^parm(4) ELSE area = 0
40 SAVE area * rate * affinity * TIME
-end

USE solution 1
KINETICS 1
Quartz_Hermanska
-formula SiO2
-parms  0  6  1  0.67
-time 0.1 10*1 year
USER_GRAPH 3
-headings H Hermanska
END

RATES
Quartz_Svd # Sverdrup, 2019
5  REM PARMS: 1 affinity, 2 m^2/mol, 3 roughness, 4 exponent
10 if parm(1) = 1 then affinity = 1 else affinity = 1 - SR("Quartz") : if affinity < parm(1) then SAVE 0 : END
20 rate = RATE_SVD("Quartz")
30 IF M > 0 THEN area = M * parm(2) * parm(3) * (M/M0)^parm(4) ELSE area = 0
40 SAVE area * rate * affinity * TIME
-end

USE solution 1
KINETICS 1
Quartz_Svd
-formula SiO2
-parms  0  6  1  0.67
-time 0.1 10*1 year
USER_GRAPH 3
-headings H Sverdup
END

RATES
Quartz_Rimstidt_Barnes
#1 rem Specific rate k = 10^-13.7 mol/m2/s (25 C), Ea = 90 kJ/mol, Rimstidt and Barnes, 1980, GCA 44, 1683
5  REM PARMS: 1 affinity, 2 m^2/mol, 3 roughness, 4 exponent
10 if parm(1) = 1 then affinity = 1 else affinity = 1 - SR("Quartz") : if affinity < parm(1) then SAVE 0 : END
20 rate = 10^-(13.7 + 4700 * (1 / 298 - 1 / TK)) * (1 + 1500*tot("Na")) # salt correction, Dove and Rimstidt, 1994, MSA Rev. 29, 259
30 IF M > 0 THEN area = M * parm(2) * parm(3) * (M/M0)^parm(4) ELSE area = 0
40 SAVE area * rate * affinity * TIME
-end

USE solution 1
KINETICS 1
Quartz_Rimstidt_Barnes
-formula SiO2
-parms  0  6  1  0.67
-time 0.1 10*1 year
USER_GRAPH 3
-headings H Rimstidt.et.al
END

SOLUTION 1
pH 7 charge; Na 30; Cl 30
KINETICS 1
Quartz_Rimstidt_Barnes
-formula SiO2
-parms  0  6  1  0.67
-time 0.1 10*1 year
USER_GRAPH 3
-headings H Rimstidt.et.al._NaCl
END


[dkonan01]

Got it Sir Parkhurst. Thank you very much