# Coupling Reaktoro into other reactive transport codes¶

In this tutorial, we show how Reaktoro can be used in other codes for reactive transport modeling. Here, you will find that we have split the mass transport and chemical reaction calculations so that you can clearly see how a dedicated and advanced transport solver could be combined with Reaktoro’s solvers for chemical reaction calculations. A basic transport solver is used here instead for the sake of software coupling demonstration.

The reactive transport problem we consider here is also relatively simple, similar to the one presented in the tutorial Reactive transport of CO2-saturated brine along a porous rock column. Unlike other tutorials, we proceed here first with a step-by-step explanation of the full script found at the end of this tutorial.

## Importing python packages¶

First, we need to import a few Python packages to enable us to perform the numerical calculations and plotting.

```
from reaktoro import *
from numpy import *
import matplotlib.pyplot as plt
from joblib import Parallel, delayed
import os, time
```

We import the **reaktoro** Python package so that we can use its classes and
methods for performing chemical reaction calculations, **numpy** for working
with arrays, **matplotlib** for plotting capabilities, **joblib** for simple
parallel computing, and, finally, **os**, to provide a portable way of using
operating system dependent functionality.

## Defining parameters for the reactive transport simulation¶

Next, we define reactive transport and numerical discretization parameters.

```
xl = 0.0 # the x-coordinate of the left boundary
xr = 1.0 # the x-coordinate of the right boundary
nsteps = 100 # the number of steps in the reactive transport simulation
ncells = 100 # the number of cells in the discretization
D = 1.0e-9 # the diffusion coefficient (in units of m2/s)
v = 1.0/day # the fluid pore velocity (in units of m/s)
dt = 10*minute # the time step (in units of s)
T = 60.0 + 273.15 # the temperature (in units of K)
P = 100 * 1e5 # the pressure (in units of Pa)
dirichlet = False # the parameter that determines whether Dirichlet BC must be used
```

We specify the considered rock domain by setting coordinates of its left and
right boundaries to 0.0 m and 100.0 m, respectively. The discretization parameters,
i.e., the number of cells and steps in time, are both set to 100. The reactive
transport modeling procedure assumes a constant fluid velocity of 1 m/day
(1.16 · 10^{-5} m/s) and the same diffusion coefficient of 10^{-9} m2/s for all
fluid species (without dispersivity). The size of the time-step is set to 10
minutes. Temperature and pressure are set to 60 °C and 100 bar, respectively,
throughout the whole tutorial.
The chemical equilibrium calculations performed in each mesh cell, at every
time step, using *Gibbs energy minimisation* algorithm (provided by the class
EquilibriumSolver).

The boolean variable `dirichlet`

is set to `True`

or `False`

depending on
which boundary condition is considered in the numerical calculation. Set to
`False`

for imposing the flux of the injected fluid, otherwise, set to
`True`

for imposing the composition of the fluid on the left boundary.

## Specifying the quantities and properties to be outputted¶

Before running the reactive transport simulations, we specify the list of
parameters we are interested in outputting. In this case, it is pH, molality of
H^{+}, Ca^{2+}, Mg^{2+}, HCO_{3}^{-}, CO_{2} (aq), as well as a phase volume of calcite
and dolomite.

```
output_quantities = """
pH
speciesMolality(H+)
speciesMolality(Ca++)
speciesMolality(Mg++)
speciesMolality(HCO3-)
speciesMolality(CO2(aq))
phaseVolume(Calcite)
phaseVolume(Dolomite)
""".split()
```

## Organizing the main function in three parts¶

Below is the main function in our script, which shows three parts, each represented by a Python function documented in the next sections:

```
if __name__ == '__main__':
# Create folders for the results
make_results_folders()
# Run the reactive transport simulations
simulate()
# Plotting the result
plot()
```

Here, we first create the required folders for the results, we then run the reactive transport simulation, and, finally, plot the outputted results.

## Creating folders for the outputted results¶

```
def make_results_folders():
os.system('mkdir -p results')
os.system('mkdir -p figures/ph')
os.system('mkdir -p figures/aqueous-species')
os.system('mkdir -p figures/calcite-dolomite')
os.system('mkdir -p videos')
```

Using **os** package, we create required folders for outputting the obtained
results and for the plot and video files later.

## Perform the reactive transport simulation¶

The reactive transport simulation is performed in the function `simulate`

(you may want to have a quick look in the code below and move to a
detailed explanation in the subsequent sections):

