A Comprehensive Modeling Environment for the Simulation of Groundwater Flow and Transport

Steven J. Owen and Norman L. Jones

Engineering Computer Graphics Laboratory, Brigham Young University, Provo, UT, 84602, USA

Jeffrey P. Holland

U.S. Army Engineer Waterways Experiment Station, Vicksburg, MS, USA

Abstract

A comprehensive graphical modeling environment has been developed to address the needs of the computer simulation of grounwater flow and transport. The Department of Defense Groundwater Modeling System (GMS), developed at the Engineering Computer Graphics Laboratory at Brigham Young University, is part of a multi-year project funded through the Department of Defense, Department of Energy and Environmental Protection Agency. GMS is a graphically based software tool providing facility through all aspects of the groundwater flow and transport modeling process. Facilities include geometric modeling of hydrostratigraphy, two and three dimensional mesh generation, graphically based model input for specific flow and transport codes, interpolation and geostatistics as well as complete three-dimensional scientific visualization.

Introduction

In recent years, there has been significant development in the area of computer modeling of groundwater flow and transport. Many powerful computer models have been written that model groundwater systems of considerable magnitude and complexity. In spite of the advances made in analyzing such problems, these software tools often remain inaccessible to the majority of the groundwater modeling community because of the overhead involved in learning a new code and writing voluminous input files requiring tedious formatting. Effective use of a groundwater modeling code also requires an in-depth understanding of the geometric relationships between hydrostratigraphic units, natural and imposed boundaries along with sources and sinks. A graphical representation of all of the components of the system is essential to effectively model these relationships.

The Department of Defense Groundwater Modeling System (GMS) is a software tool designed to address these issues. It is currently being developed by the Engineering Computer Graphics Laboratory (ECGL) at Brigham Young University. It is funded through a consortium of US government agencies, including the Department of Defense, Department of Energy and Environmental Protection Agency with principal direction provided through the U.S. Army Engineer Waterways Experiment Station (WES) at Vicksburg, Mississippi.

The main purpose of GMS is to provide a complete tool for the groundwater modeler. It is designed to provide tools throughout all aspects of the modeling process, some of which include geometric characterization of earth masses, geostatistical analysis, finite element and finite difference mesh generation, model input for specific flow and transport analysis engines as well as complete three dimensional visualization of results. This paper covers some of the main components of GMS, addressing how this new tool is applicable to the groundwater modeling process.

User Interface

GMS is currently under development in both the UNIX X-Windows and Microsoft Windows environments. A common user interface was developed on both of these platforms. The user interface is designed to take full advantage of graphical interaction and editing techniques.

The GMS user interface is divided logically according to task. Commands are divided into nine, consistently designed modules. Selecting a new module entails choosing an appropriate icon from an array of nine icons in a tool pallette. The interface changes to reflect a new set of tools and menu selections that are applicable for new module. The modules are briefly defined in Table 1. By dividing the program up into task-oriented modules, the user is effectively shielded from unnecessary complexity of other tasks in the system while focusing on specific procedures at hand.

Site Characterization

The first phase of a groundwater modeling project usually begins with site characterization. The borehole, TIN and solids modules within GMS provide tools for modeling the three dimensional stratigraphy of a region. Jones and Wright [1] describe some of the solid modeling techniques provided with GMS and how they are applicable to site characterization. These tools are primarily useful for visualization but can also be used as the basis for automated finite element meshes for groundwater numerical models, as described later.

While several CAD packages are available providing tools for solid modeling, most are geared principally for describing objects of a mechanical nature. GMS provides modeling tools that have been designed specifically for defining surfaces and solids of geologic origin. The program will accept field data in the form of borehole logs and allow the user to process the data to directly define a solid model. The approach taken is to first define surfaces separating the various stratigraphic units, followed by a series of extrusions and set operations to define a solid model of the site.

Borehole logs are typically the most basic form of information available describing subsurface characteristics. Borehole logs represent the changing material properties along strategically placed vertical traces scattered throughout a site. They can be input by means of ASCII data files. Three-dimensional oblique views of the boreholes can then be displayed (Figure 1).

