We discuss a workflow implemented for coupling arbitrary numerical simulators considering complex geological models with discrete faults. This includes grid conversion of geological model grids generated with the Petrel software package to different simulator input formats within a few minutes for multi-million element models. We introduce the conceptual workflow design and tools required for the workflow realization. In this context, different fault representations can be realized including discrete fault planes supported by a virtual element concept in the multiphase flow simulator or an equivalent porous media approach using the mechanical ubiquitous joint model. GFZ German Research Centre for Geosciences, Section 5.3 – Hydrogeology, Telefgrafenberg, 14473 Potsdam, Germany © 2014 The Authors.
Benjamin Nakaten, Thomas Kempka
We discuss a workflow implemented for coupling arbitrary numerical simulators considering complex geological models with discrete faults. This includes grid conversion of geological model grids generated with the Petrel software package to different simulator input formats within a few minutes for multi-million element models. We introduce the conceptual workflow design and tools required for the workflow realization. In this context, different fault representations can be realized including discrete fault planes supported by a virtual element concept in the multiphase flow simulator or an equivalent porous media approach using the mechanical ubiquitous joint model.
GFZ German Research Centre for Geosciences, Section 5.3 – Hydrogeology, Telefgrafenberg, 14473 Potsdam, Germany
© 2014 The Authors.
Commercial hydro-mechanical simulation workflows are provided by different software vendors. For instance, the Schlumberger Petrel-ECLIPSE-VISAGE [1-3] workflow offers the opportunity to implement geological models from seismic and/or well data, and to conduct multiphase flow simulations using available Schlumberger software products (e.g. ECLIPSE). These can then integrated into mechanical simulations using the Schlumberger VISAGE software package. Other products providing similar capabilities are e.g. the Baker & Hughes JewelSuite  combined with the commercial finite element simulators Simulia Abaqus  or DIANA  providing the option of integrating different commercial multiphase flow simulators (e.g. CMG IMEX and GEM [7,8] as well as Schlumberger ECLIPSE).
Nevertheless, the available commercial workflows show a general lack of flexibility by means of integration of scientific open source numerical simulators as well as results of specific process simulations (e.g. reactive transport simulations). For that reason, we implemented a workflow that allows for coupling of an arbitrary reservoir simulation software package, demonstrated here based on the TOUGH2-MP/ECO2N [9,10] simulator, to the commercial Itasca FLAC3D  simulator.
Using the Petrel software package, we implemented a workflow that allows for the coupling of TOUGH2- MP/ECO2N multiphase flow simulations with FLAC3D geomechanical simulations. The difference of our approach to e.g. the hydro-mechanical coupling introduced and extended by Rutqvist et al. [12,13] is that we use the Petrel software package for model gridding and intermediate data processing for data transfer between different numerical simulation grids by parameter and property upscaling methods. Furthermore, our approach allows for integration of highly complex 3D structural geological models, involving an arbitrary amount of faults and fractures.
The Petrel- to-FLAC3D workflow was implemented to use the same model geometry previously developed for multiphase flow simulations in TOUGH2-MP/ECO2N, generated via the Petrel-to-TOUGH2-MP workflow, also allowing for the application of the virtual element technique to represent discrete fault as suggested by Nakaten et al. . However, the core of the workflow is based on the equivalent porous media representation of faults introduced by exact fault geometry representation for fluid flow as previously applied in single phase fluid flow simulations .
The benefit of the introduced workflow compared to the currently available commercial solutions is its adaptability to arbitrary multiphase flow and geomechanical simulators allowing for integration of arbitrary commercial or scientific open source numerical simulator. In this context, the Petrel software package is applied for geological modelling and gridding of the numerical models as well as for result visualization. The linkage to the Petrel software package ensures a highly efficient workflow, whereby gridded Petrel models can be transferred into large-scale hydro- mechanical FLAC3D simulation models within a few minutes.
The workflow consists of two phases, the conversion and the integration phases before the coupled numerical simulation run. The conversion phase contains tree main steps (Fig. 1) comprising multiple interacting processes to convert the model grid from VIP data format  into valid TOUGH2-MP and FLAC3D input files. To increase the conversion efficiency, we developed a flexible binding in Cython  to implement modules requiring user interaction and manipulation in Python  such as input and output parsers (Fig. 1.). All other modules are implemented in C++ , since memory efficient handling of large data sets is crucial when dealing with large- scale and multi-million element models (Fig. 1.).
Fig. 1. Schematic representation of the VIP data set conversion initially exported from Petrel to valid TOUGH2-MP and FLAC3D input files in three combined modular steps. Green boxes are the input and output parser modules. Blue boxes indicate high effort computational operations and data structures.
Data processing starts in Step 1 (cf. Fig. 1) with reading the numerical grid exported from the Petrel software package and given in the VIP format. After reading the coordinates of one element, these data are integrated into the grid data structure. To reduce memory utilization for data storage, we apply the Point Cloud Library (PCL)  as main data trunk for 3D coordinate data which was extended by a mechanism for coordinate duplication reduction.
Element data are stored within an adapted grid data structure and linked to the coordinate trunk preserving all structural model features. Providing a way to traverse through the grid at any time of the conversion phase, the grid walkthrough module allows for a straightforward detection of the associated grid position of an element. This module supports following structured i-, j- and k-order paths as well as unstructured link patterns (e.g. more than six element neighbors for one element may result from a fault throw or local grid refinement). The unstructured link patterns can be supported by additional configuration files as e.g. ECLIPSE fault maps exported from the Petrel software package. After implementing the grid geometry data, the model is parameterized by reading parameter files exported from Petrel, e.g. porosity and permeability data as well as material data to identify boundary, well and/or observation elements.
Considering faults in a structured hexahedral grid, interface areas are located between two elements across a fault. In Step 2, the ECLIPSE fault map generated by the Petrel software package is used to detect all elements located at faults, to group these elements as fault-neighboring elements and to identify the fault planes related to the respective elements. Furthermore, a fault element order path for the grid walkthrough module is generated in this step. Nevertheless, an integration of virtual elements for fault representation is not mandatory, since both types of fault representation are known to provide comparable simulation results .
Step 3 is applied to prepare the conversion result output and write it to simulator-specific input files. A specific module generates all TOUGH2-MP related input data including element ids, centre coordinates and element volumes determined as discussed by Grandy  as well as element connection data (interface areas, gravity vector to cell centre connection angles and cell centre to interface distances).
Two different cases of connection data computation are applied: one is used for element connections where the four element corner coordinates at one face are equal to those at the neighbouring element. The other case is applied for element connection computation of fault inter-connections, whereby the presence of a fault throw results in elements on one fault side with corner point coordinates different to the opposite side/element (hanging nodes). These connections are addressed by polygon clipping  considering their respective interface areas.
Writing of TOUGH2-MP input files is carried out by a separate Python module which aggregates and parses the required data. Here, compared to the pyTOUGH framework , we ordinarily derive the numerical simulation grid data and all parameters from the Petrel software package instead of manually generated configuration files. Fault interface data to account for hanging nodes in the FLAC3D simulator by linking neighbouring nodes that have non- unique coordinates at a fault plane are computed by a specific FLAC3D output module written in Step 3 (Fig. 1). This Python-based module generates two different file types, one containing the grid geometry data and additional group parameters as well as another defining the interface mapping. A VTK output module (Fig. 1) can be applied to produce VTK input data files allowing for rapid visualization of assigned material and fault properties.
Phase 2 deals with the integration of Phase 1 simulator input files into TOUGH2-MP and FLAC3D. Following the completion of Phase 1, input files for the TOUGH2-MP simulator can be further parameterized (e.g. integration of well or boundary elements, etc.). At this point, the resulting input file can be directly applied for a coupled numerical simulation run. For integration of faults into the FLAC3D simulator, we implemented a specific routine written in the FLAC3D-internal programming language FISH to realize the previously computed interface locations as faults in the geomechanical model.
To determine the FLAC3D interface locations, the elements previously labeled in Phase 1 are used for orientation. Following interface generation (Fig. 2), a parameterization of the model and interfaces is carried out. Interfaces allow for a more detailed fault analysis (e.g. providing simulation results including shear and normal displacements at the fault planes), but require more computational efforts by means of their implementation in the numerical model grid. Nevertheless, fault planes cannot currently be considered for fluid flow, without application of a virtual element approach . The previously discussed option of fault integration in FLAC3D based on the equivalent porous media approach can be carried out using the FLAC3D ubiquitous joint model that may be populated by fault geometries loaded from an AutoCAD DXF file .
Both fault integration concepts are perfectly applicable for representation of discrete faults and complex fracture networks, whereas the second approach allows for a less complicated hydro-mechanical model coupling procedure. A comprehensive study using simple model geometries carried out by Cappa and Rutqvist  demonstrated that identical mechanical results can be obtained using both fault representation methods. This is mainly due to the ubiquitous joint model incorporating weak planes into the elements which allow for an exact representation of the fault dip direction and dip angle.
Fig. 2. Fault representation by discrete interfaces in a structured geological model with faults determined by fault- and face-wise labeled elements (left). Both FLAC3D interfaces are generated based on specifically labeled elements (right).
Since the applied complex numerical simulation grids are identical for the multiphase flow and geomechanical simulations, a straightforward and iteration-based parameter exchange without introduction of the upscaling or interpolation errors can be realized by application of our workflow as previously discussed for simple model geometries by Rutqvist et al. [12,13] as well as Cappa and Rutqvist . Examples of the workflow application are given e.g. in Tillner et al.  and Kempka et al. , whereas a detailed view of the fault representation strategies is given in Figs. 2 and 3.
Fig. 3. Fault representation in a structural geological model implemented and gridded with the Petrel software package (left) and in FLAC3D based on ubiquitous joint model to be applied for the equivalent porous media approach (right) [26,27].
We implemented a workflow that allows for coupling of arbitrary simulators considering complex geological models with discrete faults and/or fracture networks. The introduced workflow allows for grid conversion from geological model grids generated with the Petrel software package to different simulator formats within a few minutes for multi-million element models. Based on the scientific open source simulator TOUGH2-MP/ECO2N and the commercial geomechanical simulator FLAC3D, we demonstrated the conceptual workflow design and the tools applied for its realization. The workflow can be applied for different fault representations including discrete fault planes realized by the virtual element concept in the chosen multiphase flow simulator or an equivalent porous media approach supported by the mechanical ubiquitous joint model.
 Schlumberger. Petrel Seismic-to-Evaluation Software, Version 2011.1; 2011.
 Schlumberger. ECLIPSE Reservoir Engineering Software, Version 2011.3; 2011.
 Schlumberger. VISAGE, Technical Manual, Schlumberger Information Solutions, Version 2010; 2010.  Baker & Hughes. JewelSuite JOA, JewelSuite 2012 User Manual; 2012.
 Dassault Systèmes. F.E.A. Abaqus, D.S. Simulia Dassault Systèmes; 2012.
 DIANA. DIANA – Finite Element Analysis Release 7.2. Program and User’s Documentation. TNO Building and Construction Research.  CMG. IMEX – Three-Phase, Black-Oil Reservoir Simulator. Computer Modelling Group; 2013.
 CMG. GEM – Compositional & Unconventional Reservoir Simulator. Computer Modelling Group; 2013.
 Pruess K. ECO2N: A TOUGH2 Fluid Property Module for Mixtures of Water, NaCl, and CO2. Report LBNL-57952, Lawrence Berkeley National Laboratory, Berkeley, California, 2005.
 Zhang K, Wu YS, Pruess K. User’s Guide for TOUGH2-MP – A Massively Parallel Version of the TOUGH2 Code. Report LBNL-315E, Earth Sciences Division, Lawrence Berkeley National Laboratory, Berkeley, California; 2008.
 Itasca. FLAC3D Software Version 5.01, Advanced Three Dimensional Continuum Modelling for Geotechnical Analysis of Rock, Soil and Structural Support; 2013.
 Rutqvist J, Wu YS, Tsang CF, Bodvarsson G. A modeling approach for analysis of coupled multiphase fluid flow, heat trans fer, and deformation in fractured porous rock. Numerical Methods in Rock Mechanics 2002; 39(4): 429-442.
 Rutqvist J. Status of the TOUGH-FLAC simulator and recent applications related to coupled fluid flow and crustal deformations. TOUGH Symposium 2011; 37(6): 739-750.
 Nakaten B, Tillner E, Kempka T. Virtual Elements for Representation of Faults, Cracks and Hydraulic Fractures in Dynamic Flow Simulations. Energy Procedia 2013; 40: 447-453. doi:10.1016/j.egypro.2013.08.051.
 Zhang Y, Underschultz JR, Gartrell A, Dewhurst DN, Langhi L. Effects of regional fluid pressure gradients on strain localisation and fluid flow during extensional fault reactivation. Mar Petrol Geol 2011; 28(9):1703-1713.
 Landmark Nexus Reservoir Simulation Software.
 Wilbers I, Langtangen HP, Ødegård Å. Using Cython to Speed up Numerical Python Programs. Trondheim: Proceedings of MekIT’09 Fifth National Conference on Computational Mechanics, 2009; 495-512.
 Van Rossum G, Drake FL. An Introduction to Python. Network Theory Limited., 2011.
 Ellis MA, Stroustrup B. The Annotated C++ Reference Manual. Addison-Wesley Longman Publishing; 1990.
 Rusu R, Cousins S. 3D is here: Point Cloud Library (PCL). IEEE Robot Autom Mag 2011; 1-4. doi:10.1109/ICRA.2011.5980567.
 Cappa F, Rutqvist J. Modeling of coupled deformation and permeability evolution during fault reactivation induced by deep underground injection of CO2. Int J Greenh Gas Con 2011; 5(2): 336-346.
 Grandy J. Efficient Computation of Volume of Hexahedral Cells. DOE Scientific and Technical Information. Online:
