tutorial/rt.scaling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1c84d408",
   "metadata": {},
   "source": [
    "# One-dimensional reactive transport modeling of scaling (without oil)\n",
    "\n",
    "This tutorial demonstrates sequential reactive transport calculations of the barite\n",
    "scaling resulting from the waterflooding of the oil reservoirs.\n",
    "\n",
    "We start by importing Python packages (including **reaktoro**) to enable all the necessary numerical calculations\n",
    "as well as the analysis and visualization of the obtained results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef069a0b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from reaktoro import *\n",
    "from math import *\n",
    "import numpy as np\n",
    "from tqdm.notebook import tqdm\n",
    "import os\n",
    "import pandas as pd\n",
    "\n",
    "# Import components of bokeh library\n",
    "from bokeh.io import show, output_notebook\n",
    "from bokeh.plotting import figure\n",
    "from bokeh.models import Range1d, ColumnDataSource\n",
    "from bokeh.layouts import gridplot"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0e3572b7",
   "metadata": {},
   "source": [
    "## Initializing auxiliary time-related constants, discretization, and physical parameters\n",
    "\n",
    "We start with the initialization of the auxiliary time-related constants from seconds up to hours used in the rest of\n",
    "the code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a45e0f2a",
   "metadata": {},
   "outputs": [],
   "source": [
    "second = 1\n",
    "minute = 60\n",
    "hour = 60 * minute"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f5501fdd",
   "metadata": {},
   "source": [
    "Next, we define reactive transport and numerical discretization parameters. In particular, we specify the rock size\n",
    "by setting coordinates of its left and right boundaries to 0.0 m and 25.0 m, respectively. The discretization\n",
    "parameters, i.e., the number of cells and steps in time, are set to 243 and 900, respectively. Considering that\n",
    "the time-step `dt` is fixed to one hour, the simulation lasts 900 hours. Out of these 900 hours, the first\n",
    "45 hours the completion brine (CB) is injected, whereas the rest of the time seawater (SW) is pumped.\n",
    "The reactive transport modeling procedure assumes a constant fluid velocity of 0.8e-5 · 10<sup>-5</sup> m/s and\n",
    "the zero diffusion coefficient for all fluid species. Temperature and pressure are set to 60 &deg;C and 1 atm,\n",
    "respectively, throughout the tutorial. The porosity of the rock is assumed to be 10%. Finally, we use 1kg of water\n",
    "in all the mixtures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ba13327",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Discretization parameters\n",
    "xl = 0.0                # x-coordinate of the left boundary\n",
    "xr = 25.0               # x-coordinate of the right boundary\n",
    "ncells = 243            # number of cells in the discretization\n",
    "dx = (xr - xl) / ncells # length of the mesh cells (in units of m)\n",
    "dt = 1 * hour           # time step\n",
    "\n",
    "t_cb = 45 * hour                # duration of the completion brine (CB) injection\n",
    "t_sw = 855 * hour               # duration of the seawater (SW) injection\n",
    "nsteps_cb = 45                  # number of time steps of the completion brine (CB) injection\n",
    "nsteps_sw = 855                 # number of time steps of the seawater (SW) injection\n",
    "nsteps = nsteps_cb + nsteps_sw  # the total number of steps in the reactive transport simulation\n",
    "\n",
    "# Physical parameters\n",
    "D = 0               # diffusion coefficient (in units of m2/s)\n",
    "v = 0.8e-5          # fluid pore velocity (in units of m/s)\n",
    "T = 60.0            # temperature (in units of celsius)\n",
    "P = 200             # pressure (in units of atm)\n",
    "phi = 0.1           # the porosity\n",
    "water_kg = 1        # amount of water in the initial and injected chemical states"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87df189a",
   "metadata": {},
   "source": [
    "Next, we generate the mesh nodes (array `xcells`) by equally dividing the interval *[xr, xl]* with the number of cells\n",
    "`ncells`. The length between each consecutive mesh nodes (mesh-size) is computed and stored in `dx`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c34646bd",
   "metadata": {},
   "outputs": [],
   "source": [
    "xcells = np.linspace(xl, xr, ncells + 1)  # interval [xl, xr] split into ncells"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "50120958",
   "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 fixed fluid's composition on the left boundary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b94c93e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "dirichlet = False  # parameter that determines whether Dirichlet BC must be used"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aefcc4b1",
   "metadata": {},
   "source": [
    "To ensure that the applied finite-volume scheme is stable, we need to keep track of the Courant–Friedrichs–Lewy (CFL)\n",
    "number, which should be less than 1.0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "012732fd",
   "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": "fe67157d",
   "metadata": {},
   "source": [
    "## Specifying the quantities and properties to be outputted\n",
    "\n",
    "Before running the reactive transport simulations, we provide a list of parameters we are interested in outputting.\n",
    "In this case, it is `pH`, molality of `H+`, `Cl-`, `SO4--`, `Ba++`, `Ca++`, `Sr++`, `Na+`, as well as the\n",
    "concentration, the phase amount, and the volume of barite (BaSO<sub>4</sub>) mineral."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afdcc8a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "output_quantities = \"\"\"\n",
    "    pH\n",
    "    speciesMolality(H+)\n",
    "    speciesMolality(Cl-)\n",
    "    speciesMolality(SO4--)\n",
    "    speciesMolality(Ba++)\n",
    "    speciesMolality(Ca++)\n",
    "    speciesMolality(Sr++)\n",
    "    speciesMolality(Na+)\n",
    "    speciesMolality(Barite)\n",
    "    phaseAmount(Barite)\n",
    "    phaseVolume(Barite)\n",
    "\"\"\".split()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "83ae9de4",
   "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": "e16be17b",
   "metadata": {},
   "outputs": [],
   "source": [
    "column_quantities = \"\"\"\n",
    "    pH\n",
    "    Hcation\n",
    "    Clanion\n",
    "    SO4anion\n",
    "    Bacation\n",
    "    Cacation\n",
    "    Srcation\n",
    "    Nacation\n",
    "    Barite\n",
    "    Barite_phase_amount\n",
    "    Barite_phase_volume\n",
    "    Barite_saturation_index\n",
    "    Barite_stability_index\n",
    "\"\"\".split()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17d1f75e",
   "metadata": {},
   "source": [
    "> **Note**: All the properties mentioned above (except barite's saturation/stability index) are available using\n",
    "the class [ChemicalOutput](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalOutput.html). To retrieve properties\n",
    "that are not included in `ChemicalOutput` functionality, we add several code lines to the `outputstate_df()` function."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6a0f01b9",
   "metadata": {},
   "source": [
    "Create the list of columns stored in the dataframe structure and initialize the `DataFrame` instance with it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ac0ff82",
   "metadata": {},
   "outputs": [],
   "source": [
    "columns = ['step', 'x'] + column_quantities\n",
    "df = pd.DataFrame(columns=columns)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0298d0c0",
   "metadata": {},
   "source": [
    "Finally, we create required folders for outputting the obtained results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4e3b913",
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "folder_results = 'results-rt-scaling'\n",
    "def make_results_folders():\n",
    "    os.system('mkdir -p ' + folder_results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff35e850",
   "metadata": {},
   "source": [
    "## Performing the reactive transport simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e8a14e5b",
   "metadata": {},
   "source": [
    "Subsections below correspond to the methods responsible for each of the functional parts of the reactive transport\n",
    "simulation performed in the main part of the tutorial.\n",
    "\n",
    "### Construction of the chemical system with its phases, species, and barite reactions\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. Here,\n",
    "[supcrt07.xml](https://github.com/reaktoro/reaktoro/blob/master/databases/supcrt/supcrt07.xml) database file is used.\n",
    "In addition to that, we initialize parameters in the Debye-Huckel activity model used for aqueous mixtures. Method\n",
    "`setPHREEQC` allows setting parameters *&#229;* and *b* of the ionic species according to those used in PHREEQC v3.\n",
    "\n",
    "To specify how the chemical system should be defined, i.e., defining all phases and the chemical species they\n",
    "contain, one must use [ChemicalEditor](https://reaktoro.org/cpp/classReaktoro_1_1ChemicalEditor.html) class (see the\n",
    "method `define_chemical_system()` below). Here, we specify two phases, an *aqueous* and a *mineral*. The aqueous phase\n",
    "is defined by providing the list of elements (so that all possible combinations of these elements are considered when\n",
    "creating a set of chemical species). Function `setChemicalModelDebyeHuckel()` helps to set the chemical model of\n",
    "the phase with the Debye-Huckel equation of state, providing specific parameters `dhModel` defined earlier.\n",
    "The mineral phase corresponds to mineral barite (BaSO<sub>4</sub>)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b26d829",
   "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",
    "    # Initialize parameters in the Debye-Huckel activity model and set the corresponding parameters\n",
    "    dhModel = DebyeHuckelParams()\n",
    "    dhModel.setPHREEQC()\n",
    "\n",
    "    # Define the phases of the chemical system\n",
    "    editor = ChemicalEditor(db)\n",
    "    editor.addAqueousPhaseWithElements(\"H Cl S O Ba Ca Sr Na K Mg C Si\").\\\n",
    "        setChemicalModelDebyeHuckel(dhModel)\n",
    "    editor.addMineralPhase('Barite')\n",
    "\n",
    "    # Define the chemical system via the chemical editor\n",
    "    system = ChemicalSystem(editor)\n",
    "\n",
    "    return system"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba22b6ca",
   "metadata": {},
   "source": [
    "To evaluate the saturation indices of barite later, we define the corresponding reaction with function below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f0e0a1d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "def define_chemical_reactions():\n",
    "\n",
    "    # Define barite reaction equation\n",
    "    req_barite = ReactionEquation(\"Barite = SO4-- + Ba++\")\n",
    "\n",
    "    # Define barite chemical reaction based on the reaction equation and chemical system\n",
    "    reaction_barite = Reaction(req_barite, system)\n",
    "\n",
    "    # Create a list of reactions and return them\n",
    "    return [reaction_barite]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e03c272",
   "metadata": {},
   "source": [
    "### Initial condition (IC) of the reactive transport problem\n",
    "\n",
    "The chemical state corresponding to the **initial condition** of the reactive transport simulation is defined in\n",
    "the function `define_initial_condition_fw()`. The composition of the formation water (FW) is taken from the manuscript\n",
    "of Bethke 2008, i.e., Table 30.1 (Miller analysis). We recite it in the table below:\n",
    "\n",
    "| Aqueous species | Amount (mg / kg) |\n",
    "|-----------------|------------------|\n",
    "| Na<sup>+</sup>  | 27250            |\n",
    "| K<sup>+</sup>   | 1730             |\n",
    "| Mg<sup>2+</sup> | 110              |\n",
    "| Ca<sup>2+</sup> | 995              |\n",
    "| Sr<sup>2+</sup> | 105              |\n",
    "| Ba<sup>2+</sup> | 995              |\n",
    "| Cl<sup>-</sup>  | 45150            |\n",
    "| HCO<sub>3</sub><sup>-</sup> | 1980 |\n",
    "| SO<sub>4</sub><sup>2-</sup> | 10 · 10<sup>-3</sup>|\n",
    "\n",
    "The main characteristics of the formation water are high concentration of the Ba<sup>2+</sup> and low concentrations\n",
    "of the SO<sub>4</sub><sup>2-</sup>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "02e1995c",
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "def define_initial_condition_fw(system):\n",
    "\n",
    "    # Define composition of the initial chemical problem\n",
    "    problem_ic = EquilibriumInverseProblem(system)\n",
    "    problem_ic.setTemperature(T, \"celsius\")\n",
    "    problem_ic.setPressure(P, \"atm\")\n",
    "    problem_ic.add(\"H2O\", water_kg, \"kg\")\n",
    "    problem_ic.add(\"SO4\", 10 * water_kg, \"ug\")\n",
    "    problem_ic.add(\"Ca\", 995 * water_kg, \"mg\")\n",
    "    problem_ic.add(\"Ba\", 995 * water_kg, \"mg\")\n",
    "    problem_ic.add(\"Sr\", 105 * water_kg, \"mg\")\n",
    "    problem_ic.add(\"Na\", 27250 * water_kg, \"mg\")\n",
    "    problem_ic.add(\"K\", 1730 * water_kg, \"mg\")\n",
    "    problem_ic.add(\"Mg\", 110 * water_kg, \"mg\")\n",
    "    problem_ic.add(\"Cl\", 45150 * water_kg, \"mg\")\n",
    "    problem_ic.add(\"HCO3\", 1980 * water_kg, \"mg\")\n",
    "    problem_ic.pH(7.0, \"HCl\", \"NaOH\")\n",
    "\n",
    "    # Calculate the equilibrium states for the initial conditions\n",
    "    state_ic = equilibrate(problem_ic)\n",
    "\n",
    "    # Scale the volumes of the phases in the initial condition\n",
    "    state_ic.scalePhaseVolume('Aqueous', 0.1, '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(FW) = \", pH.val)\n",
    "\n",
    "    return state_ic"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4b11bbf",
   "metadata": {},
   "source": [
    "### Boundary condition (BC) of the reactive transport problem"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40382f55",
   "metadata": {},
   "source": [
    "Next, we define the **boundary condition** of the constructed chemical system applied the first 45 hours of\n",
    "simulations. The 7-mol sodium chloride brine below defines completion brine (CB)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fcaeebea",
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "def define_boundary_condition_cb(system):\n",
    "\n",
    "    # Define the boundary condition of the reactive transport modeling problem corresponding to completion brine\n",
    "    problem_bc = EquilibriumProblem(system)\n",
    "    problem_bc.setTemperature(T, \"celsius\")\n",
    "    problem_bc.setPressure(P, \"atm\")\n",
    "    problem_bc.add(\"H2O\", water_kg, \"kg\")\n",
    "    problem_bc.add(\"NaCl\", 7, \"mol\")\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(CB) = \", pH.val)\n",
    "\n",
    "    return state_bc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7723352d",
   "metadata": {},
   "source": [
    "For the rest of the simulation time, we inject the seawater (SW). Similarly, the seawater (SW) composition\n",
    "is taken from Bethke 2008, i.e., Table 30.1 (Seawater). Below, we list the quantities of the aqueous species:\n",
    "\n",
    "| Aqueous species | Amount (mg / kg) |\n",
    "|-----------------|------------------|\n",
    "| Na<sup>+</sup>  | 10760            |\n",
    "| K<sup>+</sup>   | 399              |\n",
    "| Mg<sup>2+</sup> | 1290             |\n",
    "| Ca<sup>2+</sup> | 411              |\n",
    "| Sr<sup>2+</sup> | 8                |\n",
    "| Ba<sup>2+</sup> | 0.01             |\n",
    "| Cl<sup>-</sup>  | 19350            |\n",
    "| HCO<sub>3</sub><sup>-</sup> | 142  |\n",
    "| SO<sub>4</sub><sup>2-</sup> | 2710 |\n",
    "\n",
    "Normally, seawater is rich in sulfate, poor in Ca<sup>2+</sup>, and nearly depleted in Sr<sup>2+</sup> and\n",
    "Ba<sup>2+</sup>. The pH of the seawater is fixed to 8.1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2e18423",
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "def define_boundary_condition_sw(system):\n",
    "\n",
    "    # Define the boundary condition of the reactive transport modeling problem corresponding to seawater\n",
    "    problem_bc = EquilibriumInverseProblem(system)\n",
    "    problem_bc.setTemperature(T, \"celsius\")\n",
    "    problem_bc.setPressure(P, \"atm\")\n",
    "    problem_bc.add(\"H2O\", water_kg, \"kg\")\n",
    "    problem_bc.add(\"SO4--\", 2710 * water_kg, \"mg\")\n",
    "    problem_bc.add(\"Ca++\", 411 * water_kg, \"mg\")\n",
    "    problem_bc.add(\"Ba++\", 0.01 * water_kg, \"mg\")\n",
    "    problem_bc.add(\"Sr++\", 8 * water_kg, \"mg\")\n",
    "    problem_bc.add(\"Na+\", 10760 * water_kg, \"mg\")\n",
    "    problem_bc.add(\"K+\", 399 * water_kg, \"mg\")\n",
    "    problem_bc.add(\"Mg++\", 1290 * water_kg, \"mg\")\n",
    "    problem_bc.add(\"Cl-\", 19350 * water_kg, \"mg\")\n",
    "    problem_bc.add(\"HCO3-\", 142 * water_kg, \"mg\")\n",
    "    problem_bc.pH(8.1, \"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",
    "    props = state_bc.properties()\n",
    "    evaluate_pH = ChemicalProperty.pH(system)\n",
    "    pH = evaluate_pH(props)\n",
    "    print(\"ph(SW) = \", pH.val)\n",
    "\n",
    "    return state_bc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "092e3ef2",
   "metadata": {},
   "source": [
    "Alternative (but very similar to the previous) definition of the SW chemical state can be found from the official\n",
    "PHREEQC tutorials, i.e., [Example 1](https://wwwbrr.cr.usgs.gov/projects/GWC_coupled/phreeqc/html/final-70.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c697fb71",
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "def define_boundary_condition_sw_phreeqc(system):\n",
    "\n",
    "    # Define the boundary condition of the reactive transport modeling problem corresponding to seawater\n",
    "    problem_bc = EquilibriumInverseProblem(system)\n",
    "    problem_bc.setTemperature(T, \"celsius\")\n",
    "    problem_bc.setPressure(P, \"atm\")\n",
    "    problem_bc.add(\"H2O\", water_kg, \"kg\")\n",
    "    problem_bc.add(\"Ca++\", 412.3 * water_kg, \"mg\")  # 412.3 mg / kg = 0.4123 kg / kg => 0.4123 * 58 = 23.9134\n",
    "    problem_bc.add(\"Mg++\", 1291.8 * water_kg, \"mg\")  # 1291.8 mg / kg = 1.2918 kg / kg => 1.2918 * 58 = 74.9244\n",
    "    problem_bc.add(\"Na+\", 10768.0 * water_kg, \"mg\")  # 10768.0  mg / kg = 10.768 kg / kg => 10.768 * 58 = 624.544\n",
    "    problem_bc.add(\"K+\", 399.1 * water_kg, \"mg\")  # 399.1 mg / kg = 0.3991 kg / kg => 0.3991 * 58 = 23.1478\n",
    "    problem_bc.add(\"Si\", 4.28 * water_kg, \"mg\")  # 4.28 mg / kg = 0.00428 kg / kg => 0.00428 * 58 = 0.24824\n",
    "    problem_bc.add(\"Cl-\", 19353 * water_kg, \"mg\")  # 19353.0 mg / kg = 19.353 kg / kg => 19.353 * 58 = 1122.474\n",
    "    problem_bc.add(\"HCO3-\", 141.682 * water_kg, \"mg\")  # 141.682 mg / kg = 0.142682 kg / kg => 0.142682 * 58 = 8.275556\n",
    "    problem_bc.add(\"SO4--\", 2712 * water_kg, \"mg\") # 2712.0 mg / kg = 2.712 kg / kg => 2.712 * 58 = 157.296\n",
    "    problem_bc.pH(8.22, \"HCl\")\n",
    "    problem_bc.pE(8.451)\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",
    "    props = state_bc.properties()\n",
    "    evaluate_pH = ChemicalProperty.pH(system)\n",
    "    pH = evaluate_pH(props)\n",
    "    print(\"ph(SW, PHREEQC) = \", pH.val)\n",
    "\n",
    "    return state_bc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d0735e68",
   "metadata": {},
   "source": [
    "### Indices of partitioning fluid and solid species\n",
    "\n",
    "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 get the indices of the\n",
    "fluid and solid species, to separate mobile from immobile ones."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "607917b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "def partition_indices(system):\n",
    "\n",
    "    # Get number of elements\n",
    "    nelems = system.numElements()\n",
    "\n",
    "    # Get indices of species in fluids and solid partitions\n",
    "    ifluid_species = system.indicesFluidSpecies()\n",
    "    isolid_species = system.indicesSolidSpecies()\n",
    "\n",
    "    return nelems, ifluid_species, isolid_species"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcef5b57",
   "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",
    "We define the arrays `b`, `bfluid`, `bsolid`, storing the concentrations (mol/m<sup>3</sup>) of each element in the\n",
    "system, and those in the fluid partition and in the solid partition at every time step.\n",
    "\n",
    "The array `b` is initialized with the concentrations of the elements at the initial chemical state, `state_ic`, using\n",
    "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 the\n",
    "concentrations of each element on the boundary in mol/m<sup>3</sup><sub>fluid</sub> and is obtained similarly from\n",
    "the `state_bc`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4f604df1",
   "metadata": {},
   "outputs": [],
   "source": [
    "def partition_elements_in_mesh_cell(ncells, nelems, state_ic, state_bc_cb, state_bc_sw):\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, while injecting completion brine\n",
    "    b_bc_cb = state_bc_cb.elementAmounts()\n",
    "\n",
    "    # Initialize the concentrations (mol/m3) of each element on the boundary, while injecting completion brine\n",
    "    b_bc_sw = state_bc_sw.elementAmounts()\n",
    "\n",
    "    return b, bfluid, bsolid, b_bc_cb, b_bc_sw"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eee0eb84",
   "metadata": {},
   "source": [
    "### Reactive transport cycle\n",
    "\n",
    "#### Transport\n",
    "\n",
    "This step updates in 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 of elements in the fluid partition. This is 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/m<sup>3</sup><sub>bulk</sub> and the\n",
    "imposed concentration `b_bc[j]` mol/m<sup>3</sup><sub>fluid</sub>, we need to scale it by the porosity `phi_bc`\n",
    "on the boundary cell m<sup>3</sup><sub>fluid</sub>/m<sup>3</sup><sub>bulk</sub>. 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",
    "m<sup>3</sup><sub>fluid</sub> and total volume m<sup>3</sup><sub>bulk</sub> in the inflow boundary cell.\n",
    "\n",
    "Finally, the updated amounts of elements in the fluid partition are summed with the amounts of elements in the\n",
    "solid partition (remained constant during the transport step), and thus updating the amounts of elements in the\n",
    "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": "011f9afa",
   "metadata": {},
   "outputs": [],
   "source": [
    "def transport(states, bfluid, bsolid, b, b_bc, nelems, ifluid_species, isolid_species):\n",
    "\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": "8da01c71",
   "metadata": {},
   "source": [
    "#### Transport calculation with the finite-volume scheme\n",
    "\n",
    "The function `transport()` expects a conservative property (argument `u`) (e.g., the concentration mol/m<sup>3</sup>\n",
    "of the *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`.\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 linear equation system is solved by the tridiagonal matrix algorithm, also known\n",
    "as the [Thomas algorithm](https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f681f721",
   "metadata": {},
   "outputs": [],
   "source": [
    "def transport_fullimplicit(u, dt, dx, v, D, ul):\n",
    "\n",
    "    # Number of DOFs\n",
    "    n = len(u)\n",
    "    # Fetch the coefficients bearing the diffusion and advection terms\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": "d57868d2",
   "metadata": {},
   "source": [
    "#### Solving the system of equations obtained from finite volume discretization\n",
    "\n",
    "The tridiagonal matrix equation is solved using the\n",
    "[Thomas algorithm](https://en.wikipedia.org/wiki/Tridiagonal_matrix_algorithm)\n",
    "(or the TriDiagonal Matrix Algorithm (TDMA))."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7c09be2",
   "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": "fd1f0a1d",
   "metadata": {},
   "source": [
    "#### Reactive chemistry\n",
    "\n",
    "The chemical equilibrium calculations performed in each mesh cell, using the *Gibbs energy minimization* algorithm\n",
    "(provided by the class [EquilibriumSolver](https://reaktoro.org/cpp/classReaktoro_1_1EquilibriumSolver.html)).\n",
    "**Note:** Before providing temperature and pressure to the `solve()` method, we need to convert Celsius and atmosphere\n",
    "to kelvin and bars, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e88662b1",
   "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",
    "        T_kelvin = T + 273.15\n",
    "        P_bar = P * 1.01325 * 1e5\n",
    "        solver.solve(states[icell], T_kelvin, P_bar, b[icell])\n",
    "    return states"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "60316c06",
   "metadata": {},
   "source": [
    "### Saving and analyzing of the results\n",
    "\n",
    "Function `outputstate_df` is the auxiliary function to add data to the `DataFrame` instance at each time step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "175d727c",
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "def outputstate_df(step, system, reactions, states):\n",
    "\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) * (100 / (1 - phi) if \"phaseVolume\" in quantity_name else 1)\n",
    "\n",
    "        # Fetch barite's stability index\n",
    "        phase_SI = state.phaseStabilityIndices()\n",
    "        barite_phase_index = system.indexPhase(\"Barite\")\n",
    "        values[-1] = phase_SI[barite_phase_index]\n",
    "\n",
    "        # Calculate barite's saturation index\n",
    "        # Fetch chemical properties\n",
    "        props = state.properties()\n",
    "        # Calculate equilibrium constant\n",
    "        lnK_barite = reactions[0].lnEquilibriumConstant(props)\n",
    "        # Calculate reaction quotient\n",
    "        lnQ_barite = reactions[0].lnReactionQuotient(props)\n",
    "        # Calculate saturation ratio\n",
    "        lnSR_barite = lnQ_barite.val - lnK_barite.val\n",
    "        # Calculate saturation index\n",
    "        SI_barite = lnSR_barite / log(10)\n",
    "        values[-2] = SI_barite\n",
    "\n",
    "        # Add values into the dataframe\n",
    "        df.loc[len(df)] = values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e2268a0",
   "metadata": {},
   "source": [
    "### Plotting of the obtained results"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e4db0efe",
   "metadata": {},
   "source": [
    "The library **bokeh** enables the plotting of the results in a Jupyter app. Below, we list auxiliary functions that\n",
    "we use in plotting. Function `titlestr` returns a string for the title of a figure in the  format Time: #h##m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3d3f2b1f",
   "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": "339e193a",
   "metadata": {},
   "source": [
    "Routines `plot_figures_ph()`, `plot_figures_barite_phase_amount()`, `plot_figures_barite_concentration()`,\n",
    "`plot_figures_barite_stability_index()`, and 'plot_figures_aqueous_species()' drawing the plots with\n",
    "chemical states or properties on the selected steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef386311",
   "metadata": {
    "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-0.1, xr+0.1)\n",
    "        p.y_range = Range1d(6.8, 10.0)\n",
    "        p.xaxis.axis_label = 'Distance [m]'\n",
    "        p.yaxis.axis_label = 'pH'\n",
    "        p.legend.location = 'top_right'\n",
    "        p.title.text = titlestr(t)\n",
    "\n",
    "        plots.append([p])\n",
    "\n",
    "    grid = gridplot(plots)\n",
    "    show(grid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3eb972a7",
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "def plot_figures_barite_phase_amount(steps):\n",
    "    plots = []\n",
    "    for i in steps:\n",
    "        print(\"On barite 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='Barite_phase_amount', color='steelblue', line_width=2, legend_label='Barite',\n",
    "               muted_color='steelblue', muted_alpha=0.2, source=source)\n",
    "        p.x_range = Range1d(xl-0.1, xr+0.1)\n",
    "        p.y_range = Range1d(-0.0001, 0.7)\n",
    "        p.xaxis.axis_label = 'Distance [m]'\n",
    "        p.yaxis.axis_label = 'Mineral Phase Amount [mol]'\n",
    "        p.legend.location = 'center_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": "code",
   "execution_count": null,
   "id": "aa05b1fb",
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "def plot_figures_barite_concentration(steps):\n",
    "    plots = []\n",
    "    for i in steps:\n",
    "        print(\"On barite 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='Barite', color='darkorchid', line_width=2, legend_label='Barite',\n",
    "               muted_color='darkorchid', muted_alpha=0.2, source=source)\n",
    "        p.x_range = Range1d(xl-0.1, xr+0.1)\n",
    "        p.y_range = Range1d(-1e-5, 8e-4)\n",
    "        p.xaxis.axis_label = 'Distance [m]'\n",
    "        p.yaxis.axis_label = 'Concentration [mol/m3]'\n",
    "        p.legend.location = 'center_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": "code",
   "execution_count": null,
   "id": "119f04c7",
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "def plot_figures_barite_stability_index(steps):\n",
    "    plots = []\n",
    "    for i in steps:\n",
    "        print(\"On barite's stability index 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='Barite_stability_index', color='indianred', line_width=2, legend_label='Barite',\n",
    "               muted_color='teal', muted_alpha=0.2, source=source)\n",
    "        p.x_range = Range1d(xl-0.1, xr+0.1)\n",
    "        p.xaxis.axis_label = 'Distance [m]'\n",
    "        p.yaxis.axis_label = 'Stability Index [-]'\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": "code",
   "execution_count": null,
   "id": "934e77b3",
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "def plot_figures_barite_saturation_index(steps):\n",
    "    plots = []\n",
    "    for i in steps:\n",
    "        print(\"On barite's saturation index 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='Barite_saturation_index', color='indianred', line_width=2, legend_label='Barite',\n",
    "               muted_color='teal', muted_alpha=0.2, source=source)\n",
    "        p.x_range = Range1d(xl-0.1, xr+0.1)\n",
    "        p.xaxis.axis_label = 'Distance [m]'\n",
    "        p.yaxis.axis_label = 'Saturation Index [-]'\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": "code",
   "execution_count": null,
   "id": "0872cd8c",
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "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=250, y_axis_type = 'log',)\n",
    "        p.line(x='x', y='Hcation', color='darkviolet', line_width=2, legend_label='H+', source=source)\n",
    "        p.line(x='x', y='Clanion', color='darkcyan', line_width=2, legend_label='Cl-', source=source)\n",
    "        p.line(x='x', y='SO4anion', color='darkorange', line_width=2, legend_label='SO4--', source=source)\n",
    "        p.line(x='x', y='Bacation', color='seagreen', line_width=2, legend_label='Ba++', source=source)\n",
    "        p.line(x='x', y='Cacation', color='indianred', line_width=2, legend_label='Ca++', source=source)\n",
    "        p.line(x='x', y='Srcation', color='darkblue', line_width=2, legend_label='Sr++', source=source)\n",
    "        p.x_range = Range1d(xl-0.1, xr+0.1)\n",
    "        p.y_range = Range1d(1e-10, 1e2)\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": "848b6454",
   "metadata": {},
   "source": [
    "# Main part of the tutorial\n",
    "\n",
    "First, we create folders for the result files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bff9e4c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "make_results_folders()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ed7d9f46",
   "metadata": {},
   "source": [
    "Construct the chemical system with its phases and species."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b2208f70",
   "metadata": {},
   "outputs": [],
   "source": [
    "system = define_chemical_system()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13c49b15",
   "metadata": {},
   "source": [
    "Define the reaction for the barite mineral."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d146a54f",
   "metadata": {},
   "outputs": [],
   "source": [
    "reactions = define_chemical_reactions()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61b0764d",
   "metadata": {},
   "source": [
    "Define the initial condition of the reactive transport modeling problem."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c2bc9b85",
   "metadata": {},
   "outputs": [],
   "source": [
    "state_ic = define_initial_condition_fw(system)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fcc455ff",
   "metadata": {},
   "source": [
    "Define the boundary condition of the reactive transport modeling problem composed of two different stages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "481b2361",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the completion brine (CB)\n",
    "state_bc_cb = define_boundary_condition_cb(system)\n",
    "\n",
    "# Define the seawater (SW)\n",
    "state_bc_sw = define_boundary_condition_sw(system)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f48c347",
   "metadata": {},
   "source": [
    "Generate indices of partitioning fluid and solid species, as well as the vectors of elements, including those in\n",
    "fluid and solid partition, and vector of element amounts corresponding to the injected completion brine and seawater."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6748253b",
   "metadata": {},
   "outputs": [],
   "source": [
    "nelems, ifluid_species, isolid_species = partition_indices(system)\n",
    "b, bfluid, bsolid, b_bc_cb, b_bc_sw \\\n",
    "    = partition_elements_in_mesh_cell(ncells, nelems, state_ic, state_bc_cb, state_bc_sw)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2eddc06",
   "metadata": {},
   "source": [
    "Create a list of chemical states for each mesh cell initialized `to state_ic`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c3c1e58",
   "metadata": {},
   "outputs": [],
   "source": [
    "states = [state_ic.clone() for _ in range(ncells + 1)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "62ae6794",
   "metadata": {},
   "source": [
    "Create the equilibrium solver object for the repeated equilibrium calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c99871cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "solver = EquilibriumSolver(system)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6572e528",
   "metadata": {},
   "source": [
    "Running the reactive transport simulation loop. We start with the completion brine injection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2dd9b9ed",
   "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, reactions, states)\n",
    "\n",
    "with tqdm(total=nsteps_cb, desc=\"45 hours of completion brine (CB) injection\") as pbar:\n",
    "    while step < nsteps_cb:\n",
    "        # Perform transport calculations\n",
    "        bfluid, bsolid, b = transport(states, bfluid, bsolid, b, b_bc_cb, nelems, ifluid_species, isolid_species)\n",
    "\n",
    "        # Perform reactive chemical calculations\n",
    "        states = reactive_chemistry(solver, states, b)\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, reactions, states)\n",
    "\n",
    "        # Update a progress bar\n",
    "        pbar.update(1)\n",
    "print(f\"time: {t / hour} hours\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39963a6f",
   "metadata": {},
   "source": [
    "Simulation continues with the seawater injection."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db2e092d",
   "metadata": {},
   "outputs": [],
   "source": [
    "with tqdm(total=nsteps_sw, desc=\"855 hours of seawater (SW) injection\") as pbar:\n",
    "    while step < nsteps_cb + nsteps_sw:\n",
    "        # Perform transport calculations\n",
    "        bfluid, bsolid, b = transport(states, bfluid, bsolid, b, b_bc_sw, nelems, ifluid_species, isolid_species)\n",
    "\n",
    "        # Perform reactive chemical calculations\n",
    "        states = reactive_chemistry(solver, states, b)\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, reactions, states)\n",
    "\n",
    "        # Update a progress bar\n",
    "        pbar.update(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f26166ed",
   "metadata": {},
   "source": [
    "To save the results in csv-format, the line below can be executed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "93356081",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "df.to_csv(folder_results + '/rt.scaling.csv', index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0abe0a7",
   "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": "11db0fca",
   "metadata": {},
   "outputs": [],
   "source": [
    "output_notebook()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d337a175",
   "metadata": {},
   "source": [
    "Select the steps on which results with pH must be output."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1383acb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_steps_to_plot = [45, 60, 300]\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}\"\n",
    "plot_figures_ph(selected_steps_to_plot)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4aab334b",
   "metadata": {},
   "source": [
    "Select the steps on which the rest of the species must be demonstrated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "33194835",
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_steps_to_plot = [120, 300, 600]\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": "99b5ed79",
   "metadata": {},
   "source": [
    "After the execution of functions `plot_figures_barite_concentration()` and `plot_figures_barite_stability_index()`,\n",
    "we see that the mineral starts to precipitate right after we initiate the injection of seawater. It happens due to\n",
    "the mixing of the formation water (FW) and seawater (SW), which have contrasting compositions. In particular, SW\n",
    "is low on Ba<sup>2+</sup> and high on SO<sub>4</sub><sup>2-</sup> concentrations, whereas FW, on the opposite, is\n",
    "high on Ba<sup>2+</sup> and low on SO<sub>4</sub><sup>2-</sup>. During the mixing, both of these ions react creating\n",
    "barite according to the reaction Ba<sup>2+</sup> + SO<sub>4</sub><sup>2-</sup> &#8594; BaSO<sub>4</sub>(s).\n",
    "Such side effect of waterflooding (as part of the oil recovery techniques) reduces the near-wellbore permeability and\n",
    "hampers well productivity/injectivity. Other mineral that have potential to scaling are CaCO<sub>3</sub> (calcite),\n",
    "CaSO<sub>4</sub> (calcium sulfate), FeCO<sub>3</sub> (siderite)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45637a02",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot barite's concentration on the selected steps:\n",
    "plot_figures_barite_concentration(selected_steps_to_plot)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0965180",
   "metadata": {},
   "source": [
    "The saturation index of the mineral is defined as a log10 of the ratio of equilibrium constant and reaction quotient.\n",
    "It is 0 for minerals that are precipitated (i.e., in equilibrium with the solution), SI > 0 for supersaturated\n",
    "minerals, and SI < 0 for undersaturated minerals."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ab962363",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot barite's saturation index on the selected steps:\n",
    "plot_figures_barite_saturation_index(selected_steps_to_plot)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c38c1708",
   "metadata": {},
   "source": [
    "In the plots of the aqueous species below, we see that we have a decrease of Ba<sup>2+</sup> concentrations during\n",
    "mixing, consumed by the barite. The concentration of SO<sub>4</sub><sup>2-</sup> in SW is considerably\n",
    "higher than in FW, which explains the monotonically decreasing curve the left to the right\n",
    "boundary. Finally, we observe a slight increase in the Cl<sup>-</sup> concentration due to the initial NaCl-brine\n",
    "injection for 45 hours."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d35ac059",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot aqueous species on the selected steps:\n",
    "plot_figures_aqueous_species(selected_steps_to_plot)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a90495c",
   "metadata": {},
   "source": [
    "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": "36323d02",
   "metadata": {
    "lines_to_next_cell": 1
   },
   "outputs": [],
   "source": [
    "step = 10\n",
    "def modify_doc(doc):\n",
    "\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-0.1, xr+0.1)\n",
    "    p1.y_range = Range1d(6.8, 10.0)\n",
    "    p1.xaxis.axis_label = 'Distance [m]'\n",
    "    p1.yaxis.axis_label = 'pH'\n",
    "    p1.legend.location = 'bottom_right'\n",
    "    p1.title.text = titlestr(0 * dt)\n",
    "\n",
    "    p2 = figure(plot_width=500, plot_height=250)\n",
    "    p2.line(x='x', y='Barite', color='darkorchid', line_width=2, legend_label='Barite',\n",
    "           muted_color='darkorchid', muted_alpha=0.2, source=source)\n",
    "    p2.x_range = Range1d(xl, xr - 1)\n",
    "    p2.y_range = Range1d(-5e-5, 8e-4)\n",
    "    p2.xaxis.axis_label = 'Distance [m]'\n",
    "    p2.yaxis.axis_label = 'Concentration [mol/m3]'\n",
    "    p2.legend.location = 'center_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='Hcation', color='darkviolet', line_width=2, legend_label='H+', source=source)\n",
    "    p3.line(x='x', y='Clanion', color='darkcyan', line_width=2, legend_label='Cl-', source=source)\n",
    "    p3.line(x='x', y='SO4anion', color='darkorange', line_width=2, legend_label='SO4--', source=source)\n",
    "    p3.line(x='x', y='Bacation', color='seagreen', line_width=2, legend_label='Ba++', source=source)\n",
    "    p3.line(x='x', y='Cacation', color='indianred', line_width=2, legend_label='Ca++', source=source)\n",
    "    p3.line(x='x', y='Srcation', color='darkblue', line_width=2, legend_label='Sr++', source=source)\n",
    "    p3.x_range = Range1d(xl-0.1, xr+0.1)\n",
    "    p3.y_range = Range1d(1e-8, 1e1)\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",
    "    p4 = figure(plot_width=500, plot_height=250)\n",
    "    p4.line(x='x', y='Barite_saturation_index', color='indianred', line_width=2, legend_label='Barite',\n",
    "           muted_color='teal', muted_alpha=0.2, source=source)\n",
    "    p4.x_range = Range1d(xl, xr - 1)\n",
    "    p4.y_range = Range1d(-1e-1, 1 + 1e-1)\n",
    "    p4.xaxis.axis_label = 'Distance [m]'\n",
    "    p4.yaxis.axis_label = 'Saturation Index [-]'\n",
    "    p4.legend.location = 'center_right'\n",
    "    p4.title.text = titlestr(0 * dt)\n",
    "    p4.legend.click_policy = 'mute'\n",
    "\n",
    "    layout =  gridplot([[p1, p2], [p3, p4]])\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",
    "                        Barite_phase_amount=new_source.data['Barite_phase_amount'],\n",
    "                        Hcation=new_source.data['Hcation'],\n",
    "                        SO4anion=new_source.data['SO4anion'],\n",
    "                        Clanion=new_source.data['Clanion'],\n",
    "                        Bacation=new_source.data['Bacation'],\n",
    "                        Cacation=new_source.data['Cacation'],\n",
    "                        Srcation=new_source.data['Srcation'],\n",
    "                        Nacation=new_source.data['Nacation'],\n",
    "                        Barite=new_source.data['Barite'],\n",
    "                        Barite_phase_volume=new_source.data['Barite_phase_volume'],\n",
    "                        Barite_saturation_index=new_source.data['Barite_saturation_index'],\n",
    "                        Barite_stability_index=new_source.data['Barite_stability_index'])\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",
    "        p4.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": "53662976",
   "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 this tutorial is executed in the *localhost*, make sure that the number provided to the variable\n",
    "`notebook_url` below coincides with the number of the localhost of the browser.\n",
    "\n",
    "In the app below, we refresh the reactive time step in a loop, automatically updating 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": "7bb46d84",
   "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
}