GMS uses triangulated irregular networks (TINs) to model terrain and the interfaces between stratigraphic units. Jones [2] presents several methods for creating and editing TINs to model surfaces of geologic origin, most of which have been incorporated into GMS. To generate a TIN, contact locations representing material interfaces are graphically or automatically selected from existing boreholes. The selected contacts serve as the vertices of a Delauney TIN. The TIN can also be extrapolated beyond the convex hull of the boreholes to a predefined site boundary. New vertices can be added and existing vertices graphically selected and dragged, allowing the user to sculpt the stratigraphy to conform to plausible conditions at the site.

The procedure of defining a TIN is repeated multiple times until a series of surfaces has been created, stacked on top of one another, representing each of the interfaces between the different material types of the hydrogeologic domain (Figure 2). Operations to fill between two adjacent TINs and extrude a TIN vertically can be used to create solids. Set operations can be used to difference or union existing solids to create other solids, allowing for complex modeling of faults or pinch-out zones within the model.

Once the solid model of the site has been generated, three dimensional visualization techniques may be employed to cut cross-sections or fence diagrams, rotate, zoom, pan, and render the solid (Figure 3).

Pre-processing

In an effort to provide an integrated environment within GMS, specific groundwater flow and transport numerical models have been incorporated into the system. Currently, several numerical models designed to simulate three dimensional flow and transport are available. While some finite element codes are in use, the finite difference approach is the most widely accepted method within the groundwater modeling community. As a result, a complete graphical user interface has been provided in GMS for the USGS, MODFLOW [3], a well-established, three-dimensional finite difference groundwater flow model.

Because of some of the inherent geometric deficiencies with finite difference modeling, a graphical interface to the three dimensional finite element model, FEMWATER has also been developed within GMS. FEMWATER is a saturated/unsaturated, density driven, coupled flow and transport model. FEMWATER was originally developed by G.T Yeh [4] and is currently being modified and maintained by the US Army Engineer Waterways Experiment Station. As time and funding permit, GMS will support additional flow and transport codes.

Semi-automated and fully-automated procedures are provided for generating two and three dimensional meshes and grids. A series of dialogs and graphical tools can then be used to define specific boundary conditions relevant to the analysis code to be used. Geometry and boundary condition files can then be saved and the solution computed using the analysis code.

Finite Elements

GMS provides several techniques for generating both two and three dimensional finite element meshes. These techniques include both manual, point-and-click methods as well as automated methods. Mesh generation for finite element analysis has been the subject of much research. It has, however, been primarily motivated by the auto and aerospace industries. While many of these techniques can be borrowed and effectively used in groundwater modeling, there are unique characteristics which can be taken advantage of when generating a mesh for a groundwater model.

two dimensional mesh generation

GMS provides a palette of tools for generating linear and quadratic two dimensional elements. Elements can be created individually by graphically selecting vertices or by filling predefined areas with a specified type of element using a Coon's patch technique.

Adaptive triangulated meshes can also be quickly defined from a polygonal boundary. This technique involves triangulating a set of adaptively defined nodal points followed by a relaxation process. A polygon can be input interactively or read from an ASCII file previously digitized from a site map (Figure 4). Feature points and lines within the boundary polygon representing well locations and rivers can also be included forcing the mesh to conform to interior characteristics.

After defining the boundary geometry, nodes are adaptively inserted into the mesh using a quadtree approach. Node density is determined from a weighting function dependant upon the edge lengths of closest feature lines and boundary edges. The root of the quadtree is defined as the bounding box of the polygon with the quadtree recursively subdivided and nodes inserted until the mesh density function is satisfied in all regions of the model (Figure 5).

The triangular elements are defined using a Delauney triangulation. Breaklines are then inserted which force the mesh to conform to the boundary polygon and feature lines. Triangles defined outside the boundary polygon are removed (Figure 6).

