Dissolution of calcite in an acidic HCl-solution

This tutorial demonstrates how Reaktoro can be used for modeling the dissolution of calcite in an acidic HCl-solution at temperature 30 °C and pressure 1 bar using chemical kinetics. A partial equilibrium assumption is considered here so that aqueous species react using a chemical equilibrium model, while calcite reacts with the aqueous solution using a chemical kinetics model.

The full script for the calculation is shown below, followed by a step-by-step explanation afterwards.

# Step 1: Import the reaktoro Python package
from reaktoro import *

# Step 2: Define a chemical system with an aqueous and a mineral phase
editor = ChemicalEditor()
editor.addAqueousPhaseWithElementsOf("H2O HCl CaCO3")
editor.addMineralPhase("Calcite")

# Step 3: Define mineral reaction for Calcite
editor.addMineralReaction("Calcite") \
    .setEquation("Calcite = Ca++ + CO3--") \
    .addMechanism("logk = -5.81 mol/(m2*s); Ea = 23.5 kJ/mol") \
    .addMechanism("logk = -0.30 mol/(m2*s); Ea = 14.4 kJ/mol; a[H+] = 1.0") \
    .setSpecificSurfaceArea(10, "cm2/g")

# Step 4: Create the ChemicalSystem and ReactionSystem instances
system = ChemicalSystem(editor)
reactions = ReactionSystem(editor)

# Step 5: Specify the equilibrium and kinetic species
partition = Partition(system)
partition.setKineticSpecies(["Calcite"])

# Step 6: Define the chemical equilibrium problem for the equilibrium partition
problem = EquilibriumProblem(system)
problem.setPartition(partition)
problem.setTemperature(30, "celsius")
problem.setPressure(1, "bar")
problem.add("H2O", 1, "kg")
problem.add("HCl", 1, "mmol")

# Step 7: Calculate the chemical equilibrium state in the equilibrium partition
state0 = equilibrate(problem)
state0.output('demo-kineticpath-calcite-hcl-before-kinetics')

# Step 8: Set the initial mass of the kinetic species
state0.setSpeciesMass("Calcite", 100, "g")

# Step 9: Define the kinetic path problem
path = KineticPath(reactions)
path.setPartition(partition)

# Step 10: Plot different properties of the chemical system during kinetics
plot1 = path.plot()
plot1.x("time(units=minute)")
plot1.y("elementMolality(Ca units=mmolal)", "Ca")
plot1.xlabel("Time [minute]")
plot1.ylabel("Concentration [mmolal]")
plot1.legend("right center")

plot2 = path.plot()
plot2.x("time(units=minute)")
plot2.y("phaseMass(Calcite units=g)", "Calcite")
plot2.xlabel("Time [minute]")
plot2.ylabel("Mass [g]")

plot3 = path.plot()
plot3.x("time(units=minute)")
plot3.y("pH")
plot3.xlabel("Time [minute]")
plot3.ylabel("pH")
plot3.legend("right center")

plot4 = path.plot()
plot4.x("time(units=minute)")
plot4.y("speciesMolality(Ca++ units=mmolal)", "Ca++")
plot4.y("speciesMolality(HCO3- units=mmolal)", "HCO3-")
plot4.xlabel("Time [minute]")
plot4.ylabel("Concentration [mmolal]")
plot4.legend("right center")

# Step 11: Perform the kinetic path calculation
t0, t1 = 0.0, 5.0
path.solve(state0, t0, t1, "minute")

# Step 12: Output the final state of the chemical system
state0.output('demo-kineticpath-calcite-hcl-after-kinetics')

Importing reaktoro Python package

# Step 1: Import the reaktoro Python package
from reaktoro import *

Here, we again import the reaktoro Python package so that we can use its classes and methods for performing the chemical reaction calculations.

Defining the phases in the chemical system

The class ChemicalEditor is used to conveniently create instances of classes ChemicalSystem and ReactionSystem.

# Step 2: Define a chemical system with an aqueous and a mineral phase
editor = ChemicalEditor()
editor.addAqueousPhaseWithElementsOf("H2O HCl CaCO3")
editor.addMineralPhase("Calcite")

In particular, we specify aqueous and mineral phases that should be considered in the chemical system. The aqueous phase is defined by the mixing of H₂O, HCl, and CaCO₃ (effectively, collecting all aqueous species in the database that contains elements H, O, C, Cl, and Ca, which are the elements in this list of compounds). There is only one pure mineral phase: the calcite phase.

Defining kinetically-controlled reactions

Next, we define a mineral reaction for calcite and set its rate parameters.

