Dissolution of carbonate minerals in a CO2-saturated brine

In this tutorial, we demonstrate how Reaktoro can be used to kinetically model the dissolution of carbonate minerals (calcite, magnesite and dolomite) in a CO2-saturated brine.

Below is the entire Python script showing the calculation, followed by a step-by-step explanation. You may want to just have a quick look at this relatively long script, and then proceed to the next sections for a more detailed discussion!

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

# Step 2: Specify the phases in the chemical system and their species
editor = ChemicalEditor()
editor.addAqueousPhaseWithElementsOf("H2O NaCl CaCO3 MgCO3")
editor.addGaseousPhase(["H2O(g)", "CO2(g)"])
editor.addMineralPhase("Calcite")
editor.addMineralPhase("Magnesite")
editor.addMineralPhase("Dolomite")
editor.addMineralPhase("Halite")

# Step 3: Define the kinetically-controlled reactions
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")

editor.addMineralReaction("Magnesite") \
    .setEquation("Magnesite = Mg++ + CO3--") \
    .addMechanism("logk = -9.34 mol/(m2*s); Ea = 23.5 kJ/mol") \
    .addMechanism("logk = -6.38 mol/(m2*s); Ea = 14.4 kJ/mol; a[H+] = 1.0") \
    .setSpecificSurfaceArea(10, "cm2/g")

editor.addMineralReaction("Dolomite") \
    .setEquation("Dolomite = Ca++ + Mg++ + 2*CO3--") \
    .addMechanism("logk = -7.53 mol/(m2*s); Ea = 52.2 kJ/mol") \
    .addMechanism("logk = -3.19 mol/(m2*s); Ea = 36.1 kJ/mol; a[H+] = 0.5") \
    .setSpecificSurfaceArea(10, "cm2/g")

# Step 4: Construct the chemical system
system = ChemicalSystem(editor)
reactions = ReactionSystem(editor)

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

# Step 6: Define the initial chemical equilibrium state
problem = EquilibriumProblem(system)
problem.setPartition(partition)
problem.setTemperature(60, "celsius")
problem.setPressure(100, "bar")
problem.add("H2O", 1, "kg")
problem.add("NaCl", 1, "mol")
problem.add("CO2", 1, "mol")

# Step 7: Calculate the initial chemical equilibrium state
state0 = equilibrate(problem)
state0.output('demo-kineticpath-carbonates-co2-before-kinetics')

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

# Step 9: Create a kinetic path solver
path = KineticPath(reactions)
path.setPartition(partition)

# Step 10: Create plots for the kinetic path calculation
plot0 = path.plot()
plot0.x("time(units=hour)")
plot0.y("pH")
plot0.xlabel("Time [hour]")
plot0.ylabel("pH")
plot0.showlegend(False)

plot1 = path.plot()
plot1.x("time(units=hour)")
plot1.y("elementMolality(Ca)", "Ca")
plot1.y("elementMolality(Mg)", "Mg")
plot1.xlabel("Time [hour]")
plot1.ylabel("Concentration [molal]")
plot1.legend("right center")

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

plot3 = path.plot()
plot3.x("time(units=hour)")
plot3.y("phaseMass(Dolomite units=grams)", "Dolomite")
plot3.xlabel("Time [hour]")
plot3.ylabel("Mass [g]")

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

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

Importing the reaktoro Python package

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

Here we import the reaktoro Python package to enable us to use all its library components (classes, methods, constants).

Defining the chemical system

In this simulation, we consider an aqueous phase (to model the brine solution), a gaseous phase (to model the CO2-rich phase with water vapour), and four pure mineral phases: halite, calcite, magnesite, and dolomite. These are the phases that will either exist initially during the simulation, or that could potentially appear later as the calculations proceed.

Attention

All potential phases that could appear in a reactive process should ideally be considered when defining the chemical system. If one or more of these phases are ignored, then the equilibrium and kinetics calculations cannot identify if they should actually exist or not with positive amounts. Unrealistic results may occur, such as, for example, an aqueous phase containing more CO2 dissolved than it could, because a gaseous phase, which should contain the excess of CO2, was not considered in the chemical system.

The code below uses class ChemicalEditor to define our chemical system with the phases of interest and their species:

# Step 2: Specify the phases in the chemical system and their species
editor = ChemicalEditor()
editor.addAqueousPhaseWithElementsOf("H2O NaCl CaCO3 MgCO3")
editor.addGaseousPhase(["H2O(g)", "CO2(g)"])
editor.addMineralPhase("Calcite")
editor.addMineralPhase("Magnesite")
editor.addMineralPhase("Dolomite")
editor.addMineralPhase("Halite")

