Exercise 4.
Groundwater Flow Between Two Rivers

CE 394K Surface Water Hydrology
University of Texas at Austin

Prepared by Kris Martinez and David R. Maidment

November 1996

Table of contents

Goals of the Exercise

Computer and Data Requirements

This exercise is intended to be performed with Geographic Information System (GIS) program ArcView, version 2.1.

In addition to ArcView, you will need some programs packaged into a project which you can load into ArcView like a file.


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 was written by Ye Zichuan of the Center for Research in Water Resources at the University of Texas at Austin. Thank you, Ye.

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. Calculate 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. Calculate 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.

To be turned in: 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. Show one example calculation for the head and one for the discharge. Don't just turn in a spreadsheet output. 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 compare our results with the model.

3. Get the necessary files.

Open 'Filemanager' and drag the directory called 'gfmodel' to your workspace. The 'gfmodel' directory is located on the G: drive under \Class\MAIDMENT\ce394k

Internet students can download the necessary files.

Unzip the file with the UNIX unzip command.

4. Start ArcView and open the project.

Once ArcView is active, open the project ('start.apr') from the file menu. The program will ask you to enter a directory, this is the drive and directory ('gfmodel') where you copied the files. After entering this information, click on the button and select the project ('gfmap.apr') from the file menu. The Start.apr project is a short project that has the function of searching through the file for the Gfmap.apr file and repathing all the file pointers to the correct location where you have the files. These pointers are programmed into the .apr files so that the next time you open the project it will look for its component files in the same places. If you have trouble operating the Start.apr project, your files will not be repathed automatically and you'll have to do that manually by finding and clicking on each file required by Gfmap.apr as Arcview prompts you to do so.

5. Run the model.

In the main window ('union') you will see a 5 x 5 grid. 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 our example are on the left and right side of the grid. The program has been set up to model the conditions we have been working with.

From the menu, select GFwModel -> GFlowSim

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.

6. Plot 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 match in each one of the subsequent menus and click the buttons.

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.

To be turned in: 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?

7. Analyze 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.

To be turned in: 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?

8. Visualize the effects of a continuous drawdown.

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

9. Run the model.

Open the project using the same procedure used in step four. We need to change the pump rate for the center square. 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 edit the attribute table. Click on the button. Go to the Table menu and select "Start Editing".

We will change the pump rate to 1 m^3/hour. Select the [Pmp0] icon and then click the button. Enter 1.0 under [Pmp0]= and click . Then go to the Table menu again and selected "Stop Editing" The table has been updated, now we can run the model.

From the menu, select GFwModel -> GFlowSim. Make sure that the fields match the menu 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.