# Step 3: Define mineral reaction for Calcite
editor.addMineralReaction("Calcite") \
    .setEquation("Calcite = Ca++ + CO3--") \
    .addMechanism("logk = -5.81 mol/(m2*s); Ea = 23.5 kJ/mol") \
    .addMechanism("logk = -0.30 mol/(m2*s); Ea = 14.4 kJ/mol; a[H+] = 1.0") \
    .setSpecificSurfaceArea(10, "cm2/g")

We set the reaction equation using the setEquation method of the class MineralReaction.

Todo

The current need for setting a reaction equation for mineral reactions will be removed in the future, which will simplify the setup of the problem.

Then we add two mineral kinetic mechanisms for the reaction: neutral and acidic. This is done with the addMechanism method, where we set, for example, logk, the kinetic rate constant of the reaction in log scale, and Ea, the Arrhenius activation energy. The values shown for logk and Ea were collected from:

Palandri, J.L., Kharaka, Y.K. (2004). A compilation of rate parameters of water-mineral interaction kinetics for application to geochemical modeling. U.S. Geological Survey Open File Report (Vol. 2004–1068). Menlo Park, California.

Note

The units of logk and Ea can be different, but convertible to mol/(m2*s) and J/mol respectively.

Finally, we provide the specific surface area of the mineral using method setSpecificSurfaceArea of class MineralReaction, which can be specified in units of m2/g or m2/m3. Compatible units are allowed, such as cm2/mg or m2/dm3, and combinations.

Creating the chemical and reaction systems

Create instances of ChemicalSystem and ReactionSystem.

# Step 4: Create the ChemicalSystem and ReactionSystem instances
system = ChemicalSystem(editor)
reactions = ReactionSystem(editor)

ChemicalSystem is a class that represents a system and its attributes and properties, such as phases (in our case aqueous and mineral ones), species, elements (5 in total, i.e., H, O, Ca, C, Cl), formula matrix, as well as chemical and thermodynamical model. Class ReactionSystem serves as a system of the chemical reaction by a collection of Reaction class instances. It provides convenient methods that calculate the equilibrium constants, reaction quotients, and rates of the reactions.

Specifying the equilibrium and kinetic species

For the partition of a chemical system into equilibrium and kinetic species, we use the class Partition. We only need to specify which species are kinetic species, and all others will be equilibrium species by default.

# Step 5: Specify the equilibrium and kinetic species
partition = Partition(system)
partition.setKineticSpecies(["Calcite"])

We set species Calcite (the only species in the mineral phase also called Calcite!) to be the only kinetic species. This will allow us to model the dissolution of calcite using chemical kinetics, while all other species (the aqueous species) are modelled using chemical equilibrium (i.e., their amounts are updated over time using chemical equilibrium calculations).

Note

In general, the chemical system can be partitioned into equilibrium, kinetic and inert species. For the equilibrium species, their amounts are governed by chemical equilibrium (calculated by solving a Gibbs energy minimization problem). The kinetic species are the species whose amounts are governed by chemical kinetics (calculated by solving a system of ordinary differential equations that model the kinetics of a system of reactions). The inert species are the species whose amounts are constant.

Attention

Aqueous species react among themselves at much faster rates than mineral dissolution reactions in general, and thus this partial equilibrium assumption is plausible and reasonably accurate.

Defining the initial state of the equilibrium species

After constructing the chemical system and specifying the partitioning of the species, we proceed to a chemical equilibrium calculation to set the initial state of the equilibrium species. For this, we use class EquilibriumProblem, as shown below:

# Step 6: Define the chemical equilibrium problem for the equilibrium partition
problem = EquilibriumProblem(system)
problem.setPartition(partition)
problem.setTemperature(30, "celsius")
problem.setPressure(1, "bar")
problem.add("H2O", 1, "kg")
problem.add("HCl", 1, "mmol")

We specify the equilibrium/kinetic partitioning of the chemical system using method setPartition of class EquilibriumProblem. We then prescribe what should be the initial state of the equilibrium species (the aqueous species in this case), before we start the chemical kinetics calculation that will simulate the dissolution of calcite in this aqueous fluid.

By mixing 1 kg of H₂O and 1 mmol of HCl at 30 °C and 1 bar, we should produce a chemical equilibrium state that corresponds to an acidic aqueous fluid. The species in this fluid will be in disequilibrium with Calcite (our single kinetic species in this setup) since only equilibrium species (i.e., the aqueous species) are considered during the next chemical equilibrium calculation.

Calculating the initial chemical equilibrium state of the fluid

We now use the function equilibrate to calculate the chemical equilibrium state of the equilibrium partition, not the entire chemical system.

# Step 7: Calculate the chemical equilibrium state in the equilibrium partition
state0 = equilibrate(problem)
state0.output('demo-kineticpath-calcite-hcl-before-kinetics')

For this calculation, Reaktoro uses an efficient Gibbs energy minimization algorithm to determine the amounts of the equilibrium species that correspond to a state of minimum Gibbs energy in the equilibrium partition only, at given conditions of temperature, pressure, and element amounts in the equilibrium partition. The result is stored in the object state0 of class ChemicalState, a computational representation of the state of a multiphase chemical system defined by its temperature (T), pressure (P), and vector of species amounts (n).

Note

Consider using class EquilibriumSolver instead of method equilibrate when a sequence of chemical equilibrium calculations are needed for increased performance. The method equilibrate is a convenient function for solving chemical equilibrium problems, but has some overhead (e.g., an object of class EquilibriumSolver is created in each call).

Setting the initial mass of calcite

To simulate the kinetic dissolution of calcite in the aqueous fluid we defined before, we need to specify its initial amount. Below, we set the initial mass of species Calcite to 100 g.

# Step 8: Set the initial mass of the kinetic species
state0.setSpeciesMass("Calcite", 100, "g")

Performing the kinetic path calculation

To be able to run the simulation of the chemical kinetic path, we use class KineticPath.

# Step 9: Define the kinetic path problem
path = KineticPath(reactions)
path.setPartition(partition)

Note that here again, we need to specify the partitioning of the chemical system into equilibrium, kinetic, and inert species.

Attention

For repeated chemical kinetics calculations (e.g., in a reactive transport simulation, where kinetics is performed for each mesh cell/node), consider using the class KineticSolver instead for avoiding some overhead of class KineticPath.

Plotting the evolution of different properties of the chemical system

Usage of the ChemicalPlot allows us to create plots that show the evolution of different properties of the chemical system during the reactive process.

# Step 10: Plot different properties of the chemical system during kinetics
plot1 = path.plot()
plot1.x("time(units=minute)")
plot1.y("elementMolality(Ca units=mmolal)", "Ca")
plot1.xlabel("Time [minute]")
plot1.ylabel("Concentration [mmolal]")
plot1.legend("right center")

plot2 = path.plot()
plot2.x("time(units=minute)")
plot2.y("phaseMass(Calcite units=g)", "Calcite")
plot2.xlabel("Time [minute]")
plot2.ylabel("Mass [g]")

plot3 = path.plot()
plot3.x("time(units=minute)")
plot3.y("pH")
plot3.xlabel("Time [minute]")
plot3.ylabel("pH")
plot3.legend("right center")

plot4 = path.plot()
plot4.x("time(units=minute)")
plot4.y("speciesMolality(Ca++ units=mmolal)", "Ca++")
plot4.y("speciesMolality(HCO3- units=mmolal)", "HCO3-")
plot4.xlabel("Time [minute]")
plot4.ylabel("Concentration [mmolal]")
plot4.legend("right center")

The code above plots four different graphics that show the concentration of calcium, mass of calcite, pH, and concentration of ions Ca²⁺ and HC0₃^- over time. In particular, plot4 shows the concentrations of ions Ca²⁺ and HC0₃^- with respect to the time (in hours). Methods plot4.x() and plot4.y() fetch properties that are meant to be plotted on x and y axes, respectively. The arguments of the function plot1.y(""speciesMolality(Ca++ units=mmolal)", "Ca++") stand for the name of the property from the chemical state we want to plot and the label to be used in the plot.

Note

A list of all possible quantities that can be plotted is shown in the class ChemicalQuantity, which provides an interface for convenient ways of their retrieval.

Solving the chemical kinetics problem

Finally, we solve the kinetic path problem.

# Step 11: Perform the kinetic path calculation
t0, t1 = 0.0, 5.0
path.solve(state0, t0, t1, "minute")

This step executes the method solve of class KineticPath, which requires the initial state of the system (100 g of calcite in disequilibrium with a 1 mmolal HCl aqueous solution at 30 °C and 1 bar, represented with the object state0), the initial and final time of the kinetic path calculation (t0 and t1, respectively), and the time unit of the specified time parameters (e.g., s, minute, day, year, etc.).

Have you got an issue?

Have you found any issue or error in this tutorial? Do you have any recommendations or you think something is not clear enough? Please, let us know by filling a new issue here:

Reaktoro’s GitHub Issues

You’ll need a GitHub account - but this is easy to sort out if you don’t have one yet!