{ "cells": [ { "cell_type": "markdown", "id": "4498c801", "metadata": {}, "source": [ "# Dissolution of carbonate minerals in a CO2-saturated brine\n", "\n", "In this tutorial, we demonstrate how Reaktoro can be used to kinetically model the dissolution of carbonate\n", "minerals (calcite, magnesite, and dolomite) in a CO2-saturated brine.\n", "\n", "### Defining the chemical system\n", "\n", "We start with importing the **reaktoro** Python package to enable us to use all its library components (classes,\n", "methods, constants)." ] }, { "cell_type": "code", "execution_count": null, "id": "67476708", "metadata": {}, "outputs": [], "source": [ "from reaktoro import *" ] }, { "cell_type": "markdown", "id": "309428e6", "metadata": {}, "source": [ "In this simulation, we consider an *aqueous phase* (to model the brine solution), a *gaseous phase* (to model the\n", "CO2-rich phase with water vapor), and four pure *mineral phases*: halite, calcite, magnesite,\n", "and dolomite. These are the phases that will either exist initially during the simulation, or that could\n", "potentially appear later as the calculations proceed.\n", "\n", "All potential phases that could appear in a reactive process should ideally be considered when defining the\n", "chemical system. If one or more of these phases are ignored, then the equilibrium and kinetics calculations cannot\n", "identify if they should be present or not with positive amounts. Unrealistic results may occur, such as,\n", "for example, an aqueous phase containing more CO2 dissolved than it could, because a gaseous phase,\n", "which should contain the excess of CO2, was not considered in the chemical system.\n", "\n", "The code below uses class [ChemicalEditor](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalEditor.html) to define\n", "our chemical system with the phases of interest and their species:" ] }, { "cell_type": "code", "execution_count": null, "id": "b133c9ba", "metadata": {}, "outputs": [], "source": [ "editor = ChemicalEditor()\n", "editor.addAqueousPhaseWithElementsOf(\"H2O NaCl CaCO3 MgCO3\")\n", "editor.addGaseousPhase([\"H2O(g)\", \"CO2(g)\"])\n", "editor.addMineralPhase(\"Calcite\")\n", "editor.addMineralPhase(\"Magnesite\")\n", "editor.addMineralPhase(\"Dolomite\")\n", "editor.addMineralPhase(\"Halite\")" ] }, { "cell_type": "markdown", "id": "fffc2ac1", "metadata": {}, "source": [ "The aqueous phase is defined by considering all aqueous species in the database that could form once the substances\n", "H2O, NaCl, CaCO3, and MgCO3 are mixed. The gaseous phase is defined\n", "so that only the gaseous species H2O(g) and CO2(g) are considered. There are four pure\n", "mineral phases: calcite (CaCO3), magnesite (MgCO3), dolomite (CaMg(CO3)2),\n", "and halite (NaCl).\n", "\n", "### Defining the kinetically-controlled reactions\n", "\n", "A *partial equilibrium assumption* is considered in this simulation. This simplifies the problem by assuming that\n", "those species that react at much faster rates do so *continually in chemical equilibrium*. They are referred to as\n", "*equilibrium species*. The remaining species (those reacting at relatively slower rates) are referred to as\n", "*kinetic species*.\n", "\n", "Because aqueous and gaseous species, as well as halite, react at relatively fast rates, they are reasonable\n", "candidates in this problem for being equilibrium species. The kinetic species are thus the carbonate minerals:\n", "calcite, magnesite, and dolomite.\n", "\n", "With this partial equilibrium assumption, there is no need to specify kinetic rate laws for the fast-reacting\n", "equilibrium species. For these species, chemical equilibrium calculations are performed to update their amounts as the\n", "amounts of each kinetic species change over time (i.e., the equilibrium species react instantaneously to a new\n", "state of equilibrium as the kinetic species react and perturb the current partial equilibrium state).\n", "\n", "Thus, we only need to define kinetic rates for the relatively slow-reacting carbonate minerals (our kinetic species\n", "in this simulation), which is shown below" ] }, { "cell_type": "code", "execution_count": null, "id": "4709cf76", "metadata": {}, "outputs": [], "source": [ "editor.addMineralReaction(\"Calcite\") \\\n", " .setEquation(\"Calcite = Ca++ + CO3--\") \\\n", " .addMechanism(\"logk = -5.81 mol/(m2*s); Ea = 23.5 kJ/mol\") \\\n", " .addMechanism(\"logk = -0.30 mol/(m2*s); Ea = 14.4 kJ/mol; a[H+] = 1.0\") \\\n", " .setSpecificSurfaceArea(10, \"cm2/g\")\n", "\n", "editor.addMineralReaction(\"Magnesite\") \\\n", " .setEquation(\"Magnesite = Mg++ + CO3--\") \\\n", " .addMechanism(\"logk = -9.34 mol/(m2*s); Ea = 23.5 kJ/mol\") \\\n", " .addMechanism(\"logk = -6.38 mol/(m2*s); Ea = 14.4 kJ/mol; a[H+] = 1.0\") \\\n", " .setSpecificSurfaceArea(10, \"cm2/g\")\n", "\n", "editor.addMineralReaction(\"Dolomite\") \\\n", " .setEquation(\"Dolomite = Ca++ + Mg++ + 2*CO3--\") \\\n", " .addMechanism(\"logk = -7.53 mol/(m2*s); Ea = 52.2 kJ/mol\") \\\n", " .addMechanism(\"logk = -3.19 mol/(m2*s); Ea = 36.1 kJ/mol; a[H+] = 0.5\") \\\n", " .setSpecificSurfaceArea(10, \"cm2/g\")" ] }, { "cell_type": "markdown", "id": "6a8dd5fe", "metadata": {}, "source": [ "We set the equation of each mineral reaction using method\n", "[MineralReaction::setEquation](https://reaktoro.org/cpp/classReaktoro_1_1MineralReaction.html#a6edcfede18eff82f74a7bd8a56f1bea8).\n", "\n", "It follows by prescribtion of neutral and acidic mechanisms for each mineral reaction using the method\n", "[MineralReaction::addMechanism](https://reaktoro.org/cpp/classReaktoro_1_1MineralReaction.html#a1bdeff51e51f42e4635208241cd54027)\n", "In particular, we set values for `logk`, the kinetic rate constant of the reaction in log scale, and `Ea`,\n", "the Arrhenius activation energy. The units of both parameters must be provided as shown in the example,\n", "and other compatible units are allowed.\n", "\n", "Finally, we define the specific surface area of the mineral using the method\n", "[MineralReaction::setSpecificSurfaceArea](https://reaktoro.org/cpp/classReaktoro_1_1MineralReaction.html#a9ea2feb68af0beddc856d6a60b863181).\n", "Any units compatible to m2/kg or m2/m3 are allowed (e.g., cm2/g,\n", "mm2/mm3).\n", "\n", "### Creating the chemical and reaction systems\n", "\n", "Create instances of\n", "[ChemicalSystem](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html) and\n", "[ReactionSystem](https://reaktoro.org/cpp/classReaktoro_1_1ReactionSystem.html)." ] }, { "cell_type": "code", "execution_count": null, "id": "50536e76", "metadata": {}, "outputs": [], "source": [ "system = ChemicalSystem(editor)\n", "reactions = ReactionSystem(editor)" ] }, { "cell_type": "markdown", "id": "d7778ce1", "metadata": {}, "source": [ "Class [ChemicalSystem](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html) is used to represent a\n", "chemical system, containing one or more phases (in this problem, aqueous, gaseous, and mineral phases). Each\n", "phase can contain one or more species (the aqueous phase includes many aqueous species, the gaseous phase two,\n", "and each mineral phase containing a single mineral species with the same name as the name of the phase). This class\n", "is also used to calculate the thermodynamic properties of the phases and species (standard thermodynamic properties,\n", "activities, chemical potentials, phase molar volume, etc.).\n", "\n", "Class [ReactionSystem](https://reaktoro.org/cpp/classReaktoro_1_1ReactionSystem.html) is a collection of\n", "[Reaction](https://reaktoro.org/cpp/classReaktoro_1_1Reaction.html) objects used to represent a system of chemical\n", "reactions that are controlled by chemical kinetics. These classes provide convenient methods for the calculation of\n", "equilibrium constants, reaction quotients, and rates of the reactions.\n", "\n", "### Specifying equilibrium and kinetic species\n", "\n", "In Reaktoro, the species in a chemical system can be partitioned into groups of:\n", "*equilibrium species*, *kinetic species*, and *inert species*.\n", "For the *equilibrium species*, their amounts are governed by chemical equilibrium (calculated via a Gibbs energy\n", "minimization). The amount of *kinetic species* are controlled by chemical kinetics (by solving a system of ordinary\n", "differential equations that models the kinetics of a system of reactions). The *inert species* maintain their\n", "amounts constant in all chemical equilibrium or kinetics calculations.\n", "\n", "This classification of species can be done using class\n", "[Partition](https://reaktoro.org/cpp/classReaktoro_1_1Partition.html). By default, all species are considered to\n", "be equilibrium species in this class. Thus, we only need to specify which are kinetic ones:" ] }, { "cell_type": "code", "execution_count": null, "id": "f8e9a652", "metadata": {}, "outputs": [], "source": [ "partition = Partition(system)\n", "partition.setKineticSpecies([\"Calcite\", \"Magnesite\", \"Dolomite\"])" ] }, { "cell_type": "markdown", "id": "7a1873b5", "metadata": {}, "source": [ "In this case, the mineral species calcite, magnesite, and dolomite are specified to be *kinetic species* using\n", "method [Partition::setKineticSpecies](https://reaktoro.org/cpp/classReaktoro_1_1Partition.html#a9691b3689b8a09280595f2b898cc5e3e).\n", "Method [Partition::setKineticPhases](\n", "https://reaktoro.org/cpp/classReaktoro_1_1Partition.html#a235c673476b0d6682148b70b04024499) could also\n", "be used here. It sets all species in the given phases to be kinetic species, and it is more convenient if\n", "a phase has many species. However, since each of the mineral phases considered here only contains a single mineral\n", "species, the method [Partition::setKineticSpecies](https://reaktoro.org/cpp/classReaktoro_1_1Partition.html#a9691b3689b8a09280595f2b898cc5e3e)\n", "is a convenient alternative.\n", "\n", "### Defining the initial state of the equilibrium species\n", "\n", "In a chemical kinetics calculation, an *initial condition* is needed for the amounts of both equilibrium and\n", "kinetic species.\n", "\n", "The equilibrium problem formulated below, using class\n", "[EquilibriumProblem](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumProblem.html), is done so that the\n", "initial condition for the amounts of each equilibrium species result from the solution of a chemical equilibrium\n", "problem in which 1 kg of water is mixed with 1 mol of NaCl and 1 mol of CO2 at 60 °C\n", "and 100 bar. This amount of CO2 is sufficient to saturate the brine solution. The excess will exist in\n", "the CO2-rich gaseous phase." ] }, { "cell_type": "code", "execution_count": null, "id": "243ad05d", "metadata": {}, "outputs": [], "source": [ "problem = EquilibriumProblem(system)\n", "problem.setPartition(partition)\n", "problem.setTemperature(60, \"celsius\")\n", "problem.setPressure(100, \"bar\")\n", "problem.add(\"H2O\", 1, \"kg\")\n", "problem.add(\"NaCl\", 1, \"mol\")\n", "problem.add(\"CO2\", 1, \"mol\")" ] }, { "cell_type": "markdown", "id": "e5bde67e", "metadata": {}, "source": [ "> **Note:** To ensure that the equilibrium calculation performed in the next step ignores the kinetic species in the\n", "> system so that we maintain a disequilibrium state between equilibrium and kinetic species, **it is important!** not to\n", "> forget to call\n", "> [EquilibriumProblem::setPartition](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumProblem.html#a53a9496c9d4ffc72a85903146b390e44)\n", "> method. Ignoring this step will produce an initial condition for the amounts of equilibrium and kinetic species\n", "> that correspond to a complete equilibrium state in the system so that no kinetics calculation makes sense afterwards.\n", "\n", "### Calculating the initial amounts of the equilibrium species\n", "\n", "We now use the convenient [equilibrate](https://reaktoro.org/cpp/namespaceReaktoro.html#af2d3b39d3e0b8f9cb5a4d9bbb06b697e)\n", "function to calculate the amounts of the equilibrium species by minimizing\n", "the Gibbs energy of the equilibrium partition only, and not of the entire system. The result is stored in the\n", "object `state0` of class [ChemicalState](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html),\n", "a computational representation of the state of a multiphase chemical system defined by its temperature (*T*),\n", "pressure (*P*), and vector of species amounts (*n*). We then output this chemical state to a file.\n", "\n", "We have now to prescribe the initial amounts of the kinetic species (i.e., the carbonate minerals). This is done\n", "below by setting the initial mass of calcite to 100 g and of dolomite to 50 g. The initial amount of magnesite is\n", "zero." ] }, { "cell_type": "code", "execution_count": null, "id": "6c31be25", "metadata": {}, "outputs": [], "source": [ "state0 = equilibrate(problem)\n", "state0.output('demo-kineticpath-carbonates-co2-before-kinetics.txt')\n", "\n", "state0.setSpeciesMass(\"Calcite\", 100, \"g\")\n", "state0.setSpeciesMass(\"Dolomite\", 50, \"g\")" ] }, { "cell_type": "markdown", "id": "9344f211", "metadata": {}, "source": [ "### Setting the kinetic path calculations\n", "\n", "To be able to simulate the kinetic process of mineral dissolution/precipitation, we introduce the\n", "instance of the class [KineticPath](https://reaktoro.org/cpp/classReaktoro_1_1KineticPath.html) that enables this\n", "functionality. The instance of kinetic path solver is provided with the partition to the equilibrium, kinetic,\n", "and inert species defined above." ] }, { "cell_type": "code", "execution_count": null, "id": "101500d8", "metadata": {}, "outputs": [], "source": [ "path = KineticPath(reactions)\n", "path.setPartition(partition)" ] }, { "cell_type": "markdown", "id": "0270f3bc", "metadata": {}, "source": [ "> **Note**: For repeated chemical kinetics calculations (e.g., in a reactive transport simulation, where kinetics is\n", "performed for each mesh cell/node), consider using the class\n", "[KineticSolver](https://reaktoro.org/cpp/classReaktoro_1_1KineticSolver.html) instead for avoiding some overhead of\n", "class [KineticPath](https://reaktoro.org/cpp/classReaktoro_1_1KineticPath.html). For equilibrium calculations,\n", "consider also the class [EquilibriumSolver](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumSolver.html),\n", "instead of method [equilibrate](https://reaktoro.org/cpp/namespaceReaktoro.html#af2d3b39d3e0b8f9cb5a4d9bbb06b697e),\n", "for similar reasons.\n", "\n", "To analyze the result of kinetic simulations, we save the evolution of different properties of the chemical system\n", "into file `result.txt`." ] }, { "cell_type": "code", "execution_count": null, "id": "35a62daf", "metadata": {}, "outputs": [], "source": [ "output = path.output()\n", "output.filename(\"results.txt\")\n", "output.add(\"time(units=minute)\")\n", "output.add(\"pH\")\n", "output.add(\"speciesMolality(Ca++ units=mmolal)\", \"Ca++ [mmolal]\")\n", "output.add(\"speciesMolality(Mg++ units=mmolal)\", \"Mg++ [mmolal]\")\n", "output.add(\"phaseMass(Calcite units=g)\", \"Calcite [units=g]\")\n", "output.add(\"phaseMass(Dolomite units=g)\", \"Dolomite [units=g]\")" ] }, { "cell_type": "markdown", "id": "bb1f5dd4", "metadata": {}, "source": [ "> **Note**: A list of all possible quantities that can be plotted is shown in the class\n", "[ChemicalQuantity](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalQuantity.html), which provides an interface\n", "for convenient ways of their retrieval.\n", "\n", "To perform the calculation of the kinetic path, we need to provide time interval, on which this path must be\n", "recovered. In this case, we choose 25 hours of simulations. The initial chemical state provided by the instance\n", "`state0`, as well as the time interval with units, in which time is measured, are given with a call of the function\n", "`solve`." ] }, { "cell_type": "code", "execution_count": null, "id": "0b5a8b01", "metadata": {}, "outputs": [], "source": [ "t0, t1 = 0.0, 25.0\n", "path.solve(state0, t0, t1, \"hours\")" ] }, { "cell_type": "markdown", "id": "08d43bd8", "metadata": {}, "source": [ "### Plotting the results of equilibrium path calculation\n", "\n", "To load results from the outputfile, we use [loadtxt](https://docs.scipy.org/doc/numpy/reference/generated/numpy.loadtxt.html)\n", "function provided by the *numpy* package:" ] }, { "cell_type": "code", "execution_count": null, "id": "33d8ade2", "metadata": {}, "outputs": [], "source": [ "import numpy\n", "\n", "filearray = numpy.loadtxt(\"results.txt\", skiprows=1) # load data from the file skipping the one row\n", "data = filearray.T # transpose the matrix with data\n", "[time_indx, ph_indx, ca_indx, mg_indx, calcite_indx, dolomite_indx] \\\n", " = numpy.arange(0, 6) # assign indices of the corresponding data" ] }, { "cell_type": "markdown", "id": "c08e1ae9", "metadata": {}, "source": [ "To visually analyze the obtained reaction path, we export\n", "[bokeh](https://docs.bokeh.org/en/latest/docs/gallery.html#standalone-examples) python plotting package." ] }, { "cell_type": "code", "execution_count": null, "id": "ff606493", "metadata": { "lines_to_next_cell": 1 }, "outputs": [], "source": [ "from bokeh.plotting import figure, show\n", "from bokeh.io import output_notebook\n", "output_notebook()" ] }, { "cell_type": "markdown", "id": "3b21b18a", "metadata": {}, "source": [ "We define a custom function that would generate figure of certain size (in this case, 600 by 300) with label `time`\n", "on the x-axis:" ] }, { "cell_type": "code", "execution_count": null, "id": "f3c81186", "metadata": {}, "outputs": [], "source": [ "def custom_figure(title, y_axis_label):\n", " return figure(plot_width=600, plot_height=300,\n", " title=title,\n", " x_axis_label='time',\n", " y_axis_label=y_axis_label)" ] }, { "cell_type": "markdown", "id": "3f9b4250", "metadata": {}, "source": [ "The plots below depict different chemical properties (x-axis) with respect to the time interval of the kinetic\n", "simulation (y-axis).\n", "In the first plot, we see the growth in the molality of Ca element. At the same time, the molality of Mg is first\n", "increasing and later decreasing. This coincides with the behavior of the mass of minerals on the plots below." ] }, { "cell_type": "code", "execution_count": null, "id": "5eda8fff", "metadata": {}, "outputs": [], "source": [ "time = data[time_indx, :] # fetch time from the data matrix\n", "fig1 = custom_figure(title=\"Ca++ and Mg++ molality w.r.t. time\", y_axis_label=\"Amount of Ca++ and Mg++ [mmolal]\")\n", "fig1.line(time, data[ca_indx], line_width=4, legend_label=\"Ca++\", color=\"orange\")\n", "fig1.line(time, data[mg_indx], line_width=4, legend_label=\"Mg++\", color=\"green\")\n", "show(fig1)" ] }, { "cell_type": "markdown", "id": "64322f36", "metadata": {}, "source": [ "We see the dissolution of the calcite along the time, which is aligned with the earlier plot, where the amount of Ca\n", "element monotonically grows (getting released)." ] }, { "cell_type": "code", "execution_count": null, "id": "b9d2f11f", "metadata": {}, "outputs": [], "source": [ "fig2 = custom_figure(title=\"Calcite dissolution w.r.t. time\", y_axis_label=\"Mass of Calcite [g]\")\n", "fig2.line(time, data[calcite_indx], line_width=4, legend_label=\"Calcite\")\n", "show(fig2)" ] }, { "cell_type": "markdown", "id": "1f661c2d", "metadata": {}, "source": [ "Similar correspondence can be seen from the dolomite plot with respect to time. We see below that amount of\n", "dolomite first decreases and then grows due to the initial dissolution and later precipitation of it in the\n", "chemical system." ] }, { "cell_type": "code", "execution_count": null, "id": "455c36a1", "metadata": {}, "outputs": [], "source": [ "fig3 = custom_figure(title=\"Dolomite behaviour w.r.t. time\", y_axis_label=\"Mass of Dolomite [g]\")\n", "fig3.line(time, data[dolomite_indx], line_width=4, legend_label=\"Dolomite\", color=\"coral\")\n", "show(fig3)" ] }, { "cell_type": "markdown", "id": "4bd241c6", "metadata": {}, "source": [ "As this happens, the chemical system becomes less acidic with pH growing with time:" ] }, { "cell_type": "code", "execution_count": null, "id": "410c20a4", "metadata": { "lines_to_next_cell": 2 }, "outputs": [], "source": [ "fig4 = custom_figure(title=\"pH behaviour w.r.t. time\", y_axis_label=\"pH [-]\")\n", "fig4.line(time, data[ph_indx], line_width=4, color=\"darkviolet\")\n", "show(fig4)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.4" } }, "nbformat": 4, "nbformat_minor": 5 }