```
def simulate():
# Step 7.1: Construct the chemical system with its phases and species
db = Database('supcrt98.xml')
editor = ChemicalEditor(db)
editor.addAqueousPhase([
'H2O(l)',
'H+',
'OH-',
'Na+',
'Cl-',
'Ca++',
'Mg++',
'HCO3-',
'CO2(aq)',
'CO3--']) # aqueous species are individually selected for performance reasons
editor.addMineralPhase('Quartz')
editor.addMineralPhase('Calcite')
editor.addMineralPhase('Dolomite')
# Step 7.2: Create the ChemicalSystem object using the configured editor
system = ChemicalSystem(editor)
# Step 7.3: Define the initial condition of the reactive transport modeling problem
problem_ic = EquilibriumProblem(system)
problem_ic.setTemperature(T)
problem_ic.setPressure(P)
problem_ic.add('H2O', 1.0, 'kg')
problem_ic.add('NaCl', 0.7, 'mol')
problem_ic.add('CaCO3', 10, 'mol')
problem_ic.add('SiO2', 10, 'mol')
# Step 7.4: Define the boundary condition of the reactive transport modeling problem
problem_bc = EquilibriumProblem(system)
problem_bc.setTemperature(T)
problem_bc.setPressure(P)
problem_bc.add('H2O', 1.0, 'kg')
problem_bc.add('NaCl', 0.90, 'mol')
problem_bc.add('MgCl2', 0.05, 'mol')
problem_bc.add('CaCl2', 0.01, 'mol')
problem_bc.add('CO2', 0.75, 'mol')
# Step 7.5: Calculate the equilibrium states for the initial and boundary conditions
state_ic = equilibrate(problem_ic)
state_bc = equilibrate(problem_bc)
# Step 7.6: Scale the volumes of the phases in the initial condition
state_ic.scalePhaseVolume('Aqueous', 0.1, 'm3')
state_ic.scalePhaseVolume('Quartz', 0.882, 'm3')
state_ic.scalePhaseVolume('Calcite', 0.018, 'm3')
# Scale the boundary condition state to 1 m3
state_bc.scaleVolume(1.0, 'm3')
# Step 7.7: Partitioning fluid and solid species
# The number of elements in the chemical system
nelems = system.numElements()
# The indices of the fluid and solid species
ifluid_species = system.indicesFluidSpecies()
isolid_species = system.indicesSolidSpecies()
# The concentrations of each element in each mesh cell (in the current time step)
b = zeros((ncells, nelems))
# The concentrations (mol/m3) of each element in the fluid partition, in each mesh cell
bfluid = zeros((ncells, nelems))
# The concentrations (mol/m3) of each element in the solid partition, in each mesh cell
bsolid = zeros((ncells, nelems))
# Initialize the concentrations (mol/m3) of the elements in each mesh cell
b[:] = state_ic.elementAmounts()
# Initialize the concentrations (mol/m3) of each element on the boundary
b_bc = state_bc.elementAmounts()
# Step 7.8: Create a list of chemical states for the mesh cells
# The list of chemical states, one for each cell, initialized to state_ic
states = [state_ic.clone() for _ in range(ncells)]
# Step 7.9: Create the equilibrium solver object for the repeated equilibrium calculation
# The interval [xl, xr] split into ncells
x = linspace(xl, xr, ncells + 1)
# The length of each mesh cell (in units of m)
dx = (xr - xl)/ncells
# Step 7.10: Create the equilibrium solver object for the repeated equilibrium calculation
solver = EquilibriumSolver(system)
# Step 7.11: The auxiliary function to create an output file each time step
def outputstate():
# Create the instance of ChemicalOutput class
output = ChemicalOutput(system)
# The number of digits in number of steps (e.g., 100 has 3 digits)
ndigits = len(str(nsteps))
# Provide the output file name, which will correspond
output.filename('results/{}.txt'.format(str(step).zfill(ndigits)))
# We define the columns' tags filled with the name of the quantities
# The first column has a tag 'x' (which corresponds to the center coordinates of the cells )
output.add('tag', 'x') # The value of the center coordinates of the cells
# The rest of the columns correspond to the requested properties
for quantity in output_quantities:
output.add(quantity)
# We update the file with states that correspond to the cells' coordinates stored in x
output.open()
for state, tag in zip(states, x):
output.update(state, tag)
output.close()
# Step 7.12: Running the reactive transport simulation loop
step = 0 # the current step number
t = 0.0 # the current time (in seconds)
while step <= nsteps:
# Print the progress of the simulation
print("Progress: {}/{} steps, {} min".format(step, nsteps, t/minute))
# Output the current state of the reactive transport calculation
outputstate()
# Collect the amounts of elements from fluid and solid partitions
for icell in range(ncells):
bfluid[icell] = states[icell].elementAmountsInSpecies(ifluid_species)
bsolid[icell] = states[icell].elementAmountsInSpecies(isolid_species)
# Transport each element in the fluid phase
for j in range(nelems):
transport(bfluid[:, j], dt, dx, v, D, b_bc[j])
# Update the amounts of elements in both fluid and solid partitions
b[:] = bsolid + bfluid
# Equilibrating all cells with the updated element amounts
for icell in range(ncells):
solver.solve(states[icell], T, P, b[icell])
# Increment time step and number of time steps
t += dt
step += 1
print("Finished!")
```

