David R. Maidment and Kris Martinez

Center for Research in Water Resources

University of Texas at Austin

November 1996

- Goals of the Exercise
- Computer and Data Requirements
- Introduction
- Mathematical Model
- Procedure
- Start ArcView and open the project.
- Running the Script for the Steady-State Condition.
- Running the Script for a Continuous-Drawdown (Unsteady-State) Condition.
- Technical Appendix Defining the Model Parameters

To develop a better understanding of groundwater hydrology, by visualizing flow direction and magnitude, and the piezometric head changes in an unsteady flow situation.

This exercise requires the ArcView GIS, version2.1 or later, developed
by the Environmental Systems Research Institute (ESRI). Additionally, you
will need the ArcView projects **start.apr** and **gfmap.apr** which
include the scripts that model the groundwater flow. The files required
for this exercise have already been loaded onto your computer and can be
found in the \gisfiles\ex5af\gfmap directory.

Internet students can download the necessary files.

- gfmap.zip (binary)

Site: **ftp.crwr.utexas.edu
**Login:

In this exercise, piezometric heads and flows will be calculated using a map-based groundwater flow model that applies simultaneously Darcy's law and the continuity equation. The program also provides graphic display of the calculated values in the form of arrows that represent flow direction and magnitude, and of numbers that represent piezometric heads. Calculated and observed values are then compared.

The case to be studied (Bear, 1979) involves a 50 km wide unconfined aquifer, that lays over an impermeable bed at an elevation of 0 m, bounded by two parallel rivers . Initially, the river levels and the aquifer water level are at an elevation of 50 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, and it is equal to

q = -K(dh/dx) (1)

where:

q is the volumetric flux [L/T]

h is the piezometric head [L]

x is the distance in the x-direction [L]

K is the hydraulic conductivity [L/T]

For a unit width, the volumetric flow rate (Q) = q*h [L^2/T].

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)

where:

Q(in) is the incoming volumetric flow rate [L^2/T]

Q(out) is the outgoing volumetric flow rate [L^2/T]

N is the uniform recharge rate [L/T]

x is the incremental length [L]

- For a differential control volume we get the equation,

d[K*h(dh/dx)]/dx + N = dh/dt (4)

For the steady-state conditions (dh/dt=0), the equation can be solved to obtain a relationship describing piezometric head 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)

where:

h is the piezometric head elevation [L]

h(0) is the elevation of the river to the west [L]

h(L) is the elevation of the river to the east [L]

x is the distance along the aquifer from the west river [L]

L is the total length of the aquifer [L]

Differentiating equation (5) with respect to x, and combining it with equation (2), we get,

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

in which a positive Q implies that the is flow in the positive x direction.

The location of the highest point of the water table 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 system.

Using equation (7), calculate the water levels at points 15 km, 25 km and 35 km along the length of the aquifer. Recall that h(0) = 50 m, h(L) = 50 m, K = 1 m/hr, N = 0.5 mm/day and L = 50 Km. 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.

Using equation (6), calculate Q(x) at points 15 km, 25 km and 35 km along the length of the aquifer.Note that, because of symmetry, 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.

Once ArcView is active, click **File/Open Project** and open the
project **start.apr** located in the **\gisfiles\ex5af\gfmap** directory.You
will be prompted to enter the directory where the files are located. This
directory is also **\gisclass\ex5af\gfmap**. After entering this information,
click on the button and select the
project **gfmap.apr** in the Demo Startup window.

In the view window **Union**, you will see a 5 x 5 grid. Each of
the squares represents a 100 km^2 area, and the entire grid represents
an area 50 km from south to north and 50 km from west to east (2500 km^2).
In this example, the rivers are on the left and right sides of the grid.
The project has already been set up with the hydrologic conditions described
above, in which the water level in both rivers is 50 m.

From the menu, select **GFwModel/GFlowSim**. Make sure the parameters
in the "Parameters for Gflow" window match those shown below
and click .

Set the number of iterations to 400 and click again. Just press OK for the memory messages that follow. The script will take a few minutes to run.

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 cell border. The numbers are the piezometric head levels in the center of each of the cells. As you can see, the water levels according to the ArcView calculations are pretty close to those obtained with equations (5) and (6).

Close the project and exit ArcView without saving the changes.

To model the condition in which the water level in the rivers is not
the same, the boundary conditions have to be changed. To do this, make
the project main window **Gfmap.apr** active, click on the **Tables**
icon, select **Gfmap.lin** from the table list, and click the **Open**
button. Once the **Gfmap.lin** attribute table is open, make the **Union**
view window active, and the **Gfmap.lin** theme active. Next, select
the lines that correspond to the river its water level you want to change.
Make the **Gfmap.lin** attribute table active again, and select the
**Bhead** field by clicking on its header. Go to **Tables/Start Editing**
to make the table editable, and then click on the calculator icon. You
will see a dialog window that prompts you for a value. Enter the water
level of the river you selected and click OK. You may want to follow the
same procedure to change the water level of the other river. After having
changed the water levels in the rivers, run the script as you did before.

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.

Open the project using the same procedure used before. 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** field
by clicking on its header and then click the
button. Enter 1.0 under [Pmp0]= and click .
The table has been updated.

To run the script, go to **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. A
tabulated version of the results can be found in the table **headtb.dbf**.
To see them, make the project main window **Gfmap.apr** active, click
on the **Tables** icon, select **headtb.dbf** from the table list,
and click the **Open** button. In this table, the headers of the columns
refer to the terrain cell, being G14 the cell where the pump is located.

After the simulation finishes, close the project and exit ArcView without saving the changes. You are now finished!

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

**5**the 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.

In theme Gfmap.ply the attributes are:

**KV** hydraulic conductivity in m/s

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

**SV** 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.

**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

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.

Go to Index of Exercises

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