The aqueous phase is defined by considering all aqueous species in the database that could form once the substances H2O, NaCl, CaCO3, and MgCO3 are mixed. The gaseous phase is defined so that only the gaseous species H2O(g) and CO2(g) are considered. There are four pure mineral phases: calcite (CaCO3), magnesite (MgCO3), dolomite (CaMg(CO3)2), and halite (NaCl).

Defining the kinetically-controlled reactions

A partial equilibrium assumption is considered in this simulation. This simplifies the problem by assuming that those species that react at much faster rates are continually in chemical equilibrium. They are referred to as equilibrium species. The remaining species (those reacting at relatively slower rates) are referred to as kinetic species.

Because aqueous and gaseous species, as well as halite, react at relatively fast rates, they are reasonable candidates in this problem for being equilibrium species. The kinetic species are thus the carbonate minerals: calcite, magnesite, and dolomite.

With this partial equilibrium assumption, there is no need to specify kinetic rate laws for the fast reacting equilibrium species. For these species, chemical equilibrium calculations are performed to update their amounts as the amounts of each kinetic species change over time (i.e., the equilibrium species react instantaneously to a new state of equilibrium as the kinetic species react and perturb the current partial equilibrium state).

Thus, we only need to define kinetic rates for the relatively slow reacting carbonate minerals (our kinetic species in this simulation), which is shown below:

# Step 3: Define the kinetically-controlled reactions
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")

editor.addMineralReaction("Magnesite") \
    .setEquation("Magnesite = Mg++ + CO3--") \
    .addMechanism("logk = -9.34 mol/(m2*s); Ea = 23.5 kJ/mol") \
    .addMechanism("logk = -6.38 mol/(m2*s); Ea = 14.4 kJ/mol; a[H+] = 1.0") \
    .setSpecificSurfaceArea(10, "cm2/g")

editor.addMineralReaction("Dolomite") \
    .setEquation("Dolomite = Ca++ + Mg++ + 2*CO3--") \
    .addMechanism("logk = -7.53 mol/(m2*s); Ea = 52.2 kJ/mol") \
    .addMechanism("logk = -3.19 mol/(m2*s); Ea = 36.1 kJ/mol; a[H+] = 0.5") \
    .setSpecificSurfaceArea(10, "cm2/g")

We set the equation of each mineral reaction using method setEquation of 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.

We then prescribe neutral and acidic mechanisms for each mineral reaction using the method addMechanism of class MineralMechanism, using values for logk, the kinetic rate constant of the reaction in log scale, and Ea, the Arrhenius activation energy. The units of both parameters must be provided as shown in the example, and other compatible units are allowed.

Finally, we define the specific surface area of the mineral using method setSpecificSurfaceArea of class MineralReaction. Any units compatible to m2/kg or m2/m3 are allowed (e.g., cm2/g, mm2/mm3).

Note

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.

Creating the chemical and reaction systems

Create instances of ChemicalSystem and ReactionSystem.

# Step 4: Construct the chemical system
system = ChemicalSystem(editor)
reactions = ReactionSystem(editor)

Class ChemicalSystem is used to represent a chemical system, containing one or more phases (in this problem, aqueous, gaseous, and mineral phases), each phase containing one or more species (the aqueous phase containing many aqueous species, the gaseous phase containing two, and each mineral phase containing a single mineral species with the same name as the name of the phase). This class is also used to calculate thermodynamic properties of the phases and species (standard thermodynamic properties, activities, chemical potentials, phase molar volume, etc.).

Class ReactionSystem is a collection of Reaction objects used to represent a system of chemical reactions that are controlled by chemical kinetics. These classes provide convenient methods for the calculation of equilibrium constants, reaction quotients, and rates of the reactions.

Specifying equilibrium and kinetic species

In Reaktoro, the species in a chemical system can be partitioned into groups of:

  • equilibrium species,
  • kinetic species, and
  • inert species.

For the equilibrium species, their amounts are governed by chemical equilibrium (calculated via a Gibbs energy minimization). The amount of kinetic species are controlled by chemical kinetics (by solving a system of ordinary differential equations that models the kinetics of a system of reactions). The inert species maintain their amounts constant in all chemical equilibrium or kinetics calculations.

This classification of species can be done using class Partition. By default, all species are considered to be equilibrium species in this class. Thus, we only need to specify which ones are kinetic species:

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

In this case, the mineral species calcite, magnesite and dolomite are specified to be kinetic species using method setKineticSpecies of class Partition.

Note

