PHREEQC can handle kinetic reactions with complicated rates in a general way with the BASIC compiler that is part of the program, but writing the rates may be cumbersome and prone to errors. However, for rate equations with the same form, a subroutine can be defined that only needs numbers for the individual minerals for doing the calculations. The numbers are entered with -parms in the KINETICS keyword, which hands them over to the rate that stores the numbers in memory and calls the subroutine. Rates have been compiled by Sverdrup and Warfvinge (1995) and Sverdrup et al. (2019) for practical calculation of weathering rates in soils, and a.o. by Palandri and Kharaka (2004), Zhang et al. (2019) and Hermanská et al. (2022, 2023), and an example is given with data from a table.
Kinetic rates are derived from the moles of mineral that reacts with time in an experiment:
where r is the rate (mol/m2/s), A is the surface area of the mineral, m is moles of mineral, t is time, and SR the saturation ratio. The factor (1 - SR) is the affinity that lets the rate go to zero when the solution is in equilibrium with the mineral.
The rates are a complex function of the concentrations of protons, hydroxyls, and other ions, and of course, of temperature and pressure. For example, the pH dependence of Na- and K-feldspar dissolution rates shown in Figure 1a indicates that the rate is smallest at neutral pH and increases towards lower and higher pH. We can imagine that it is related to protons entering the mineral at low pH, displacing the metal ions, and to detachment of the structural [Al,Si3]O8 tetrahedra into solute complexes of Al(OH)4- and H3SiO4- with OH- at high pH. To describe this behavior, the equation for the rate must have the form:
where k and n are coefficients that may depend on temperature, pressure and other solutes, and the brackets indicate activities of the protons and the hydroxyl ions.
How is the surface area of the mineral defined? It can be estimated assuming that the mineral grains are spheres or cubes with a smooth surface (Ageom), or by measurement with gas adsorption (ABET). Generally, the measured surface is larger than calculated, since grinding increases the surface roughness λ:
The surface roughness also increases with time because crystal defects and the more soluble members of twins and zoned crystals tend to dissolve faster, creating pits and dissolution valleys. Figure 1b illustrates the appearance of solution valleys in an etching experiment with an albite/orthoclase twin in which albite (NaAlSi3O8) dissolved faster than orthoclase (KAlSi3O8), in agreement with Figure 1a. Figure 1c shows that the surface roughness of freshly ground minerals is about 10, and that it can increase to over 500 in 1 Ma. Clearly, the surface area is an uncertain factor and its definition requires some expertise or local knowledge for field application.
Figure 1. (a) The measured dissolution rate of albite- and orthoclase-feldspars as a function of pH (Blum, 1994);
(b) solution valleys developing in albite/orthoclase twins (Lee et al., 1998);
(c) surface roughness of primary minerals in soils as a function of soil age, and of freshly ground minerals with error bars on the left (White and Brantley, 2003).
So, what is needed in a general model for calculating the dissolution rates is the specific surface area, the roughness, a parameter that allows to modify the affinity, and the numbers that define the rate as a function of the physical and chemical conditions. For example, the dissolution of albite with the numbers from Table 13 of Palandri and Kharaka:
KINETICS 1
Albite_PK
-formula NaAlSi3O8
# parms affinity_factor m^2/mol roughness, lgkH e_H nH, lgkH2O e_H2O, lgkOH e_OH nOH
# parm number 1 2 3, 4 5 6, 7 8, 9 10 11
-parms 0 1 1, -10.16 65.0 0.457, -12.56 69.8, -15.60 71.0 -0.572 # parms 4-11 from TABLE 13
-m0 1; -time 1
RATES
Albite_PK # name of the reactant, rate from Palandri and Kharaka, 2004
10 if parm(1) = 1 then affinity = 1 else affinity = 1 - SR("Albite") :\
if affinity < parm(1) then SAVE 0 : END # the rate is for dissolution only
20 put(affinity, -99, 1) # store number in memory
30
40 for i = 2 to 11 : put(parm(i), -99, i) : next i
50
60 SAVE calc_value("Palandri_rate")
-end
SOLUTION 1
The KINETICS block starts with the rate name 'Albite_PK', and the formula of the mineral since the rate name is different from the name in PHASES in PHREEQC.DAT. The first three 'parms' give the affinity factor, the specific surface area (m2/ mol) and the surface roughness, respectively. Then follow the numbers for the 3 pH domains, each with log k, the reaction enthalpy for calculating the temperature dependence, and the exponent for the H+ dependence. We will compare the calculated rate with experimental data from Chou and Wollast, and define m0 = 1 mol and time = 1 s.
In the rate, the parameters are stored in memory, and the rate is calculated by the subroutine calc_value("Palandri_rate"). Numbers for other minerals can be copied from the other tables; when unavailable enter -30 for log k, 0 for exponents and 1 for other parameters.
The subroutine is defined with keyword CALCULATE_VALUES:
CALCULATE_VALUES
Palandri_rate
10 affinity = get(-99, 1) # retrieve number from memory
20
30 REM # specific area m2/mol, surface roughness
40 sp_area = get(-99, 2) : roughness = get(-99, 3)
50
60 REM # temperature, gas constant
70 dif_temp = 1 / TK - 1 / 298 : R = 2.303 * 8.314e-3 : dT_R = dif_temp / R
80
90 REM # rate by H+
100 lgk_H = get(-99, 4) : e_H = get(-99, 5) : nH = get(-99, 6)
110 rate_H = 10^(lgk_H - e_H * dT_R) * ACT("H+")^nH
120
130 REM # rate by hydrolysis
140 lgk_H2O = get(-99, 7) : e_H2O = get(-99, 8)
150 rate_H2O = 10^(lgk_H2O - e_H2O * dT_R)
160
170 REM # rate by OH-
180 lgk_OH = get(-99, 9) : e_OH = get(-99, 10) : nOH = get(-99, 11)
190 rate_OH = 10^(lgk_OH - e_OH * dT_R) * ACT("H+")^nOH
200
210 rate = rate_H + rate_H2O + rate_OH
220 area = sp_area * M
230
240 rate = area * roughness * rate * affinity
250 SAVE rate * TIME
-end
The subroutine is available in the PHREEQC file kinetic_rates.dat. It can be used in other PHREEQC input files with INCLUDE$ kinetic_rates.dat, as in PHREEQC file kinetic_rates.phr.
Download the 2 files, install them in the same directory, and run the file kinetic_rates.phr to produce the graph shown below.
The circles are measured rates by Chou and Wollast, the red line gives the pH dependent rate according to Palandri, and the green and blue lines result when the compilations of Sverdrup et al. (2019) and Hermanská et al. (2022) are used (and for which other subroutines were written in kinetic_rates.dat).
Unfortunately, the results are not altogether satisfactory: the Sverdrup rate is an order of magnitude smaller than the data over the whole pH range, and the Palandri rate is an order too high at high pH. However, Figure 1a shows a similar divergence in the experimental data, probably related to uncertainty in defining the reactive surface area and the naturally varying properties of the minerals. The smaller rate in the Sverdrup and Warfvinge model is the outcome of many detailed field studies in Sweden, and the empirical connection makes this database useful for calculating silicate weathering in similar climates and in soils with an age of about 104 years. However, we expect that the surface roughness of the minerals is then about 20 (figure 1c), which would increase the logarithmic rate of dissolution by 1.3 rather than decrease it by that factor as the figure indicates.
Can you show the effect of changing the surface roughness from 1 to 20 in the rate for albite according to Sverdrup?
References
Blum, A.E., 1994. Feldspars in Weathering, in: Parsons, I. (Ed.), Feldspars and Their Reactions. Springer Netherlands, Dordrecht, pp. 595-630.
Chou, L. and Wollast, R., 1985. Steady-state kinetics and dissolution mechanisms of albite. Am. J. Sci. 285, 963-993.
Hermanská, M., Voigt, M.J., Marieni, C., Declercq, J. and Oelkers, E.H., 2022. A comprehensive and internally consistent mineral dissolution rate database: Part I: Primary silicate minerals and glasses. Chemical Geology, 597, p.120807
Hermanská, M., Voigt, M.J., Marieni, C., Declercq, J. and Oelkers, E.H., 2023. A comprehensive and consistent mineral dissolution rate database: Part II: Secondary silicate minerals. Chemical Geology, p.121632.
Lee, M.R., Hodson, M.E. and Parsons, I., 1998. The role of intragranular microtextures and microstructures in chemical and mechanical weathering: direct comparisons of experimentally and naturally weathered alkali feldspars. Geochim. Cosmochim. Acta 62, 2771-2788. And pers. comm.
Palandri, J.L. and Kharaka, J.K., 2004. A compilation of rate parameters of water-mineral interaction kinetics for application to geochemical modeling. USGS Open-File Report 2004-1068, 71 p.
Sverdrup, H. and Warfvinge, P., 1995. Estimating field weathering rates using laboratory kinetics. Rev. Mineral. Geochem. 31, 485-541.
Sverdrup, H., Oelkers, E., Erlandsson Lampa, M., Belyazid, S., Kurz, D. and Akselsson, C., 2019. Reviews and syntheses: Weathering of silicate minerals in soils and watersheds: Parameterization of the weathering kinetics module in the PROFILE and ForSAFE models. Biogeosciences Discuss. 1-58.
White, A.F. and Brantley, S.L., 2003. The effect of time on the weathering of silicate minerals: why do weathering rates differ in the laboratory and field? Chem. Geol. 202, 479-506.
Zhang, Y., Hu, B., Teng, Y., Tu, K., Zhu, C., 2019. A library of BASIC scripts of reaction rates for geochemical modeling using phreeqc. Comput. Geosci. 133, 104316.