Environ Geol (2005) 47: 1017–1027 DOI 10.1007/s00254-005-1236-z

Thomas Kalbacher Wenqing Wang Chris McDermott Olaf Kolditz Takeo Taniguchi

Received: 19 August 2004 Accepted: 7 January 2005 Published online: 5 March 2005
Springer-Verlag 2005

T. Kalbacher (&) Æ W. Wang C. McDermott Æ O. Kolditz Center for Applied Geoscience, University of Tuebingen, Tuebingen, Germany E-mail: [email protected] T. Taniguchi Faculty of Environmental Science and Technology, Okayama University, Okayama, Japan

ORIGINAL ARTICLE

Development and application of a CAD interface for fractured rock

Abstract The implementation of numerical methods involving geometrical description of geosystems with the application of numerical meshes is of vital importance for investigation and modelling of complex, coupled, natural, or man-induced processes. Innovative application of CAD-Systems in combination with computer aided modeling of such systems provides a user-friendly powerful tool for the often complex geometrical descriptions necessary for comprehensive modeling of the dominant in situ processes. Here the development of such a CAD interface is presented, based on GOCAD, coupled with the software engineering of the finite element program RockFlow/RockMech and the implementation of polylines to describe complex geometrical structures and surfaces. This program allows the calculation of flow, transport, and deformation

Introduction The implementation of numerical methods involving geometrical description of geosystems with the application of numerical meshes is of vital importance for investigation and modeling of complex, coupled, natural or man-induced processes. Innovative application of CAD Systems in combination with computer aided modeling of such systems provides a user-friendly powerful tool for the often complex geometrical descriptions necessary for comprehensive modeling of

processes in one or more fluid phases in fractured and porous media. Such simulation is becoming increasingly important, especially in the fields of waste disposal, geothermal energy, and groundwater supply management. Here an example of the application of the CAD interface to RockFlow is presented for the modeling of flow during artificial stimulation (a 300-day long pump test) of a potential geothermal reservoir in the deep continental borehole (KTB) program at Windischeschenbach, Germany. Keywords CAD interface Æ RockFlow Æ Polyline and surface objects Æ Fractured rock

the dominant processes. In this paper, the development of such a CAD interface is presented in Sect 2, coupled with the software engineering of the finite-element program RockFlow/RockMech. This program allows the calculation of complex flow, transport, and deformation processes of one and more fluid phases in fractured and porous media, Kolditz (1999), Kolditz et al. (2003). Such simulation is becoming increasingly important, especially in the fields of waste disposal, geothermal energy and groundwater supply management. Section 3 deals with the conversion of complex geometries into polylines

1018

and surfaces and the software engineering implementation. Correspondingly, Sect 4 presents an example of the application of the CAD interface to RockFlow for modeling flow during artificial stimulation (a 300-day long pump test) of a potential geothermal reservoir comprising several 3D fracture elements in the deep continental borehole (KTB) program at Windischeschenbach, Germany.

CAD interface The development of a long-term concept for a windoworiented user interface for RockFlow is presented. In terms of preprocessing, an environment similar to CAD (computer aided design) has been created for RockFlow. Similarly, the creation of interfaces between existing CAD programs and RockFlow is currently under development. Data and data preparation CAD systems are suitable for many scientific applications as a preprocessor to the data preparation before mesh generation. A significant advantage of preprocessing by CAD is the ability to blend and combine different datasets with one an other for the model generation. This allows the preparation of well matching geometry for the model. In addition, visual control of the data and their spatial distribution, the possibility to add, remove, and modify geometric data, and the use of standard tools for the geometric and topological processing (e.g. create new polyline from points or close polyline) may also be easily applied. The application of CAD Systems in combination with computer aided modelling of fractured rock is becoming increasingly important for the simulation of the flow, and transport processes, especially in the fields of waste disposal, geothermal energy, and groundwater supply management. Fractured rock systems in subsurface are very complex systems described by many parameters with spatial data distributions, and their realistic modelling is based on precise calculations, which require precisely located data. The advantage of CAD Systems is that they allow the editing and management of such definite geographically located data. Data for fractured rock models include the geometry of the domain to be considered and its fractures. Boundary conditions, initial conditions and material groups (specifying the material properties of rock) also need to be taken into account as additional geometric information. CAD is used to prepare the fractured-rock model geometry for the mesh generation according to

