Thomas Heinzer, Michael Sebhat, William Greer, Dave Hansen, and Charles B. Johnson
Abstract
Geographic Information Systems are the future of pre and post processors for many
types of modeling. This paper examines the use of ArcInfo menu driven GRID
applications to produce a robust interface to the MODFLOW finite difference model.
The interface enables the setup of the finite difference mesh, the generation and raster
editing of model layers, the use of map algebra GUI equation mbuilders for rapid
population of data layers, and provides automated generation of MODFLOW ready
files. The model runs directly from the interface, and provides for viewing, contouring,
and map production of model run results. This paper further details the AML/Unix
programming approaches used to create the model links, the application of GRID 7.0
functionality, and the advantages and difficulties of this type of interface.
Introduction
The integration of Geographic Information Systems and modeling has long been a
topic of discussion and research. Indeed, the two disciplines have been somewhat
segregated from each other, as they both tend to be quite complicated by nature.
However, it is clear that a marriage of the two disciplines would have wide
reaching ramifications. This paper is one such attempt at unification.
A model possesses a mathematical component, along with certain assumptions, that
define its function. Another component exists, however, which we term the
spatial component: the data sets that the model accesses. We have found through
experience that this second component is of equal importance to the first; if
the data structures and/or the georeferencing of the model itself are in error,
the overall results are inevitably erroneous, no matter how elegant the
mathematics.
Our interest is to examine methods by which we can talk to models through GIS
systems. In other words, data integration and manipulations are performed using
a GIS based Graphical User Interface (GUI) which additionally automates any tasks
that are required for the model to access the data. If done properly, this type
of approach can eliminate many sources of error and reduce run preparations
significantly.
Originally, we incorporated the finite difference equations directly into
Grid map algebra expressions. The model mathematics were performed directly
on the Grid layers, using Gauss-Seidel interation techniques within the map algebra.
This worked, and although the concept was esthetically pleasing, disk I/O slowed
run times beyond what is generally acceptable.
We decided therefore to interface with a proven finite difference engine, which
was the USGS MODFLOW package (1). Considerable work has been done in this area,
most notably that of Leonard Orzol(2), who was the first to link ArcInfo vector
coverages to MODFLOW source code. Our interest was to interface the raster
component of ArcInfo, Grid, with MODFLOW. The goal was to create a GIS
interface that would allow one to:
1) Generate the finite difference mesh. This would enable one to specify
cell size, and allow mesh rotation and moving in reference to
other geo-relational data.
2) Automatically generate model layers and provide for overall data
management. Enable graphic layer addition, deletion and editing
using ArcInfo Grid 7.0 GridEdit functionality. Also develop an
'equation builder' to facilitate rapid layer setup and manipulation.
3) Visualize Model Output to enable one to view and analyze model
output, usually heads and drawdowns using shade ramping and surface
contouring.
There are two main concerns raised when discussing the use of ArcInfo Grid to
interface with a finite difference ground-water flow model. The first is that
if there is significant horizontal anisotropy, it is necessary to rotate the
mesh to align it with the principal directions of transmissivity.
The problem arises in that ArcInfo grids are always upright and cannot be rotated.
We felt that the advantages of using Grid far outweighed the ramifications of this
problem, and sought to overcome it. If the mesh had to be rotated, then the most
logical solution was to automatically rotate the geodata sets behind the mesh,
enabling viewing and better decision making. So, if the interface detects mesh rotation,
it automatically rotates the geodata sets so they will be properly oriented in
relation to the mesh.
The second concern is the fact that ArcInfo Grid cells cannot be irregular in
shape in the horizontal plane. In areas where hydraulic gradients and/or aquifer
properties are rapidly changing spatially, it has been the practice to refine the mesh in
these areas to 'more accurately' capture these effects. Additionally, sometimes outlying
cells are made larger to reduce processing times. However, the authors contend that being
restricted to a uniform cell size (in the horizontal plane) is not a serious limitation,
especially in view of rapidly increasing computer speeds.
Conceptual Approach
The goal was to develop a GIS Graphical User Interface to MODFLOW that would be
completely 'point and click', providing the ability to generate a model from
scratch. It soon became clear that an additional scenario had to be addressed:
we also needed the ability to interface with pre-existing models. A current
application of this technique will be discussed in the Conclusion section later
in this paper.
The basic functions of the interface are: Store the ASCII array data in the form
of ArcInfo grids; provide menu driven editing and viewing capabilities; use
GRIDASCII to drop the ASCII array components of the grids; automatically
generate MODFLOW ready files; run MODFLOW; view and analyze the model results,
generate contours, and produce maps or displays.
The interface is broken into three main components:
1. The Mesh Generator
2. The Setup of the Model DataBase Structures
3. The DataBase Populator, Editor, and Viewer.
Components 1 and 2 are utilized once and then locked to prevent inadvertant over-
writing. If one wants to start a new model, a small set of programs or 'seed'
directory is copied elsewhere and a new model is spawned. Component 3 can be invoked
as many times as desired, as this is where changes are made and the executable program
is launched.
The directory structure is straightforward. A directory exists for each MODFLOW
module, and underneath each one are the associated grids and MODFLOW data files.
For example, in the [pumps] directory one would find a number of grids, each
pertaining to well layers that exist in the model, with the possibility of
additional grids representing stress period changes.
Whatever Grid layers exist in a given directory, the interface uses them to
create the MODFLOW ready file. Therefore, modifying a model amounts to adding,
deleting, or editing Grid layers. What you see, MODFLOW gets. The grids have
the naming convention layer(layer)_(period). Hence, a layer denoting layer 6 at
stress period 10 would be named layer6_10.
Below is an example of a possible [pumps] directory:
Here pumping starts at time step 1, stress period 1 in layers 5 and 6. At stress
period 10, however, the pumping rates in layer 6 are changed. The two programs,
reset.aml and pumpfilegen, access the grid layers and build the MODFLOW ready
"pumps.data" file.
After the model run has executed, the interface scans MODFLOW's output file and
breaks the results (drawdowns and heads at specified times) into corresponding
grids, which can be viewed and analyzed. Popup menus allow the setting of
contouring parameters (base contour and interval). A command entry area also
facilitates advanced ArcInfo manipulations if needed.
To start the interface, the aml program 'GMI.aml' (which stands for Grid-MODFLOW
Interface) is executed in Arc. This calls the interface menu below which contains
three buttons pertaining to three categories, which are discussed below in some
detail.
The Mesh Generator
This module facilitates generation of the finite difference mesh. A simple menu
(above) and an ArcInfo canvas appear. At this point any ArcInfo data coverages
we want to be able to view should have already been placed in the [coverages]
directory. The first button is "view GIS data". Here we can call ArcInfo geo-
data sets in different colors as an aid in generating and placing our mesh.
The next button, "Set Origin Point", prompts us to place the lower-left coor-
dinate of the mesh. Below it, 3 slider bars allow us to (1) set the number of
rows, (2) set the number of columns, and (3) set the cellsize of our model. At
this point, we can view our settings, and change them if we so choose. The
"Generate" button then generates the mesh, which appears on the screen.
Usually we want to move and rotate it a number of times, and the next button
"Move Grid" and subsequent button and slider "Rotate Grid" provide for this.
When we are satisfied with the mesh placement, the "Save" button saves the
current mesh for use with the model.
It should be noted that when the mesh is first generated, the specified settings
are passed to the generate command, which creates the mesh as a coverage. The
mesh coverage is then displayed in Arcedit, with the other geodata sets shown as
Arcplot background calls. The rotations and moves are accomplished in Arcedit,
and then saved out as a coverage. The specified move and rotation parameters
are used to generate a corresponding Grid layer. If the mesh is rotated, all of
the geodata sets in the directory [coverages] are rotated accordingly and placed
in the directory [geodata]. This allows viewing of coverage data at any future
time. If the mesh is not rotated, the coverage mesh and the Grid are coincident
and the coverages are merely copied to the directory [geodata] for future visual
use.
In addition to rapidly generating the mesh in conjuction with real spatial data,
we have ensured the availablity of our geospatial data at any time during the
modeling process, rotation or no rotation. Furthermore, the pan/zoom viewing
capabilities of the ArcInfo canvas may also be utilized.
The Model Setup
This module sets up the database and directory structure that the interface will
access. A menu pops up (above) with slider bars which set the number of model
layers and number of stress periods. The interface then creates all the MODFLOW
module directories (ibound, storage, tops, etc.) with corresponding empty grids (we
call them protogrids; item value = 0) needed for a base model run- ie. those
layers that are mandatory for MODFLOW to execute.
The optional package module directories are set up (drains, wells, recharge...)
and can be implemented, if desired, by a form menu which sets up the bas.data file.
Two programs are also copied to each module directory: reset.aml and modgenfile.
Reset.aml drops all of the ASCII arrays down from whatever Grid layers exist in
that directory. For example, if the [pumps] directory contains layer5_1 and
layer6_10, two ASCII arrays would be dropped, representing pumping at layer 5
period 1 and layer 6 period 10. The modgenfile program is a shell script that
assimilates these ASCII arrays into the file that MODFLOW ultimately reads.
This is done using UNIX filters, pipes, and shell programming, the details of
which will be presented later in this discussion. These programs vary for
different modules, as the module files take on differing formats. In general,
when the Grid layers are added, deleted, or changed, these programs reset the
MODFLOW files.
The Database Populator, Editor, and Viewer.
This is the main module of the interface. Below is an example of a typical
session.
At the top of the interface menu (above), all of MODFLOW's modules appear. Upon the
selecting a module, a pop-up menu presents all of the current layers in that
module. After selecting a layer, it is graphically displayed in the ArcInfo
canvas to the left. The layer is identified in the small box to the right,
along with the min, max, mean, and standard deviation of the cell values in the
layer. The functions that follow will be discussed separately.
1. The "ColorMap" button threads the above menu that allows control of the layer
shading or coloring. Shade color ramps can be used to ramp continuous data
sets such as starting heads. Here, two colors are selected which are linked
to the minimum and maximum values in the layer, and a 'linear stretch'
produces gradient effects. The 'identity' remap is useful when visualizing
categorical data, which may occur, for example, in pumping layers. Custom
remap tables are also available to color nonstandard data types. The goal
is to provide the user enough 'point and click' tools to reasonably visualize
any data layer that one might encounter.
2. The "GeoData" button facilitates viewing of vector ArcInfo geodata sets in
conjunction with the model layers, enabling rapid and accurate decision
making. The above menu appears, and the user may select from various line
weights and colors, and then select a data set for viewing. Attribute data may
be identified if relational database information is needed. This functionality
has proved to be extremely useful, not only for populating a model, but also in
applications such as the USGS program ZONEBUDGET (5), which performs water-budget
analysis on specified model subregions by utilizing MODFLOW.
3. The "LayerEdit" utility uses the ArcInfo 7.0 Grid raster editing tools to
update model layers. The above menu appears. To populate or change layer
cell values, the new value is entered, and the user can change one cell at a
time, or select a box or polygon and change a group of cells all at once. An
'oops' button permits backing up to erase mistakes. At any point the layer can be
saved, and the changes written to disk in the form of a grid. The interface
detects if changes have occurred in the grid layers and updates the MODFLOW files
accordingly; if a user changes a data layer graphically, the MODFLOW program
will see it.
4. The "LayerBuilder" uses the map algebra functionality to manipulate model layer
data (see above). The first equation builder allows the user to set an existing
or new layer equal to a single numeric value, which is useful for setting a
regional storage coefficient, for example. The second equation builder allows the
setting of a layer as a function of itself or another layer. This is useful, for
example, when setting up the layer dimensions (tops and bottoms), or varying
model stresses, for example doubling the recharge. The third equation
builder invokes the 'con' functionality of the map algebra. This enables
model layer cells to be set to algebraic combinations of other layers if they
meet a certain condition. For example, we may want to double the pumping
rates only in a certain area or 'zone' of the model. If we had this 'zone'
layer, we could build the expression:
layer5_1 = con( zonelayer = 6, layer5_1 * 2, layer5_1 )
which means: on the condition that a model cell lies in zone six, double its
pumping rate, otherwise, keep it at its current value. This procedure can be
extremely useful when performing sensitivity analyses in conjunction with zone
budgeting operations. We have found this functionality to drastically speed
certain aspects of layer manipulations.
For ibound array setups, a 'make all layers like layer1' button has been
added, and layers can be subsequently edited if need be.
Required by MODFLOW is the bas.data (Basic Package) file, which defines boundary
conditions, starting head conditions, stress period durations and other model
variables. A menu provides for the generation of this file. Also required by
MODFLOW is the bcf.data (Block Centered Flow Package) file, the foc.data (Output
Control) file, and the sip/sor.data (solution algorithm) file. These files
can be generated via a menu interface, which will not be addressed in detail
here.
Viewing Output
MODFLOW's foc.data file provides model output control, enabling the user to set
various points during a run that head and/or drawdown values are saved to disk.
A menu interface is provided to accomplish this. When a run is complete, the
specified output can subsequently be viewed. This is accomplished by a set of
UNIX filters acting on MODFLOW's output file, which breaks out and loads the data
into their corresponding grids. At that point, the data can be viewed by selecting
an output layer and using the ColorMap utility (see above).
Generating Surface Contours.
It is usually helpful to utilize contouring algorithms to depict model output
head or drawdown distributions. The interface provides a menu (see above) that
allows the user to set a base contour and contour interval, and subsequently
interpolate the given surface, using the standard ArcInfo SURFACECONTOUR command.
We have found this to be sufficient in many applications, especially when dealing
with model output, where many times it becomes futile to use extravagant
interpolators on model generated data sets. However, other interpolators are
available such as Kriging and Inverse Distance Weighting since the interface
executes out of the ArcInfo environment.
The Map Recorder
At certain points in the modeling procedure it is desirable to obtain maps or
hardcopys of the interface canvas. At any time the user may 'make a map' by
invoking the 'start record' button (see above). Any layers, contours, and vector
geodata sets available to the interface may then be called. After recording is
stopped, the map is saved as an ArcInfo map composition, which can be further
devoloped, or sent directly to an electrostatic plotter or postscript output device.
This feature provides for quick generation of maps, from working plots up to
presentation quality displays.
Technical Details
Originally, the supposition was that it would be nice if MODFLOW could read the
arrays directly from the ArcInfo grids. While this seemed appealing, there are
certain problems with such an approach. The first is that special modification
of the MODFLOW code is necessary, which we felt was unneeded, since in the Grid
system, ASCII arrays can be rapidly dropped from grids. Conversely, the ASCII
arrays can readily be loaded into grids. Using this method, the interface
becomes relatively independent of software enhancements on either the MODFLOW or
ArcInfo platforms.
Secondly, we do not expect modelers be required to blindly trust this interface.
This interface simply generates the files that the modeler would generate anyway
except much more rapidly and reliably. We feel this is a sound approach, since
the files that the interface generates can be inspected manually for integrity.
When a Grid layer is edited and saved, the command is used to drop
its ASCII array. Some of the Basic Submodule files are ASCII arrays of set
format (e.g. ibound); these are piped through awk filters containing programs
with fprint statements to format the file appropriately, and then all arrays are
concatenated to form the MODFLOW ready file. Other modules are not that simple.
While it would be too tedious to examine all modules, it is appropriate to show
how one, the pumps, was programmed.
This program (aml) is named reset.aml, which detects model layers and drops the
ASCII grids, and processes them into sub-files (modarray*_*). A script called
pumpfilegen is then invoked, which sets up the final form of the MODFLOW ready
pump.data file.
/*******************************************************
/* This program is reset.aml:
/* Note: The embedded awk program, pumps.awk, is displayed below.
/* Note: The embedded UNIX script, pumpfilegen, is displayed below.
/*******************************************************
&sys \rm array*
&sys \rm modarray*
&do per = 1 &to %.nper%
&do lay = 1 &to %.nlay%
&if [exists layer%lay%_%per% -Grid] &then &do
&type Resetting Layer %lay% Period %per% ...
arc GridASCII layer%lay%_%per% array%lay%_%per%
&sys awk -v R=`cat ../amls/numrows` -v C=`cat ../amls/numcols` ~
-f pumps.awk `tail +7 array%lay%_%per%` | grep -v ' 0' |~
awk '{printf "%10s%10s%10s%10s\n",'%lay%',$1,$2,$3}' > modarray%lay%_%per%
&end
&end
&end
&type Array Dumps Complete...
&sys ./pumpfilegen /* Note This script is displayed below.
&type MODFLOW pump.data file complete
/*******************************************************
/* The awk program pumps.awk is defined by:
/*******************************************************
BEGIN {
for (j = 1; j < R+1; j++)
for (i = 1; i < C+1; i++)
printf j" " i" %s\n", ARGV[(C * j) - C + i]
exit }
/*******************************************************
/* The pumpfilegen script is defined by:
/*******************************************************
# This script sets up MODFLOW-ready pump files
if [ -f pump.data ]; then rm pump.data;fi
if [ -f pumparray ]; then rm pumparray;fi
cat /dev/null > pumparray
# For all timesteps, check for modarray(lay)_(per). If it exists,
# write to file pumparray. If not, send -10 to first column
# to signify 'use last pump layout'.
# The maximum wc for modarray(lay)_(per) is calculated and
# put at top of file. Final output file is pump.data.
MAX=0
for VAR in `cat ../amls/periodallfile`
do
echo $VAR
if [ -f modarray*_$VAR ] ; then
cat /dev/null > temp
for MODFILE in `ls modarray*_$VAR`
do
cat $MODFILE >> temp
done
wc temp | awk '{printf "%10s%55s\n",$1,"Stress'$VAR'"}' >> pumparray
cat temp >> pumparray
COUNT=`cat temp | wc -l`
if [ $COUNT -gt $MAX ]; then MAX=$COUNT;fi
rm temp
else
printf "%10s%55s\n" -10 Stress$VAR >> pumparray
fi
done
echo $MAX `cat ../amls/var/iwelcb` | awk '{printf "%10s%10s\n",$1,$2}' > hd
cat hd pumparray > pump.data
rm pumparray hd
Similar programs exist for the other modules. One problem of interest is that of
format. The GRIDASCII command drops arrays consisting of numbers separated by
one space. MODFLOW, on the other hand, expects to see a specified format for
loading. In practice, some of the modules could contain numerics ranging from
large integers down to small, negative exponents. We had to develop a method
that would permit these various possibilities. ArcInfo Grid stores numbers in
a four byte input, 12 byte output field format; the maximum bytes a number
coming out of an ArcInfo Grid could assume is twelve (including exponential).
We found that for the general case, if we passed the ASCII arrays through an
fprint %12s format specifier, and the MODFLOW fortran format read statements
were set at F12.0, data could be passed, as far as we could tell, regardless of
values. This eliminates the need for a free format version of MODFLOW, however
when it is released, we will probably eliminate this additional filter, as it
will not be needed.
Discussion and Conclusion
The goal of this research was to develop a menu driven GUI to the ArcInfo Grid
model that interfaces with MODFLOW software. Many hours have been dedicated to
this endeavor, from hydrologists to modelers to GIS programmers to database,
system, and program managers. It truly is an interdisciplinary project. It is
our hope that this attempt to interface GIS with modeling will encourage others
to develop new and better ways to integrate new technologies into their modeling
projects. Some sites will be beta testing this software, and when it is found
stable enough, it will be available at no cost to anybody who wants to use it
(ArcInfo 7.0 is required). Comments or suggestions regarding this project are
certainly welcome.
This interface has been utilized within the Mid Pacific Region in conjunction
with a West San Joaquin Valley MODFLOW simulation (Belitz et al, 1993).
The model is a rotated, 37 by 20 mile mesh consisting of 1 mile square cells.
It is calibrated over a 20 year period, and Reclamation has been using it to
evaluate different pumping and recharge scenarios in the area. The interface
implemented numerous scenarios on both the base case and a 5 year extrapolation
into the future. This application enabled us to debug and enhance the interface
program's capabilities. One result was clear: Once set up, scenario evaluations
could be rapidly performed.
It should also be noted that we have not discussed all of the ArcInfo
functionality that can be used in this application. We felt it illogical to
incorporate all of ArcInfo into the menu systems; we placed the functionality
into the interface that is required to get the job done. We wanted to avoid
forcing the user into the intricacies of ArcInfo, and present a more intuitive
menu system that users could readily start to work with. However, we did not
want to preclude advanced GIS users from performing other functions, so a
command line entry point is available at all times. In addition, for example,
simple AML routines can be developed to present head variations in a given
layer for a model run (see below).
When modeling, we have found over the years that using GIS to manage, view and
analyze model data is extremely useful. As new data sources become available,
implementation of these data is certainly necessary, however managing the data
becomes overwhelming without GIS. The integration of GIS and modeling has been
and will continue to be a major topic of research. We hope this paper has shown
how useful a GIS interface can be in assisting the ground-water flow modeler.
References:
(1) Micheal G. McDonald and Arlen W. Harbaugh, 1988, A Modular Three-Dimensional
Finite-Difference Ground-Water Flow Model: U.S.G.S. Techniques of Water-
Resources Investigations, Book 6.
(2) Leonard L. Orzol and Timothy S. McGrath, Summary or Modifications of the
U.S. Geological Survey Modular, Finite-Difference, Ground-Water Flow Model
to Read and Write Geographic Information System Files.
(3) David R. Maidment, 1994, Hydrologic Modeling Using ArcInfo, Seminar Presented
at the 14th Annual Esri User Conference.
(4) Esri, 1991, Cell Based Modeling With Grid, ArcInfo User Manual.
(5) Arlen W. Harbaugh, 1990, A Computer Program for Calculating Subregional Water
Budgets Using Results from the U.S. Geological Survey Modular Three-Dimensional
Finite-Difference Ground-Water Flow Model: Open-File Report 90-392
(6) Kenneth Belitz, Steven Phillips, and J.M. Gronberg, 1993, Numerical Simulation
or Ground-Water Flow in the Central Part of the Western San Joaquin Valley,
California: United States Geological Survey Water-Supply Paper 2396
Author information:
Thomas Heinzer
Consulting GIS Analyst
Chemical Engineer
MPGIS Service Center
U.S. Bureau of Reclamation
2800 Cottage Way
Sacramento CA., 95825
Telephone: (916)979-2442
E-mail: theinzer@mpgis2.mp.usbr.gov
Michael Sebhat
Consulting GIS Analyst
Electrical Engineer/Computer Science
Project Manager, Systems Administrator
MPGIS Service Center
U.S. Bureau of Reclamation
2800 Cottage Way
Sacramento CA., 95825
Telephone: (916)979-2442
E-mail: msebhat@mpgis1.mp.usbr.gov
William Greer
Hydrologist
2800 Cottage Way
Sacramento CA., 95825
Telephone: (916)979-2322
E-mail: bgreer@mp710b.mp.usbr.gov
David Hansen
GIS Specialist
Soil Scientist
2800 Cottage Way
Sacramento CA., 95825
Telephone: (916)979-2419
E-mail: dhansen@mpgis7.mp.usbr.gov
Charles B. Johnson Jr.
Regional GIS Program Manager, MPGIS
Regional Soil Scientist
U.S. Bureau of Reclamation
2800 Cottage Way
Sacramento CA., 95825
Telephone: (916)979-2419
E-mail: cjohnson@mpgis6.mp.usbr.gov