Chemical kinetics for mineral reactions using Palandri-Kharaka model
Contents
Chemical kinetics for mineral reactions using Palandri-Kharaka model#
Written by Allan Leal (ETH Zurich) on Nov 30th, 2022
Attention
Always make sure you are using the latest version of Reaktoro. Otherwise, some new features documented on this website will not work on your machine and you may receive unintuitive errors. Follow these update instructions to get the latest version of Reaktoro!
In the previous tutorial, we learned how to define our own reaction rate model to calculate the dissolution of halite (NaCl) over time in pure water. In this tutorial, we’ll instead use a reaction rate model already available in Reaktoro based on the model of Palandri and Kharaka [2004].
Palandri-Kharaka reaction rate model#
Consider a mineral reacting with an aqueous electrolyte solution according to the reaction:
where \(\alpha_{i}\) is the \(i\)th species in the aqueous solution and \(\nu_{i}\) is its stoichiometric coefficient in the reaction. For example, for the mineral reaction:
\(\nu_i=1\) for both \(\mathrm{Ca}^{2+}\) and \(\mathrm{Mg}^{2+}\); \(\nu_i=2\) for \(\mathrm{CO_{3}^{2-}}\); and \(\nu_i=0\) for all other species.
The Palandri-Kharaka rate model for mineral dissolution and precipitation reactions is mathematically formulated as follows:
where \(r_m\) is the mineral reaction rate (in mol/s), \(\mathrm{SA}\) is the current surface area of the mineral in m2 and \(\Omega\) is its current saturation ratio:
where \(K\) is the equilibrium constant of the mineral reaction; \(a_{i}\) is the activity of the \(i\)th aqueous species; and \(\mathrm{IAP}\) is an acronym for ion activity product.
Note
Reaktoro adopts the following sign convention for reaction rates: positive when the reaction proceeds from left to right (i.e., reactants forming products), negative otherwise. Thus, for a mineral reaction, with the mineral considered as a reactant (which is the convention in Reaktoro), the reaction rate is positive when the mineral is dissolving and negative when precipitating.
Note
Reaktoro calculates the mineral saturation ratio, \(\Omega\), using dual chemical potentials of elements and charge determined from the current chemical potentials of the aqueous species (see Smith and Missen [1982], Kulik [2006] and Leal et al. [2017]). As long as the mineral reactions proceed under the assumption that all aqueous species are in equilibrium at all times (i.e., a partial equilibrium assumption in which the aqueous species instantaneously reequilibrate as soon as the mineral reacts slightly), then the mineral saturation ratio calculated using this methodology is identical to using equation (4). If this assumption is not applicable, equation (4) should then be used. In geochemical modeling, this partial equilibrium assumption is generally considered, since aqueous species react at much faster rates than that of mineral dissolution and precipitation.
Note
The Palandri–Kharaka reaction rate model (3) considers several reaction mechanisms (e.g., acidic, basic, neutral, carbonate mechanisms) and lumps them all in a single reaction rate model (that’s why there is a sum in its definition). There is actually no need to specify the aqueous product species of the reaction, as we’ll see later, because these are assumed to be instantaneously equilibrated as the mineral reaction proceeds.
As mentioned before, the Palandri-Kharaka rate model (3) lumps different reaction mechanisms that control the reactivity of the mineral and this is reflected by the sum term on the right-hand side. Each reaction mechanism depend on the following parameters:
\(k_{i}^{\circ}\), a reaction rate constant parameter (in mol/(m2·s));
\(E_{i}\), an activation energy parameter (in J/mol);
\(w_{i,j}\), an exponent parameter on the activity of species \(j\) (dimensionless);
\(p_{i}\), an exponent parameter on the saturation ratio \(\Omega\) of the mineral (dimensionless);
\(q_{i}\), an exponent parameter on the term \((1 - \Omega)^{p_i}\) that measures disequilibrium (dimensionless).
For mineral gibbsite (Al(OH)3), the Palandri-Kharaka rate model contains three reaction mechanisms: acid, neutral, and base, and it depends on the activity of species H+. Below are the values of the above parameters for each mechanism, taken from Palandri and Kharaka [2004]:
Mechanism |
\(\log k^\circ\) [log mol/(m2·s)] |
\(E\) [kJ/mol] |
\(w_{\mathrm{H^+}}\) |
|---|---|---|---|
acid |
-7.65 |
47.5 |
0.992 |
neutral |
-11.50 |
61.2 |
0.0 |
base |
-16.65 |
80.1 |
-0.784 |
Let’s next take advantage of Palandri–Kharaka rate model already implemented in Reaktoro to simulate the dissolution/precipitation of carbonate minerals.
Dissolving and precipitating carbonate minerals#
Let’s now simulate the simultaneous dissolution and precipitation of carbonate minerals: calcite (CaCO3) and magnesite (MgCO3) dissolve while dolomite (CaMg(CO3)2) precipitates. We consider the following configuration: 1 mol of calcite and 1 mol of magnesite are immersed in a CO2 saturated water solution under ambient conditions. The solution is prepared by mixing 1 kg of water with 5 moles of CO2(g), which is sufficient to saturate the water.
The reaction rates of minerals depend on their surface area exposed to the aqueous solution, which varies over time. Let’s assume that the surface area of all minerals changes linearly with their volume. At any point in the simulation, the surface area of a mineral is thus its current volume multiplied by 6 cm²/cm³ (as if it maintained a cubic geometry at all times).
Attention
This cubic geometric consideration is completely arbitrary. Mineral surface areas often change in unpredictable ways. Reaktoro, however, does provide the means to attach a custom surface area model (which you can make as complicated as you like!). This will be covered in future tutorials.
Defining the chemical system#
We define below a chemical system representative of the reacting materials involved, containing the following phases:
an aqueous phase, with constituent chemical species automatically selected based on the chemical elements present in the minerals and in the gaseous phase below
a gaseous phase, with only CO2(g) species (necessary because we want to saturate the aqueous phase with CO22)
mineral phases calcite (CaCO3), magnesite (MgCO3) and dolomite (CaMg(CO3)2).
In terms of activity models assigned to each phase, we select:
for the aqueous phase, a combination of Davies activity model (for solvent and solutes) and Drummond [1981] model for solute species CO2(aq)
for the gas phase, an activity model based on the Peng-Robinson equation of state [Peng and Robinson, 1976]
for mineral phases, these are pure minerals with unit activities, so the default ideal activity model for solids covers this automatically.
The code below creates the ChemicalSystem object with the specifications described above, with the mineral reactions defined with the Palandri–Kharaka model presented earlier:
from reaktoro import *
db = SupcrtDatabase("supcrtbl") # let's use the SUPCRTBL database for this simulation
params = Params.embedded("PalandriKharaka.yaml")
system = ChemicalSystem(db,
AqueousPhase().set(ActivityModelPitzer()),
GaseousPhase("CO2(g)").set(ActivityModelPengRobinsonPhreeqcOriginal()),
MineralPhases("Calcite Magnesite Dolomite"),
MineralReaction("Calcite").setRateModel(ReactionRateModelPalandriKharaka(params)),
MineralReaction("Magnesite").setRateModel(ReactionRateModelPalandriKharaka(params)),
MineralReaction("Dolomite").setRateModel(ReactionRateModelPalandriKharaka(params)),
MineralSurface("Calcite", 6.0, "cm2/cm3"),
MineralSurface("Magnesite", 6.0, "cm2/cm3"),
MineralSurface("Dolomite", 6.0, "cm2/cm3"),
)
For your convenience, please note above the usage of the embedded database file PalandriKharaka.yaml containing all parameters of mineral reactions from Palandri and Kharaka [2004]. Note also the class Params, which was specifically designed to make it easier to initialize models (e.g., reaction rate models, activity models) when many parameters are needed. In case you want to use a customized version of this embedded database file (e.g., with different parameter values, or extended for new minerals), you’ll need to use the Params.local method instead of Params.embedded as shown below:
params = Params.local("/home/username/path/to/CustomPalandriKharaka.yaml")
params = Params.local("C:\\Users\\Username\\Path\\To\\CustomPalandriKharaka.yaml")
Performing chemical kinetics steps#
We’ll now initiate the chemical kinetics computation, by letting the system to react for 1,500 days in 500 time steps. First, however, we need to create a ChemicalState object that will represent the initial state of the system, as detailed before:
state = ChemicalState(system)
state.set("H2O(aq)" , 1.0, "kg")
state.set("CO2(g)" , 5.0, "mol")
state.set("Calcite" , 1.0, "mol")
state.set("Magnesite", 1.0, "mol")
state.set("O2(aq)" , 1.0, "umol") # Important: when O2(aq) or H2(aq) are in the system, add an insignificant tiny amount of one of them to avoid numerical problems due to floating-point rounding errors!
We can now apply chemical kinetics on this chemical state multiple times and determine its composition as the reactions take place, which is done next:
solver = KineticsSolver(system) # the chemical kinetics solver to simulate the dissolution/precipitation of minerals
day = units.seconds(1.0, "day") # 1 day in seconds
duration = 1500 * day # simulate kinetics for 1500 days
steps = 500 # consider 500 time steps
dt = duration/steps # compute time step
# Initiate the chemical kinetics simulation
table = Table() # used to create a table of collected data during the simulation
for i in range(steps + 1):
result = solver.solve(state, dt) # compute the chemical state of the system after it reacted for given time length
assert result.succeeded(), f"Calculation failed at time step #{i}!"
props = state.props() # get the current thermodynamic and chemical properties of the system
aprops = AqueousProps.compute(props) # compute the current aqueous properties of the system
# At every time step, collect chemical state data to be used later for plotting!
table.column("Time") << i*dt / day # convert time from seconds to days
table.column("Calcite") << props.speciesAmount("Calcite") # in mol
table.column("Magnesite") << props.speciesAmount("Magnesite") # in mol
table.column("Dolomite") << props.speciesAmount("Dolomite") # in mol
table.column("Ca+2") << props.speciesAmount("Ca+2") # in mol
table.column("Mg+2") << props.speciesAmount("Mg+2") # in mol
table.column("RateCalcite") << props.reactionRate("Calcite") # in mol/s
table.column("RateMagnesite") << props.reactionRate("Magnesite") # in mol/s
table.column("RateDolomite") << props.reactionRate("Dolomite") # in mol/s
table.column("pH") << aprops.pH()
table.column("OmegaCalcite") << aprops.saturationRatio("Calcite")
table.column("OmegaMagnesite") << aprops.saturationRatio("Magnesite")
table.column("OmegaDolomite") << aprops.saturationRatio("Dolomite")
# Once time steps have finished, you may want to save collected data to a file (e.g., to use in a spreadsheet software)
table.save("data.txt")
Plotting collected chemical properties#
Let’s now create plots for the results we collected in each time step:
from reaktplot import *
fig = Figure()
fig.title("AQUEOUS SPECIES AMOUNTS OVER TIME")
fig.xaxisTitle("Time [day]")
fig.yaxisTitle("Amount [mol]")
fig.drawLine(table["Time"], table["Ca+2"], "Ca<sup>2+</sup>")
fig.drawLine(table["Time"], table["Mg+2"], "Mg<sup>2+</sup>")
fig.show()
fig = Figure()
fig.title("MINERAL AMOUNTS OVER TIME")
fig.xaxisTitle("Time [day]")
fig.yaxisTitle("Amount [mol]")
fig.drawLine(table["Time"], table["Calcite"], "Calcite")
fig.drawLine(table["Time"], table["Magnesite"], "Magnesite")
fig.drawLine(table["Time"], table["Dolomite"], "Dolomite")
fig.show()
fig = Figure()
fig.title("PH OVER TIME")
fig.xaxisTitle("Time [day]")
fig.yaxisTitle("pH")
fig.drawLine(table["Time"], table["pH"], "pH")
fig.show()
fig = Figure()
fig.title("REACTION RATES OVER TIME")
fig.xaxisTitle("Time [day]")
fig.yaxisTitle("Rate [mol/s]")
fig.drawLine(table["Time"], table["RateCalcite"], "Calcite")
fig.drawLine(table["Time"], table["RateMagnesite"], "Magnesite")
fig.drawLine(table["Time"], table["RateDolomite"], "Dolomite")
fig.show()
And that’s it! If you have suggestions on how to improve this tutorial, please contact us. Read on to learn more advanced chemical kinetics modeling capabilities.