.
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.
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: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.
q = -K(dh/dx) (1)
where:
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)
Q(in) + N*x - Q(out) = (dh/dt)*x (3)
where:
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]
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)
where:
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.
Using equation (5), calculate the water levels at points 5 km, intervals across the aquifer. In the case previously described,
h(0) = 50 m
h(L) = 50 m
K = 1 m/hr
N = 0.5 mm/day
L = 50000 m
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.
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.
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.
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.
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?
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?
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!
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:
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.
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.