The interior nodes, not belonging to a feature line or point are then relaxed. The relaxation method used involves iteratively repositioning nodes to the centroid of their adjacent triangles (Figure 7). If quadrilateral elements are desired, adjacent triangles whose interior angles are less than a specified threshold value are merged and the relaxation process repeated.

Two dimensional groundwater models often require the surface elevation to be defined at each node of the mesh. A TIN or set of sparse data points scattered throughout the site can be used as the basis to interpolate the elevations of any nodes added to the mesh during the adaptive meshing procedure. A C2 continuous interpolation scheme is used to interpolate elevations at each node during the insertion process and during the relaxation procedure.

Several tools for graphically editing an existing mesh are also provided. Individual or groups of nodes and elements can be selected and deleted. Nodes can be dragged, while the mesh is dynamically updated. Node elevations can be edited by dragging the node in an oblique view or adjusting a vertical scroll bar while contours are displayed and updated dynamically. Quadrilateral elements can also be selected to be split and adjacent triangular elements selected to be merged. Breaklines can also be forced through the existing mesh by graphically defining a poly-line.

Future versions of GMS will also include Blacker's [5] paving algorithm. Figure 8 shows an all-quadrilateral mesh of the same domain used for the adaptive meshing technique. The same approach to defining exterior boundaries and feature points is used. Feature line capability has not yet been developed for use with paving inside GMS. Mesh density can also be controlled by the boundary polygon edge lengths.

Presently, GMS does not support a two-dimensional groundwater numerical model. The primary purpose for the two-dimensional mesh generation tools are as a preliminary step to defining a three-dimensional model.

three dimensional mesh generation

While three-dimensional mesh generation for groundwater models is still an area of on-going research, GMS currently provides a limited number of automated methods. These methods are motivated by the type of data typically available from a hydrogeologic study and by the requirements of the numerical model to be used.

One technique, developed by Davis [6], allows the user to generate the mesh directly from a set of boreholes obtained from a site. The first step in this procedure is to develop a two dimensional mesh of the site in plan view. The purpose of the two-dimensional mesh is to provide a template from which to extrude a three dimensional mesh. Using the tools available for site characterization, a series of TINs defining the interfaces between the various hydrostratigraphic units can be created. This procedure is the same as that used for developing a solid model of the site.

Once the interfaces between layers have been defined, two adjacent TINs can be selected and a mesh created between them. The two dimensional mesh is extruded vertically to create any number of three dimensional element layers between an upper and lower TIN. The result is a regular three-dimensional mesh defined from either hexahedral or wedge shaped elements (Figure 9). The elevation of the nodes at the top and bottom of the mesh are interpolated from the bounding TINs using a C2 continuous interpolation method. The actual vertices of the TINs are not necessarily used in the final mesh, but rather used to define a continuous surface for the finite element mesh boundary. Since most geological surfaces are statistically determined from sparse data collected at a site, interpolation methods are quite appropriate for defining the boundaries of the mesh. Elevations of intermediate elements within the mesh are determined by linearly interpolating along a vertical line between the top and bottom TINs for each node of the two-dimensional extrusion mesh.

A simplified form of this procedure has also been defined in which the user selects a material region to be meshed. After specifying the number of vertical layers of elements to be included within the region, the program automatically selects all boreholes containing this material and extrudes the pre-defined two-dimensional mesh through the region. The process of defining top and bottom bounding TINs is automatically performed resulting in a well defined three-dimensional mesh. This procedure can be repeated multiple times through each different material type defined by the boreholes. This technique is, however limited to well-defined stratigraphic layers extending throughout the horizontal extent of the domain.

Once the mesh has been defined, three-dimensional editing of nodes and elements may be performed. Dynamic checking for valid element shapes is also incorporated, as nodes are dragged in 3D.

Additional automated three dimensional meshing techniques are currently being researched for future inclusion into GMS. These methods include both tetrahedral and hexahedral approaches. Work is currently underway to implement the three dimensional counterpart to the two-dimensional adaptive meshing algorithm as well as the three dimensional counterpart to paving, known as plastering [7]. The ultimate goal of some of these more advanced techniques is to be able to define the mesh directly from the solid model of the hydrostratigraphy with minimal interaction from the user.

