FAO/UNESCO Water Balance of Africa

Groundwater Flow Between Two Rivers

Prepared by
Kris Martinez , David R. Maidment and Kwabena Asante
Center for Research in Water Resources
The University of Texas at Austin
November 1996

Table of contents

Goals of the Exercise

This exercise explores the use of a map-based simulation model in determining water levels and flow rates in an unconfined aquifer located between two rivers. It enables the user to develop a better understanding of possible applications of GIS in Groundwater hydrology. It also illustrates the use of a map-based approach to quantifying the interaction between surface water and groundwater in a shallow aquifer.

Computer and Data Requirements

This exercise was developed on a customized version of the Geographic Information System (GIS) software, ArcView Version 2.1. The installed ground water module was developed by Zichuan Ye at the Center for Research in Water Resources. The coverages and output files required for the model have already been preprocessed and are available in the directory, rivers/ex4/gisfiles.

Coverages include: Tables required by the model include:


In this exercise, water levels in an unconfined aquifer will be calculated using Darcy's law for a given set of conditions. The conditions will then be simulated using a map-based groundwater flow model. The program creates a visual representation of the water level in the aquifer and its movement. Calculated and model results will then be compared. The program used for this exercise

The case to be studied involves a strip of land 50 km wide that is bounded by two parallel rivers (Bear, 1979). The rivers partially penetrate an unconfined aquifer. Initially, the river levels and the aquifer water level are at an elevation of 50 m. There is an impermeable bed at an elevation of 0 m. A uniformly distributed recharge of 0.5 mm/day occurs on the surface.

The volumetric flux in a porous medium can be found using Darcy's law,
q = -K(dh/dx)   (1)

q is the volumetric flux [m/s]
h is the elevation of the flow [m]
x is the length in the x-direction [m]
K is the hydraulic conductivity [m/hr]

For a unit width, the volumetric flow rate (Q) is q*h [m^2/s].
Incorporating Darcy's law,

Q = -K(dh/dx) * h   (2)
A mass balance for an incremental volume in the x-direction yields,
Q(in) + N*x - Q(out) = (dh/dt)*x   (3)

Q(in) is the incoming volumetric flow rate [m^2/s]
Q(out) is the outgoing volumetric flow rate [m^2/s]
N is the uniform recharge rate [mm/day]
x is the incremental length [m]

For a differential control volume we get the equation,
d/dx(K*h(dh/dx)) + N = dh/dt   (4)

For steady-state conditions (dh/dt=0), the equation can be solved to obtain a relationship describing water elevation along the length of the aquifer,

K(h^2-h(0)^2) - N*x(L-x) + K*x/L(h(0)^2-h(L)^2) = 0   (5)

h is the water elevation [m]
h(0) is the elevation of the river to the west [m]
h(L) is the elevation of the river to the east [m]
x is the distance along the aquifer from the west river [m]
L is the total length of the aquifer [50000 m]

Differentiating with respect to x, we get,

-Q(x) = K(dh/dx) * h = N(L/2-x) - K/2L(h(0)^2-h(L)^2)   (6)

A positive Q means that there is flow in the +x direction (towards the East river).

The location of the peak water level in the aquifer can be found by setting dh/dx equal to zero and solving for x,

x = L/2 - K/(2*N*L) (h(0)^2-h(L)^2)   (7)

We will now use these equations to predict the response of the groundwater flow model for the given conditions.


1. Calculating the steady-state water levels

Using equation (5), calculate the water levels at points 5 km, intervals across the aquifer. In the case previously described,

The recharge value can be divided by (1000 mm/m) and (24 hr/day) to convert into (m/hr).

Note that the aquifer water level is the same at 15 km and 35 km. This is a result of having both rivers at the same level. We will look at the effects of different river levels a little later.

2. Calculating the steady-state flow rates

Using equation (6), calculate Q(x) at 5 km intervals across the aquifer.

The flow rates at 15 km and 35 km are of equal magnitude but opposite direction. The rate increases as the water moves out towards the rivers. In this case, the aquifer is at a higher water level and acts as a source.

Things To Do: Prepare a table showing the piezometric head and the discharge flowing toward the river at 5 km intervals (i.e. 0, 5, 10 15, 20, 25, 30, 35, 40, 45, 50 km) across the aquifer. Points 0 and 50 km are at the two bounding rivers. Perform one sample calculation for the head and one for the discharge. Prepare a graph showing the piezometric head elevation as a function of distance between the rivers and one showing the magnitude of the discharge as a function of distance between the rivers.

Now, let's use the model to perform the same computation and compare the results of the two approaches.

3. Setting up the model

Create a workspace on your local computer by making a new directory, say gfmodel, where you will store your files and open a MS-DOS command prompt window. Copy the files from the CD into the workspace by typing the following command at the DOS command prompt. (Substitute your CD drive for D in the command below if the CD drive on your computer is not D).

xcopy D:\rivers\ex4\gisfiles C:\gfmodel\gisfiles /e

where the D drive is the directory containing the CD-ROM and C drive is the drive containing the newly created gfmodel directory. Copying files in this manner ensures that the permissions on the files are changed from the Read Only (automatically assigned to files on the CD-ROM) to Read-Write. If the file manager or the copy command is used to copy these files, the permissions on the individual files must be changed to allow Write access.

4. Running the model

Start ArcView and open the project ('gfmap.apr') from the file menu. A view called Union containing a 5 x 5 grid of the site comes up. Each of the squares represents a 100 km^2 area of land. The entire grid represents an area 50 km from south to north and 50 km from west to east (2500 km^2). The rivers in the example are represented by the outer left and right boundaries of the grid. The program has been set up to model the conditions we have been working with.

Select GFwModel from the GFlowSim menu of the View Menu bar.

Make sure that the fields match the menu below and click .

The meaning of these parameters is described in the Appendix at the end of this exercise.

Set the number of iterations to 400 and click again. Just press return for the memory messages that follow. The model will take a few minutes to run. At the bottom of the screen several diagnostic numbers are displayed, whose values are described in the Appendix.

5. Plotting the results.

Once the model has finished running, click the button to plot the steady-state flow rates and water levels in the aquifer. Make sure the fields in each of the subsequent menus match the parameters shown below and click .

Click on anywhere on the grid.

The results are shown below.

The red vectors represent the magnitude and direction of flow rates across the lines. The numbers are the water levels in the center of each of the 100 km^2 areas of land.

As you can see, the water levels from the model are pretty close to what we calculated. The values of the discharge on the lines can be obtained by highlighting the theme Gfmap.lin, and using the information tool to click on each line. In the table displayed the discharge towards the river is given by xflux in m3/day for the line. Since each line is 10km long, the discharge for a 1m wide section of aquifer can be obtained by dividing the value you obtain from the query by 10,000. This is the discharge which is comparable to the values obtained in the previous part of the exercise.

After you have looked at the graph, close the project and exit ArcView without saving the changes.

Things To Do: Prepare a table showing the comparison between the theoretical values of the heads and discharges and those determined by the Arcview model. Note that the heads are available only at the cell centers, at 5, 15, 25, 35, 45 km while the discharges are available only on the lines at 0, 10, 20, 30, 40, 50 km. Find the % difference between the computed and observed values. Where is the difference the greatest? Why is that so?

6. Analyzing the effects of varying water levels

Create a spreadsheet to calculate the steady-state water levels at various points along the aquifer using equation (5). Calculate the levels for three different cases, H(L) = 50 m, 25 m and 0 m. A graph showing results is shown below.

Notice that the peak water level decreases as the river levels decrease.

Now, calculate the location of the peak water levels for each of the three cases. The peak location moves west as the level in the east river decreases. These effects can be seen in the graph below. The markers represent the position of the peak water level.

Things To Do: Prepare a graph of the piezometric head profile in the aquifer like the one shown above for the three cases for H(L). What proportion of the flow goes to the left as compared to the right hand river for each of the cases?

7. Visualizing the effects of a continuous drawdown

Pumping from an aquifer will cause water levels to drop over time. We will now use the model to simulate a continuous drawdown for the aquifer previously described.

Open the project using the same procedure as in step four. The pumping rate in the cell at the center of the grid is to be changed. Click on the blue square and the inner square of the icon. Then click on to switch to selection mode. Select the middle square in the grid and it will turn yellow.

To change the conditions for the selected area, we need to edit the attribute table. Click on the button to open the table, Attributes of Gfmap.ply. Notice that the record corresponding to the selected square is also highlighted in the table. Go to the Table menu and select "Start Editing".

Lets assume the rate of pumping has now changed to 1 m^3/hour. Select the [Pmp0] icon and then click the field calculator button, . Type 1.0 in the input box under [Pmp0]= and click . Then go to the Table menu again and selected "Stop Editing". Select YES when prompted on whether or not to save edits. The table has been updated, now we can run the model.

Select GFlowSim from the GFwModel menu as before. Make sure the input parameters match the parameters given below and press enter.

Set the number of iterations to 50 and press enter for the memory messages that follow.

