{
"cells": [
{
"cell_type": "markdown",
"id": "f5f30665",
"metadata": {},
"source": [
"# Reactive transport modeling of the H2S scavenging (with Hematite and Pyrite)\n",
"\n",
"In this tutorial is dedicated to sequential reactive transport calculations of the injected\n",
"H2S-rich brine into a siderite- and hematite-bearing reservoir with subsequent pyrite precipitation.\n",
"All minerals are assumed to be handle by the local equilibrium.\n",
"\n",
"> **Note**: Similar example of H2S scavenging (with dissolving siderite and precipitating pyrrhotite)\n",
"> can also be found in the tutorial\n",
"> [Reactive transport modeling of the H2S scavenging process along a rock core](rt.scavenging.ipynb).\n",
"\n",
"## Importing python packages\n",
"\n",
"First, we need to import a few Python packages to enable numerical calculations and plotting."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a39cf3a2",
"metadata": {},
"outputs": [],
"source": [
"from reaktoro import *\n",
"import numpy as np\n",
"from tqdm.notebook import tqdm\n",
"import os\n",
"\n",
"# Import components of bokeh library\n",
"from bokeh.io import show, output_notebook\n",
"from bokeh.layouts import column\n",
"from bokeh.plotting import figure\n",
"from bokeh.models import Range1d, ColumnDataSource\n",
"from bokeh.layouts import gridplot"
]
},
{
"cell_type": "markdown",
"id": "6d8841f8",
"metadata": {},
"source": [
"In particular, we import the **reaktoro** package to use Reaktoro's classes and methods for performing chemical\n",
"reaction calculations, **numpy** for working with arrays, **tqdm** for the progress bar functionality and **os**,\n",
"to provide a portable way of using operating system dependent functionality. For plotting capabilities,\n",
"we use **bokeh** library.\n",
"\n",
"## Initializing auxiliary time-related constants\n",
"In this step, we initialize auxiliary time-related constants from seconds up to years used in the rest of the code."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1a7ce44d",
"metadata": {},
"outputs": [],
"source": [
"second = 1\n",
"minute = 60\n",
"hour = 60 * minute\n",
"day = 24 * hour\n",
"week = 7 * day\n",
"year = 365 * day"
]
},
{
"cell_type": "markdown",
"id": "937035aa",
"metadata": {},
"source": [
"## Defining parameters for the reactive transport simulation\n",
"\n",
"Next, we define reactive transport and numerical discretization parameters. In particular, we specify the considered\n",
"rock domain by setting coordinates of its left and right boundaries to 0.0 m and 100.0 m, respectively. The\n",
"discretization parameters, i.e., the number of cells and steps in time, are set to 100 and 2500, respectively.\n",
"The reactive transport modeling procedure assumes a constant fluid velocity of 3 · 10-5 m/s and\n",
"the zero diffusion coefficient for all fluid species. The size of the time-step is set to 0.1 days (2.4 hours).\n",
"Temperature and pressure are set to 60 °C (333.15 K) and 200 atm (200 * 1.01325 * 1e5), respectively, throughout\n",
"the tutorial. The porosity of the rock is set to 10%."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9c0cee68",
"metadata": {},
"outputs": [],
"source": [
"# Discretization parameters\n",
"xl = 0.0 # x-coordinate of the left boundary\n",
"xr = 100.0 # x-coordinate of the right boundary\n",
"ncells = 100 # number of cells in the discretization\n",
"nsteps = 1000 # number of steps in the reactive transport simulation\n",
"dx = (xr - xl) / ncells # length of the mesh cells (in units of m)\n",
"dt = 0.1*day # time step\n",
"\n",
"# Physical parameters\n",
"D = 0 # diffusion coefficient (in units of m2/s)\n",
"v = 3e-5 # fluid pore velocity (in units of m/s)\n",
"T = 60.0 + 273.15 # temperature (in units of K)\n",
"P = 200 * 1.01325 * 1e5 # pressure (in units of Pa)\n",
"phi = 0.1 # the porosity"
]
},
{
"cell_type": "markdown",
"id": "f3d3836f",
"metadata": {},
"source": [
"Next, we generate the coordinates of the mesh nodes (array `xcells`) by equally dividing the interval *[xr, xl]* with\n",
"the number of cells `ncells`. The length of the mesh-node is denoted by `dx`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "32352b6d",
"metadata": {},
"outputs": [],
"source": [
"xcells = np.linspace(xl, xr, ncells + 1) # interval [xl, xr] split into ncells"
]
},
{
"cell_type": "markdown",
"id": "ccdf7a14",
"metadata": {},
"source": [
"The boolean variable `dirichlet` is set to `True` or `False`, depending on which boundary condition is considered in\n",
"the numerical calculation. `False` corresponds to imposing the flux of the injected fluid; otherwise, `True` means\n",
"imposing the composition of the fluid on the left boundary."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "394a3f89",
"metadata": {},
"outputs": [],
"source": [
"dirichlet = False # parameter that determines whether Dirichlet BC must be used"
]
},
{
"cell_type": "markdown",
"id": "c8c7a4fb",
"metadata": {},
"source": [
"To make sure that the applied finite-volume scheme is stable, we need to keep track of Courant–Friedrichs–Lewy (CFL)\n",
"number, which should be less than 1.0."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6d94423c",
"metadata": {},
"outputs": [],
"source": [
"CFL = v * dt / dx\n",
"print(\"CFL = \", CFL)\n",
"assert CFL <= 1.0, f\"Make sure that CFL = {CFL} is less that 1.0\""
]
},
{
"cell_type": "markdown",
"id": "eaa4951e",
"metadata": {},
"source": [
"## Specifying the quantities and properties to be outputted\n",
"\n",
"Before running the reactive transport simulations, we specify the list of parameters to be output, i.e.,\n",
"`pH`, molality of `H+`, `HS-`, `S2--`, `CO3--`, `HSO4-`, `H2S(aq)`, `Fe++`, as well as the concentrations\n",
"of siderite, pyrite, and hematite."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1b10b6ba",
"metadata": {},
"outputs": [],
"source": [
"output_quantities = \"\"\"\n",
" pH\n",
" speciesMolality(H+)\n",
" speciesMolality(HS-)\n",
" speciesMolality(S2--)\n",
" speciesMolality(CO3--)\n",
" speciesMolality(HSO4-)\n",
" speciesMolality(H2S(aq))\n",
" speciesMolality(Fe++)\n",
" speciesMolality(Pyrite)\n",
" speciesMolality(Hematite)\n",
"\"\"\".split()"
]
},
{
"cell_type": "markdown",
"id": "84e1a31f",
"metadata": {},
"source": [
"Then, we define the list of names for the DataFrame columns. Note that they must correspond to the order of the\n",
"properties defined in the `output_quantities` list."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e7fda35d",
"metadata": {},
"outputs": [],
"source": [
"column_quantities = \"\"\"\n",
" pH\n",
" Hcation\n",
" HSanion\n",
" S2anion\n",
" CO3anion\n",
" HSO4anion\n",
" H2Saq\n",
" Fe2cation\n",
" pyrite\n",
" hematite\n",
"\"\"\".split()"
]
},
{
"cell_type": "markdown",
"id": "5b1947d9",
"metadata": {},
"source": [
"Create the list of columns and initialize `DataFrame` instance with it."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d69def85",
"metadata": {},
"outputs": [],
"source": [
"columns = ['step', 'x'] + column_quantities\n",
"import pandas as pd\n",
"df = pd.DataFrame(columns=columns)"
]
},
{
"cell_type": "markdown",
"id": "b86d071c",
"metadata": {},
"source": [
"## Organization of the program\n",
"\n",
"The main part of the tutorial consists of three parts documented in the following sections:\n",
"* creation of folders for the results,\n",
"* the main part of the reactive transport simulations, and\n",
"* plotting of the obtained results.\n",
"\n",
"## Creating folders for the outputted results\n",
"\n",
"Using **os** package, we create required folders for outputting the obtained results and for the plot and video\n",
"files later."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e19bd45a",
"metadata": {
"lines_to_next_cell": 1
},
"outputs": [],
"source": [
"folder_results = 'results-rt-scavenging-with-hematite'\n",
"def make_results_folders():\n",
" os.system('mkdir -p ' + folder_results)"
]
},
{
"cell_type": "markdown",
"id": "e2dd99e6",
"metadata": {},
"source": [
"## Performing the reactive transport simulation\n",
"\n",
"The reactive transport simulation is performed below (after all the functions are defined) and consists of several\n",
"stages:\n",
"* initialization of the reactive transport problem and\n",
"* performing the reactive transport simulation along with a defined time interval.\n",
"\n",
"The preparatory initialization step consists of the following sub-steps:\n",
"* definition of the chemical system with its phases and species using `define_chemical_system()`,\n",
"* definition of the initial condition of the reactive transport problem in `define_initial_condition()`,\n",
"* definition of the boundary condition of the reactive transport problem in `define_initial_condition()`,\n",
"* generation of auxiliary indices to partition elements using `partition_indices()` and elements' partitioning\n",
"corresponding to fluid and solid species with function `partition_elements_in_mesh_cell()`, and, finally,\n",
"* definition of the instance of [EquilibriumSolver](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumSolver.html).\n",
"\n",
"The simulation of the reactive transport problem is represented by the loop over a discretized time interval until\n",
"the final time is reached. On each step of this loop, the following functionality of performed:\n",
"* transport calculations using method `transport()`,\n",
"* reactive chemical calculations with `reactive_chemistry()` function, and\n",
"* saving the obtained results using `outputstate()`.\n",
"\n",
"Performing of the transport and reactive chemistry sequentially is possible due to the *operator splitting procedure*.\n",
"First, we update the amount of elements `b` with the transport procedure. Then, `b` are used to evaluate its new\n",
"chemical equilibrium state, producing new amounts of the species in both the fluid and solid phases (available in\n",
"the list `states` of [ChemicalState](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html) objects).\n",
"This chemical reaction equilibrium calculation step, at each mesh cell, permits aqueous species and minerals to react,\n",
"causing mineral dissolution or precipitation, depending on how much the amount of mineral species changes.\n",
"It can then be used, for example, to compute new porosity value for the cell.\n",
"\n",
"Subsections below correspond to the methods responsible for each of the functional parts listed above.\n",
"\n",
"### Construction of the chemical system with its phases and species\n",
"\n",
"To define the chemical system, we need to initialize the class\n",
"[Database](https://reaktoro.org/cpp/classReaktoro_1_1Database.html)\n",
"that provides operations to retrieve physical and thermodynamic data of chemical species. To achieve that, we use\n",
"[supcrt07.xml](https://github.com/reaktoro/reaktoro/blob/master/databases/supcrt/supcrt07.xml) database file.\n",
"\n",
"In addition to the database, we also need to initialize parameters in the Debye-Huckel activity model used for aqueous\n",
"mixtures. Method `setPHREEQC` allows setting parameters *å* and *b* of the ionic species according to those used\n",
"in PHREEQC v3.\n",
"\n",
"Reaktoro is a general-purpose chemical solver that avoids as much as possible presuming specific assumptions about\n",
"considered problems. Thus, one needs to specify how the chemical system of interest, which encompasses the\n",
"specification of all phases in the system as well as the chemical species that compose each phase. Using the\n",
"[ChemicalEditor](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalEditor.html) class, one can conveniently achieve\n",
"this, as shown below in method `define_chemical_system()`.\n",
"\n",
"In this step, we create an object of class [ChemicalEditor](\n",
"https://reaktoro.org/cpp/classReaktoro_1_1ChemicalEditor.html) and specify four phases, an *aqueous* and three\n",
"*mineral* phases that should be considered in the chemical system. The aqueous phase is defined by specifying the\n",
"list of chemical species. Function `setChemicalModelDebyeHuckel()` helps to set the chemical model of the phase\n",
"with the Debye-Huckel equation of state, providing specific parameters `dhModel` defined earlier. The mineral\n",
"phases are defined as two mineral species: siderite (FeCO3), hematite (Fe2O3),\n",
"and pyrite (FeS2).\n",
"\n",
"Finally, we create an object of class [ChemicalSystem](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html)\n",
"using the chemical system definition details stored in the object `editor`.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "61036353",
"metadata": {
"lines_to_next_cell": 1
},
"outputs": [],
"source": [
"def define_chemical_system():\n",
"\n",
" # Construct the chemical system with its phases and species\n",
" db = Database('supcrt07.xml')\n",
"\n",
" # Set PHREEQC parameters to the model\n",
" dhModel = DebyeHuckelParams()\n",
" dhModel.setPHREEQC()\n",
"\n",
" # Define the editor\n",
" editor = ChemicalEditor(db)\n",
" # Select species of the of the aqueous phase by providing the list of elements\n",
" editor.addAqueousPhaseWithElements(\"C Ca Cl Fe H K Mg Na O S\").\\\n",
" setChemicalModelDebyeHuckel(dhModel)\n",
" # Select mineral phases\n",
" editor.addMineralPhase('Pyrite')\n",
" editor.addMineralPhase('Hematite')\n",
" editor.addMineralPhase('Quartz')\n",
"\n",
" # Define chemical system\n",
" system = ChemicalSystem(editor)\n",
"\n",
" return system"
]
},
{
"cell_type": "markdown",
"id": "3c4c0125",
"metadata": {},
"source": [
"### Initial condition (IC) of the reactive transport problem\n",
"\n",
"We have defined and constructed the chemical system of interest, enabling us to move on to the next step in\n",
"Reaktoro's modeling workflow: *defining our chemical reaction problems*. Below, we define its **initial condition**\n",
"with already prescribed equilibrium conditions for *temperature*, *pressure*, and *amounts of elements* that are\n",
"consistent to model reactive transport of injected hydrogen sulfide brine into the rock-fluid composition of siderite\n",
"at 25 °C and 1.01325 bar. The resident fluid in the rock is obtained by the mixture of the aqueous species\n",
"summarized in the following table:\n",
"\n",
"| Aqueous species | Amount (kg) |\n",
"|-----------------|--------------------------------------|\n",
"| H2 | 58.0 |\n",
"| Cl- | 1122.3 · 10-3 |\n",
"| Na+ | 624.08 · 10-3 |\n",
"| SO42- | 157.18 · 10-3 |\n",
"| Mg2+ | 74.820 · 10-3 |\n",
"| Ca2+ | 23.838 · 10-3 |\n",
"| K+ | 23.142 · 10-3 |\n",
"| HCO3- | 8.236 · 10-3 |\n",
"| O2(aq) | 58 · 10-12 |\n",
"\n",
"Both siderite and hematite are added in quantities of 0.5 mol to the reservoir. For equilibrating, the class\n",
"[EquilibriumInverseProblem](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumInverseProblem.html) is used,\n",
"where specific fixed pH and pE are be prescribed to 8.951 and 8.676, respectively."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "484a843a",
"metadata": {
"lines_to_next_cell": 1
},
"outputs": [],
"source": [
"def define_initial_condition(system):\n",
" problem_ic = EquilibriumInverseProblem(system)\n",
" problem_ic.setTemperature(T)\n",
" problem_ic.setPressure(P)\n",
" problem_ic.add(\"H2O\", 58.0, \"kg\")\n",
" problem_ic.add(\"Cl-\", 1122.3e-3, \"kg\")\n",
" problem_ic.add(\"Na+\", 624.08e-3, \"kg\")\n",
" problem_ic.add(\"SO4--\", 157.18e-3, \"kg\")\n",
" problem_ic.add(\"Mg++\", 74.820e-3, \"kg\")\n",
" problem_ic.add(\"Ca++\", 23.838e-3, \"kg\")\n",
" problem_ic.add(\"K+\", 23.142e-3, \"kg\")\n",
" problem_ic.add(\"HCO3-\", 8.236e-3, \"kg\")\n",
" problem_ic.add(\"O2(aq)\", 58e-12, \"kg\")\n",
" problem_ic.add(\"Pyrite\", 0.0, \"mol\")\n",
" problem_ic.add(\"Hematite\", 0.5, \"mol\") # 10 % of all minerals\n",
" problem_ic.add(\"Quartz\", 0.5 * 9, \"mol\") # 90 % of all minerals\n",
" problem_ic.pH(8.951, \"HCl\", \"NaOH\")\n",
"\n",
" # Calculate the equilibrium states for the initial conditions\n",
" state_ic = equilibrate(problem_ic)\n",
" state_ic.scalePhaseVolume('Aqueous', phi, 'm3') # 10% of porosity\n",
" state_ic.scaleVolume(1.0, 'm3')\n",
"\n",
" # Fetch teh value of the ph in the initial chemical state\n",
" props = state_ic.properties()\n",
" evaluate_pH = ChemicalProperty.pH(system)\n",
" pH = evaluate_pH(props)\n",
" print(\"ph(IC) = \", pH.val)\n",
"\n",
" return state_ic"
]
},
{
"cell_type": "markdown",
"id": "221bd5bd",
"metadata": {},
"source": [
"To calculate the system's chemical equilibrium state with the given initial conditions, we use the method\n",
"[equilibrate](https://reaktoro.org/cpp/namespaceReaktoro.html#af2d3b39d3e0b8f9cb5a4d9bbb06b697e), the numerical\n",
"solution of which is written in the objects `state_ic`. It is an instance of the class\n",
"[ChemicalState](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html) that stores the temperature,\n",
"pressure, and the amounts of every species in the system.\n",
"For this calculation, Reaktoro uses an efficient **Gibbs energy minimization** computation to determine the species\n",
"amounts that correspond to a state of minimum Gibbs energy in the system, while satisfying the prescribed amount\n",
"conditions for temperature, pressure, and element amounts. In an inverse equilibrium problem, however, not all\n",
"elements have known molar amounts. Their amount's constraints are replaced by other equilibrium constraints such as\n",
"fixed pH and pE.\n",
"\n",
"The function ends with scaling the volume to 1 m3. Moreover, we specify the 10% porosity of the rock\n",
"by calling `state_ic.scalePhaseVolume('Aqueous', 0.1, 'm3')`.\n",
"\n",
"### Boundary condition (BC) of the reactive transport problem\n",
"\n",
"Next, we define the **boundary condition** of the constructed chemical system with its *temperature*, *pressure*,\n",
"and *amounts of elements*. We prescribe the amount of injected hydrogen sulfide brine, in particular,\n",
"0.0196504 mol of hydrosulfide ion (HS-) and 0.167794 mol of aqueous hydrogen sulfide (H2S(aq)).\n",
"None of the minerals are presented in the injected fluid. Here, the ph is lowered in comparison to the initial\n",
"state to 5.726.\n",
"\n",
"After equilibration, the obtained chemical state representing the boundary condition for the injected fluid\n",
"composition, we scale its volume to 1 m3. It is done so that the amounts of the species in the fluid are\n",
"consistent with a mol/m3 scale."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7c5018ce",
"metadata": {},
"outputs": [],
"source": [
"def define_boundary_condition(system):\n",
"\n",
" # Define the boundary condition of the reactive transport modeling problem\n",
" problem_bc = EquilibriumInverseProblem(system)\n",
" problem_bc.setTemperature(T)\n",
" problem_bc.setPressure(P)\n",
" problem_bc.add(\"H2O\", 58.0, \"kg\")\n",
" problem_bc.add(\"Cl-\", 1122.3e-3, \"kg\")\n",
" problem_bc.add(\"Na+\", 624.08e-3, \"kg\")\n",
" problem_bc.add(\"SO4--\", 157.18e-3, \"kg\")\n",
" problem_bc.add(\"Mg++\", 74.820e-3, \"kg\")\n",
" problem_bc.add(\"Ca++\", 23.838e-3, \"kg\")\n",
" problem_bc.add(\"K+\", 23.142e-3, \"kg\")\n",
" problem_bc.add(\"HCO3-\", 8.236e-3, \"kg\")\n",
" problem_bc.add(\"O2(aq)\", 58e-12, \"kg\")\n",
" problem_bc.add(\"Pyrite\", 0.0, \"mol\")\n",
" problem_bc.add(\"Hematite\", 0.0, \"mol\")\n",
" problem_bc.add(\"HS-\", 0.0196504, \"mol\")\n",
" problem_bc.add(\"H2S(aq)\", 0.167794, \"mol\")\n",
" problem_bc.pH(5.726, \"HCl\", \"NaOH\")\n",
"\n",
" # Calculate the equilibrium states for the boundary conditions\n",
" state_bc = equilibrate(problem_bc)\n",
" # Scale the boundary condition state to 1 m3\n",
" state_bc.scaleVolume(1.0, 'm3')\n",
"\n",
" # Fetch ph of the evaluated chemical state\n",
" props = state_bc.properties()\n",
" evaluate_pH = ChemicalProperty.pH(system)\n",
" pH = evaluate_pH(props)\n",
" print(\"ph(BC) = \", pH.val)\n",
"\n",
" return state_bc"
]
},
{
"cell_type": "markdown",
"id": "8896d0e6",
"metadata": {},
"source": [
"### Indices of partitioning fluid and solid species\n",
"\n",
"Only species in fluid phases are mobile and transported by advection and diffusion mechanisms. The solid phases are\n",
"immobile. The code below identifies the indices of the fluid and solid species. We use methods\n",
"[indicesFluidSpecies](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html#ac2a8b713f46f7a66b2731ba63faa95ad)\n",
"and\n",
"[indicesSolidSpecies](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html#a8b0c237fff1d827f7bf2dbc911fa5bbf)\n",
"of class [ChemicalSystem](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalSystem.html) to store corresponding data\n",
"in the lists `ifluid_species` and `isolid_species`, respectively."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "12cfc848",
"metadata": {},
"outputs": [],
"source": [
"def partition_indices(system):\n",
" nelems = system.numElements()\n",
"\n",
" ifluid_species = system.indicesFluidSpecies()\n",
" isolid_species = system.indicesSolidSpecies()\n",
"\n",
" return nelems, ifluid_species, isolid_species"
]
},
{
"cell_type": "markdown",
"id": "e8f3408d",
"metadata": {},
"source": [
"### Partitioning fluid and solid species\n",
"\n",
"Next, we create arrays to track the amounts of elements in the fluid and solid partition\n",
"(i.e., the amounts of elements among all fluid phases, here only an aqueous phase, and the amounts of elements among\n",
"all solid phases, here the mineral phases). Arrays `b`, `bfluid`, and `bsolid` will store these concentrations\n",
"(mol/m3) at every time step.\n",
"\n",
"The array `b` is initialized with the concentrations of the elements at the initial chemical state, `state_ic`,\n",
"using method\n",
"[elementAmounts](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html#a827457e68a90f89920c13f0cc06fda78)\n",
"of class [ChemicalState](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html). The array `b_bc` stores\n",
"the concentrations of each element on the boundary in mol/m3fluid."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b77cae6e",
"metadata": {},
"outputs": [],
"source": [
"def partition_elements_in_mesh_cell(ncells, nelems, state_ic, state_bc):\n",
"\n",
" # The concentrations of each element in each mesh cell (in the current time step)\n",
" b = np.zeros((ncells, nelems))\n",
"\n",
" # Initialize the concentrations (mol/m3) of the elements in each mesh cell\n",
" b[:] = state_ic.elementAmounts()\n",
"\n",
" # The concentrations (mol/m3) of each element in the fluid partition, in each mesh cell\n",
" bfluid = np.zeros((ncells, nelems))\n",
"\n",
" # The concentrations (mol/m3) of each element in the solid partition, in each mesh cell\n",
" bsolid = np.zeros((ncells, nelems))\n",
"\n",
" # Initialize the concentrations (mol/m3) of each element on the boundary\n",
" b_bc = state_bc.elementAmounts()\n",
"\n",
" return b, bfluid, bsolid, b_bc"
]
},
{
"cell_type": "markdown",
"id": "0690e83e",
"metadata": {},
"source": [
"### Reactive transport cycle\n",
"\n",
"#### Transport\n",
"\n",
"This step updates the fluid partition `bfluid` using the transport equations (without reactions).\n",
"The `transport_fullimplicit()` function below is responsible for solving an advection-diffusion equation, that is\n",
"later applied to transport the concentrations (mol/m3) of elements in the fluid partition (*a\n",
"simplification that is possible because of common diffusion coefficients and velocities of the fluid species,\n",
"otherwise, the transport of individual fluid species would be needed*).\n",
"\n",
"To match the units of concentrations of the elements in the fluid measure in mol/m3bulk and the\n",
"imposed concentration `b_bc[j]` mol/m3fluid, e need to multiply it by the porosity `phi_bc`\n",
"on the boundary cell m3fluid/m3bulk. We use function\n",
"[properties](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html#ad3fa8fd9e1b948da7a698eb020513f3d)\n",
"of the class [ChemicalState](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalState.html) to retrieve fluid volume\n",
"m3fluid and total volume m3bulk in the inflow boundary cell.\n",
"\n",
"The updated amounts of elements in the fluid partition are then summed with the amounts of elements in the solid\n",
"partition `bsolid`, which remained constant during the transport step), and thus updating the amounts of elements\n",
"in the chemical system `b`. Reactive transport calculations involve the solution of a system of\n",
"advection-diffusion-reaction equations."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "58481bc2",
"metadata": {},
"outputs": [],
"source": [
"def transport(states, bfluid, bsolid, b, b_bc, nelems, ifluid_species, isolid_species):\n",
" # Collect the amounts of elements from fluid and solid partitions\n",
" for icell in range(ncells):\n",
" bfluid[icell] = states[icell].elementAmountsInSpecies(ifluid_species)\n",
" bsolid[icell] = states[icell].elementAmountsInSpecies(isolid_species)\n",
"\n",
" # Get the porosity of the boundary cell\n",
" bc_cell = 0\n",
" phi_bc = states[bc_cell].properties().fluidVolume().val / states[bc_cell].properties().volume().val\n",
"\n",
" # Transport each element in the fluid phase\n",
" for j in range(nelems):\n",
" transport_fullimplicit(bfluid[:, j], dt, dx, v, D, phi_bc * b_bc[j])\n",
"\n",
" # Update the amounts of elements in both fluid and solid partitions\n",
" b[:] = bsolid + bfluid\n",
"\n",
" return bfluid, bsolid, b"
]
},
{
"cell_type": "markdown",
"id": "3aef1cf6",
"metadata": {},
"source": [
"##### Transport calculation with finite-volume scheme\n",
"\n",
"The function `transport()` expects a conservative property (argument `u`) (e.g., the concentration mol/m3\n",
"of *j*th element in the fluid given by `bfluid[j]`), the time step (`dt`), the mesh cell length (`dx`),\n",
"the fluid velocity (`v`), the diffusion coefficient (`D`), and the boundary condition of the conservative property\n",
"(`g`) (e.g., the concentration of the *j*th element in the fluid on the left boundary).\n",
"\n",
"The transport equations are solved with a finite volume method, where diffusion and convection are treated implicitly.\n",
"Its discretization in space and time (implicit) results in the constants `alpha` and `beta`. These correspond to\n",
"the diffusion and advection terms in the equation: `D*dt/dx**2` and `v*dt/dx`, respectively.\n",
"\n",
"Arrays `a`, `b`, `c` are the diagonals in the tridiagonal matrix that results by writing all discretized equations\n",
"in a matrix equation. This system of linear equations is solved by the tridiagonal matrix algorithm, also known\n",
"as the Thomas algorithm."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9bb37edd",
"metadata": {},
"outputs": [],
"source": [
"def transport_fullimplicit(u, dt, dx, v, D, ul):\n",
"\n",
" # Number of DOFs\n",
" n = len(u)\n",
" alpha = D * dt / dx ** 2\n",
" beta = v * dt / dx\n",
"\n",
" # Upwind finite volume scheme\n",
" a = np.full(n, -beta - alpha)\n",
" b = np.full(n, 1 + beta + 2 * alpha)\n",
" c = np.full(n, -alpha)\n",
"\n",
" # Set the boundary condition on the left cell\n",
" if dirichlet:\n",
" # Use Dirichlet BC boundary conditions\n",
" b[0] = 1.0\n",
" c[0] = 0.0\n",
" u[0] = ul\n",
"\n",
" else:\n",
" # Flux boundary conditions (implicit scheme for the advection)\n",
" # Left boundary\n",
" b[0] = 1 + alpha + beta\n",
" c[0] = -alpha # stays the same as it is defined -alpha\n",
" u[0] += beta * ul # = dt/dx * v * g, flux that we prescribe is equal v * ul\n",
"\n",
" # Right boundary is free\n",
" a[-1] = - beta\n",
" b[-1] = 1 + beta\n",
"\n",
" # Solve a tridiagonal matrix equation\n",
" thomas(a, b, c, u)"
]
},
{
"cell_type": "markdown",
"id": "7ae51061",
"metadata": {},
"source": [
"##### Solving the system of equations obtained from finite volume discretization\n",
"\n",
"The tridiagonal matrix equation is solved using the Thomas algorithm (or the TriDiagonal Matrix Algorithm (TDMA)).\n",
"It is a simplified form of Gaussian elimination that can be used to solve tridiagonal systems of equations."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b50bfc29",
"metadata": {},
"outputs": [],
"source": [
"def thomas(a, b, c, d):\n",
" n = len(d)\n",
" c[0] /= b[0]\n",
" for i in range(1, n - 1):\n",
" c[i] /= b[i] - a[i] * c[i - 1]\n",
" d[0] /= b[0]\n",
" for i in range(1, n):\n",
" d[i] = (d[i] - a[i] * d[i - 1]) / (b[i] - a[i] * c[i - 1])\n",
" x = d\n",
" for i in reversed(range(0, n - 1)):\n",
" x[i] -= c[i] * x[i + 1]\n",
" return x"
]
},
{
"cell_type": "markdown",
"id": "4f008995",
"metadata": {},
"source": [
"#### Reactive chemistry\n",
"\n",
"The chemical equilibrium calculations performed in each mesh cell, using *Gibbs energy minimization* algorithm (\n",
"provided by the class [EquilibriumSolver](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumSolver.html))."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "245a5467",
"metadata": {},
"outputs": [],
"source": [
"def reactive_chemistry(solver, states, b):\n",
" # Equilibrating all cells with the updated element amounts\n",
" for icell in range(ncells):\n",
" solver.solve(states[icell], T, P, b[icell])\n",
" return states"
]
},
{
"cell_type": "markdown",
"id": "25c1858c",
"metadata": {},
"source": [
"### Results saving and analyzing\n",
"\n",
"Function `outputstate_df` is the auxiliary function to add data to the DataFrame at each time step."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c2c67ed6",
"metadata": {
"lines_to_next_cell": 1
},
"outputs": [],
"source": [
"def outputstate_df(step, system, states):\n",
" # Define the instance of ChemicalQuantity class\n",
" quantity = ChemicalQuantity(system)\n",
"\n",
" # Create the list with empty values to populate with chemical properties\n",
" values = [None] * len(columns)\n",
" for state, x in zip(states, xcells):\n",
"\n",
" # Populate values with number of reactive transport step and spacial coordinates\n",
" values[0] = step\n",
" values[1] = x\n",
"\n",
" # Update the\n",
" quantity.update(state)\n",
" for quantity_name, i in zip(output_quantities, range(2, len(states))):\n",
" values[i] = quantity.value(quantity_name)\n",
" df.loc[len(df)] = values"
]
},
{
"cell_type": "markdown",
"id": "4cfb34bd",
"metadata": {},
"source": [
"# Main parts of the tutorial\n",
"\n",
"First, we create folders for the results."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "20d8b967",
"metadata": {},
"outputs": [],
"source": [
"make_results_folders()"
]
},
{
"cell_type": "markdown",
"id": "618d8f11",
"metadata": {},
"source": [
"Construct the chemical system with its phases and species."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8dc029e7",
"metadata": {},
"outputs": [],
"source": [
"system = define_chemical_system()"
]
},
{
"cell_type": "markdown",
"id": "8e116289",
"metadata": {},
"source": [
"Define the initial condition of the reactive transport modeling problem."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "fb596aeb",
"metadata": {},
"outputs": [],
"source": [
"state_ic = define_initial_condition(system)"
]
},
{
"cell_type": "markdown",
"id": "8533c53c",
"metadata": {},
"source": [
"Define the boundary condition of the reactive transport modeling problem."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "9699512a",
"metadata": {},
"outputs": [],
"source": [
"state_bc = define_boundary_condition(system)"
]
},
{
"cell_type": "markdown",
"id": "9b42ba7f",
"metadata": {},
"source": [
"Generate indices of partitioning fluid and solid species."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "15ce5b6e",
"metadata": {},
"outputs": [],
"source": [
"nelems, ifluid_species, isolid_species = partition_indices(system)"
]
},
{
"cell_type": "markdown",
"id": "312d9cf8",
"metadata": {},
"source": [
"Partitioning fluid and solid species."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "eb40b94e",
"metadata": {},
"outputs": [],
"source": [
"b, bfluid, bsolid, b_bc = partition_elements_in_mesh_cell(ncells, nelems, state_ic, state_bc)"
]
},
{
"cell_type": "markdown",
"id": "368dab64",
"metadata": {},
"source": [
"Create a list of chemical states for the mesh cells (one for each cell, initialized to state_ic)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c202a651",
"metadata": {},
"outputs": [],
"source": [
"states = [state_ic.clone() for _ in range(ncells + 1)]"
]
},
{
"cell_type": "markdown",
"id": "b5509612",
"metadata": {},
"source": [
"Create the equilibrium solver object for the repeated equilibrium calculation."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "354cd0f3",
"metadata": {},
"outputs": [],
"source": [
"solver = EquilibriumSolver(system)"
]
},
{
"cell_type": "markdown",
"id": "c917949c",
"metadata": {},
"source": [
"Running the reactive transport simulation loop."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "7d9ff799",
"metadata": {},
"outputs": [],
"source": [
"step = 0 # the current step number\n",
"t = 0.0 # the current time (in seconds)\n",
"\n",
"# Output the initial state of the reactive transport calculation\n",
"outputstate_df(step, system, states)\n",
"\n",
"with tqdm(total=nsteps, desc=\"Reactive transport simulations\") as pbar:\n",
" while step < nsteps:\n",
" # Perform transport calculations\n",
" bfluid, bsolid, b = transport(states, bfluid, bsolid, b, b_bc, nelems, ifluid_species, isolid_species)\n",
"\n",
" # Perform reactive chemical calculations\n",
" states = reactive_chemistry(solver, states, b)\n",
"\n",
"\n",
" # Increment time step and number of time steps\n",
" t += dt\n",
" step += 1\n",
"\n",
" # Output the current state of the reactive transport calculation\n",
" outputstate_df(step, system, states)\n",
"\n",
" # Update a progress bar\n",
" pbar.update(1)"
]
},
{
"cell_type": "markdown",
"id": "8d8a4195",
"metadata": {},
"source": [
"To inspect the collected data, one can run:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "59a5fe07",
"metadata": {},
"outputs": [],
"source": [
"df"
]
},
{
"cell_type": "markdown",
"id": "71bcdbe2",
"metadata": {},
"source": [
"To save the results in csv-format, please execute:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "b6510c28",
"metadata": {},
"outputs": [],
"source": [
"df.to_csv(folder_results + '/rt.scavenging-with-hematite.csv', index=False)"
]
},
{
"cell_type": "markdown",
"id": "b0b229e1",
"metadata": {},
"source": [
"### Plotting of the obtained results\n",
"\n",
"The last block of the main routine is dedicated to the plotting of the results in a Jupyter app generated by the\n",
"library **bokeh**. It is an interactive visualization library that provides elegant, concise construction of\n",
"versatile graphics, and affords high-performance interactivity over large or streaming datasets.\n",
"\n",
"Below, we list auxiliary functions that we use in plotting. Function `titlestr` returns a string for the title\n",
"of a figure in the format Time: #h##m"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c767a41a",
"metadata": {
"lines_to_next_cell": 1
},
"outputs": [],
"source": [
"def titlestr(t):\n",
" t = t / minute # Convert from seconds to minutes\n",
" h = int(t) / 60 # The number of hours\n",
" m = int(t) % 60 # The number of remaining minutes\n",
" return 'Time: %2dh %2dm' % (h, m)"
]
},
{
"cell_type": "markdown",
"id": "d2ab8333",
"metadata": {},
"source": [
"Routines `plot_figures_ph()`, `plot_figures_minerals()`, and `plot_figures_aqueous_species()' are\n",
"dedicated to producing the plots with chemical properties on the selected steps that are specified by the user below."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "98eccde0",
"metadata": {
"lines_to_end_of_cell_marker": 0,
"lines_to_next_cell": 1
},
"outputs": [],
"source": [
"def plot_figures_ph(steps):\n",
"\n",
" plots = []\n",
" for i in steps:\n",
" print(\"On pH figure at time step: {}\".format(i))\n",
" t = i * dt\n",
" source = ColumnDataSource(df[df['step'] == i])\n",
"\n",
" p = figure(plot_width=600, plot_height=250)\n",
" p.line(x='x', y='pH', color='teal', line_width=2, legend_label='pH', source=source)\n",
" p.x_range = Range1d(xl - 1, xr - 1)\n",
" p.y_range = Range1d(5.0, 10.0)\n",
" p.xaxis.axis_label = 'Distance [m]'\n",
" p.yaxis.axis_label = 'pH'\n",
" p.legend.location = 'center_right'\n",
" p.title.text = titlestr(t)\n",
"\n",
" plots.append([p])\n",
"\n",
" grid = gridplot(plots)\n",
" show(grid)\n",
"\n",
"def plot_figures_minerals(steps):\n",
" plots = []\n",
" for i in steps:\n",
" print(\"On pyrite-hematite-siderite figure at time step: {}\".format(i))\n",
" t = i * dt\n",
" source = ColumnDataSource(df[df['step'] == i])\n",
"\n",
" p = figure(plot_width=600, plot_height=250)\n",
" p.line(x='x', y='pyrite', color='blue', line_width=2, legend_label='Pyrite',\n",
" muted_color='blue', muted_alpha=0.2, source=source)\n",
" p.line(x='x', y='hematite', color='orange', line_width=2, legend_label='Hematite',\n",
" muted_color='orange', muted_alpha=0.2, source=source)\n",
" p.x_range = Range1d(xl - 1, xr - 1)\n",
" p.y_range = Range1d(-1e-3, 0.012)\n",
" p.xaxis.axis_label = 'Distance [m]'\n",
" p.yaxis.axis_label = 'Concentration [mol/m3]'\n",
" p.legend.location = 'top_right'\n",
" p.title.text = titlestr(t)\n",
" p.legend.click_policy = 'mute'\n",
" plots.append([p])\n",
"\n",
" grid = gridplot(plots)\n",
" show(grid)\n",
"\n",
"def plot_figures_pyrite(steps):\n",
" plots = []\n",
" for i in steps:\n",
" print(\"On pyrite figure at time step: {}\".format(i))\n",
" t = i * dt\n",
" source = ColumnDataSource(df[df['step'] == i])\n",
"\n",
" p = figure(plot_width=600, plot_height=250)\n",
" p.line(x='x', y='pyrite', color='blue', line_width=2, legend_label='Pyrite',\n",
" muted_color='blue', muted_alpha=0.2, source=source)\n",
" p.x_range = Range1d(xl - 1, xr + 1)\n",
" p.xaxis.axis_label = 'Distance [m]'\n",
" p.yaxis.axis_label = 'Concentration [mol/m3]'\n",
" p.legend.location = 'top_right'\n",
" p.title.text = titlestr(t)\n",
" p.legend.click_policy = 'mute'\n",
" plots.append([p])\n",
"\n",
" grid = gridplot(plots)\n",
" show(grid)\n",
"\n",
"def plot_figures_hematite(steps):\n",
" plots = []\n",
" for i in steps:\n",
" print(\"On hematite figure at time step: {}\".format(i))\n",
" t = i * dt\n",
" source = ColumnDataSource(df[df['step'] == i])\n",
"\n",
" p = figure(plot_width=600, plot_height=250)\n",
" p.line(x='x', y='hematite', color='orange', line_width=2, legend_label='Hematite',\n",
" muted_color='orange', muted_alpha=0.2, source=source)\n",
"\n",
"\n",
" p.x_range = Range1d(xl - 1, xr + 1)\n",
" p.xaxis.axis_label = 'Distance [m]'\n",
" p.yaxis.axis_label = 'Concentration [mol/m3]'\n",
" p.legend.location = 'top_right'\n",
" p.title.text = titlestr(t)\n",
" p.legend.click_policy = 'mute'\n",
" plots.append([p])\n",
"\n",
" grid = gridplot(plots)\n",
" show(grid)\n",
"\n",
"def plot_figures_aqueous_species(steps):\n",
" plots = []\n",
" for i in steps:\n",
" print(\"On aqueous species figure at time step: {}\".format(i))\n",
" source = ColumnDataSource(df[df['step'] == i])\n",
" t = dt * i\n",
"\n",
" p = figure(plot_width=600, plot_height=300, y_axis_type = 'log',)\n",
" p.line(x='x', y='HSanion', color='darkcyan', line_width=2, legend_label='HS-', source=source)\n",
" p.line(x='x', y='S2anion', color='darkorange', line_width=2, legend_label='S2--', source=source)\n",
" p.line(x='x', y='CO3anion', color='seagreen', line_width=2, legend_label='CO3--', source=source)\n",
" p.line(x='x', y='HSO4anion', color='indianred', line_width=2, legend_label='HSO4-', source=source)\n",
" p.line(x='x', y='H2Saq', color='gray', line_width=2, legend_label='H2S(aq)', source=source)\n",
" p.line(x='x', y='Hcation', color='darkviolet', line_width=2, legend_label='H+', source=source)\n",
" p.line(x='x', y='Fe2cation', color='darkblue', line_width=2, legend_label='Fe++', source=source)\n",
" p.x_range = Range1d(xl - 1, xr - 1)\n",
" p.y_range = Range1d(1e-16, 1e-2)\n",
" p.xaxis.axis_label = 'Distance [m]'\n",
" p.yaxis.axis_label = 'Concentration [molal]'\n",
" p.legend.location = 'top_right'\n",
" p.title.text = titlestr(t)\n",
" p.legend.click_policy = 'mute'\n",
" plots.append([p])\n",
"\n",
" grid = gridplot(plots)\n",
" show(grid)"
]
},
{
"cell_type": "markdown",
"id": "89435cf9",
"metadata": {},
"source": [
"Select the steps, on which results must plotted:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2c7121bf",
"metadata": {},
"outputs": [],
"source": [
"selected_steps_to_plot = [120, 480, 960]\n",
"assert all(step <= nsteps for step in selected_steps_to_plot), f\"Make sure that selected steps are less than \" \\\n",
" f\"total amount of steps {nsteps}\""
]
},
{
"cell_type": "markdown",
"id": "a5e2b5bf",
"metadata": {},
"source": [
"Outputting the plots to the notebook requires the call of `output_notebook()` that specifies outputting the plot\n",
"inline in the Jupyter notebook:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a621b978",
"metadata": {},
"outputs": [],
"source": [
"output_notebook()"
]
},
{
"cell_type": "markdown",
"id": "0eb01e8c",
"metadata": {},
"source": [
"Plot ph on the selected steps:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "45b0ec1c",
"metadata": {},
"outputs": [],
"source": [
"plot_figures_ph(selected_steps_to_plot)"
]
},
{
"cell_type": "markdown",
"id": "bc6681d7",
"metadata": {},
"source": [
"Plot pyrite, hematite, and siderite on the selected steps:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "792bd63d",
"metadata": {},
"outputs": [],
"source": [
"plot_figures_minerals(selected_steps_to_plot)"
]
},
{
"cell_type": "markdown",
"id": "d0bc8fac",
"metadata": {},
"source": [
"Plot aqueous species on the selected steps:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "d3c47f62",
"metadata": {},
"outputs": [],
"source": [
"plot_figures_aqueous_species(selected_steps_to_plot)"
]
},
{
"cell_type": "markdown",
"id": "5241710a",
"metadata": {},
"source": [
"We see on the plots above that the main chemical reactions can be divided into three parts.\n",
"On the first step, the iron ions Fe2+ are being released by the siderite (FeCO3)\n",
"reacting with the injected brine. On the second one, some of the free Fe2+ ions are reacting\n",
"to slightly precipitate hematite, which also tend to react with injected water and generate Fe2+. On the third\n",
"stage, we observe all Fe2+ is consumed by the forming pyrite:\n",
"Fe2+ + 2HS- → FeS2 + 2H+.\n",
"\n",
"The minerals' dissolution and precipitation are accompanied by the formation and dilution of aqueous species.\n",
"Both curves representing HS- and H2S(aq) have two points of the sharp decrease.\n",
"The first and second one is where both species are involved in the dissolution of FeCO3 and\n",
"Fe3O2, whereas the second one is where they are being exhausted by the reaction\n",
"with iron ions to form iron sulfide.\n",
"\n",
"To study the time-dependent behavior of the chemical properties, we create a Bokeh application using the function\n",
"`modify_doc(doc)`. It creates Bokeh content and adds it to the app. The speed of streaming of reactive transport\n",
"data can be controlled by the parameter `step` defined below (bigger the step, faster we will run through available\n",
"data set):"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c548d01e",
"metadata": {
"lines_to_next_cell": 1
},
"outputs": [],
"source": [
"step = 25"
]
},
{
"cell_type": "markdown",
"id": "02524b51",
"metadata": {},
"source": [
"The data streaming is looped, i.e., we will return to the initial time step when reaching the end of the reactive\n",
"transport simulations."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a8a059e2",
"metadata": {
"lines_to_next_cell": 1
},
"outputs": [],
"source": [
"def modify_doc(doc):\n",
" # Initialize the data by the initial chemical state\n",
" source = ColumnDataSource(df[df['step'] == 0])\n",
"\n",
" # Auxiliary function that returns a string for the title of a figure in the format Time: #h##m\n",
" def titlestr(t):\n",
" t = t / minute # Convert from seconds to minutes\n",
" h = int(t) / 60 # The number of hours\n",
" m = int(t) % 60 # The number of remaining minutes\n",
" return 'Time: %2dh %2dm' % (h, m)\n",
"\n",
" # Plot for ph\n",
" p1 = figure(plot_width=500, plot_height=250)\n",
" p1.line(x='x', y='pH', color='teal', line_width=2, legend_label='pH', source=source)\n",
" p1.x_range = Range1d(xl - 1, xr - 1)\n",
" p1.y_range = Range1d(5.0, 10.0)\n",
" p1.xaxis.axis_label = 'Distance [m]'\n",
" p1.yaxis.axis_label = 'pH'\n",
" p1.legend.location = 'center_right'\n",
" p1.title.text = titlestr(0 * dt)\n",
"\n",
" # Plot for calcite and dolomite\n",
"\n",
" p2 = figure(plot_width=500, plot_height=250)\n",
" p2.line(x='x', y='pyrite', color='blue', line_width=2, legend_label='Pyrite', muted_color='blue',\n",
" muted_alpha=0.2, source=source)\n",
" p2.line(x='x', y='hematite', color='orange', line_width=2, legend_label='Hematite', muted_color='orange',\n",
" muted_alpha=0.2, source=source)\n",
" p2.x_range = Range1d(xl - 1, xr - 1)\n",
" p2.y_range = Range1d(-1e-3, 0.012)\n",
" p2.xaxis.axis_label = 'Distance [m]'\n",
" p2.yaxis.axis_label = 'Concentrations [mol/m3]'\n",
" p2.legend.location = 'top_right'\n",
" p2.title.text = titlestr(0 * dt)\n",
" p2.legend.click_policy = 'mute'\n",
"\n",
" p3 = figure(plot_width=500, plot_height=250, y_axis_type='log')\n",
" p3.line(x='x', y='HSanion', color='darkcyan', line_width=2, legend_label='HS-', source=source)\n",
" p3.line(x='x', y='S2anion', color='darkorange', line_width=2, legend_label='S2--', source=source)\n",
" p3.line(x='x', y='CO3anion', color='seagreen', line_width=2, legend_label='CO3--', source=source)\n",
" p3.line(x='x', y='HSO4anion', color='indianred', line_width=2, legend_label='HSO4-', source=source)\n",
" p3.line(x='x', y='H2Saq', color='gray', line_width=2, legend_label='H2S(aq)', source=source)\n",
" p3.line(x='x', y='Hcation', color='darkviolet', line_width=2, legend_label='H+', source=source)\n",
" p3.line(x='x', y='Fe2cation', color='darkblue', line_width=2, legend_label='Fe++', source=source)\n",
" p3.x_range = Range1d(xl - 1, xr - 1)\n",
" p3.y_range = Range1d(1e-16, 1e-2)\n",
" p3.xaxis.axis_label = 'Distance [m]'\n",
" p3.yaxis.axis_label = 'Concentration [molal]'\n",
" p3.legend.location = 'top_right'\n",
" p3.title.text = titlestr(0 * dt)\n",
" p3.legend.click_policy = 'mute'\n",
"\n",
" layout = column(p1, p2, p3)\n",
"\n",
" # Function that return the data dictionary with provided index of the file\n",
" def update():\n",
"\n",
" if source.data['step'][0] + 1 <= nsteps:\n",
" step_number = source.data['step'][0] + step\n",
" else:\n",
" step_number = 0\n",
"\n",
" new_source = ColumnDataSource(df[df['step'] == step_number])\n",
" new_data = dict(index=np.linspace(0, ncells, ncells + 1, dtype=int),\n",
" step=new_source.data['step'],\n",
" x=new_source.data['x'],\n",
" pH=new_source.data['pH'],\n",
" pyrite=new_source.data['pyrite'],\n",
" hematite=new_source.data['hematite'],\n",
" HSanion=new_source.data['HSanion'],\n",
" S2anion=new_source.data['S2anion'],\n",
" CO3anion=new_source.data['CO3anion'],\n",
" HSO4anion=new_source.data['HSO4anion'],\n",
" H2Saq=new_source.data['H2Saq'],\n",
" Hcation=new_source.data['Hcation'],\n",
" Fe2cation=new_source.data['Fe2cation'])\n",
"\n",
" p1.title.text = titlestr(step_number * dt)\n",
" p2.title.text = titlestr(step_number * dt)\n",
" p3.title.text = titlestr(step_number * dt)\n",
"\n",
" source.stream(new_data, rollover=ncells+1)\n",
"\n",
" doc.add_periodic_callback(update, 500)\n",
" doc.add_root(layout)"
]
},
{
"cell_type": "markdown",
"id": "904c4773",
"metadata": {},
"source": [
"Outputting the plots to the notebook requires the call of `output_notebook()` that specifies outputting the plot\n",
"inline in the Jupyter notebook. Finally, the function `modify_doc()` must be passed to `show` so that the app defined\n",
"by it is displayed inline.\n",
"\n",
"> **Important:** If you run this tutorial in the *localhost*, make sure that number provided to the variable\n",
"`notebook_url` below coincides with the number of the localhost you have in your browser.\n",
"\n",
"In the app below, we refresh the reactive time step in a loop, which automatically updates the data source for the\n",
"plots for ph, volume phases of calcite and dolomite, and mollalities of aqueous species (in logarithmic scale)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a9751163",
"metadata": {},
"outputs": [],
"source": [
"output_notebook()\n",
"show(modify_doc, notebook_url=\"http://localhost:8888\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
}
},
"nbformat": 4,
"nbformat_minor": 5
}