the domain, boundary conditions, initial conditions, and material groups. Such a preparation step could include the division of a closed polygon into two polygons in accordance with the necessary boundary conditions. These preparations are upgraded and also restricted by the input prerequisites of the mesh generator. An upgrade of the mesh generator for example is the possibility to define a point density along polylines. Restrictions arise with the convention of the data input to the mesh generator, such as that the generator may only handle closed polylines for the mesh generation. Therefore, the interaction of CAD and mesh generator decides which and how many geometric/topological methods need to be developed after mesh generation to optimally prepare the fractured-rock model for the process simulation. GOCAD The visualization software Geological Object Computer Aided Design(GOCAD) is a product of the GOCAD consortium (http://www.ensg.inpl-nancy.fr/GOCAD/). GOCAD is a geologically object -oriented CAD software that provides various tools to construct a range of models. These models are particularly for application in geology, geophysics, and reservoir engineering. GOCAD also offers tools to generate and modify triangulated meshes. It can also handle hexahedra meshes. The following part of this section illustrates the method of data preparation using GOCAD for a simplified 2.5D fractured rock (Moenickes 2004). An example of the application of this method is given for the KTB program at Windischeschenbach, Germany, which provides datasets from a pilot borehole and a

Fig. 1 Borehole data of a fractured rock model

1019

main borehole situated some 200 m from the pilot borehole. The first step is to input and visualize the data sets to be used. In our example, (Fig. 1) those datasets include borehole coordinates (x, y, z), dip and strike of fractures, and a hydrogeological/geophysical interpretation of the fracture system, allowing the dominant features to be identified. Once the model space has been created, in this case a simple cuboid, the fractures are introduced as plain surfaces. The example in Fig. 2 shows seven simple, flat and complete surfaces, which would further be subdivided by a higher data density (e.g. additional drill holes). As more information is added to the model so the complexity of the model increases, i.e. the complexity of the model is a function of the information density within the model. With the help of these surfaces, GOCAD is able to produce lines defining the intersection of the various surfaces, the ‘‘intersection lines’’. These intersection lines are later applied in the generation of finite elements used to solve the processes being investigated. During the mesh generation, the intersection lines are applied to ensure an exact fit of the finite elements of one surface with the next. Figures 3 and 4 show the difference between a method not considering the intersection lines and the method we apply which considers the intersection lines. Notice particularly the difference between the linkage of the various fracture planes; in Fig. 3 it can be seen that there is a poor matching with the intersecting planes regarding the elements, whereas in Fig. 4 the match is perfect. GOCAD creates the surface using variations of a Delaunay-triangulation algorithm, which enforce adding points or homogeneous triangles. An example of such a procedure is given in Fig. 5; here the triangula-

Fig. 2 Measured fractures which are introduced as surfaces inside the model area

Fig. 3 Triangulation without consideration of intersection lines

tion method has been used to create a new surface from a closed polyline. To enable this approach, the polylines of the surface edges and the intersection lines must be combined to the smallest possible closed polylines (Fig. 5). This approach functions where the entire surface domain has been cut by an intersection line, but requires a manual postprocessing at intersection lines which only partly cut the domains. Using this method, the original 7 planes were subdivided into 88 sub-planes (closed polylines) as illustrated

Fig. 4 Triangulation with consideration of intersection lines

1020

Fig. 7 Schema GOCAD-GOCAD2RF- RockFlow (Geometry Objects, Mesh Objects)

Fig. 5 2.5D fractured rock model with a triangulated mesh generated within GOCAD

in Fig. 6. A point density for defining the mesh generation can be specified along the individual lines (intersection lines, boundary lines) within GOCAD prior to mesh generation. Shared polylines are described repeatedly, and exactly those parts need identical inter-points to enable the sharing of the finite element nodes along the intersection lines. If such points do not correspond, then internal boundaries arise inside the model area, which can cause errors in process simulation later. Due to the finite element model requirements (i.e. maximum number of elements, mesh density, local mesh refinements, mesh quality, etc.), such nets are sometimes not completely suitable for the process simulation and

Fig. 6 The 88 sub-planes of the 2.5D fractured rock model

require a revision. This revision will be described in the following sections.

GOCAD2RF—Interface The GOCAD2RF is an interface for RockFlow into the world of CAD (Fig. 7). The object-oriented data concept of RockFlow is generally designed for numerical simulation of multiphase and multicomponent processes in porous-fractured media. The GOCAD2RF converter uses the data structures of RockFlow for generating the output (used RockFlow-Objects: ELEMENTS, NODES, MATERIALS). This converter allows the transformation of input data comprising polylines (curves), triangulated surfaces, and tetrahedra meshes (solid). The polylines are written as RFM files (Sect. 2.4) and RFI files are produced from the triangle meshes or the tetraheda meshes (Sect. 2.4). The interface is in the development stage and works at the moment with a self-sufficient text-based GUI of its own. The advantage of the development of an interface instead of own CAD tools simply lies in the shorter development time for the use of professional CAD tools. The possibility to manage geological, geophysical, and other spatial allocated properties already in the CAD system is an additional advantage of GOCAD. However, it is also necessary to handle the weaknesses of the selected CAD program. To date, the difficulties we have had in the application of the CAD systems to the finite element process simulations include: – double nodes (neighboring edges that are supposed to be joined by sharing a common node, but have separate defined, discrete nodes) – distorted elements with no area inside (all nodes along a line)

1021

It is therefore necessary to clean-up the meshes within the converter. The GOCAD2RF converter looks for double nodes and unnecessary distorted elements, deletes them, and rebuilds the meshes. The polylines of the domain are transformed unchanged.

methods, which can be applied before, during, and after the mesh generation. For the fractured-rock model GSMH is particularly interesting because it offers the possibility of a point-obtained mesh refinement or mesh coarsening (Fig. 14, 17, 18).

RockFlow RFI and RFM file

GeoObj

RFI files contain the information of the geometry, topology, and material groups of the already meshed model area. With the application of the material groups to define different material zones in the volume to be modeled, the material properties are connected directly with the elements. An additional advantage of GOCAD is that it can assign properties to the geometric objects and therefore it enables the allocation of material properties to elements. However, the mesh generation within GOCAD and the transformation of already meshed domains also have some disadvantages: particularly some complex challenges can cause numerical errors during the mesh generation or other processes inside the CAD. These effects are hard to locate, analyze, comprehend, and difficult to remove afterwards. Such numerical errors can be cleared directly by direct access to the source code of the mesh generator or additional specific geometric methods, which are better suited for the corresponding applications can be developed. For this reason the GOCAD2RF interface can convert pure polyline information. The RFM file is the RockFlow input file for the polylines. RFM files contain information of the geometry and topology on the unmeshed model area and are suitable for the mesh generator GSMH (http:// www.geuz.org/gmsh/) as an input file (Sect 3.22). Another advantage of the conversion of pure polylines information together with the mesh generation by using GSMH is the possibility to develop geometric

Fig. 8 Geometric objects

The geometric objects (GeoObj) show a hierarchical structure, which orientates itself at the topology. The domain can be described by points, polylines, surfaces, and volumes (Fig. 8). Polylines Polylines are universal geometric objects, which can be applied for several purposes, such as domain boundaries in 2D, line structures in 2D (e.g. fractures), line structures in 3D (e.g. boreholes, flow channels), and contour lines. Polylines can be used for very different line-related information, such as initial (IC) and boundary conditions (BC), source terms (ST), and output data (OUT). Material groups can be specified by closed polylines. Additionally, we can relate polylines to geometric and topologic objects such as nodes (NOD) and finite elements (ELE) (Fig. 9). Polylines are based on points and build the basis for surfaces (Fig. 9). As stated above, we use polylines for meshing fracture networks (MshObj).Due to the universality of polylines, the design and implementation of an object (PLY) is given below. Data constructor A polyline consists of an ordered array of nodes (*point array) defined by their points (x, y, z), which are the basic geometric data of PLY. A tolerance (eps) can be specified to define a geometric environment along a polyline. This is useful

Fig. 9 Object relationships of polylines

1022

identify geometric objects (such as points), which are near a polyline. Additionally, polyline points can posses properties (values) to specify boundary conditions etc. Closed polylines are automatically recognized if the first and last points are identical. Using the polygons defined in this manner, the material groups can be determined. Every polyline has a name to refer to for access. The corresponding data construct (ANSI-C/C++) of the PLY object is given below:

typedef struct { char *name; POINT_ARRAY *point_array; long number_of_points; double eps; int closed; int number_this; int type; int subdivision; int material_group; int Index; int selected; }POLYLINE;

Construction and destruction of polyline instances is done by following functions using names.

Fig. 10 Data organization of polylines - object list

Polyline applications Here two applications for polylines are considered, (1). identifying points near a polyline and (2) meshing of polygons. Identifying points near a polyline

Data organization

We consider an example of a relationship between the PLY object and geometric objects such as nodes NOD. Here the grid points of an existing mesh are searched for

We use a user-defined double-chained list for the organization of polylines.

typedef struct { char *list_name; List *list_items; long number_of_polylines; } LIST_POLYLINES;

A graphical representation of an object list is given in Fig. 10. Each item of such a list (i.e. polyline) can be identified by a polyline name. A typical set of functions of lists for construction/destruction, insertion and removal of items from the list as well as stepping through the list is given below.

LIST_POLYLINES *create_list_polylines (char *name) void destroy_list_polylines(void) long insert_polyline_list(POLYLINE *polyline) long remove_polyline_list (POLYLINE *polyline) void *get_list_polylines_next (void)

Fig. 11 Points close to a polyline (=50 m)

1023

After removing double nodes from the composite grid (see above), a finite element mesh is obtained. The polyline-based meshing strategy allows easy mesh manipulations (refinements/coarsening) around selected polyline nodes. In the present study, the surface is set as the last elementary entity in the priori work to mesh a fracture network. It is governed by a polyline entity that is defined in the last section. The definition of the surface object is given below:

Fig. 12 Intersection polylines of a fracture surface (2D view)

those which are close to a polyline (Fig. 11). This can be used for instance to assign source terms to a borehole. long *GEOGetPointsClose2Polyline(char *name,long *number_of_nodes) Meshing of polygons The second polyline application deals with mesh generation for fracture networks. The networks of intersecting planes can be represented by sets of polylines. Based on the CAD system introduced in Sect. 1, intersection polylines of all fractures can be produced. Figure 12 illustrates these intersection traces for one selected fracture. These polylines are building surfaces, which are geometric objects described in the next section. There is a polyline editor (GS PLY) available in RockFlow for polyline handling (creation, drawing, editing,...).

class Surface {public: Surface(const int Number Of Polylines);  Surface(); void SetPolylineName(const char* PlyName, const int PlyIndex, const int Size); char (??) GetPolylineName(const int PlyIndex); const int GetNumberofPly() {return number_of_polylines;} const int GetSizeOfName(const int Indx) {return SizeOfName[Indx];} void FreeMem(); private: char **NameOfPolyLines;//Handle of polylines int number of polylines; int* SizeOfName; };

The handles of the polylines are defined by polyline names. Therefore, the surface object contains the message of polyline names and methods to manipulate the names. To illustrate this, two polylines, A and B, are depicted in Fig. 13 to define the surfaces. If A and B are selected together they form the surface object with their names, as shown in Fig. 13a. Otherwise, if A and B were selected to form the surface object and then A was selected to form the object,

Surfaces Now, sub-surfaces enclosed by these polylines can be triangulated using the meshing algorithm (GMSH).

Fig. 13 Ruled surface by polylines

Fig. 14 Meshing of polylines for a fracture surface (2D view)

1024

Fig. 15 2.5D fractured rock network with an increasing mesh density to the drill hole

we would have two surfaces as shown in Fig. 13b. In the former case, the names of two polylines, A and B, are accepted by member NameOfPolyLines as an instance of the surface object. In the latter case, we have another instance of the surface object, member NameOfPolyLines, which keeps the name polyline A, besides the instance created by the names of polyline A and polyline B. The polylines that include interface lines could be the boundaries of any patches in the facets of fracture network. The instances of the surface object are put in a list to be used by the mesh procedure (see Figs. 14, 15)

Application to the modelling of fractured rock (stimulation of a potential geothermal reservoir) Quantification of thermal, hydraulic, and geomechanic coupled processes is very important for the understanding and characterization of natural as well as man-

induced processes in geosystems. The implementation of numerical methods involving geometrical description of in situ fracture systems on hand in0 numerical meshes is an important step in the investigation of coupled processes in deep crystalline rock. Experimental methods for the tomographical investigation and three- dimensional characterization of fractured material are presented in McDermott (2003). Here an example of the application of geoscientific understanding of fracture networks and the implementation of numerical meshes as described above is provided for the deep continental borehole (KTB) program at Windischeschenbach, Germany, comprising a pilot hole to a depth of 4 km and a main hole with a depth of 9.1 km. Modelling techniques are being used for the general analysis of the flow, transport, and deformation processes in the fracture systems and for the predictive experimental design of investigation programs to further define the characteristics of the system. The information and experience gained will then be of importance for the long-term predictive behaviour of deep fracture systems, for instance, in the field of geothermal reservoirs or waste depositories. Accurate mesh representations of the fracture networks coupled with finite element modelling techniques are necessary for prognostic modelling of the in situ natural conditions and possible applications later. Concerning the KTB program, a great deal of literature is available on the geological setting of the KTB boreholes and the fracture geometry, type, and nature. Emmermann and Lauterjung (1997) provided a lead paper to a special edition in the Journal of Geophysical Research, concerned with research at the KTB site. The work of Hirschmann and others in several KTB reports and particularly Hirschmann et al. (1997) coupled with geophysical investigations of the borehole site by Harjes et. al. (1997) provide a necessary basis for the development of a regional hydraulic model comprising principally a fracture network to describe flow and transport. Figure 16 presents the structural geology in terms of the dominant fracture zones of the borehole site taken from Hirschmann et al. The fracture zones presented can be taken as being fluid bearing providing the necessary contrast for the geophysical detection and potential fluid pathways for flow and transport, forming the basis of the fracture network model. The accuracy of the location of the fracture zones, particularly in the vertical orientation, is of the order of a few tens of meters, totally acceptable given the scale of the entire fracture network we are dealing with (here at least 20 km · 20 km · 10 km depth). In addition, as pointed out and modelled by Kessels (2000), the fracture zones comprise several interacting smaller fractures and are

1025

Fig. 16 Model fractures superimposed on three-dimensional Geological Model after Hirschmann 1996

not localized breaks in the rock, rather shear zones which can be of the order of 100 m+ in thickness. As a first approximation of the hydraulic system, to form a basis for the later investigation of the geomechanical processes and transport behaviour, the structural information was combined to provide a threedimensional fracture network model of the dominant features within the fracture system. Numerically the fracture systems are represented by an equivalent med-

Fig. 17 Three-dimensional fracture network

Fig. 18 Meshing and head distribution on a single fracture plane within the network intersected by the borehole

Fig. 19 Meshing and head distribution on the entire fracture network, head distribution for steady state conditions

1026

Fig. 20 Modelling of field pump test

ium combining the properties of the fractures and the surrounding matrix. Details of the steps included in the modelling of the system are presented in (McDermott and Kolditz 2004) The boundary conditions of the model are hydrostatically defined constant pressure boundaries, taking into account the density of the saline fluid found at such depths in the KTB site. Figures 16 and 17 provide a comparison of the structural geology and the features considered to be dominant fluid pathways (Fig. 16), and the fracture network model (Fig. 17).

For the purpose of a more detailed modelling approach, the entire fracture zone was reduced to an area of 4 km · 4 km · 6 km deep, and the meshing density was increased around the vicinity of the drill hole to enable a more accurate numerical representation of the system. An extract of one of the fractures intersected by the borehole under steady-state flow conditions is presented in Fig. 18, illustrating the increase in density of the meshing around the areas of interest. Since June 2002, a pump test has been carried out at the KTB site. This pump test involves the extraction of geothermal fluids from the fracture network, and the measurement of the drawdown in the borehole coupled with the rate of the extraction. Modelling this pump test allows the determination of relevant flow parameters important for latter predictive modelling of the hydraulics of the system. Figure 19 presents the threedimensional meshed fracture network of the natural system, here with steady state conditions modelled and the head distribution within the network, and Fig. 20 illustrates the comparison of the modelled and measured data for transient conditions. It can be seen from the accurate fitting of these curves that the three-dimensional fracture network model provides an important tool for the predictive modelling of the fracture system. Acknowledgments The authors wish to thank particularly the Deutsche Forschungsgemeinschaft (DFG) and the GeoForschungsZentrum, Potsdam (GFZ). The application of RockFlow to modelling the pump data from the KTB site was made possible through the Deutsche Forschungsgemeinschaft (DFG) in the framework of the SPP-ICDP program. Significant financial support for carrying out the pump test was also received from the GeoForschungsZentrum, Potsdam (GFZ). In addition, the data from the KTB field pump test was provided by the KTB science team.

References Emmermann R, Lauterjung J (1997) The German Continental Deep Drilling Program KTB: overview and major results. J Geophys Res 102(B8):18179– 18201 Harjes HP et al (1997) Origin and nature of crustal reflections: results from integrated seismic measurements at the KTB superdeep drilling site. J Geophys Res 102(B8):18267–18288 Hirschmann G (1996) The structure of a Variscan terrane boundary: seismic investigation-drilling-models. Tectonophysics 264:327–339

Hirschmann G, Duyster J, Harms U, Konty A, Lapp M, de Wall H, Zulauf G (1997) The KTB superdeep borehole: petrology and structure of a 9-km-deep crustal section. Geol Rundsch 86(Suppl):S3–S14 Kessels W (2000) Elastische Eigenschaften von Klu¨ften und Interpretation hydraulischer Tests- Untersuchungen im Rahmen des KTB- Projektes, 3. Workshop Kluft- Aquifere ,, Gekoppelte Prozesse in Geosystemen‘‘,Institute fu¨r Stro¨mungsmechanik und Elektron. Rechnen im Bauwesen der Universita¨t Hannover

Kolditz O (1999) Computational methods in environmental fluid mechanics, ISBN 3-540-42895-x Springer-Verlag, Berlin Heidelberg New York, ss378 Kolditz O, de Jonge J, Beinhorn M, Xie M, Kalbacher T, Wang W, Bauer S, McDermott CI, Kaiser R, and Kohlmeier M (2003) RockFlow manual, RFD input description, version 3.9. Center for Applied Geosciences, University of Tu¨bingen and Institute of Fluid Mechanics, University of Hannover

1027

McDermott CI, Kolditz O (2004) Hydraulic-geomechanical effective stress model: determination of some discrete fracture network parameters from a pump test and application to geothermal reservoir modelling., Twenty-ninth workshop on geothermal reservoir engineering, Stanford University, Stanford

McDermott CI, Sauter M, Liedl R (2003) New experimental techniques for pneumatic tomographical determination of the flow and transport parameters of highly fractured porous rock samples. J Hydrology (in print)

Moenickes S (2004) Grid generation for simulation of flow and transport processes in fractures-porous media. Dissertation, Institut fu¨r Stro¨mungsmechanik, Universita¨t Hannover, Bericht Nr. 68/2004