Reactive transport of CO2-saturated brine along a porous rock column¶
In this tutorial, we show how Reaktoro can be used for one-dimensional reactive transport calculations for modeling the geochemical reactions that occur along a porous rock column as an aqueous fluid is continuously injected on its left side.
The injected fluid is a brine with 0.9 molal NaCl, 0.05 molal MgCl2, 0.01 molal CaCl2 and almost CO2-saturated, with 0.75 molal of CO2 dissolved.
The porous rock is initially composed of minerals quartz (SiO2) and calcite (CaCO3). The initial porosity is 10 %, and the initial volume percentages of the minerals are 98 %vol of quartz, 2 %vol calcite. The initial conditions for the fluid in the rock is a 0.7 molal NaCl brine in equilibrium with the existing rock minerals calcite and quartz. These initial fluid and rock composition conditions are uniform throughout the rock core. We assume a rock column length of 100 m at temperature 60 °C and 100 bar throughout.
Assumptions
To simplify this tutorial, the following assumptions are made:
Chemical equilibrium is used for modeling the chemical reactions in this problem, not only for reactions between aqueous-aqueous species but also for those between mineral and aqueous species.
A uniform constant velocity field is imposed and it is not updated by solving, for example, Darcy equation.
Both temperature and pressure are also kept constant along the rock.
Todo
Create other more complicated tutorials that do not consider the previous assumptions so that chemical kinetics is considered for mineral dissolution and precipitation reactions and the partial chemical equilibrium assumption is assumed for aqueous and gaseous species. Also, create tutorials in which velocity and pressure fields are physically consistent with each other and computed via the solution of Darcy equation for porous media flow coupled with mass transport equations.
More details about the problem setup can be found in the script below, and also in the step-by-step explanation that follows it.
# Step 1: Import the reaktoro Python package (and other packages)
from reaktoro import *
import os
# Step 2: Initialise auxiliary time-related constants
second = 1
minute = 60
hour = 60 * minute
day = 24 * hour
year = 365 * day
# Step 3: Define parameters for the reactive transport simulation
xl = 0.0 # the x-coordinate of the left boundary
xr = 100.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 (1 m/day in units of m/s)
dt = 0.5*day # the time step (0.5 day in units of s)
T = 60.0 # the temperature (in units of degC)
P = 100 # the pressure (in units of bar)
# Step 4: Construct the chemical system with its phases and species
editor = ChemicalEditor()
editor.addAqueousPhaseWithElementsOf('H2O NaCl CaCl2 MgCl2 CO2')
editor.addMineralPhase('Quartz')
editor.addMineralPhase('Calcite')
editor.addMineralPhase('Dolomite')
# Step 5: Create the ChemicalSystem object using the configured editor
system = ChemicalSystem(editor)
# Step 6: Define the initial condition of the reactive transport modeling problem
problem_ic = EquilibriumProblem(system)
problem_ic.setTemperature(T, 'celsius')
problem_ic.setPressure(P, 'bar')
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: Define the boundary condition of the reactive transport modeling problem
problem_bc = EquilibriumProblem(system)
problem_bc.setTemperature(T, 'celsius')
problem_bc.setPressure(P, 'bar')
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 8: Calculate the equilibrium states for the initial and boundary conditions
state_ic = equilibrate(problem_ic)
state_bc = equilibrate(problem_bc)
# Step 9: Scale the volumes of the phases in the initial condition
state_ic.scalePhaseVolume('Aqueous', 0.1, 'm3') # corresponds to the initial porosity of 10%.
state_ic.scalePhaseVolume('Quartz', 0.882, 'm3')
state_ic.scalePhaseVolume('Calcite', 0.018, 'm3')
# Step 10: Scale the boundary condition state
state_bc.scaleVolume(1.0, 'm3')
# Step 11: Create the mesh for the column
mesh = Mesh(ncells, xl, xr)
# Step 12: Create a chemical field object with every cell having state given by state_ic
field = ChemicalField(mesh.numCells(), state_ic)
# Step 13: Initialize the reactive transport solver
rt = ReactiveTransportSolver(system)
rt.setMesh(mesh)
rt.setVelocity(v)
rt.setDiffusionCoeff(D)
rt.setBoundaryState(state_bc)
rt.setTimeStep(dt)
rt.initialize(field)
# Step 14: Set the output of the reactive transport simulation
output = rt.output()
output.add('pH')
output.add('speciesMolality(H+)')
output.add('speciesMolality(Ca++)')
output.add('speciesMolality(Mg++)')
output.add('speciesMolality(HCO3-)')
output.add('speciesMolality(CO2(aq))')
output.add('phaseVolume(Calcite)')
output.add('phaseVolume(Dolomite)')
output.filename('results/reactive-transport.txt') # Set the name of the output files
os.system('mkdir -p results') # Ensure a results folder exist
# Step 15: Perform given number of reactive tranport steps
t = 0.0 # current time variable
step = 0 # current number of steps
while step <= nsteps: # step until the number of steps are achieved
# Print the progress of the simulation
print("Progress: {}/{} steps, {} min".format(step, nsteps, t/minute))
# Perform one reactive transport time step
rt.step(field)
# Increment time step and number of time steps
t += dt
step += 1
Importing the reaktoro package¶
# Step 1: Import the reaktoro Python package (and other packages)
from reaktoro import *
import os
First, we import the reaktoro Python package so that we can use its classes and methods for performing the chemical reaction calculations.
Defining parameters for the reactive transport simulation¶
Next, we define reactive transport and numerical discretization parameters.
# Step 3: Define parameters for the reactive transport simulation
xl = 0.0 # the x-coordinate of the left boundary
xr = 100.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 (1 m/day in units of m/s)
dt = 0.5*day # the time step (0.5 day in units of s)
T = 60.0 # the temperature (in units of degC)
P = 100 # the pressure (in units of bar)
First, we define the considered rock domain by setting coordinates of its left and right boundaries to 0.0 m and 100.0 m, respectively. The number of cells and steps in time are both set to 100. The reactive transport modeling procedure assumes a constant fluid velocity of v = 1 m/day (1.16 · 10-5 m/s) and the same diffusion coefficient D = 10-9 m2/s for all fluid species. The size of the time-step is set to be half of the day. As highlighted earlier, the temperature and pressure of the fluids are selected to be 60 °C and 100 bar, respectively.
Defining the chemical system¶
We need to define a chemical system that can represent both our fluid and rock. We use class ChemicalEditor below to define a system with an aqueous phase and three mineral phases: quartz, calcite and dolomite. Initially, our rock has no dolomite (CaMg(CO3)2), but since this is a mineral that could potentially precipitate given the fluid composition injected (containing CaCl2 and MgCl2 dissolved), we add it here in the chemical system to ensure that the calculations are able to model dolomite precipitation.
# Step 4: Construct the chemical system with its phases and species
editor = ChemicalEditor()
editor.addAqueousPhaseWithElementsOf('H2O NaCl CaCl2 MgCl2 CO2')
editor.addMineralPhase('Quartz')
editor.addMineralPhase('Calcite')
editor.addMineralPhase('Dolomite')
Note
The aqueous phase is defined above by using a list of compounds, which is then broken automatically by Reaktoro into a list of element names. These element names are then used to find in the database all the aqueous species that could be formed out of them.
Constructing 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.
# Step 5: Create the ChemicalSystem object using the configured editor
system = ChemicalSystem(editor)
Initial condition (IC) for the fluid composition¶
Below, we define the initial condition for the fluid composition in the rock. We want an aqueous fluid that is 0.7 molal of NaCl and in equilibrium with calcite (CaCO3) and quartz (SiO2). To achieve this, we mix 1 kg of H2O, 0.7 mol of NaCl, and plenty of calcite and quartz (10 mol each) to ensure that the aqueous solution is saturated with respect to these minerals.
# Step 6: Define the initial condition of the reactive transport modeling problem
problem_ic = EquilibriumProblem(system)
problem_ic.setTemperature(T, 'celsius')
problem_ic.setPressure(P, 'bar')
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')
Boundary condition (BC) for the fluid composition¶
Next, we define the boundary condition for the fluid composition on the left side of the rock, which should be the one that represents the fluid being continuously injected: 0.9 molal NaCl, 0.05 molal MgCl2, 0.01 molal CaCl2 and almost CO2-saturated, with 0.75 molal of CO2 dissolved. To provide that, we mix 1 kg of H2O with 0.9 mol of NaCl, 0.05 mol of MgCl2, 0.01 mol of CaCl2, and 0.75 mol of CO2.
# Step 7: Define the boundary condition of the reactive transport modeling problem
problem_bc = EquilibriumProblem(system)
problem_bc.setTemperature(T, 'celsius')
problem_bc.setPressure(P, 'bar')
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')
Calculating the IC and BC fluid compositions¶
# Step 8: 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 the equilibrate function to calculate the chemical
equilibrium state of the system with the given initial and boundary equilibrium
conditions stored in the object problem_ic and problem_bc, respectively.
The result is stored in the corresponding instances of the class ChemicalState,
i.e., state_ic and state_bc.
Scaling the phases in the initial condition¶
The initial chemical state state_ic computed before has, at this point,
phases with volumes that do not correspond to our desired porosity of 10 % and
rock mineral composition of 98 %vol of quartz and 2 %vol of calcite.
To obtain this, we scale the volumes of the aqueous and mineral phases as shown below:
# Step 9: Scale the volumes of the phases in the initial condition
state_ic.scalePhaseVolume('Aqueous', 0.1, 'm3') # corresponds to the initial porosity of 10%.
state_ic.scalePhaseVolume('Quartz', 0.882, 'm3')
state_ic.scalePhaseVolume('Calcite', 0.018, 'm3')
Note
After this scaling step, the sum of the phase volumes in state_ic is 1
m3. This also ensures that the amounts of the species in the chemical
system are normalized by m3, and thus they can be regarded as
concentrations in a unit of mol/m3 (bulk volume, not fluid volume!).
Scaling the boundary condition state¶
Next, we scale the boundary condition state to 1 m3, so that we have the
amounts of fluid species in state_bc also normalized by m3.
# Step 10: Scale the boundary condition state
state_bc.scaleVolume(1.0, 'm3')
Note
The chemical state represented by state_bc has no other stable phase
than the aqueous phase (i.e., all mineral phases have zero or negligible
amounts such as 10-21 mol).
Creating the mesh¶
We define the mesh with the class Mesh in order to use in the initialization of class ReactiveTransportSolver later.
# Step 11: Create the mesh for the column
mesh = Mesh(ncells, xl, xr)
Here, we specify the number of cells in the mesh and the x-coordinates of the left and right boundaries (in m).
Creating a chemical field object¶
We have been using class ChemicalState to represent an individual chemical state. We will now use class ChemicalField to represent a collection of chemical states: one for each mesh cell.
# Step 12: Create a chemical field object with every cell having state given by state_ic
field = ChemicalField(mesh.numCells(), state_ic)
Note
Different initial conditions across the mesh cells are possible by assigning
different chemical states to each mesh cell. Here, the same chemical state
in state_ic is used for all cells.
Initializing the reactive transport solver¶
At last, we define the object responsible for solving the reactive transport problem, which is handled by the class ReactiveTransportSolver.
# Step 13: Initialize the reactive transport solver
rt = ReactiveTransportSolver(system)
rt.setMesh(mesh)
rt.setVelocity(v)
rt.setDiffusionCoeff(D)
rt.setBoundaryState(state_bc)
rt.setTimeStep(dt)
rt.initialize(field)
Here, we set the mesh and problem parameters such as velocity, diffusion coefficient, the chemical state representing the boundary condition, and the time step. We also initialize the reactive solver object with the chemical field object specified on the previous step, at this point containing the initial condition for the chemical state of each mesh cell.
Defining the output quantities¶
Before starting the reactive transport calculations, we define the quantities that will be output for every mesh cell, at every time step. For this, we use an object of the class ChemicalOutput:
# Step 14: Set the output of the reactive transport simulation
output = rt.output()
output.add('pH')
output.add('speciesMolality(H+)')
output.add('speciesMolality(Ca++)')
output.add('speciesMolality(Mg++)')
output.add('speciesMolality(HCO3-)')
output.add('speciesMolality(CO2(aq))')
output.add('phaseVolume(Calcite)')
output.add('phaseVolume(Dolomite)')
output.filename('results/reactive-transport.txt') # Set the name of the output files
os.system('mkdir -p results') # Ensure a results folder exist
The name of the output file is to reactive-transport.txt. We specify the
parameters that we are interested in outputting. In this case, it is pH,
molality of H+, Ca2+, Mg2+, HCO3-, CO2 (aq), as well as a phase volume
of calcite and dolomite.
Running the reactive transport simulation¶
As shown below, we perform a sequence of reactive transport calculations, one for each time step, during which the chemical state of each mesh cell is updated. The iterations continue until the maximum number of steps is achieved.
# Step 15: Perform given number of reactive tranport steps
t = 0.0 # current time variable
step = 0 # current number of steps
while step <= nsteps: # step until the number of steps are achieved
# Print the progress of the simulation
print("Progress: {}/{} steps, {} min".format(step, nsteps, t/minute))
# Perform one reactive transport time step
rt.step(field)
# Increment time step and number of time steps
t += dt
step += 1
At each time step, we print the progress of the simulation. We then use the method
step of class ReactiveTransportSolver to perform a single reactive
transport time-stepping. This method also produces a new output file containing
the requested output properties for every mesh cell. In each such file, rows
correspond to cells, whereas the columns correspond to the requested output
properties, i.e., pH, molality of H+, Ca2+, Mg2+, HCO3-, CO2 (aq), as
well as the phase volume of calcite and dolomite.
Todo
Use the result files to general plots and videos, and add them here in the 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!