Method setKineticPhases of class Partition could also be used here. This method sets all species in the given phases to be kinetic species, and it is more convenient if a phase has many species. However, since each of the mineral phases considered here only contains a single mineral species, method setKineticSpecies is a convenient alternative.

Defining the initial state of the equilibrium species

In a chemical kinetics calculation, an initial condition is needed for the amounts of both equilibrium and kinetic species.

The equilibrium problem formulated below, using class EquilibriumProblem, is done so that the initial condition for the amounts of each equilibrium species result from the solution of a chemical equilibrium problem in which 1 kg of water is mixed with 1 mol of NaCl and 1 mol of CO2 at 60 °C and 100 bar. This amount of CO2 is sufficient to saturate the brine solution. The excess will exist in the CO2-rich gaseous phase.

# Step 6: Define the initial chemical equilibrium state
problem = EquilibriumProblem(system)
problem.setPartition(partition)
problem.setTemperature(60, "celsius")
problem.setPressure(100, "bar")
problem.add("H2O", 1, "kg")
problem.add("NaCl", 1, "mol")
problem.add("CO2", 1, "mol")

Attention

To ensure that that the equilibrium calculation performed in the next step ignores the kinetic species in the system so that we maintain a disequilibrium state between equilibrium and kinetic species, it is important not to forget the following command:

problem.setPartition(partition)

Ignoring this step will produce an initial condition for the amounts of equilibrium and kinetic species that correspond to a complete equilibrium state in the system so that no kinetics calculation makes sense afterwards!

Calculating the initial amounts of the equilibrium species

We now use the convenient equilibrate function to calculate the amounts of the equilibrium species by minimizing the Gibbs energy of the equilibrium partition only, and not of the entire system. 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). We then output this chemical state to a file.

# Step 7: Calculate the initial chemical equilibrium state
state0 = equilibrate(problem)
state0.output('demo-kineticpath-carbonates-co2-before-kinetics')

Setting the initial mass of minerals

We have now to prescribe the initial amounts of the kinetic species (i.e., the carbonate minerals). This is done below by setting the initial mass of calcite to 100 g and of dolomite to 50 g. The initial amount of magnesite is zero.

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

Setting the kinetic path calculations

To be able to run the simulation of the kinetic process of mineral dissolution/precipitation, we introduce the instance of the class KineticPath that enables this functionality.

# Step 9: Create a kinetic path solver
path = KineticPath(reactions)
path.setPartition(partition)

The instance of kinetic path solver is provided with the partition to the equilibrium, kinetic, and inert species defined above.

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. Consider also the class EquilibriumSolver, instead of method equilibrate, for similar reasons.

Plotting the evolution of thermochemical properties

Usage of the ChemicalPlot allows us to create plots illustrating the evolution of different properties of the chemical system over the time interval in which kinetics is calculated.

# Step 10: Create plots for the kinetic path calculation
plot0 = path.plot()
plot0.x("time(units=hour)")
plot0.y("pH")
plot0.xlabel("Time [hour]")
plot0.ylabel("pH")
plot0.showlegend(False)

plot1 = path.plot()
plot1.x("time(units=hour)")
plot1.y("elementMolality(Ca)", "Ca")
plot1.y("elementMolality(Mg)", "Mg")
plot1.xlabel("Time [hour]")
plot1.ylabel("Concentration [molal]")
plot1.legend("right center")

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

plot3 = path.plot()
plot3.x("time(units=hour)")
plot3.y("phaseMass(Dolomite units=grams)", "Dolomite")
plot3.xlabel("Time [hour]")
plot3.ylabel("Mass [g]")

The code above plots four different figures with the evolution of pH, the concentration of calcium and magnesium, and the mass of calcite dolomite. In particular, plot1 shows how two properties, i.e., calcium and magnesium concentrations, with respect to the time (in hours). Methods plot1.x() and plot1.y() fetch properties that are meant to be plotted on x and y axises, respectively. The arguments of the function plot1.y("elementMolality(Ca)", "Ca") stand for the name of the property we want to fetch from the chemical state and the tag, we want to illustrate it with.

Note

A list of all possible quantities that can be plotted is shown in the class ChemicalQuantity.

Solving the chemical kinetics problem

Finally, we calculate the transient state of the entire system comprised of the carbonate minerals, the CO2-saturated brine fluid, and the CO2-rich gaseous phase. This is performed by using the method solve of the KineticPath class, which expects the initial state (state0), the initial and final times (t0 and t1 respectively) of the kinetic path, and time units of the specified time (e.g., s, minute, hour, day, year, etc.).

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

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!