Author manuscript, published in "Computers and Geosciences 37, 4 (2011) 610-620"
DOI : 10.1016/j.cageo.2010.11.020
A workflow for handling heterogeneous 3D models with the TOUGH2 family of codes: Applications to
numerical modeling of CO2 geological storage
Pascal Audigane, Christophe Chiaberge, Frédéric Mathurin, Julie Lions, Géraldine Picot-Colbeaux
Corresponding author email adresse: p.audigane@brgm.fr
BRGM: Water Division,
3 Av. Claude Guillemin, BP 36009
45060 ORLEANS Cedex 2 (FRANCE)
tel:+33(0)2 38 64 33 72
fax:+33(0)2 38 64 3719
Abstract
This paper is addressed to the TOUGH2 user community. It presents a new tool for handling simulations run
hal-00614438, version 1 - 11 Aug 2011
with the TOUGH2 code with specific application to CO2 geological storage. This tool is composed of separate
FORTRAN subroutines (or modules) that can be run independently using input and output files in ASCII format
for TOUGH2. These modules have been developed specifically for modeling of carbon dioxide geological
storage and their use with TOUGH2 and the Equation of State module ECO2N, dedicated to CO2-water-salt
mixture systems, with TOUGHREACT, which is an adaptation of TOUGH2 with ECO2N and geochemical
fluid-rock interactions, and with TOUGH2 and the EOS7C module dedicated to CO2-CH4 gas mixture is
described. The objective is to save time for the pre-processing, execution and visualization of complex geometry
for geological system representation. The workflow is rapid and user-friendly and future implementation to other
TOUGH2 EOS modules for other contexts (e.g. nuclear waste disposal, geothermal production) is
straightforward. Three examples are shown for validation: i) leakage of CO2 up through an abandoned well, ii)
3D reactive transport modeling of CO2 in a sandy aquifer formation in the Sleipner gas Field, (North Sea,
Norway), and iii) an estimation of enhanced gas recovery technology using CO2 as the injected and stored gas to
produce methane in the K12B Gas Field (North Sea, Denmark).
Keywords: TOUGH2, pre and post processing, CO2 geological storage, 3D visualization, multiphase flow,
reactive transport modeling
1. Introduction
The capture and geological storage of carbon dioxide (CCS) has been under consideration as one technical
option for greenhouse gas emission mitigation since the Intergovernmental Panel on Climate Change published
its Special Report on Carbon Dioxide Capture and Storage (CCS) (IPCC, 2005). The elaboration of a legal
1
framework for CCS requires pertinent estimates of storage feasibility, efficiency and safety during and after
injection, and numerical modeling is increasingly being used for a better understanding in the preparation,
approbation and elaboration stages of pilot sites for CCS. Extrapolation of laboratory results from core to
reservoir scale can be done by numerical modeling with correct implementation of the identified relevant
physical processes. Modeling of carbon dioxide geological storage requires at least the implementation of
multiphase flow and phase partitioning (dissolution of supercritical CO2 in the brine and vice versa). The CO2 is
to be injected deep enough (> 800 m) so that it will stabilize in its supercritical state with a high density and
therefore with a better storage efficiency. This injection will have a huge impact on local and regional
hydrogeology and hydrogeochemistry. Since the CO2 will be injected for several years, a large volume of water
will be displaced and the integrity of the storage zone must be assessed in order to avoid any upward migration
hal-00614438, version 1 - 11 Aug 2011
of the supercritical CO2 due to over pressure or buoyancy.
The TOUGH2 family of codes is one of several numerical models dedicated to CCS that is capable of handling
identified trapping mechanisms (Bachu et al., 1994) when CO2 is injected and stored in saline aquifers (Pruess et
al., 1999, Pruess and García 2002, Pruess and Spycher, 2007). Implementation of geochemical fluid-rock
interactions in the code extended its modeling applications to mineral trapping characterization (Xu et al., 2003;
Xu et al., 2006). As a start, the 1D radial demonstration problem (one layer system centered around the injection
well) provides a simplified but complete picture of gas bubble evolution, storage capacity, pressure build up,
dissolution and mineral trapping (Xu et al., 2003, André et al., 2007). The 2D radial model better describes the
gravity impact on pressure build up, CO2 migration, induced convection effects and mineral trapping efficiency
(Audigane et al., 2007; Zhou et al. 2008, Zhang et al., 2009). The weakness of the radial models is that they
consider a flat geological formation, disregarding any dip or regional hydrogeology induced by the geometry of
the system. Simulations of long term effects would require a better representation of layer horizons as well as the
heterogeneity of the cap rock formation overlying the storage zone afforded by 3D models.
3D modeling and meshing of geological formations is a vast subject in itself, and is not the aim of the present
paper. Our study deals solely with the use of TOUGH2 (and associated modules) to develop a 3D model of a
geometrically complex system. The mesh maker (“MESHMAKER”) embedded in TOUGH2 generates numerical
grids with Cartesian or radial geometrical configurations. Although the TOUGH family of codes has great
flexibility in handling the wide variety of grid information required to describe complex objects, designing and
generating a suitable irregular grid can be tedious and error-prone. For this reason, graphical user interfaces
(GUI) have been developed specifically for TOUGH2 (O‟sullivan and Bullivant, 2005; Sato et al., 2003 ;
2
Burnell et al., 2003; Swenson et al., 2003; Alcott et al., 2006, Pan, 2003 and 2008). These tools aim to provide
efficient and user-friendly integrated gridding software with graphic interfaces and easy viewing. Though this
creates a more convenient working environment for the TOUGH2 family of codes, various constraints remain.
First, updating the TOUGH2 module requires a complete recompilation of the tool, making the user highly
dependent on the GUI developers' technical support teams. Second, the option for editing results may be limited
in comparison with the specific finite element mesh visualization tool. Finally, the GUI may have to be
purchased.
Our objective was, therefore, to provide a convenient, user-friendly tool for modeling CCS with TOUGH2 on a
geometrically complex geological system while preserving a maximum flexibility for updating and rapid release.
Our strategy involved developing a set of separate modules (hereafter called a “workflow”) written as an
hal-00614438, version 1 - 11 Aug 2011
executable FORTRAN subroutine aimed at providing specific tasks using only input and output files from a
“standard” TOUGH2 run (Chiaberge and Audigane, 2009). The TOUGH2 code can then be easily modified or
updated in one of the subroutines. The workflow outputs ASCII files in specific formats for two classical FEM
visualization tools: Tecplot (www.tecplot.com), a licensed computational fluid dynamics post-processing tool
and ParaView (www.paraview.org), an open-source, multi-platform application designed to visualize small to
very large data sets.
After briefly describing the TOUGH2 family of codes, focusing on the CCS application, we describe here in
detail the structure of the input and output files needed for a “standard” run, and associated equation-of-state
(EOS) routines. Mesh generation and boundary condition management is also addressed.
Thereafter, we
describe the workflow structure and its connection with TOUGH2 input and output files. Several other tools are
also documented (Local Grid Refinement implementation, import from geo-modeler, streamline visualization).
Three applications are presented. The first involves using the workflow in a benchmarking trial on a CCS model
recently proposed by the University of Stuttgart (Class et al., 2009). The second application is a case study of
reactive transport modeling of CO2 injection in the Utsira sand aquifer at the Sleipner Gas Field in the North Sea
(Picot-Colbeaux et al., 2009). The third concerns a pilot site for enhanced methane recovery by carbon dioxide
injection in the K12B gas field, also in the North Sea (Audigane et al., 2008).
2. General description of TOUGH2
2.1 TOUGH2 and CCS
3
TOUGH2 is a non-isothermal multiphase fluid flow simulator used extensively for modeling heat and mass
transfer in various situations such as geothermal systems, nuclear waste disposal and basin hydrogeology. Recent
developments enable its use in assessing CO2 geological storage. The code‟s versatility is due to the
development of a set of equation-of-state (EOS) modules adapted to various gas- and liquid-mixture behaviors.
An EOS routine called ECO2N has recently been developed specifically for treating CO2-water-salt systems
typically encountered in the case of CO2 geological storage in saline aquifers (Pruess and Spycher, 2007).
TOUGHREACT is an adaptation of TOUGH2 to reactive transport modeling considering geochemical fluidrock interactions during multiphase fluid flow (Xu et al., 2006). TOUGHREACT is designed to run with several
EOS routines, including ECO2N.
hal-00614438, version 1 - 11 Aug 2011
Numerous publications (not all of which are cited here) describe studies of CO2 geological storage done using
TOUGH2/ECO2N and TOUGHREACT/ECO2N codes. One of these, involving a 1D model using radial
geometry centered on the injection well in a single-layer porous geological formation, has been proposed as a
code intercomparison study (Pruess et al., 2004). In a similar fashion, Xu et al. (2003) developed a 1D reactive
transport model of CO2 geological storage in the Gulf Coast sediment using TOUGHREACT. Audigane et al.
(2007) have studied a 2D model of short- and long-term evolution of a CO2 gas bubble and the related
geochemical reactions of the Utsira sandy aquifer in the North Sea. TOUGH2 and TOUGHREACT have also
been used to predict near-well dry out induced by several years of CO2 injection (Pruess and Muller 2009), along
with salt precipitation (André et al., 2007). Oldenburg et al. (2004) developed a specific EOS called EOS7C to
study the injection of CO2 in a depleted methane reservoir. Audigane et al. (2008) have used TOUGHREACT
and TOUGH2/EOS7C to predict induced geochemical interactions during Enhanced Gas Recovery at the K12B
site in the North Sea.
2.2 TOUGH2/ECO2N Input/Output files
The TOUGH2 and TOUGHREACT User‟s Guides (Pruess et al., 1999; Xu et al., 2005) describe the TOUGH2
code structure in detail. We focus our attention, therefore, on the input and output files needed to run the
workflow presented here, from a run using TOUGH2 and associated components.
The numerical scheme used in TOUGH2 is based on the Integral Finite Difference (IFD) method (Narasimhan
and Witherspoon, 1976) with a fully implicit time discretization. The residual form of the system of equations is
4
solved using the Newton-Raphson iteration method. The IFD method offers the advantage of being applicable to
regular or irregular discretization in three dimensions. In our study we will focus on the use of a hexahedral grid
block commonly used in reservoir fluid flow simulators.
A set of primary variables, which differs depending on the EOS routine used, must be defined for each gridblock in order to initiate the simulation run. When TOUGH2 is run with the ECO2N module, four variables are
needed: pressure, CO2 amount (either as gas saturation in two-phase conditions or mass fraction in single-phase
conditions), salinity (either as solid saturation or mass fraction of the dissolved part), and temperature. When the
EOS7C module is used, the amount of methane (either as gas saturation in two-phase conditions or mass fraction
hal-00614438, version 1 - 11 Aug 2011
in single-phase conditions) must be provided.
2.2.1 Input files
The input files containing the parameters needed to run TOUGH2 and TOUGHREACT are (Figure 1):
-
“flow.inp” This file mainly includes rock properties, time stepping and multi-phase fluid and heat flow
simulation parameters. Each of these parameters is defined for a domain which can contain more than
one grid block. This file also allows date selection for data editing.
-
“INCON” This file includes initial condition for porosity, permeability, and the four primary variables
(pressure, gas saturation (or mass fraction of the dissolved gas), salinity and temperature) for each gridblock.
-
“MESH” This file includes all grid block geometric properties (volume, exchange area, distance
between grid block centers, coordinates, definition of rock properties domain).
The structure of the input and output files in TOUGH2 (and associated EOS) is flexible. For example,
information from the “INCON” and “MESH” files can be incorporated in the “flow.inp” file. We use the
“INCON” and “MESH” files below as a basis in order to simplify the structure of the workflow.
The following files are needed for TOUGHREACT only.
-
“inchem” This file includes initial condition for the geochemical system (aqueous species, solid
mineral phases, etc.).
5
-
“chemical.inp” This file is used for defining the general parameters of the geochemical system (kinetic
rate of mineral species, secondary solid and aqueous phases, sorbed species, etc.).
-
“solute.inp” This file includes transport and other run parameters such as diffusion coefficients,
tolerance limits for convergence, printout flags, grid block attribution of a domain with a different
chemical composition (in agreement with the definition supplied in the “chemical.inp” file), etc.
The structure of the “inchem” file depends directly on the geochemical system defined in the “chemical.inp” file.
The assumed number of chemical species and mineral phases in the geochemical system will determine the
number of variables to be calculated and stored in the “inchem” file. Although TOUGHREACT can run without
the “inchem” file, in this study we use separate files for pre- and post-processing data treatment. In addition, it is
hal-00614438, version 1 - 11 Aug 2011
very easy to generate the “INCON”, “inchem” and “MESH” files with TOUGH2 and TOUGHREACT (see
User‟s Guide).
2.2.2 Output files
The output files are:
-
“flow.out” This file includes, for each date selected in the “flow.inp” file, the parameters calculated
during the simulations for each grid block (pressure, temperature, gas saturation, salinity, gas and liquid
density, gas and liquid velocity vectors, etc.).
-
“SAVE” This file gives the values calculated for porosity, permeability, and the four primary variables
(pressure, gas saturation, salinity, temperature) for each grid-block, at the end of the simulation. It can
be used as a restart file by being renamed “INCON”. There is some redundancy in data editing between
this file and the “flow.out” file.
The following output files are provided for the TOUGHREACT run only:
-
“savechem” This file provides the result of the calculated variables initially in the “inchem” file at the
end of the simulation. It can be used as a restart file by being renamed “inchem”.
-
“conc.dat” This file includes all aqueous chemical species concentrations for each grid block for every
selected date in the “flow.inp” file.
-
“min.dat” This file includes all mineral phase concentrations for each grid block for every selected date
in the “flow.inp” file.
6
Other files are produced during a TOUGHREACT/ECO2N run (“chemical.out”, “runlogout.out” etc.), but these
files are not necessary, in our case, for a complete editing of the simulated variables. On the other hand they
provide additional information very helpful for a better understanding of the simulation (run convergence, mass
balances etc.).
3. The workflow tool
Pre- and post-processing of data using TOUGH2 and TOUGHREACT are challenging for complex 3D models
including geological heterogeneities and multi-domain geological facies with associated mineralogy. We,
therefore, developed a workflow module as a set of Fortran subroutines able to handle (i) visualization of
hal-00614438, version 1 - 11 Aug 2011
TOUGHREACT file results, (ii) grid block attribution of heterogeneous properties, (iii) boundary condition
attribution, (iv) refinement and coarsening of 3D grid mesh. These modules can be broken down into three main
categories (Figure 2):
-
building or modification of existing mesh (GRID)
-
physical property attribution to grid blocks (PROPERTY)
-
3D Visualization of the simulation results (RESULT)
The 15 workflow modules, along with their input and output files, are listed below by category. The flow-chart
in Figure 3 describing the entire procedure can serve as a guide for using the workflow.
3.1 GRID category
“Mesh2visu”
This module enables the visualization of the TOUGH2 mesh generated by MESHMAKER. It enables the
conversion of “MESH” files to “grid.dat” files that can be opened with Tecplot or Paraview software. This file
will be used by many other modules. Visualization of a mesh grid containing local refinement is also possible.
Input file: “MESH”
Output file: “grid.dat”
“Petrel2Tough”
7
This module converts a mesh generated by the Petrel geomodeler (Schlumberger) to a “MESH” file, ready to be
run with TOUGH2. The Petrel mesh file must be in ASCII format, readable by the ECLIPSE “black oil”
simulator (Schlumberger). The “Petrel2Tough” module is programmed to identify ECLIPSE keywords such as
„SPECGRID‟, „COORD‟ and „ZCORN‟, to recognize grid block coordinates, as well as „PERMX‟, „PORO‟,
„ACTNUM‟ to identify active grid block properties.
It is important to notice that Petrel software can produce several types of mesh according the topological feature
of the generated grid. For structured grid cartesian block centered or rectangular corner point geometry is
required while unstructured grid models can be built with control volume finite elements, composite tetrahedral
grids or PEBI (Perpendicular bisector or voronoï) element techniques (Palagi and Aziz, 1994). The IFD method
assumed that node connections are orthogonal to interface surface which makes the PEBI method perfectly
hal-00614438, version 1 - 11 Aug 2011
suited for the use of TOUGH2 In the present study we used grids provided with a cartesian block centered
approach. Therefore, the mesh generated with Petrel2Tough module contain elements with identical volumes,
identical surface area and identical inter-nodal distances with respect to the original grid, but the orthogonality
between inter-nodal distance and surface area might not be preserved in agreement with the dip angle of the
modeled geological formations.
Input file: ASCII Petrel Mesh file
Output file: “MESH”
“Coarser”
This module enables the coarsening of a Petrel grid mesh in order to reduce the total number of grid blocks. It is
often needed when importing large grid meshes from geological models that can contain up to several million
grid blocks. Up-scaling of the reservoir properties is needed for assigning average property values to the new
grid blocks of the coarse mesh. A separate input file is needed to define the selected zone of the mesh to be
coarsened, and the degree of coarsening.
Input files: “Coarser_Change.itt”, “MESH”
Output file: “MESH” (modified)
“LGR”
This module is dedicated to “Local Grid Refinement” technique in mesh generation. Nested gridding is done in
order to refine one part of the domain. This is often needed when addressing near well zone issues like scaling,
8
dry out or salt precipitation. An input file is needed to describe the domain to be refined, and also the degree of
refinement. Several refinement levels are supplied with this option.
It is important to specify that for the LGR technique proposed in this study the node connections (parameters D1
and D2 in the CONNE block of the MESH file) of gridblocks generated for the sub-gridding decomposition is
kept orthogonal to the interface surface. Therefore the real distance “D” between the nodes of two adjacent grid
blocks (one in the sub-grid part and one from the initial mesh) is not preserved (D is not equal to D1+D2). The
correct implementation of sub-gridding elements would require introduction of composite elements containing
interpolated information and modification of the TOUGH2 source code files (as proposed by Garcià and Pruess,
2000). Based on several application case study recommendation of using the “LGR” module with a smooth
transition for sub-gridding decomposition of a factor of 5 is reported (Wertz and Audigane, 2010).
hal-00614438, version 1 - 11 Aug 2011
Input files: “Gigogne.itt”, “MESH”
Output file: “MESH” (modified)
“ActiveCell”
This module enables the inactivation of mesh elements for the simulation. It is very efficient for reducing the
total number of grid blocks. An input file is needed to select the domain where inactive elements should be
removed from the original mesh file.
Input files: “ActiveCell.itt“, “MESH”
Output file: “MESH” (modified)
3.2 PROPERTY category
“Rocktype”
This module enables the attribution of „ROCK‟ domains as defined in the TOUGH2 “flow.inp” input file. In
TOUGH2, several intrinsic rock properties, such as relative permeability function, capillary pressure, rock grain
density, pore compressibility and specific heat, are defined for each domain. The attribution of a domain to each
grid block defined in the “MESH” file becomes difficult when working with 3D mesh. This module requires an
input file to select the geometrical coordinates of the domain under consideration.
Input files: “Rocktype.itt“, “MESH”
Output file: “MESH”
9
“Chemtype”
In much the same way, “Chemtype” assigns domains representing (i) initial water formation, (ii) boundary water
formation and (iii) mineralogical assemblage. This is input in TOUGHREACT by the “solute.inp” input file.
This domain attribution is independent of the „ROCK‟ domain from a numerical point of view, but the “MESH”
file is needed as an input in order to recover the grid block coordinates.
Input files: “Chemtype.itt”, “solute.inp”, “MESH”
Output file: “solute.inp” (modified)
“Init_P_t”
This is a tool for initializing a pressure and temperature field according to a hydrostatic and a geothermal profile
hal-00614438, version 1 - 11 Aug 2011
along the depth. This is done by modifying the TOUGH2 “INCON” input file, which contains all pressure and
temperature data for each grid block before the simulation begins. The “MESH” file provides grid block
coordinates. A separate input file is needed for defining the imposed pressure and temperature gradient (linear
according to depth).
Input files: “Init_p_t.itt”, “MESH” , “INCON”
Output file: “INCON” (modified)
“Condinit”
This is a tool for initializing all parameters in the “INCON” file (porosity, permeability, pressure, salinity, gas
saturation field and temperature). This is done by entering values for each parameter by keyword (e.g. „PORO‟
for porosity) for each grid block in the order defined in the “MESH” file. A separate include file is, therefore,
needed, in addition to the “MESH” and “INCON” files, as an input.
Input files: “Condinit.itt”, “MESH”, “INCON”
Output file: “INCON” (modified)
“Volinf”
This module enables the implementation of Dirichlet boundary conditions on either existing grid blocks or new
grid blocks. In TOUGH2, assignment of such conditions is implemented by considering large volume grid
blocks, which implicitly keeps all grid block properties constant during the simulation run. We kept the same
approach by imposing an arbitrary volume of 1.0 E+29 m3 (easily modifiable). A separate input file makes it
10
possible to define the domain that contains either grid blocks for assignment, or grid blocks to which new grid
blocks need to be generated and connected – in which case, these new grid blocks are not used for the final
visualization of the simulation results. The “MESH” and “INCON” files are therefore needed for modifications.
Input files: “Volinf.itt”, “MESH”, “INCON”
Output files: “MESH” and “INCON” (modified)
3.3 RESULTS category
“Visuflow”
This module was specifically developed to visualize simulation results from the TOUGH2 “flow.out” file. The
hal-00614438, version 1 - 11 Aug 2011
input files, “flow.out” and “grid.dat” (generated by the “Mesh2visu” module described in the GRID category
above) are, therefore, needed. The output file, “gridflow.dat”, is in a Tecplot format and contains all properties
calculated in “flow.out” for each date selected.
Input files: “flow.out”, “grid.dat”
Output file: “gridflow.dat”
“Visuchem”
In much the same way as for the “Visuflow” module, this module was developed to visualize simulations
resulting from a TOUGHREACT run. Therefore, the “grid.dat” file (generated with the “Mesh2visu” module)
must be supplied along with the files containing aqueous, mineral and gas phase concentrations – "conc.dat”,
“min.dat” and “gas.dat”, respectively (default names). An output file, “gridchimie.dat”, in Tecplot format, is
produced, which contains all chosen parameters (chemical species, minerals) for editing and for each date
selected.
Input files: “grid.dat” and “conc.dat” and/or “min.dat” and/or “gas.dat”
Output file: “gridchimie.dat”
“Visuvectflow”
This module is specific to the edition of flow vectors which can be (optional) edited in the “flow.out” file. Flow
through connected grid blocks is provided for liquid, gas and heat and for each selected date. Vector coordinates
are calculated at the center of each grid-block averaging the vector components provided by TOUGH2 for each
11
grid-block connection (interface). As for the “Visuflow” module, the “grid.dat” and “gridflow.dat” input files are
needed to generate an output file, “gridvect.dat”, in Tecplot format.
Input files: “grid.dat” and “gridflow.dat”
Output file: “gridvect.dat”
“Mass-balance”
This module provides a mass balance of the CO2, salt and water system under aqueous, gaseous and solid phases
for the entire model. This information is provided during a TOUGH2 run in the “flow.out” file for each selected
date and for each domain defined in the ROCK block of the “flow.inp” file. The output is a text file easily
manageable using any table editor for plotting history evolution of such mass balances.
hal-00614438, version 1 - 11 Aug 2011
Input file: “flow.out”
Output files: A list of text files for each data field (gaseous CO2, aqueous/liquid CO2, solid salt etc.) and for
each rock domain.
“Tec2Para”
This module converts files in ASCII Tecplot format to ASCII Paraview format, because with Paraview software
it is necessary to define a separate file for each parameter (pressure, temperature etc.) whereas Tecplot handles
all parameters in one large file.
Input file: “gridflow.dat”
Output files: “geometry.insight”, “ensight.case”, list of files for each parameter
4. Application and case study
The workflow has been tested with a 3D case of (i) Local Grid Refinement (“LGR”) and (ii) inactive grid block
(“Activecell”) implementation for hydrodynamic multiphase modeling of a generic synthetic case, and two case
studies with realistic data using reactive transport data modeling from real pilot sites.
4.1 Stuttgart Benchmark
The University of Stuttgart recently published a code intercomparison study on CCS modeling (Class et al.,
2009). Three exercises were proposed and we participated in exercise 1 involving the leakage of CO2 from a
storage zone to an overlying aquifer by way of a leaky well (Figure 4). A semi-analytical solution has been
12
provided, to be compared with numerical modeling results (Nordbotten et al., 2005), the rate of leakage being
calculated in the middle section of the leaky well.
To develop the model, we first created a 3D grid with the TOUGH2 “MESHMAKER” tool. This 3D grid
includes the storage aquifer, the overlying aquifer and the aquitard. Using the “LGR” module, we then refined
the grid around the leaky well for the entire domain, and around the injection well only within the storage
aquifer. Using the “Activecell” module, we removed all grid blocks defining the aquitard domain except those
located within 1 m of the leaky well. “Mesh2Visu” shows of the final 3D grid mesh (Figure 5).
CO2 is injected for 1,000 days at a constant rate of 8.87 kg/s (1,600 m3/day). Using the “Init_p_t” module, the
pressure is initialized for the whole system (the temperature is assumed to be constant). Constant hydrostatic
pressure is imposed on all sides of the grid with the “Volinf” module. The simulated leakage rate is then
hal-00614438, version 1 - 11 Aug 2011
compared with the analytical solution (Figure 6). The results differ mainly because of the grid mesh resolution
around the wells, although they are satisfactory in comparison with other simulation results (see Figure 5 in
Class et al., 2009). Visualization of the CO2 gas bubble is possible with the “Visuflow” module (Figure 7). Using
“Visuvectflow”, a fluid flow vector field was built, enabling the visualization of the liquid and gas phase in the
system (Figure 8).
4.2 Case study 1: Injection of CO2 in the Sleipner field (Saline aquifer)
One of the first pilot sites for CO2 injection into a saline aquifer is the Sleipner gas field operated by Statoil in
the North Sea (Torp and Gale, 2004). The gas produced via several wells contains a relatively high CO2
concentration ranging from 4.5 to 9% (Korbøl and Kaddour, 1995). As a result of financial deterrents (a
Norwegian tax of $55 per ton of CO2 emitted), Statoil decided to store carbon dioxide underground in the Utsira
formation located above the gas reservoir (Kongsjorden et al., 1997). Injection started in 1996 at an approximate
average rate of 1 million tonnes per year. The Utsira formation is a large sandy aquifer extending over 26,100
km² at a depth of between 700 and 1000 m. Several seismic surveys have been conducted since the start of the
injection in order to visualize the evolution of the gas bubble (Torp and Gale, 2004; Arts et al., 2004). Several
modeling studies have also been done for history matching and for long-term predictions of the hydrodynamics
and geochemical impacts of this injection (Lindeberg and Bergmo, 2002; Johnson et al., 2001; Audigane et al.,
2007). BRGM has recently conducted a study to predict the multiphase flow and associated fluid-rock
interactions induced by CO2 injection at Sleipner (Picot-Colbeaux et al., 2009). The geometry of the reservoir
13
was built from a geological 3D model elaborated by the Sintef Institute with a geo-modeler Petrel (courtesy of
Bergmo from Lindeberg and Bergmo, 2002).
Using the “Petrel2Tough” module, the Petrel mesh deduced from the geological model is converted into a
"MESH" file readable by TOUGH2 and TOUGHREACT (Figure 9a). However, this resulted in 428,400 grid
blocks, which is definitely too many for running a reactive transport model with TOUGHREACT in a reasonable
length of time. The number of blocks was, therefore, reduced by coarsening the mesh using the “Coarser” tool
(Figure 9b).
Intra-shale layers were identified by seismic interpretation. This geological heterogeneity enables the local
accumulation of CO2 during its upward migration due to buoyancy. The mineralogical assemblage of the model
was divided into two domains in order to consider both sand and shale formations in the system. Heterogeneous
hal-00614438, version 1 - 11 Aug 2011
hydrogeological and geochemical properties were attributed to these two geological formations using the
“Rocktype” and “Chemtype” modules.
Simulations were done with the coarser mesh (41,400 grid blocks) for a 3-year injection period. Visualization of
both multiphase flow and fluid-rock interactions in a Tecplot format was done with the “Visuflow” and
“Visuchem” modules. Results predict a rather low geochemical impact, as seen on a 2D cross-section of the
change in porosity induced by CO2 injection (Figure 9c). After three years of injection, the primary trapping
process is the geological structure of the system – mainly the cap rock formation overlying the reservoir aquifer,
which prevents the upward migration of the free CO2 due to buoyancy. Solubilization has started but does not
represent more than 10% of the injected mass of carbon dioxide.
4.3 Case study 2: Injection of CO2 in the K12B Field (Enhanced Gas Recovery)
As part of the Dutch Climate Policy to implement the Kyoto protocol, the nearly depleted K12-B gas reservoir in
the North Sea operated by Gaz de France Netherlands was chosen as a demonstration site for the Offshore Reinjection of CO2 project (ORC project). The K12-B field is located in the Dutch sector of the North Sea (150 km
from Amsterdam). Gaz de France Production Nederland B.V. has been producing gas with four wells (K12-B1,
K12-B2, K12-B5 and K12-B7) since 1985 and the field is now almost depleted. Located at a depth of about
3,800 m, the gas field consists of a number of tilted fault blocks that are not, or only barely, in pressure
communication. Fluid flow is strongly influenced by the presence of wrench fault tectonics resulting in some
further compartmentalization of individual fault blocks (van der Meer et al., 2004). The present study addresses
CO2 injection into compartment 3 of the gas field during phase 3 of the ORC project. We used the geological
14
model and the associated 3D grid mesh generated by Geel et al. (2006) in agreement with previous studies
(Figure 10). This 3D model was built with Petrel and compartment 3 initially contained 27,747 grid blocks. A
coarser 3D grid in a TOUGH2 "MESH" file format was built using the “Petrel2Tough” and “Coarser” tools in
order to be able to run both hydrodynamic and reactive simulations with TOUGH2/EOS7C and
TOUGHREACT, respectively. Four different geological formations have been identified in compartment 3. A
clay cap rock formation (Ten Boer) with two sandstone formations (Upper and Lower Slochteren) separated by a
clay layer near the cap rock unit (Ameland). We used the “Rocktype” and “Chemtype” modules to implement
the heterogeneities on the final 3D grid containing 3,171 grid blocks (Figure 11).
The principal aim of the operator was to assess the efficiency of the Enhanced Gas Recovery technique using
CO2 injection. Carbon dioxide was injected in order to both reduce CO2 emissions and recover CH4. This
hal-00614438, version 1 - 11 Aug 2011
technology is often envisaged at the final stage of field production, which is the case in the K12B gas field. We
used a specific EOS module for TOUGH2 called EOS7C (Oldenburg et al., 2004), which can handle a CO2-CH4
gas mixture. The injection rate of 10 kg/s corresponds to phase 3 of the research project considering a full-scale
injection of about 310-475 Kton per year. K12-B1 and K12-B5 are production wells, while CO2 is injected in
K12-B6. We chose a production rate of 1 kg/s for each producer to limit the CO2 breakthrough time. Enhanced
Gas Recovery could be greater if the production rate was greater but this simulation was carried out as if the
injection of large quantities of CO2 was a primary objective.
Gas recovery and geochemical impacts were assessed and visualized with the “Visuflow” and “Visuchem”
modules (Figure 12). Since the reservoir originally contained 13% of CO2 in the gas phase, the aqueous phase
was already in equilibrium with high carbon dioxide fugacities and therefore the initial brine was very acidic
with a pH of around 4.7. Simulations predict, after 10 years of injection, a very small modification of pH in the
entire reservoir causing very little geochemical fluid-rock interaction (Figure 13). Structural and solubility
trapping are the two main carbon dioxide storage processes, as shown by the “MassBalance” module (Figure
14). Since the reservoir is a closed system, a large pressure increase of 50 bar is predicted, which is in agreement
with the very great under-pressurization (from about 350 bar to 40 bar) previously induced by methane
production. These results are discussed in greater detail in Audigane et al. (2008).
5. Conclusion
Modeling of carbon dioxide geological storage requires a precise description of the geological formation
envisaged as a storage target in order to improve storage capacity prediction and risk assessment of such a new
15
technology. 3D geological modeling tools are now able to reproduce such a detailed geometry and pertinent geostatistical simulations for modeling hydrogeological reservoir property fields. Multiphase flow modeling done
with the TOUGH2 simulator, developed by the Lawrence Berkeley National Laboratory (LBNL) and recently
used, along with other numerical tools, for CCS modeling, is poorly suited to handling complex mesh geometry
and 3D heterogeneous porous media. Pre and post processing of such a tool is still lacking, though several
attempts have been made to provide graphic-interface software.
The workflow proposed in this study is a set of Fortran subroutines that are easy to use and adapt for any
TOUGH2 code up-date. This set of modules can generate TOUGH mesh and corresponding property fields to
assess liquid and gas hydrodynamics as well as fluid-rock geochemical interactions. Post processing is provided
as export in two classical mesh visualizer formats: Tecplot and ParaView. Additional options are also provided
hal-00614438, version 1 - 11 Aug 2011
to handle boundary condition attribution and mesh refinement and coarsening. A converter of geo-modeler Petrel
files in ASCII format into TOUGH mesh files has also been developed.
The workflow has been applied to different test cases in order to validate its use in a carbon dioxide geological
storage context. An unstructured mesh geometry has been modeled to simulate an abandoned well linking a
saline storage reservoir with an overlying aquifer. Import of the Petrel format and coarsening of the mesh in
order to reduce the number of grid blocks enables the simulation of induced geochemical fluid-rock interaction
associated with CO2 injection in the sandy aquifer of the Utsira formation in the North Sea, Norway. Another
context of carbon dioxide storage was tested in a context of methane production using reservoir data from the
K12b gas field, also in the North Sea.
The flexibility of the workflow enables it to be easily extrapolated to other contexts such as nuclear waste
disposal or geothermal production by adapting the workflow to specific equation of state modules (EOS1,
EOS3…) developed by LBNL for TOUGH2.
References
Alcott, A., Swenson, D., Hardeman, B., 2006. Using Petrasim to create, execute, and post-process TOUGH2
models. In: Proceedings TOUGH Symposium, Lawrence Berkeley National Laboratory, Berkeley, California,
May 15-17.
16
André, L; Audigane, P; Azaroual, M, Menjoz, A., 2007. Numerical modeling of fluid-rock chemical interactions
at the supercritical CO2-liquid interface during CO2 injection into a carbonate reservoir, the Dogger aquifer
(Paris Basin, France). Energy Conversion and Management 48(6), 1782-1797.
Arts. R., Eiken. O., Chadwick. A., Zweigel. P., Van Der Meer. B., Kirby., G., 2004. Seismic monitoring at the
Sleipner underground CO2 storage site (North Sea), In: Baines, S.J & Worden. R. H. (Eds.) Geological storage of
carbon dioxide, Geological Society of London, Special publications, 233, pp. 181-191.
Audigane, P; Gaus, I; Czernichowski-Lauriol, I, Pruess, K., Xu, T., 2007. Two-dimensional reactive transport
modeling of CO2 injection in a saline Aquifer at the Sleipner site, North Sea. American Journal of Science 307,
hal-00614438, version 1 - 11 Aug 2011
974-1008.
Audigane, P., Oldenburg, C.M., van der Meer, B., Geel, K., Lions, J., Gaus, I., Robelin, Ch., Durst, P., Xu, T.,
2008. Geochemical modeling of the CO2 injection into a methane gas reservoir at the K12-B Field, North Sea. In
M. Grobe, J.C. Pashin, and R.L. Dodge, (Eds.), Carbon dioxide sequestration in geological media – State of the
science, AAPG Studies, pp. 1-20.
Bachu, S., Gunter, W.D., Perkins, E.H., 1994. Aquifer disposal of CO2: hydrodynamic and mineral trapping.
Energy Conversion Management 35(4), 269-279.
Burnell, J.G., White, S.P., Osato, K., Sato, T., 2003. GeoCad, A pre and postprocessor for TOUH2. In:
Proceedings TOUGH Symposium, Lawrence Berkeley National Laboratory, Berkeley, California, May 12-14.
Chiaberge, Ch., Audigane, P., 2009. Manuel d‟utilisation des outils de pre et post processing pour la
modélisation de stockage géologique de CO2 avec les logiciels TOUGH. Brgm Report RP57457, 50 pp.
Class, H., Ebigbo, E., Helmig, R., Dahle, H.K., Nordbotten, J.M., Celia, M.A., Audigane, P., Darcis, M., EnnisKing, J., Fan, Y., Flemisch, B., Gasda, S.E., Jin, M., Krug, S., Labregere, D., Naderi, A., Pawar, R.J., Sbai, A.,
Thomas, S.G., Trenty, L., Wei, L., 2009. A benchmark study on problems related to CO2 storage in geologic
formations, Computational Geosciences, onlinefirst. doi 10.1007/s10596-009-9146-x
17
Garcià, J., Pruess, K. (2000). Local grid refinement for multi-scale geothermal reservoir simulation with
TOUGH2. Lawrence Berkeley National Laboratory, Berkeley, CA, Technical Report LBNL-45646.)
Geel, C.R., Arts, R.J., van Eijs, R.M.H.E., Kreft, E., Hartman, J., D‟Hoore, D., 2006. Geological site
characterisation of the nearly depleted K12-B gas field, offshore the Netherlands. In: Proceedings CO2SC
Symposium Lawrence Berkeley National Laboratory, Berkeley, California, March 20-22, pp. 124-127.
Intergovernmental Panel on Climate Change, 2005. IPCC Special Report on Carbon Dioxide Capture and
hal-00614438, version 1 - 11 Aug 2011
Storage. Cambridge University Press, New York, USA, 431 p.
Johnson, J.W., Nitao, J.J., Steefel, C., Knaus, K.G., 2001. Reactive transport modeling of CO2 storage in saline
aquifers to elucidate fundamental processes; trapping mechanisms, and sequestration partionning, In: Geological
Society of London, Special publication on carbon sequestration technologies, UCRL-JRNL-205627, 45 p.
Kongsjorden. H., Karstad. O., Torp. T. A., 1997. Saline aquifer storage of carbon dioxide in the Sleipner project.
Waste Management 17(5/6), 303-308.
Korbøl. R., Kaddour. A., 1995. Sleipner Vest CO2 disposal injection of removed CO2 into the Utsira Formation.
Energy Conversion Management 36(6-9), 509-512.
Lindeberg, E., Bergmo, P., 2002. The long term fate of CO2 injected into an aquifer, In: Gale J.and Kaya Y.
(Eds.), Proceedings of the 6th International Conference on Greenhouse Gas Control Technology, Kyoto, Japan,
1-4 October, pp 489.
Narasimhan, T.N., and Witherspoon, P.A., 1976, An integrated finite difference method for analyzing fluid flow
in porous media: Water Resources Research, v. 12, p. 57-64.
Nordbotten, J.M., Celia, M.A., Bachu, S., Dahle, H., 2005. Semi-analytical solution for CO2 leakage through an
abandoned well. Environmental Science and Technology 39(2), 602-611.
18
Oldenburg, M.C., Moridis, G.J., Spycher, N., Pruess, K., 2004. EOS7C Version 1.0: TOUGH2 module for
carbon dioxide or nitrogen in natural gas (methane) reservoirs, Lawrence Berkeley National Laboratory, Paper
LBNL-56589, Berkeley, USA, 53 p.
O‟Sullivan, M.J., Bullivant, D.P., 1995. A Graphical Interface for the TOUGH Family of Flow Simulators,
Proceedings of the TOUGH Workshop 1995, Lawrence Berkeley Laboratory Report LBL-37200, Berkeley, CA,
USA, 90-95.
Palagi, C. L., Aziz, K., 1994. Use of Voronoi Grid in Reservoir Simulation. SPE Advanced Technology Series,
Vol. 2, No 2, 69-77.
hal-00614438, version 1 - 11 Aug 2011
Pan, L., 2003. WinGridder-An interactive grid generator for TOUGH2. In: Proceedings TOUGH Symposium,
Lawrence Berkeley National Laboratory, Berkeley, California, May 12-14.
Pan, L., 2008. User information for WinGridder Version 3.0, Lawrence Berkeley National Laboratory. Report
LBNL-273E, Berkeley, California, USA, 76 p.
Picot-Colbeaux G., Mathurin F., Audigane P., 2009. Short term reactive multiphase flow modeling for
performance assessment of CO2 storage at Sleipner, North Sea. BRGM, Orléans, France, BRGM Report
RP57267-FR, 2009, 39 pp.
Pruess, K., Oldenburg, C., Moridis, G., 1999. TOUGH2 user‟s guide, version 2.0, Lawrence Berkeley National
Laboratory, Berkeley, California, USA, Report LBNL-43134, 198 pp.
Pruess, K., García, J., 2002. Multiphase flow dynamics during CO2 injection in saline aquifers. Environmental
Geology 42, 282-295. doi:10.1007/s00254-001-0498-3.
Pruess, K, Garcia J, Kovscek T, Oldenburg C, Rutqvist J, Steefel C, Xu T, 2004. Code intercomparison builds
confidence in numerical simulation models for geologic disposal Of CO2, In: J. Gale and Y. Kaya, (Eds), 6th
International Conference on Greenhouse Gas Control Technologies, 1-4 October, 2002, Kyoto, Japan. Energy
29(9-10), pp 1431-1444.
19
Pruess, K., Spycher, N., 2007. ECO2N – A fluid property module for the TOUGH2 code for studies of CO2
storage
in
saline
aquifers,
Energy
Conversion
Management
48(6),
1761-1767.
doi:10.1016/j.enconman.2007.01.016
Pruess, K., Muller, N., 2009. Formation dry-out from CO2 injection into saline aquifers: 1. Effects of solids
precipitation and their mitigation. Water Resources Research 45, W03402. doi:10.1029/2008WR007101.
Sato, T., Ohsato, K., Shiga, T., 2003. G-Star-Base (G*Base)-A data base system for underground information
and post-processing for TOUGH2, In: Proceedings TOUGH Symposium, Lawrence Berkeley National
hal-00614438, version 1 - 11 Aug 2011
Laboratory, Berkeley, California, May 12-14.
Swenson, D., Hardeman, B., Persson, C., Thornton, Ch., 2003. Using Petrasim to create, execute, and postprocess TOUGH2 models, In: Proceedings TOUGH Symposium, Lawrence Berkeley National Laboratory,
Berkeley, California, May 12-14.
Torp, T.A., Gale, J., 2004. Demonstrating storage of CO2 in geological reservoirs: the Sleipner and SACS
projects, In: J. Gale and Y. Kaya, (Eds), 6th International Conference on Greenhouse Gas Control Technology,
Kyoto, Japan, 1-4 October, 2002. Energy 29(9-10), pp 1361-1369.
van der Meer, L.G.H., Kreft, E., Geel, C., Hartman, J., 2004. K12-B Test site for CO2 storage and enhanced gas
recovery, In: 14th Europec Biennal Conference, Madrid, Spain, 13-16 June, paper SPE/IADC 94128.
Werz, F. Audigane, P., 2010. Maillage non-voronoï et erreure induite : application au code de calcul TOUGH2
pour le CO2. Brgm Report RP 58594, 37 pp.
Xu. T., Apps, J.A., Pruess, K., 2003. Reactive geochemical transport simulation to study mineral trapping for
CO2 disposal in deep arenaceous formations. Journal of Geophysical Research 108(B2), 2071,
doi:10.1029/2002JB001979
20
Xu, T; Sonnenthal, E; Spycher, N, Pruess, K., 2005. TOUGHREACT User‟s guide: A simulation program for
non isothermal multiphase reactive geochemical transport in variably saturated geologic media, Lawrence
Berkeley national Laboratory, Berkeley, California, USA, Report LBNL-55460, 192 p.
Xu, T; Sonnenthal, E; Spycher, N, Pruess, K., 2006. TOUGHREACT-A simulation program for non-isothermal
multiphase reactive geochemical transport in variably saturated geologic media: Applications to geothermal
injectivity and CO2 geological sequestration, Computers & Geosciences 32(2), 145-165.
Zhou., Q., Birkholzer J.T., Tsang C-F., Rutqvist. J., 2008. A method for quick assessment of CO2 storage
capacity in closed and semi-closed saline formations, International Journal of Greenhouse Gas Control 2(4), 626-
hal-00614438, version 1 - 11 Aug 2011
639.
Zhang, W., Li Y., Xu T., Cheng HL, Zheng Y, Xiong P., 2009. Long-term variations of CO2 trapped in different
mechanisms in deep saline formations: A case study of the Songliao Basin, China. International Journal of
Greenhouse Gas Control 3(2), 161-180.
21
hal-00614438, version 1 - 11 Aug 2011
Figure 1: Input and output files for TOUGH2 and TOUGHREACT. Perm=permeability, Poro=porosity,
P=Pressure, T=temperature, Sg=gas saturation, NaCl=salt mass concentration.
22
hal-00614438, version 1 - 11 Aug 2011
Figure 2: Modules developed for the workflow. Each module is a separate FORTRAN executable. Associated
input and output files are listed under the module name.
23
hal-00614438, version 1 - 11 Aug 2011
Figure 3: Workflow descriptive flow chart including modules and links with TOUGH2 and TOUGHREACT
input and output files.
Figure 4: Schematic drawing of leakage scenario proposed for code benchmarking by Class et al. (2009).
24
hal-00614438, version 1 - 11 Aug 2011
Figure 5: Mesh generation using “LGR” and “ActiveCell” modules applied to problem used for code
benchmarking (Class et al., 2009).
Figure 6: Results of code benchmarking (Class et al., 2009) and comparison with semi-analytical solution.
25
hal-00614438, version 1 - 11 Aug 2011
Figure 7: 3D visualization of gas bubble simulated with TOUGH2 for code benchmarking (Class et al., 2009)
using the “Visuflow” module.
Figure 8: 3D visualization of gas flow vector field simulated with TOUGH2 for code benchmarking (Class et
al., 2009) using the “Visuvectflow” module.
26
hal-00614438, version 1 - 11 Aug 2011
Figure 9: (a) Petrel mesh of the Utsira geological formation, (b) converted to a TOUGH2 mesh, and (c) reactive
transport simulation with TOUGHREACT using mesh from (b)
Figure 10: K12B-field model built with Petrel geo-modeler by SINTEF Institute (Geel et al., 2006). CO2
injection and CH4 production is simulated in compartment 3 of the field.
27
hal-00614438, version 1 - 11 Aug 2011
Figure 11: TOUGH2 mesh of compartment 3 for the K12-B gas field using “Petrel2Tough” and “Coarser”
modules.
Figure 12: 3D visualization of multiphase flow simulation of Enhanced Methane production injecting carbon
dioxide on the K12-B field using the TOUGH2-EOS7C simulator and the “Visuflow” module. “XCO2” and
“XCH4” are the mass fraction of CO2 and CH4 in the gas phase, respectively.
28
hal-00614438, version 1 - 11 Aug 2011
Figure 13: 3D visualization of reactive flow simulation of CO2 injection in the K12-B field using the
TOUGHREACT simulator and the “Visuchem” module.
Figure 14: Mass balance between injected CO2 and produced CH4 simulated at the K12-B gas field using the
“Massbalance” module.
29