We now explain step-by-step the function `simulate`

. We begin by defining the
chemical system using the ChemicalEditor class:

```
db = Database('supcrt98.xml')
editor = ChemicalEditor(db)
editor.addAqueousPhase([
'H2O(l)',
'H+',
'OH-',
'Na+',
'Cl-',
'Ca++',
'Mg++',
'HCO3-',
'CO2(aq)',
'CO3--']) # aqueous species are individually selected for performance reasons
editor.addMineralPhase('Quartz')
editor.addMineralPhase('Calcite')
editor.addMineralPhase('Dolomite')
```

We specify aqueous and mineral phases that should be considered in the chemical
system. For performance reasons, the aqueous phase is defined by manually
specifying the chemical species, instead of automatic collection of species
from database. There are three pure mineral phase considered: quartz (SiO_{2}),
calcite (CaCO_{3}), and dolomite (CaMg(CO_{3})_{2}).

### Creating the chemical system¶

This step is where we create an object of class ChemicalSystem using the
chemical system definition details stored in the object `editor`

.

```
system = ChemicalSystem(editor)
```

### Initial condition (IC) of the reactive transport problem¶

After constructing the chemical system of interest, we can proceed to the
definition of a chemical equilibrium problem to set the **initial condition**
for the fluid and rock composition that is consistent with our intention of
modeling reactive transport in a porous rock column composed of quartz and
calcite, at 60 °C and 100 bar, with a 0.7 NaCl molal brine in equilibrium
with the rock minerals. To ensure equilibrium with both quartz and calcite, the
equilibrium problem setup below considers a relatively large amount for these
minerals (10 moles each).

```
problem_ic = EquilibriumProblem(system)
problem_ic.setTemperature(T)
problem_ic.setPressure(P)
problem_ic.add('H2O', 1.0, 'kg')
problem_ic.add('NaCl', 0.7, 'mol')
problem_ic.add('CaCO3', 10, 'mol')
problem_ic.add('SiO2', 10, 'mol')
```

We will later scale the volumes of the aqueous and mineral phases so that they
are consistent with a 10 % porosity and the required volume percentages of the
rock minerals (98 %_{vol} of quartz and 2 %_{vol} of calcite).

### Boundary condition (BC) of the reactive transport problem¶

For the boundary condition, we need to specify the composition of the fluid
that is injected into the rock. This is done below, by defining an equilibrium
problem that will later be solved to produce an aqueous fluid with 0.9 molal
NaCl, 0.05 molal MgCl_{2}, 0.01 CaCl_{2}, and 0.75 molal CO_{2}, in a state very
close to CO_{2} saturation.

```
problem_bc = EquilibriumProblem(system)
problem_bc.setTemperature(T)
problem_bc.setPressure(P)
problem_bc.add('H2O', 1.0, 'kg')
problem_bc.add('NaCl', 0.90, 'mol')
problem_bc.add('MgCl2', 0.05, 'mol')
problem_bc.add('CaCl2', 0.01, 'mol')
problem_bc.add('CO2', 0.75, 'mol')
```

The composition of the injected aqueous fluid results from the solution of the
above equilibrium problem, in which 1 kg of H_{2}O is mixed with 0.90 moles of
NaCl, 0.05 moles of MgCl_{2}, 0.01 moles of CaCl_{2}, and 0.75 moles of CO_{2}.

### Calculate the equilibrium states for the initial and boundary conditions¶

```
state_ic = equilibrate(problem_ic)
state_bc = equilibrate(problem_bc)
```

In this step, we use method `equilibrate`

to calculate the chemical
equilibrium state of the system with the given initial and boundary equilibrium
conditions stored in the objects `problem_ic`

and `problem_bc`

. The
numerical solution of each problem results in the objects `state_ic`

and `state_bc`

, respectively, of class ChemicalState, which stores the
temperature, pressure, and the amounts of every species in the system.

### Scaling the volumes of the phases¶

We scale the volumes of the phases in the initial condition according to the
condition that porosity is 10 % (ratio of fluid volume to total volume), and
the volume fractions of minerals are 98 %_{vol} of quartz (SiO_{2}) and 2 %_{vol}
of calcite (ratio of mineral volume to solid volume).

For the chemical state representing the boundary condition for the injected
fluid composition, we scale its volume to 1 m^{3}. This is done so that the
amounts of the species in the fluid is consistent with a mol/m^{3} scale.

```
state_ic.scalePhaseVolume('Aqueous', 0.1, 'm3')
state_ic.scalePhaseVolume('Quartz', 0.882, 'm3')
state_ic.scalePhaseVolume('Calcite', 0.018, 'm3')
# Scale the boundary condition state to 1 m3
state_bc.scaleVolume(1.0, 'm3')
```

### Partitioning fluid and solid species¶

Only species in fluid phases are mobile and transported by advection and diffusion mechanisms. The solid phases are immobile.

The code below identify the indices of the fluid and solid species, and also create arrays to keep track of the amounts of elements in the fluid and solid partition (i.e., the amounts of elements among all fluid phases, here only an aqueous phase, and the amounts of elements among all solid phases, here the mineral phases).

```
# The number of elements in the chemical system
nelems = system.numElements()
# The indices of the fluid and solid species
ifluid_species = system.indicesFluidSpecies()
isolid_species = system.indicesSolidSpecies()
# The concentrations of each element in each mesh cell (in the current time step)
b = zeros((ncells, nelems))
# The concentrations (mol/m3) of each element in the fluid partition, in each mesh cell
bfluid = zeros((ncells, nelems))
# The concentrations (mol/m3) of each element in the solid partition, in each mesh cell
bsolid = zeros((ncells, nelems))
# Initialize the concentrations (mol/m3) of the elements in each mesh cell
b[:] = state_ic.elementAmounts()
# Initialize the concentrations (mol/m3) of each element on the boundary
b_bc = state_bc.elementAmounts()
```

We use methods `indicesFluidSpecies`

and `indicesSolidSpecies`

of class
ChemicalSystem to get the indices of the fluid and solid species, which are
then stored in the lists `ifluid_species`

and `isolid_species`

,
respectively. Then, we define the arrays `b`

, `bfluid`

, `bsolid`

, that
will store, respectively, the concentrations (mol/m^{3}) of each element in the
system, in the fluid partition, and in the solid partition at every time step.

The array `b`

is initialized with the concentrations of the elements at the
initial chemical state, `state_ic`

, using method `elementAmounts`

of class
ChemicalState. The array `b_bc`

stores the concentrations of each element
on the boundary.

### Create a list of chemical states for the mesh cells¶

Every mesh cell needs a ChemicalState object, which is done below with the
creation of a list `states`

of length `ncells`

, initialized with `state_ic`

for each cell (i.e., all cells have initially the same fluid and rock composition).

```
# The list of chemical states, one for each cell, initialized to state_ic
states = [state_ic.clone() for _ in range(ncells)]
```

### Cell coordinates and lengths¶

Next, we generate the coordinates of the mesh nodes (array `x`

) by equally
dividing the interval *[xr, xl]* into `ncells`

. The length between each
consecutive mesh nodes is computed and stored in `dx`

(the length of the mesh
cells).

```
# The interval [xl, xr] split into ncells
x = linspace(xl, xr, ncells + 1)
# The length of each mesh cell (in units of m)
dx = (xr - xl)/ncells
```

### Creating the equilibrium solver¶

For the repeated equilibrium calculation, we define an equilibrium solver object using EquilibriumSolver class.

```
solver = EquilibriumSolver(system)
```

### Define an auxiliary function for creating the output file¶

The auxiliary function `outputstate`

uses the class ChemicalOutput to
output quantities and properties at each mesh cell, at every time step.

