Carl Gable, Harold Trease and Terry Cherry

Automated Grid Generation From Models of Complex Geologic Structure and Stratigraphy

Abstract

The construction of computational grids which accurately reflect complex geologic structure and stratigraphy for flow and transport models poses a formidable task. With an understanding of stratigraphy, material properties and boundary and initial conditions, the task of incorporating this data into a numerical model can be difficult and time consuming. Most GIS tools for representing complex geologic volumes and surfaces are not designed for producing optimal grids for flow and transport computation. We have developed a tool, GEOMESH, for generating finite element grids that maintain the geometric integrity of input volumes, surfaces, and geologic data and produce an optimal (Delaunay) tetrahedral grid that can be used for flow and transport computations. GEOMESH also satisfies the constraint that the geometric coupling coefficients of the grid are positive for all elements.

GEOMESH generates grids for two dimensional cross sections, three dimensional regional models, represents faults and fractures, and has the capability of including finer grids representing tunnels and well bores into grids. GEOMESH also permits adaptive grid refinement in three dimensions. The tools to glue, merge and insert grids together demonstrate how complex grids can be built from simpler pieces. The resulting grid can be utilized by unstructured finite element or integrated finite difference computational physics codes.

Grid Generation and GIS

Grid generation is a broad field with applications in aerodynamics, material science, biology, earth science, chemistry and physics, for example. For any problem where there is a need to represent continuous functions with discrete points, or a system of equations needs to be solved, the first step is often to construct a discrete tessellation of the region. In applications where complex geometries do not need to be represented, regular orthogonal grids are simply constructed. When it is important to represent complex geometries, multiple materials and properties, and optimize the shape of the elements used for tessellation, automatic grid generation tools are necessary.

Grid generation tools must provide the means to build and represent volumes from surfaces, distribute nodes within a volume, define a connectivity of nodes that is compatible with the computational tools being utilized for modeling, and provide a means of assigning initial and boundary conditions. In addition, grid generation must insure the quality of the grid to insure the accuracy and stability of the physics being modeled.

There are a wide variety of geologic applications where accurate representation of complex engineering systems (wells, tunnels, reservoirs) and geologic structure and stratigraphy (layering, folding, domes, faulting) is critical to producing accurate numerical models of fluid flow and mass transport. Oil and gas reservoir production, groundwater resource development, hazardous waste site characterization and remediation, and nuclear waste disposal in a geologic repository are examples of the areas where modeling must be used to predict the long term behavior of a system. In all the systems, grid generation is a key link between the geoscientific information systems (GSIS, Turner, 1991) and the numerical models. They must capture complex geometry and insure the computational grid is optimized to produce accurate and stable solutions.

Issues dealing with the general approach and methods of integrating GIS and various process modeling studies are presented in Fisher (1993) and Moore (1993). Although many of the GSIS modeling methods and tools (Fisher, 1993) claim to be developed as solid geometry modeling tools they often lack open architecture or the ability to export the models they create. Consequently, their use as a preprocessor for physics packages is limited.

This work identifies the data needs for a particular grid generation tool, GEOMESH (Gable, 1995), discusses data that is required for GSIS to be useful as a preprocessor for grid generation and physical modeling, and presents examples of grid generation projects utilizing GEOMESH. This discussion concentrates on integrating GSIS data into deterministic numerical models and discusses issues in producing optimized grids which accurately maintain geometries defined in GSIS. Deterministic 3D geometry modeling assumes that the goal is to accurately represent the GSIS model. It does not insure the quality of input data or how errors in input data affect the GSIS model.

Structured and Unstructured Grids

Grid generation is far more than the tessellation of space with polyhedra. There are two major classes of grids, structured and unstructured. In the simplest case, a structured grid is a logically rectangular array of points with a well-defined and regular relationship of each point with its neighbor. For example, a typical 2D finite difference grid (Figure 1) can represent a function f(i,j) in a regular pattern where the neighbor to the right is f(i+1,j) and the neighbor above is f(i,j+1). A 2D rectangular grid will have four connections to any interior node, however they are not necessarily orthogonal. Irregular structured grids can be thought of as regular orthogonal grids but with all of their connections made of rubber bands so the grid can be stretched to a complicated shape as long as none of the connections are broken and none of the grid elements are turned inside out. This type of gridding scheme has many limitations, primarily with the complexity of geometries that can be represented. Numerical schemes for solving equations on these grid often have problems with accuracy, as the stretched grid becomes nonorthogonal.

Unstructured grids, often used in finite element computations, (Figure 1) offer far more flexibility in representing complex geometries. Grids can be constructed from a variety of elements (triangles, quadrilaterals, tetrahedra, hexahedra). This flexibility comes at a higher cost. The spacial relationship of one element to another must be explicitly stated since there is no implied or logical relationship of one element to its neighbor. The number of connections a particular node has can vary greatly.

GEOMESH Grid Generation

GEOMESH is a software tool for automatically producing unstructured finite element grids tuned to the special needs of geologic and geo-engineering applications. The code provides for 2D and 3D grids and includes elements that are triangles, quadrilateral, hexahedral and tetrahedral. Capabilities have been expanded to improve grid generation for complex geometries and for providing more choices in designing a grid. At the same time, old routines have been tested and solutions found for degenerate cases. The design of the code is modular, allowing for flexibility and consistency.

The core functions of GEOMESH utilize the X3D grid generation package developed at Los Alamos National Laboratory. GEOMESH developed out of a need for accurate and automated grid generation for 3D modeling of subsurface porous flow and transport. Since the grids represent the geology being modeled, the accuracy of the grid directly affects the accuracy of the model. It was also found that grid generation was tedious, time consuming, and prone to errors, especially for models with complex structures such as faults and stratigraphy such as pinch outs and layer truncations.

Automated grid generation algorithms not only streamline grid construction, they offer more flexibility and consistency. As input constraints change, and as better data sets become available, these data are easily incorporated into the computational mesh. Automated grid generation can also be used to produce coarsened grids for preliminary calculations and refined grids for final, high resolution calculations. The grids produced by GEOMESH are widely applicable and can be used by any numerical algorithm that can utilize unstructured grids. They are not specific to any particular computational code

Defining a geometry is the first step in grid generation and is where the grid generation tool must interface with GSIS systems. The GEOMESH program can input geometry definitions in terms of their bounding surfaces. The bounding surfaces may be defined by triangular irregular networks (TINS), multi-dimensional analytic function representations such as nonuniform rational B-splines (NURBS), or by regular grids f(i,j) where f is the grid elevation and (i,j) are regular spaced x and y coordinates. Regardless of the means of defining surfaces, it is important that the geometry defined by the surfaces is unambiguous. Surfaces must form a closed volume, the normal to each surface must be consistently defined, the line defined by the intersection of two planes must be unique and there cannot be gaps where the interface of one surface joins another surface.

An alternative to surfaces are representations that already define volumes. For example, there are some GSIS packages that produce a mesh as output. However, the mesh may not be usable for computations. In this case the grid generation tool must fix the problems to allow computations while maintaining the geometry of the GSIS model. This utilizes grid generation tools for insuring grid quality.

GEOMESH grid generation uses three criteria to insure grid quality. They are: (1) the final grid preserves the input geometry model, (2) the grid is Delaunay (Figure 2) and (3) all coupling coefficients are positive. It is beyond the scope of this work to give the mathematical details of Delaunay triangulation. However, an intuitive definition of a Delaunay triangulation is that the circumscribed circle of any triangle contain the three points of the triangle and no other points. The extension to three dimensions is that the circumscribed sphere of any tetrahedra contains the four points of the tetrahedra and no other points. There are special cases and degeneracies when the above description is not absolutely correct, but these cases are not important to this discussion. The third constraint, that the coupling coefficients related to grid geometry produce a semi-positive definite matrix (Trease and Dean, 1990, Gable et al. 1995), insures that the solutions to nonlinear flow and transport will be stable and accurate. In addition, there are no negative transmissibilities when solving porous flow problems.

Many of the problems associated with importing data into grid generation tools results from insufficient GSIS output to completely characterize the solid geometry model. This may be a consequence of ambiguities in the output geometry model. Some GSIS packages produce beautiful graphical representations of the geometry model, but do not provide utilities or an open architecture for exporting the solid geometry model.

Examples

The first example, (Figure 3) shows a 3D grid produced from a GSIS 3D geometry model. The GSIS model is composed of hexahedral (8 node, 6 sided elements). However there are many degenerate elements that must be eliminated. For example, when a layer pinches out the GSIS model retaines zero volume elements as part of the output. While this produces a correct representation of the geometry, physics codes cannot compute on zero volume elements. Furthermore, the GSIS output produces elements without regard element shape or grid quality.

In order to produce a grid for computation from the GSIS output the following steps are taken: (1) extract the 3D geometry model, (2) convert the hexahedral elements to tetrahedral elements by subdividing, (3) eliminate the zero volume and degenerate elements, (4) reconnect the grid to insure a Delaunay grid while never disrupting material interfaces, (5) subdivide elements to insure positive coupling coefficients and(6) output node coordinates, connectivity, coupling coefficients and a list of boundary and material nodes. This grid is then used for computation of unsaturated porous flow and contaminant transport (Robinson et al. 1995).