http://www.osti.gov/bridge/purl.cover.jsp?purl=/632793-4p2OLa/webviewable/632793.pdf, 1997. doi:10.2172/632793.  Simonson L, Suto G. Geometry Template Library for STL-like 2D Operations. Colorado: GTL Boostcon 2009.
 Croucher A. PyTOUGH: a Python Scripting Library for Automating TOUGH2 Simulations. New Zealand: Geothermal Workshop. 2011.
 AutoCAD 2012 DXF Reference. Autodesk Inc. 2011. Online: http://images.autodesk.com/adsk/files/autocad_2012_pdf_dxf- reference_enu.pdf, 2011.
 Tillner E, Shi JQ, Bacci G, Nielsen CM, Frykman P, Dalhoff F, Kempka T. Coupled dynamic flow and geomechanical simulations for an integrated assessment of CO2 storage impacts in a saline aquifer. Austria: General Assembly European Geosciences Union; 2013.
 Kempka T, Nielsen CM, Frykman F, Shi JQ, Bacci G, Dalhoff F. Coupled Hydro-Mechanical Simulations of CO2 Storage Supported by Pressure Management Demonstrate Synergy Benefits from Simultaneous Formation Fluid Extraction. Oil Gas Sci Technol 2014. doi:10.10501/ogst/2014029. (in press).
Corresponding author Benjamin Nakaten:
Tel.: +49-331-288-28725; fax: +49-331-288-1529.
E-mail address: firstname.lastname@example.org
This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/3.0/).
Peer-review under responsibility of the Organizing Committee of GHGT-12