```
def outputstate():
# Create the instance of ChemicalOutput class
output = ChemicalOutput(system)
# The number of digits in number of steps (e.g., 100 has 3 digits)
ndigits = len(str(nsteps))
# Provide the output file name, which will correspond
output.filename('results/{}.txt'.format(str(step).zfill(ndigits)))
# We define the columns' tags filled with the name of the quantities
# The first column has a tag 'x' (which corresponds to the center coordinates of the cells )
output.add('tag', 'x') # The value of the center coordinates of the cells
# The rest of the columns correspond to the requested properties
for quantity in output_quantities:
output.add(quantity)
# We update the file with states that correspond to the cells' coordinates stored in x
output.open()
for state, tag in zip(states, x):
output.update(state, tag)
output.close()
```

### Running the reactive transport simulation¶

Before proceeding to the reactive transport calculations, we set the initial time and a counter for the number of time steps.

```
step = 0 # the current step number
t = 0.0 # the current time (in seconds)
while step <= nsteps:
# Print the progress of the simulation
print("Progress: {}/{} steps, {} min".format(step, nsteps, t/minute))
# Output the current state of the reactive transport calculation
outputstate()
# Collect the amounts of elements from fluid and solid partitions
for icell in range(ncells):
bfluid[icell] = states[icell].elementAmountsInSpecies(ifluid_species)
bsolid[icell] = states[icell].elementAmountsInSpecies(isolid_species)
# Transport each element in the fluid phase
for j in range(nelems):
transport(bfluid[:, j], dt, dx, v, D, b_bc[j])
# Update the amounts of elements in both fluid and solid partitions
b[:] = bsolid + bfluid
# Equilibrating all cells with the updated element amounts
for icell in range(ncells):
solver.solve(states[icell], T, P, b[icell])
# Increment time step and number of time steps
t += dt
step += 1
print("Finished!")
```

The reactive transport simulation continues until the maximum number of time
steps is achieved. At each time step, we output the progress of the simulation
using function `outputstate`

. Then, we collect the amounts of elements from
fluid and solid partition. We now apply an *operator splitting procedure*, in
which we first update the amounts of elements in the fluid partition
(`bfluid`

) using the transport equations (without reactions). These updated
amounts of elements in the fluid partition are now summed with the amounts of
elements in the solid partition (`bsolid`

, which remained constant during the
transport step), and thus updating the amounts of elements in the chemical
system (`b`

). We now use these updated amounts of elements in the cell
(`b`

) to evaluate its new chemical equilibrium state, thus producing new
amounts of the species in both the fluid and solid phases (available in the
list `states`

of ChemicalState objects). This chemical reaction
equilibrium calculation step, at each mesh cell, permits, for example,
aqueous species and minerals to react, and thus causes mineral dissolution or
precipitation, depending on how much the amount of mineral species changes. This
can then be used, for example, to compute new porosity value for the cell.

### Solving the transport equations¶

Reactive transport calculations involve the solution of a system of advection-diffusion-reaction equations. Due to the use of operator splitting scheme in this tutorial, we solve the transport equations and then perform the chemical reaction calculations for each cell before proceeding to the next time step.

The `transport`

function below is responsible for solving an
advection-diffusion equation, that is later applied to transport the
concentrations (mol/m^{3}) of elements in the fluid partition (*a simplification
that is possible because of common diffusion coefficients and velocities of the
fluid species, otherwise the transport of individual fluid species would be
needed*).

```
def transport(u, dt, dx, v, D, g):
# Number of DOFs
n = len(u)
alpha = D*dt/dx**2
beta = v*dt/dx
a = full(n, -beta - alpha)
b = full(n, 1 + beta + 2*alpha)
c = full(n, -alpha)
# Set the boundary condition
if dirichlet:
# Use Dirichlet BC boundary conditions
b[0] = 1.0
c[0] = 0.0
u[0] = g
else:
# Use flux boundary conditions
u[0] += beta*g
b[0] = 1 + beta + alpha
b[-1] = 1 + beta + alpha
# Solve a tridiagonal matrix equation
thomas(a, b, c, u)
```

The function `transport`

expects a conservative property (argument `u`

)
(e.g., the concentration (mol/m^{3}) of *j*-th element in the fluid given by
`bfluid[j]`

), the time step (`dt`

), the mesh cell length (`dx`

), the
fluid velocity (`v`

), the diffusion coefficient (`D`

), and the boundary
condition of the conservative property (`g`

) (e.g., the concentration of the
*j*-th element in the fluid on the left boundary).

The transport equations are solved with a finite volume method. Its
discretization in space and time (implicit) results in the constants `alpha`

and `beta`

. These correspond to the diffusion and advection terms in the
equation: `D*dt/dx**2`

and `v*dt/dx`

, respectively.

Arrays `a`

, `b`

, `c`

are the diagonals in the tridiagonal matrix that
results by writing all discretized equations in a matrix equation. This system
of linear equations is solved by the tridiagonal matrix algorithm, also known
as the Thomas algorithm:

```
def thomas(a, b, c, d):
n = len(d)
c[0] /= b[0]
for i in range(1, n - 1):
c[i] /= b[i] - a[i]*c[i - 1]
d[0] /= b[0]
for i in range(1, n):
d[i] = (d[i] - a[i]*d[i - 1])/(b[i] - a[i]*c[i - 1])
x = d
for i in reversed(range(0, n - 1)):
x[i] -= c[i]*x[i + 1]
return x
```

## Plotting of the results¶

The last block of the main routine is plotting of the results.

```
def plot():
# Plot all result files
files = sorted(os.listdir('results'))
Parallel(n_jobs=16)(delayed(plotfile)(file) for file in files)
# Create videos for the figures
ffmpegstr = 'ffmpeg -y -r 30 -i figures/{0}/%03d.png -codec:v mpeg4 -flags:v +qscale -global_quality:v 0 videos/{0}.mp4'
os.system(ffmpegstr.format('calcite-dolomite'))
os.system(ffmpegstr.format('aqueous-species'))
os.system(ffmpegstr.format('ph'))
```

The full script for the reactive transport calculation is shown below.