The second example, Figure 4, is a 2D cross section extracted from the 3D model shown in Figure 3. All of the material interfaces in the geometry model are maintained, and the computational grid quality has been insured by producing a Delaunay grid with all positive coupling coefficients.

These cross sections have a vertical extent of approximately 700 meters with very thin layers of a few meters. The automatic grid quality algorithms made it simple to produce these grids with very large differences in length scale from one layer to another.

A final example demonstrate the ability to add or merge one grid with another. This provides a powerful tool for building complex grids from simpler building blocks. In this example (Figure 5), a mesh representing a curved well is inserted into a background grid. This mesh merge utility can be used in a wide variety of other applications. For example, instead of using grid refinement to produce a high resolution grid, one can produce two grids, a course grid that covers a large volume and a high resolution grid that covers a smaller volume and merge them together to produce a single grid with high resolution in a specific region.

Other features not demonstrated by these examples are interpolation functions and grid smoothing. Interpolation of node attributes (i.e. temperature, pressure) from one grid to another is useful when applying grid refinement algorithms so that the refined grid also has estimates of attribute values. Grid smoothing algorithms allow grid optimization by letting nodes move rather than by changing connections. This method preserves material interfaces fixing the position of interface nodes while allowing nodes in the interior to move until grid quality constraints are achieved. This method can also be used to increase grid resolution where necessary by allowing the grid to migrate to zones where finer gridding is needed.

Future Outlook

There are a number of GSIS software packages available whose market is mainly oil and gas reservoir modeling. However, most of these packages suffer from a lack an open architecture. These tools often provide means of building complex representations of complex geology, however they do not provide open, accurate and flexible means of exporting 3D geometry models once they are created. Software packages that produce impressive 3D visualizations are a dead end for getting the information into physics modeling packages. This severely limits their utility as a tool in a fully integrated approach to computational physics modeling.

As tools for various aspects of an integrated modeling approach mature and become more sophisticated, the problems that the user often faces are issues of data compatibility. Well logging software may not talk to the GSIS software. The GSIS software may not interface with the grid generation tools. The grids produced are not tuned to the special needs of a physics code. The results of physics modeling cannot be easily used in reinterpretation and modification of the GSIS model. As each of the separate pieces mature, more effort will need to be spent creating seamless connections amongst them. The ability to incorporate and utilize a wide variety of data into modeling and interpretation will have a strong impact towards achieving the goal of fully integrated modeling. Automated grid generation is one piece of this puzzle.

References

Abdou, M. K., H. D. Pham and A. S. Al-Aqeeli, (1993), Impact of grid selection of reservoir simulation, Journal of Petroleum Technology, 45(7), pp 664-669.

Fisher, T., (1993), Use of 3D geographic information systems in hazardous waste site investigations, Environmental modeling with GIS, ed. Michael F. Goodchild, Bradley O. Parks and Louis T. Steyaert, Oxford University Press, pp 238-247.

Gable, C. W., T. A. Cherry and H. E. Trease and G. A. Zyvoloski, (1995), GEOMESH grid generation, LA-UR-95-4143.

Gable, C.W., George Zyvoloski, (1994), Site Scale Modeling of Radionuclide Transport At Yucca Mountain, NV: Grid Generation and Reactive Tracers, LA-UR-94-1041.

Khamayseh, A. and Andrew Kuprat, (1995), Anisotropic Smoothing and Solution Adaption for Unstructured Grids, LA-UR-95-2205, International Journal for Numerical Methods in Engineering, (submitted).

Moore, I. D., A. K. Turner, J. P. Wilson, S. K. Jenson and L. E. Band (1993), GIS and Land-Surface-Subsurface Process Modeling, Environmental modeling with GIS, ed. Michael F. Goodchild, Bradley O. Parks and Louis T. Steyaert, (Oxford University Press, pp 196-230.

Robinson, B. A., A. V. Wolfsberg, G. A. Zyvoloski and C. W. Gable, (1995), An unsaturated zone flow and transport model of Yucca Mountain, in press.

Trease, H. E. and S. H. Dean, (1990), Thermal Diffusion in the X-7 Three-Dimensional Code, Proceeding of the Next Free-Lagrange Conference, Jackson Lake Lodge, Wyoming, June 3-7, Springer-Verlag Press, Vol. 395, pp. 193-202.

Turner, A. K., (ed.) (1991), Three Dimensional Modeling With Geoscientific Information Systems, Nato ASI Series C: Mathematical and Physical Sciences, Vol. 354, Dordrecht: Kluwer Academic Publishers.

Author Information

Carl W. Gable

mail MS F665, Los Alamos National Laboratory, Los Alamos NM 87545

email gable@lanl.gov

voice 505-665-3533

fax 505-665-3687

Last Modified: 03:15pm MST, January 05, 1996