Bruce Rindahl

Groundwater Modelling and Analysis

Using the USGS MODFLOW program
and ArcView

Abstract

The City of Aurora's Utilities Department is utilizing ArcView 2.1 to analyze groundwater modelling and aquifer response the Cherry Creek alluvium. ERSI's ArcView 2.1 was chosen for displaying drawdowns, stream flow and aquifer elevations simulated from the USGS MODFLOW program because of the clear display of the complex model simulation results. Utilizing ArcView's ability to join GIS coverages and text tables, an easy to use interface for displaying MODFLOW results was developed. This application is the basis for a water rights application and a basin-wide management plan for the upper Cherry Creek basin. This paper will lay out the steps in the development and use of this tool including Avenue script customization.


Introduction

The City of Aurora is a community of approximately 247,000 located in the southeast Denver metro area. Aurora is the third largest city in the State of Colorado. The service area of the Utilities Department is about 135 square miles. Figure 1 shows the general location of the model grid.

Since 1956, some of the water supply of Aurora has come from a seven wells in the Cherry Creek aluvium located just southwest of the city. As a result of State of Colorado regulations enacted since the wells were drilled, an analysis of the potential impacts to the Cherry Creek and South Platte stream system was required for continued use of the wells. In 1991 the City of Aurora hired Bishop-Brogden Associates, Inc. to analize the well field and the effected reach of Cherry Creek.

Bishop-Brogden Associates (BBA) selected the MODFLOW computer program with an additional streamflow package developed by the USGS. The purpose of this paper is to describe the application of ArcView to display and analyze the results of this modelling effort.


Location Map

Figure 1 - Location of the City of Aurora


Model Description

The model developed by BBA consisted of a 72 column by 244 row grid of 200 foot squares aligned along the axis of Cherry Creek upstream of Cherry Creek Reservoir. Due to the limited scope of this paper, a detailed description of the parameters of the model is not included but can be optained from the City of Aurora. Input parameters to the model for each cell are tabulated in input files read by the MODFLOW program. Additional parameters such as streamflows, well pumping, and time periods are also input to the program. Outputs are typically groundwater elevation, drawdowns from the starting conditions, and streamflows in either printed form or binary files. Printed output is quite voluminous and requires a through knowledge of both the particular model and groundwater hydraulics.

GIS Coverages

In order to evaluate the model results using ArcView, the 72 column by 244 row grid was created as a polygon coverage by ArcInfo using the GENERATE and FISHNET commands. Additional attributes of Row, Column, and Loc-Tag were added to correspond to the model. Loc-tag is simply a character string consisting of row number and the column number seperated by a "-". The Loc-Tag is the vital attribute through which model results will be joined to the coverage (or Theme in ArcView). Additional fixed attributes such as bottom elevation and transmisivity are added to the coverage from the model input files. The basic grid can then be displayed using ArcView as a standard theme. Finally, the coverage was translated to the correct coordinates of the actual location of the model in order to utilize the existing GIS coverages of the City of Aurora. Figures 2 and 3 show the model grid displaying the bottom elevation and transmisivity respectively. Figure 4 shows a portion of the attribute table for the coverage.

A copy of the model grid coverage was created using only those cells where Cherry Creek streamflow is modeled. This coverage is used to analyze the surface flow in the model. In addition to the model grid, coverages were developed for other features in the model such as well locations and monitoring holes to display these features on the same reference as the model.

Bottom Elevation

Figure 2 - Bottom Elevation of Cherry Creek Aquifer


Transmisivity

Figure 3 - Transmisivity of Cherry Creek Aquifer


Attribute Table

Figure 4 - Attribute Table of Cherry Creek Model

Displaying Model Results

As stated before, the MODFLOW program can output results in both printed output and binary files. To utilize the GIS coverages described above, the binary files were translated into ASCI text tables. Several simple programs were developed to read the binary output files and output the text tables. The program first skips the header of the binary file then reads and stores the value for every cell in the first time step (or stress period as it is called in MODFLOW). This is then repeated for every time step in the file. Finally, the output file is created first by writing the header, then one line for each cell consisting of the Loc-Tag and the values for the drawdown (or head elevation or streamflow) for each time step of the model run. As an added bonus, this procedure can be used to compare two model runs. The first run is read in and saved in memory, then the second is read and the output file is the cell by cell difference between the two files. In this manner the impacts of well pumping can be analyzed by running two models - one with and one without pumping - and comparing the results. Copies of the program codes are included in the appendix.

The final step to displaying the model is deceptively simple - just join the attribute table of the model grid with the output text file using the Loc-Tag field. The model can then be displayed using any time step using the legend editor of ArcView. In addition, the identify tool can be pointed at any cell and the complete drawdown pattern for that cell is displayed. Maximum and average drawdowns can be computed for each time step using the field menu on the attribute table. Drawdowns greater than a specified criteria can be displayed using the query tool. Thus a complete model simulation can be displayed and analyzed in one View in ArcView. Figure 5 shows a typical model run displayed on ArcView.