```
from reaktoro import *
from numpy import *
import matplotlib.pyplot as plt
from joblib import Parallel, delayed
import os, time
# Step 2: Auxiliary time related constants
second = 1
minute = 60
hour = 60 * minute
day = 24 * hour
year = 365 * day
# Step 3: Parameters for the reactive transport simulation
xl = 0.0 # the x-coordinate of the left boundary
xr = 1.0 # the x-coordinate of the right boundary
nsteps = 100 # the number of steps in the reactive transport simulation
ncells = 100 # the number of cells in the discretization
D = 1.0e-9 # the diffusion coefficient (in units of m2/s)
v = 1.0/day # the fluid pore velocity (in units of m/s)
dt = 10*minute # the time step (in units of s)
T = 60.0 + 273.15 # the temperature (in units of K)
P = 100 * 1e5 # the pressure (in units of Pa)
dirichlet = False # the parameter that determines whether Dirichlet BC must be used
# Step 4: The list of quantities to be output for each mesh cell, each time step
output_quantities = """
pH
speciesMolality(H+)
speciesMolality(Ca++)
speciesMolality(Mg++)
speciesMolality(HCO3-)
speciesMolality(CO2(aq))
phaseVolume(Calcite)
phaseVolume(Dolomite)
""".split()
# Step 7: Perform the reactive transport simulation
def simulate():
# Step 7.1: Construct the chemical system with its phases and species
db = Database('supcrt98.xml')
editor = ChemicalEditor(db)
editor.addAqueousPhase([
'H2O(l)',
'H+',
'OH-',
'Na+',
'Cl-',
'Ca++',
'Mg++',
'HCO3-',
'CO2(aq)',
'CO3--']) # aqueous species are individually selected for performance reasons
editor.addMineralPhase('Quartz')
editor.addMineralPhase('Calcite')
editor.addMineralPhase('Dolomite')
# Step 7.2: Create the ChemicalSystem object using the configured editor
system = ChemicalSystem(editor)
# Step 7.3: Define the initial condition of the reactive transport modeling problem
problem_ic = EquilibriumProblem(system)
problem_ic.setTemperature(T)
problem_ic.setPressure(P)
problem_ic.add('H2O', 1.0, 'kg')
problem_ic.add('NaCl', 0.7, 'mol')
problem_ic.add('CaCO3', 10, 'mol')
problem_ic.add('SiO2', 10, 'mol')
# Step 7.4: Define the boundary condition of the reactive transport modeling problem
problem_bc = EquilibriumProblem(system)
problem_bc.setTemperature(T)
problem_bc.setPressure(P)
problem_bc.add('H2O', 1.0, 'kg')
problem_bc.add('NaCl', 0.90, 'mol')
problem_bc.add('MgCl2', 0.05, 'mol')
problem_bc.add('CaCl2', 0.01, 'mol')
problem_bc.add('CO2', 0.75, 'mol')
# Step 7.5: Calculate the equilibrium states for the initial and boundary conditions
state_ic = equilibrate(problem_ic)
state_bc = equilibrate(problem_bc)
# Step 7.6: Scale the volumes of the phases in the initial condition
state_ic.scalePhaseVolume('Aqueous', 0.1, 'm3')
state_ic.scalePhaseVolume('Quartz', 0.882, 'm3')
state_ic.scalePhaseVolume('Calcite', 0.018, 'm3')
# Scale the boundary condition state to 1 m3
state_bc.scaleVolume(1.0, 'm3')
# Step 7.7: Partitioning fluid and solid species
# The number of elements in the chemical system
nelems = system.numElements()
# The indices of the fluid and solid species
ifluid_species = system.indicesFluidSpecies()
isolid_species = system.indicesSolidSpecies()
# The concentrations of each element in each mesh cell (in the current time step)
b = zeros((ncells, nelems))
# The concentrations (mol/m3) of each element in the fluid partition, in each mesh cell
bfluid = zeros((ncells, nelems))
# The concentrations (mol/m3) of each element in the solid partition, in each mesh cell
bsolid = zeros((ncells, nelems))
# Initialize the concentrations (mol/m3) of the elements in each mesh cell
b[:] = state_ic.elementAmounts()
# Initialize the concentrations (mol/m3) of each element on the boundary
b_bc = state_bc.elementAmounts()
# Step 7.8: Create a list of chemical states for the mesh cells
# The list of chemical states, one for each cell, initialized to state_ic
states = [state_ic.clone() for _ in range(ncells)]
# Step 7.9: Create the equilibrium solver object for the repeated equilibrium calculation
# The interval [xl, xr] split into ncells
x = linspace(xl, xr, ncells + 1)
# The length of each mesh cell (in units of m)
dx = (xr - xl)/ncells
# Step 7.10: Create the equilibrium solver object for the repeated equilibrium calculation
solver = EquilibriumSolver(system)
# Step 7.11: The auxiliary function to create an output file each time step
def outputstate():
# Create the instance of ChemicalOutput class
output = ChemicalOutput(system)
# The number of digits in number of steps (e.g., 100 has 3 digits)
ndigits = len(str(nsteps))
# Provide the output file name, which will correspond
output.filename('results/{}.txt'.format(str(step).zfill(ndigits)))
# We define the columns' tags filled with the name of the quantities
# The first column has a tag 'x' (which corresponds to the center coordinates of the cells )
output.add('tag', 'x') # The value of the center coordinates of the cells
# The rest of the columns correspond to the requested properties
for quantity in output_quantities:
output.add(quantity)
# We update the file with states that correspond to the cells' coordinates stored in x
output.open()
for state, tag in zip(states, x):
output.update(state, tag)
output.close()
# Step 7.12: Running the reactive transport simulation loop
step = 0 # the current step number
t = 0.0 # the current time (in seconds)
while step <= nsteps:
# Print the progress of the simulation
print("Progress: {}/{} steps, {} min".format(step, nsteps, t/minute))
# Output the current state of the reactive transport calculation
outputstate()
# Collect the amounts of elements from fluid and solid partitions
for icell in range(ncells):
bfluid[icell] = states[icell].elementAmountsInSpecies(ifluid_species)
bsolid[icell] = states[icell].elementAmountsInSpecies(isolid_species)
# Transport each element in the fluid phase
for j in range(nelems):
transport(bfluid[:, j], dt, dx, v, D, b_bc[j])
# Update the amounts of elements in both fluid and solid partitions
b[:] = bsolid + bfluid
# Equilibrating all cells with the updated element amounts
for icell in range(ncells):
solver.solve(states[icell], T, P, b[icell])
# Increment time step and number of time steps
t += dt
step += 1
print("Finished!")
# Step 10: Return a string for the title of a figure in the format Time: #h##m
def titlestr(t):
t = t / minute # Convert from seconds to minutes
h = int(t) / 60 # The number of hours
m = int(t) % 60 # The number of remaining minutes
return 'Time: {:>3}h{:>2}m'.format(h, str(m).zfill(2))
# Step 9: Generate figures for a result file
def plotfile(file):
step = int(file.split('.')[0])
print('Plotting figure', step, '...')
t = step * dt
filearray = loadtxt('results/' + file, skiprows=1)
data = filearray.T
ndigits = len(str(nsteps))
plt.figure()
plt.xlim(left=-0.02, right=0.52)
plt.ylim(bottom=2.5, top=10.5)
plt.title(titlestr(t))
plt.xlabel('Distance [m]')
plt.ylabel('pH')
plt.plot(data[0], data[1])
plt.tight_layout()
plt.savefig('figures/ph/{}.png'.format(str(step).zfill(ndigits)))
plt.figure()
plt.xlim(left=-0.02, right=0.52)
plt.ylim(bottom=-0.1, top=2.1)
plt.ticklabel_format(style='sci', axis='y', scilimits=(0,0))
plt.title(titlestr(t))
plt.xlabel('Distance [m]')
plt.ylabel('Mineral Volume [%$_{\mathsf{vol}}$]')
plt.plot(data[0], data[7] * 100, label='Calcite')
plt.plot(data[0], data[8] * 100, label='Dolomite')
plt.legend(loc='center right')
plt.tight_layout()
plt.savefig('figures/calcite-dolomite/{}.png'.format(str(step).zfill(ndigits)))
plt.figure()
plt.yscale('log')
plt.xlim(left=-0.02, right=0.52)
plt.ylim(bottom=0.5e-5, top=2)
plt.title(titlestr(t))
plt.xlabel('Distance [m]')
plt.ylabel('Concentration [molal]')
plt.plot(data[0], data[3], label='Ca++')
plt.plot(data[0], data[4], label='Mg++')
plt.plot(data[0], data[5], label='HCO3-')
plt.plot(data[0], data[6], label='CO2(aq)')
plt.plot(data[0], data[2], label='H+')
plt.legend(loc='lower right')
plt.tight_layout()
plt.savefig('figures/aqueous-species/{}.png'.format(str(step).zfill(ndigits)))
plt.close('all')
# Step 8: Plot all result files and generate a video
def plot():
# Plot all result files
files = sorted(os.listdir('results'))
Parallel(n_jobs=16)(delayed(plotfile)(file) for file in files)
# Create videos for the figures
ffmpegstr = 'ffmpeg -y -r 30 -i figures/{0}/%03d.png -codec:v mpeg4 -flags:v +qscale -global_quality:v 0 videos/{0}.mp4'
os.system(ffmpegstr.format('calcite-dolomite'))
os.system(ffmpegstr.format('aqueous-species'))
os.system(ffmpegstr.format('ph'))
# Step 7.10.2: Solve a tridiagonal matrix equation using Thomas algorithm (or TriDiagonal Matrix Algorithm (TDMA))
def thomas(a, b, c, d):
n = len(d)
c[0] /= b[0]
for i in range(1, n - 1):
c[i] /= b[i] - a[i]*c[i - 1]
d[0] /= b[0]
for i in range(1, n):
d[i] = (d[i] - a[i]*d[i - 1])/(b[i] - a[i]*c[i - 1])
x = d
for i in reversed(range(0, n - 1)):
x[i] -= c[i]*x[i + 1]
return x
# Step 7.10.1: Perform a transport step
def transport(u, dt, dx, v, D, g):
# Number of DOFs
n = len(u)
alpha = D*dt/dx**2
beta = v*dt/dx
a = full(n, -beta - alpha)
b = full(n, 1 + beta + 2*alpha)
c = full(n, -alpha)
# Set the boundary condition
if dirichlet:
# Use Dirichlet BC boundary conditions
b[0] = 1.0
c[0] = 0.0
u[0] = g
else:
# Use flux boundary conditions
u[0] += beta*g
b[0] = 1 + beta + alpha
b[-1] = 1 + beta + alpha
# Solve a tridiagonal matrix equation
thomas(a, b, c, u)
# Step 6: Creating folders for the results' output
def make_results_folders():
os.system('mkdir -p results')
os.system('mkdir -p figures/ph')
os.system('mkdir -p figures/aqueous-species')
os.system('mkdir -p figures/calcite-dolomite')
os.system('mkdir -p videos')
# Step 5: Define the main function
if __name__ == '__main__':
# Create folders for the results
make_results_folders()
# Run the reactive transport simulations
simulate()
# Plotting the result
plot()
```

Todo

Add figures and videos here in this tutorial.

## Have you got an issue?¶

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

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