button and select the
project gfmap.apr in the Demo Startup window.
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 \exercise\ex5\gisfiles directory.
Internet students can download the necessary files.
Site: ftp.crwr.utexas.edu
Login: anonymous
Password: your e-mail address
Directory: /pub/crwr/gishydro/africa/ex5
File(s):
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.
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)
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]
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
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.
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.
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.