Sample Display

Figure 5 - Display of Drawdown Impacts of Aurora Pumping

Aurora Pumping 4000 AF/Year - Drawdowns are in feet


Enhancements to ArcView using Avenue

In addition to the standard features of ArcView, several functions were created using Avenue for this project. Three scripts were developed to control the legend editor. One allowed the user to specify specific steps in the legend, i.e. show the drawdown in even one foot increments. Another allowed ramping of colors between specific values in the legend so negative drawdowns (or higher elevation than the starting elevation) could be graduated shades of one color while positive drawdowns could be graduated shades of another color. Finally, a third script allowed grouping several categories together into one legend (grouping all negative drawdowns together as the positive drawdown is of more concern). All thes scripts were useful in this project, however they will all be standard features in ArcView 3.0
Animation is also possible through the use of Avenue. By sequentially displaying the groundwater model one stress period at a time an animated picture of the effects of pumping can be displayed. This allows an clear visualization of the reaction of the aquifer during the model simulation.
One enhancement worth noting specifically is the use of "gnuplot" to graph both model results and monitoring hole data on the same graph. The program "gnuplot" is a copyrighted but free command-driven interactive function plotting program which runs on most UNIX platforms and Windows machines. The tremendous benefit of "gnuplot" for this application instead of the charting functions built into ArcView is "gnuplot" allows time-series data and data from separate tables to be plotted on one graph. Some of the model runs in the analysis of Cherry Creek are calibration runs tied to specific dates and streamflow and pumping runs. By using "gnuplot" the data from both of these sources can be graphed together and compared for calibration purposes. In the ArcView application a tool is created where the active themes are selected and the data is read from the respective tables. The data is then written to an ASCI file in the command format of "gnuplot". Finally, "gnuplot" is launched using the data file as input. Figure 6 shows a sample of "gnuplot" displaying the data from a model run. Copies of the Avenue Scripts are included in the appendix.


gnuplot demo

Figure 6 - gnuplot of Modelled Results vs Monitoring Hole Data


Summary and Conclusions

Through the use of ArcView 2.1, the City of Aurora can now analyze and display complex groundwater model results from the USGS's MODFLOW program. The ease of use and clear display of the effects on the Cherry Creek aquifer allow the ability to analyze more simulations and "What-If" scenarios than would be possible through conventional analysis. The areas of most significant impact are visually obvious even of highly complex models. This application is the basis for a water rights application and a basin-wide management plan for the upper Cherry Creek basin.

Aknowledgements:

ArcView and ArcInfo were developed by Environmental Systems Research Institute, Inc. (Esri)

MODFLOW was developed by Michael G. McDonald and Arlen W. Harbaugh for the USGS

The streamflow package for MODFLOW was developed by David E. Purdic for the USGS

The paticular model described in this paper was developed by Bishop-Brogden and Associates, Lakewood, Colorado under contract with the City of Aurora.

gnuplot was developed by Thomas Williams and Collin Kelley
David Denholm is the coordinator of gnuplot 3.6
Major contributors to gnuplot (alphabetic order):

  • John Campbell
  • Robert Cunningham
  • David Denholm
  • Gershon Elber
  • Roger Fearick
  • Carsten Grammes
  • Thomas Koenig
  • David Kotz
  • Ed Kubaitis
  • Russell Lang
  • Alexander Lehmann
  • Carsten Steger
  • Tom Tkacik
  • Jos Van der Woude
  • Alex Woo

    Appendix

    Sample Source Code to translate MODFLOW binary files to test tables:
    
    #include 
    #include 
    
    main(argc,argv)
    int argc;
    char *argv[];
    {
    
    	int i,j,c,d;
    	int sign,k;
    	float z;
    	float store[72][244][26];
    	FILE *fp;
    
    	fp = fopen(argv[1],"r");
    	for(i=0;i<56;i++)   /* 56 bytes in the header - machine dependent */
    		c=getc(fp);
    	printf("\"Loc_Tag\"");
    	for(k=0;k<26;k++) {
    	printf(",\"%d\"",k+1);
    	for(j=0;j<244;j++) {     /* 244 rows by 72 columns - model specific */
    		for(i=0;i<72;i++) {
    			fread(&z,sizeof(float),1,fp);
    			store[i][j][k] = z;
    		}
    	}    
    	for(i=0;i<60;i++)   /* 60 bytes between stress periods - machine dependent */
    		c=getc(fp);
    	}
    	fclose(fp);
    	printf("\n");
    	for(j=0;j<244;j++) {
    		for(i=0;i<72;i++) {
    			if(store[i][j][0] < 999.9) {  /* 999.9 is a flag in MODFLOW that this cell */
    			printf("%d-%d",j+1,i+1);      /* is not active - this reduces the file size */
    			for(k=0;k<26;k++)
    				printf(",%.3f",store[i][j][k]);
    			printf("\n");
    			}
    		}
    	}
    }
    
    Sample Avenue Script to animate groundwater model:
    
    if(av.GetActiveDoc.GetActiveThemes.Count > 1) then
      MsgBox.Info("Select only one Theme!","ERROR")
      return -1
    end
    theView = av.GetActiveDoc
    theChart = av.GetProject.FindDoc("Discharge along Stream")
    theCFList = theChart.GetFields
    theCVTab = theChart.GetVTab
    cList = theCVTab.GetFields
    theTheme = theView.GetActiveThemes.Get(0)
    theTName = theTheme.GetName
    aLegend = theTheme.GetLegend
    theFTab = theTheme.GetFtab
    aList = theFTab.GetFields
    aField = MsgBox.List(aList,"Select the First Field","Legend")
    theFF = aList.Find(aField)
    aField = MsgBox.List(aList,"Select the Last Field","Legend")
    theLF = aList.Find(aField)
    if(theCFList.Count > 1) then
      for each i in 1..(cList.Count-1)
        theChart.GetFields.Remove(i)
      end
    end
    for each i in theFF..theLF
      aLegend.SetField(aList.Get(i))
      aLegend.Load("/resource/aview/legends/drawdown.avl".asFileName)
      theTheme.SetName(aList.Get(i).GetAlias)
      theTheme.SetVisible(False)
      theTheme.UpdateLegend
      theTheme.SetVisible(True)
      theView.Draw(theView.GetDisplay)
      newField = aList.Get(i)
      CVTindex = cList.Find(theCVTab.FindField(newField.GetAlias))
      theChart.GetFields.Set(0,cList.Get(CVTindex))
      theChart.UndoErase
    end
    theTheme.SetName(theTName)
    
    Sample Avenue Script to launch gnuplot:
    
    theBitmap0 = _theVTab0.GetSelection  'Globals are set in another script
    theBitmap1 = _theVTab1.GetSelection
    lf = LineFile.Make("/tmp/tmp.gnu".AsFileName,#FILE_PERM_WRITE)
    
    lf.WriteElt("set title 'Cherry Creek Groundwater Model - 1995 Results'")
    lf.WriteElt("set xlabel 'Date'")
    lf.WriteElt("set ylabel 'Elevation'")
    lf.WriteElt("set timefmt '%m/%d/%y'")
    lf.WriteElt("set xdata time")
    lf.WriteElt("set format x '%m/%d/%y'")
    lf.WriteElt("set grid")
    lf.WriteElt("set key left")
    lf.WriteElt("plot \")
    for each r in theBitmap0
      lf.WriteElt("'-' using 1:2 t '"+_theVTab0.ReturnValue(_theFieldList0.Get(_theNF0),r).AsString+"' w l,\")
      lf.WriteElt("'-' using 1:2 t '' w p,\")
    end
    for each r in theBitmap1
      lf.WriteElt("'-' using 1:2 t '"+_theVTab1.ReturnValue(_theFieldList1.Get(_theNF1),r).AsString+"' w l,\")
      lf.WriteElt("'-' using 1:2 t '' w p,\")
    end
    lf.WriteElt("'-' using 1:2 t '' w l")
    for each r in theBitmap0
      for each i in _theFF0.._theLF0
        aField = _theFieldList0.Get(i)
        x = aField.GetAlias
        y = _theVTab0.ReturnValue(aField,r)
        if(y = 0) then
          lf.WriteElt(x)
        else
          lf.WriteElt(x++y.AsString)
        end
      end
      lf.WriteElt("end")
      for each i in _theFF0.._theLF0
        aField = _theFieldList0.Get(i)
        x = aField.GetAlias
        y = _theVTab0.ReturnValue(aField,r)
        if(y = 0) then
          lf.WriteElt(x)
        else
          lf.WriteElt(x++y.AsString)
        end
      end
      lf.WriteElt("end")
    end
    for each r in theBitmap1
      for each i in _theFF1.._theLF1
        aField = _theFieldList1.Get(i)
        x = aField.GetAlias
        y = _theVTab1.ReturnValue(aField,r)
        if(y = 0) then
          lf.WriteElt(x)
        else
          lf.WriteElt(x++(y+5000).AsString)
        end
      end
      lf.WriteElt("end")
      for each i in _theFF1.._theLF1
        aField = _theFieldList1.Get(i)
        x = aField.GetAlias
        y = _theVTab1.ReturnValue(aField,r)
        if(y = 0) then
          lf.WriteElt(x)
        else
          lf.WriteElt(x++(y+5000).AsString)
        end
      end
      lf.WriteElt("end")
    end
    lf.WriteElt("end")
    lf.Close
    System.Execute("/usr/local/bin/gnuplot -persist /tmp/tmp.gnu")
    

    Author

    Bruce Rindahl, Senior Water Resources Engineer (Primary & Presenting Author)
    City of Aurora, Utilities Department
    1470 South Havana Street
    Aurora, CO 80012
    Phone: (303) 695-7383
    Fax: (303) 695-7491
    E-mail: brucer@dilbert.ci.aurora.co.us