The model should now be plotting the flow rate vectors and water levels for each time step. Notice that water is flowing into the center square to meet the pumping demands while the water level continues to drop. Pretty cool!! If you want to stop the computations at any time just hit the "Stop" button in the bottom left hand corner of the screen. After the simulation finishes, close the project and exit ArcView without saving the changes. You are now finished!

Technical Appendix Defining the Model parameters

1. Initial Dialog Box

The parameters in this dialog box represent the following:

fluxplot: the ratio of the groundwater flow rate to the length of the flow arrow on the plot. The default value is 0.0001. Increasing this number makes the arrows longer, decreasing it makes them shorter.

bndfact: this factor adjusts the line parameter Slength for the flux computation at the boundary line to allow for the fact that the boundary has no width while all the interior cells of the grid have finite width. The factor 1.21 was obtained so as to match the numerical solution obtained for this problem to its theoretical solution.

Sctrl: this is a global value for specific storage of the groundwater system. This value is multiplied by the value (0.1) of SV in the Attributes of Gfmap.ply table to determine the specific storage of the unconfined flow system. Hence in this case, the specific storage = 0.01515 * 0.1 = 0.001515.

Yescolor this is a code that governs what is plotted on the View.

0 no items are plotted on the screen

1 the boundary line is highlighted when the flow across that line is being computed

2 the polygon is highlighted as the flow balance is done in that cell

3 both the lines and polygons are highlighted

5the flow vector and piezometric head levels are plotted

ISSteady a flag that signals steady or unsteady flow computation:

true a steady state computation is done and the governing values are taken from the attribute tables of the coverages.

false an unsteady state computation is done and the governing values for some of the variables are taken from time series tables whose table names are the same as the attributes of Gfmap.ply defined below. For example headtb.dbf stores time series of computed heads in m.

2. Attributes of Gfmap.ply

In theme Gfmap.ply the attributes are:

KV1000 hydraulic conductivity in m/hr

Head0 initial head in the cells (m)

Rch0 recharge rate in mm/DT, where DT is the time step of the computation. DT is 1 day in this computation.

Spr0 spring flow rate in m^3/s

Pmp0 pumping rate in m^3/s

gbh0 a value for head in a general head boundary cell. When the value is zero, the head can fluctuate according to flow conditions. When the value is nonzero, the head in this cell is held fixed at the specified level.

evt0 is the evaporation rate in mm/DT where Dt is 1 day in this computation

Bot elevation of the bottom of the aquifer (m)

Top elevation of the top of the aquifer (m).

Cnfd is a flag for confined and unconfined flow. When Cnfd = 0 the flow is considered unconfined and the level Top is not considered. When Cnfd = 1 the flow is considered confined and the piezometric head level is checked against Top to verify that the flow is confined.

SV1000 is a default value that is multiplied by Sctrl to give the specific storage.

headn is the piezometric head in the aquifer as computed in the last program iteration (m)

dvol is net outflow from the cell across its four boundary lines at the last time step. (m^3/DT) where DT = 1 day in this case.

3. Attributes of Gfmap.lin

ldx length of the line in the x-direction

ldy length of the line in the y-direction

fcosx the cosine of the angle formed by the flux across this line and the x-axis

fcosy the cosine of the angle formed by the flux across this line and the y-axis

CCX, CCY are the coordinates of the center point of this line

Slength is the length to be applied when computing fluxes across this line. Slength = cell size when the line is interior to the mesh and Slength = 0.5 * cell size when the line is on the boundary.

isbnd is a flag for boundary lines.

1 means the line is a boundary,and the interior of the domain is to the right of the line, where the direction of the line is from its from-node to its to-node

0 means the line is not a boundary

-1 means the line is a boundary and the interior of the domain is to the left of the line, where the direction of the line is from its from node to its to-node.

bndtp is a boundary type

0 means nonflow boundary

1 means constant head boundary

Bhead is the head on a constant head boundary

xflux, yflux are components of the flow across the line in m^3/DT where DT = 1 day in this computation

4. Program Bottom Bar

When the program is running the bar along the bottom of the screen indicates the degree of progress of the computations and reports some summary numbers:

T = the number of time steps (days)

MCT = the maximum Courant number in the cells (should be less than 1)

MX/m = the maximum and minimum change in head that has occured in any time step so far in the computations, in m e.g 0.66/0 means that the maximum change in any time step was 0.66m and the minimum was 0.

maxhead|minhead the maximum and minimum head at any cell center in the mesh in m. This number is followed by the maximum difference between the heads at the previous time step and those at the current time step.

These materials may be used for study, research, and education, but please credit the authors and the Center for Research in Water Resources, The University of Texas at Austin. All commercial rights reserved. Copyright 1997 Center for Research in Water Resources.