Finite Difference

Finite difference codes require that a regular Cartesian grid be defined over the domain of the region to be modeled. Since finite difference grids are limited to rectangular or regular hexahedral cells, the difficulties in generating the grid are greatly reduced from that of finite element meshes. GMS provides several graphical techniques for building and editing a finite difference grid in both two and three dimensions. Active and inactive cells may be automatically defined using a polygon either input interactively or read from an ASCII file representing the boundary of an existing hydrogeologic domain. Rows, columns and layers may be easily inserted into the model. Refining and editing operations are accomplished by graphically pointing and dragging with the mouse. Automated methods for adaptively biasing the grid to refine near specific points have also been included (Figure 10).

Interpolation of scattered data

One of the principal problems involved with working with field data, is the difficulty of transferring the data into a useful form that can be viewed or assigned as input to a model. Since only a limited number of scattered data points are typically available, the data must be interpolated to an underlying mesh or grid to be of use. GMS provides a variety of interpolation and geostatistical schemes in both two and three dimensions for characterizing hydrogeologic data. For example, borehole samples of contaminant concentrations scattered throughout a site may be interpolated to a dense grid and visualized in the form of contours or iso-surfaces (Figure 11). Another application of geostatistics is for defining hydrogeologic properties for input to a numerical model. For instance, tests measuring hydraulic conductivity may be performed at scattered locations throughout the site. Results may then be interpolated to each cell of an underlying finite difference grid.

Interpolation has been the topic of much research producing many different schemes possessing a variety of different characteristics. No technique is guaranteed to give the best results under all circumstances. Because of this, GMS provides a number of widely accepted methods, leaving it to the user to decide which is the most appropriate for their application. The primary interpolation schemes supported include Inverse Distance Weighted [8], Clough Tocher [9], Natural Neighbor [10,11] and Kriging [12].

Post-Processing

In keeping with the objective of providing a fully integrated environment for groundwater modeling, significant effort has been made to provide state of the art tools for scientific visualization. Once the analysis engine has been run, GMS will directly read the analysis results and display them in a variety of formats. Methods for visualizing two and three dimensional data include contours and fringes, vector arrows, isosurfaces, cutting planes as well as several animation techniques.

Often, the results from a groundwater simulation are in the form of transient data. In order to visualize this type of data effectively, animation techniques must be used. GMS provides the capability to generate a sequence of images representing different time steps in the simulation. For example, each frame may represent the changing head of a water-table over time or the changing position and geometry of a contaminant plume.

GMS also utilizes animation for visualizing steady state data. Cutting planes can be swept through a volume while fringes of contaminant concentration or pressure head are mapped onto its surface. The value of an iso-surface may also be changed with each frame of an animation sequence to represent how the contaminant concentration may vary spatially within the subsurface.

Conclusion

While the state-of-the art in analysis of groundwater flow and transport has advanced rapidly, the tools needed to effectively take advantage of these codes have not advanced at the same rate. This has occurred during a time when many other scientific and engineering disciplines have taken full advantage of new computer graphics technologies. In spite of this, groundwater modeling has continued to use older, less effective methods for processing and visualizing data. GMS has been designed specifically for the discipline of groundwater modeling incorporating many of the latest state-of-the art techniques in user interfaces and computer graphics.

The synthesis of any number of manual and automated methods is often required for a Hydrogeologist to complete a successful groundwater modeling project. An extensive knowledge of specific file formats, operating systems, computer languages and a variety of software tools is often a necessity. While most Hydrogeologists have an in-depth understanding of the principles behind groundwater flow and transport, many will avoid highly effective numerical analysis tools because of the substantial overhead of learning the computer technology required. GMS is an attempt to integrate and simplify the process of groundwater flow and transport modeling by bringing together many new and existing technologies and applying them in a consistent and user friendly environment.

References

Figures