Chapter 8 Two demensional (Matrix) groundwater model

Matrix model

8.1

Chapter 8

This chapter deals with a simple groundwater head model that can approximate head changes in an infinite confined aquifer. The changes in head are assumed to be due to discharges or recharges from a number of wells who's locations have been given by the user. Also the duration of operation of all the wells must be specified. In this model, although initial heads may be specified, flows due to head differences in the aquifer are not considered; changes in head are calculated only on a summation of drawdowns The model should be useful in due to the discharging/recharging wells. getting a first approximation of the changes in an existing piezometric surface, due to a superimposed pattern of drawdown cones. 1. GROUNDWATER MATRIX MODELS There are two major classes of models dealing with a matrix of piezometric heads, depending on the method used to obtain a solution. 1/ Analytical, where an exact solution is at least theoretically obtainable, given that the aquifer fits the conceptual model exactly. (In practice it never does.) 2/ Numerical, where in effect, numbers are allowed to flow through the model similarly (it is hoped) to the way water flows through the real aquifer, or system of aquifers. Here the aquifer is divided up into smaller units (called elements) which are assumed to have constant hydraulic properties within themselves. An equation approximates the flow between each element and it's neighbours over the given time span, and for the given initial conditions. The equation constants vary from one element junction to the next depending on the values given to transmissivity, and storage coefficient for the elements on either side of the junction. The equations used for flow to and from elements along the boundary of the model are slightly different to those relating to the interior of the model. All of the equations must be solved simultaneously. This book will not give any examples of numerical models. If the reader wishes to pursue the subject he might try Kinzelbach's Groundwater Modelling, (1986) which includes programs in Basic. As mentioned in the introduction, in all real groundwater problems where numbers are given to aquifer properties, we are dealing only with approximations to reality. This becomes especially apparent in groundwater matrix models. The art is to fit the model to reality sufficiently closely so that the projections it gives will be similar to the behaviour of the real system, and, at the same time, to keep it as simple as possible. (A generous application of Ochams razor is called for.) Matrix models can also be divided into two classes depending on the behaviour of the flow system in time. 1/ In the steady state class, it is assumed that all changes to the flow system have occurred sufficiently far into the past for all heads aoo

8.2

Matrix model

flows to reach an equilibrium. That is, not only does flow into the aquifer system equal flow out of the system, but flows into and out of any definable part of the system are also equal (this excepts points of recharge or discharge) . The steady state class of model has a limited application because in practice it is changes in an aquifer system that cause problems, and therefore it is a changing system that must be modelled. 2/ So we come to the unsteady, or transient, state. In models of this type time is a variable that must be considered. Depending on the model, it may be the time from the beginning of discharge from a well or series of wells, or it may be the time over which groundwater has been flowing through an aquifer under the influence of unequal piezometric heads. Other examples are the duration of the leakage of a contaminant from a land fill, or of saline groundwater from the raised water table beneath an irrigation area. The model described in this chapter is an analytical model of transient state discharge or recharge. Since flows due to variations in the preexisting piezometric surface are not considered, the model is neither fully transient state, nor fully steady state, but is a hybrid type. 2. PROGRAM ANMODL The basis of this program is in the fact that the change in piezometric surface at a point in an aquifer under the influence of several discharging or recharging wells, is the sum of the drawdown due to each well individually. Puting it another way, suppose at a given time, t, the water level in a piezometer has changed by A units due to discharge from well A1. At some other time drawdown in the piezometer under the influence of another well, B1 is measured as B units, and on a third occasion, drawdown due to discharge from well C1 is found to be C. (In all three cases the pumps have operated for the same length of time before the measurement of the drawdown.) If all three pumps start now at the same time, and go on discharging for the same length of time as before, the drawdown in the piezometer should be A+B+C. 2.1. Program verification Load and run the program. You will first be asked if you want to use the printer. I suggest answering in the affirmative for this run. The computer will next display "Size of the array", and then ask "No. of nodes in x direction?" Answer with 8. Answer with 10 to the question "No. of nodes in the y direction?" The program will use this information to dimension three 8 by 10 matrices, one for heads, one for discharges/recharges, and the third for times; ie. duration of discharge/recharge. Give 100 (metres) as the spacing between (There is no need for these two nodes in both the x and y direction. distances to be the same in other cases.) Transmissivity and storage coefficient will be asked for; reply With 40 cubic metres per day per metre and 0.0005 respectively. The time unit to use is days, so press d or D at the prompt.

Matrix model

8.3

Now we have come to the part of the program that allows the setting up of the matrices. This part will require some care, although you can recover from any mistakes by causing the program to write over the erroneous data. In this example run you will be given data to enter, but in practical applications it would be all but essential to have a drawing of the matrices that you want to set up laid out in front of you. You will notice that you are given the option of setting up either the initial head, or discharge/recharge matrix first. For the example run leave the initial heads, they will default to zero right across the matrix. Press 2 to set up the discharge/recharge matrix. Menu 2 will now be displayed, giving the options of several groupings for the setting up of parts of the matrix. Choose the column option by pressing 2. Now you must define which column it is that you are interested in, the portion of that column, and the discharge rate for each node within that portion. Type 2,4,6,300. This indicates to the computer that you want to set discharge rates in column 2, which is the third column from the left as numbering starts at zero. The 4 and 6 indicate that within column 2 you want to set up three nodes, numbers 4, 5, and 6. ie. You will put discharge points at elements (2,4), (2,5), and (2,6). Each of the discharges will be (Note that even if the time unit, minute was 300 metres cubed per day. chosen, discharge rates would still be based on days. The minute would only apply to the unit used for duration of discharge/recharge. Press enter. The program checks for errors in entries at this stage, but may not pick up all operator mistakes. For example, entry of a negative number for duration of discharge will probably lead to some Basic error message. No difficulty should be experienced if a little care is used. After accepting the entry of the discharge rate, the program will request the duration of the pumping. Reply with 10 (days). Before requesting a printout, enter a last single recharge node by pressing 1, and then typing 7,8,-500. [ie. A recharge rate of 500 cubic metres per day at element (7,8).] Give a duration for this recharge of 5 days. Now, if all has gone well, you have set up the matrices as shown below. Get a printout and check by pressing 5 in menu 2. Table 8.1; 0.(\0 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

The discharge rate matrix of the example run. 0.00 0.00 0.00 300.00 300.00 300.00 0.00 0.00 0.00 0.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.00 0.00 R 9 0.00 -500.00 R 8 0.00 0.00 R 7 0.00 0.00 R 6 0.00 0.00 R 5 0.00 0.00 R 4 0.00 0.00 R 3 0.00 0.00 R 2 0.00 0.00 R 1 0.00 0.00 R 0

8.4

Matrix model Table 8.2;

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

The duration of operation of discharge/recharge (days). 0.00 0.00 0.00 10.00 10.00 10.00 0.00 0.00 0.00 0.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

0.00 5.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00

R9 R8 R1 R6 R5 R4 R3 R2 R 1 R0

After printing the tables above, you will be returned to the main menu. Choose option 3, "finished" (with setting up). The calculations of drawdowns due to each discharging node, and the summation of all such drawdowns at each node will now commence. The total drawdowns will appear on the screen, and on your printer as they are calculated. Depending on what printer you have, it may wait to build up a line of type before printing anything, or it may print as it receives data. There is a special case in the drawdown calculations which should be mentioned. In the example data given above, for the calculation of drawdown at node (2,4) due to discharge from node (2,4), a distance of zero will be calculated. If this is fed into the Theis equation subroutine so that the drawdown may be calculated, an error will result. Whenever a distance of zero is detected then, it is changed to a quarter of the distance between nodes in the x direction, in this case it is changed to 25m (in program line 1080). This only applies to the distance used in calculation of the drawdown at a node due to the discharge from the same node, drawdowns due to discharges at other nodes will be dependent on the distances calculated. The values for change in piezometric head for the example data are given below. Table 8.3; -3·21 -3.92 -4.14 -5.41 -5.81 -5.61 -5.01 -4.29 -3.64 -3·09 C0

-3.28 -4.21 -5.48 -6.90 -1.53 -7.09 -5.85 -4.71 -3.86 -3.22 C 1

Calculated heads. -3.11 -4.13 -5.19 -8.83 -9.19 -9.11 -6.31 -4.82 -3.89 -3.23 C2

-2.62 -3.51 -4.82 -6.32 -1.05 -6.72 -5.56 -4.49 -3.69 -3·09 C 3

-1.83 -2.41 -3.31 -4.30 -4.86 -4.88 -4.44 -3.86 -3·31 -2.84 C4

-0.18 -1.10 -1.84 -2.68 -3·29 -3.52 -3.44 -3.18 -2.86 -2.53 C5

0.53 0.82 -0.21 -1·35 -2.11 -2.50 -2.62 -2.56 -2.40 -2.19 C 6

1.62 4.08 1.08 -0.45 -1.28 -1.75 -1.98 -2.03 -1.98 -1.86 C7

R9 R8 R1 R6 R5 R4 R3 R2 R 1 R0

Matrix model

8.5

If you set up the model for more than 8 or 9 columns, the printout to the screen and/or the printer will become less tidy because of insufficient space on one row for all the columns of data. Little can be done about this problem on the screen, the limit of 80 character columns is fixed. With the printer, however, there is usually an answer. Most readers will be using dot matrix printers which can usually be set up to print a character every 1.45mm instead of the default of a character every 2.54mm. This will give you about 132 characters on a 'standard' 210mm wide page, and much more on a wide printer. Unfortunately each printer has it's own set of commands, and you will have to check your manual for the command to set your printer to condensed print. A second problem may arise if you attempt to send more than 80 characters from your computer to your printer without an instruction to change to the next line. Several computers I've used (including an IBM PC clone) count the characters sent to the printer, and if the count gets to 80 before a command to roll on to a new line (line feed, or just LF) is sent, the computer will take matters into it's own hands and send a line feed (with usually a carriage return as well. ie. an instruction to go back to the start of the line. Abbreviated to CR.) This problem is fixed by an instruction to the computer to send longer lines, usually with the WIDTH command. As the way the command is used will depend on your computer, you will need to refer to your manuals. 2.2. The method used by the program Initialization and the loading of constants are done first, similarly to many other of the programs in this book. Then the first part of interactive input, up to and including the entry of the time unit for discharges. Line 60 dimensions arrays before the subroutine to set-up those arrays is called. This subroutine uses variables LX, LY, UX, UY, and V to handle the values entered for lower x, lower y, upper x, upper y, and value, respectively. In the case of entry of a single point, or a single column, UX is made to equal LX, and similarly in the case of a single row, UY is made to equal LY. A flag, PA, keeps track of the parameter that is being entered. Lines 7216 to 7228 enter the changes into the matrices. The core of the program is in the double double loop in lines 1000 to 1240. The outer double loop sets the element for which drawdown is to be calculated, and the inner one sums drawdowns at the current node due to all discharging nodes. Line 1060 tests each element in the inner double loop for a discharge value, if it has no discharge the next element is tested, and so on until a discharging node is found. Line 1070 uses Pythagoras's theorem to calculate the distance between the current element for which the drawdown is required, and the current discharge element (or more correctly, between the two nodes). As mentioned above, line 1080 tests for a distance of zero, and if found, alters it. Lines 1090 to 1130 use the polynomial approximations to the Theis equation to calculate drawdown due to the current discharge element, and then line 1140 sums that discharge with any calculated for the current element previously. Back in the outer double loop, line 1160 sums

8.6

Matrix model

the calculated stage.

total

drawdown with

that

entered by the user at the setup

2.3. Usage of variables in ANMODL used in the polynomial approximation of the Well CO - cg Constants Equation. Total of all calculated drawdowns for the current node. rn The storage coefficient. CS The drawdown calculated for the current node, due to the current DD discharge/recharge well. A string for numerical formatting. G1 H(x,y) The matrix of piezometric heads. The primary (integer) counter. I The second level primary counter. 10 The secondary (integer) counter. J The second level secondary counter. JO A flag which is set when the line printer is to be used. LP Lower x value. LX Lower y value. LY The number of columns in the matrices (and also the model). NC The number of rows in the matrices. NR Parameter currently being entered, or printed. PA Has the value 4 times Pi times transmissivity. PT Distance from the current discharge well to the current node being R calculated. Note that the R printed at the end of each row of output is an indicator only, and not a variable used in the program. A string used for handling yes/no questions. SY Transmissivity. T Duration of current discharge/recharge. T1 Time factor, used to convert times to days. TF TI(x,y) The matrix of durations of discharges. Q(x,y) The matrix of discharge rates. U The variable u, as used in the Theis Equation. U2 Has the value u squared. UX Upper x value. Upper y value. UY Value to be entered into current matrix (as indicated by PAl at the V nodes indicated by LX, LY, UX, and UY. The distance between nodes in the x direction. X1 The distance between nodes in the y direction. Y1 ZX, Z$, and ZZ$ are used in program flow control. REFERENCE Kinzelbach, W., 1986. Programs in Basic.

Groundwater Modelling, An Introduction with Sample Elsevier. Amsterdam-Oxford-New York-Tokyo. 333PP.

Matrix model 3. 3.1.

PROGRAM LISTING. An analytical matrix model; program ANMODL

1 REM ••••• ANMODL/BAS ••••• 2 REM • Last modified 18th. May, 1986 5 CLS: PRINT TAB( 10) i" ANMODL/BAS" 10 DEFINT I-L:DEFSTR G,S 20 PRINT"Transient state. Homogeneous, isotropic aquifer." 30 GOSUB 8500:REM U Load constants 50 GOSUB 7000:REM U First part interactive input 60 DIM H(NC,NR),Q(NC,NR),TI(NC,NR) 70 GOSUB 7064:REM U Second part interactive input 1000 REM === Calculation section === 1010 CLS 1020 PRINT"Calculated heads":IF LP=1 THEN LPRINT:LPRINT"Calculated heads" 1030 FOR J=NR TO 0 STEP-1:FOR 1=0 TO NC:REM # Element being calculated 1040 CD=O:DD=O 1050 FOR JO=O TO NR:FOR 10=0 TO NC:REM U Element causing head change 1060 IF Q(IO,JO)=O GOTO 1150 1070 R=SQR«IO-I).(IO-I).X1.X1+(JO-J)·(JO-J).Y1'Y1) 1080 IF R=O THEN R=X1/4 1090 U=R.R'CS/(4'T'TI(IO,JO)):IF U>80 THEN 1150 1100 U2=U'U:IF U>1 GOTO 1130 1110 DD=(Q(IO,JO)/PT)'(-LOG(U)+CO+C1'U+C2'U2+C3'U'U2+ C4'U2'U2+C5'U'U2'U2) 1120 GOTO 1140 1130 DD=(Q(IO,JO)/PT)'(1/(U'EXP(U)))'(C6+C7'U+U2)/(C8+C9'U+U2) 1140 CD=CD+DD 1150 NEXT IO:NEXT JO 1160 H(I,J)=H(I,J)-CD 1170 PRINT USING G1iH(I,J)i:IF LP=1 THEN LPRINT USING G1iH(I,J); 1180 NEXT I:PRINT" R"iJ:IF LP=1 THEN LPRINT" R"iJ 1190 NEXT J 1200 FOR 1=0 TO NC 1210 PRINT USING" CU# "iIi 1220 IF LP=1 THEN LPRINT USING" C## ";Ii 1230 NEXT I 1240 PRINT:GOTO 9999 7000 REM === Keyboard input === 7004 PRINT"Do you want to use the printer?" 7008 Z$=SY:GOSUB 9800:IF ZX<3 THEN LP=1 ELSE LP=O 7012 PRINT:PRINT"Size of array" 7016 INPUT"No. of nodes in x direction"iNC:NC=NC-1 7020 INPUT"No. of nodes in y direction";NR:NR=NR-1 7024 INPUT"Spacing of nodes in the x direction (m)";X1 7028 INPUT"Spacing of nodes in the y direction (m)"iY1

8.7

8.8

Matrix model

7032 7040 7044 7048 7052 7056 7060 7064 7068 7072 7076 7080 7084 7088 7092 7096 7100 7104 7108 7112 7116 7120 7124 7128 7132 7136 7140 7144 7148 7152 7156 7160 7164 7168 7172 7176 7180 7184 7188 7192 7196 7200 7204 7208 7212 7216 7220

INPUT"Transmissivity (m*m/Day)";T INPUT"Storage coefficient";CS PRINT"What time unit (for duration of discharge/recharge) ." PRINT"inutes or ays." Z$="MmDd":GOSUB 9800 IF ZX<3 THEN TF=1440 ELSE TF=1 RETURN REM -- Second part of keyboard input CLS PRINT"Note that the first column and row are numbered zero." PRINT"== Setting up ==" PRINT" Menu 1" PRINT"To enter" PRINT"Initial head, press 1" PRINT"Discharge/Recharge, 2" PRINT"or i f finished; 3" Z$="123":GOSUB 9800:IF ZX=3 THEN RETURN PA=ZX PRINT:PRINT" Menu 2. Type of 'block'" PRINT"Single node, 1" PRINT" Single column, 2" PRINT"Single row, 3" PRINT"Block, 4" PRINT"Get print-out, 5" PRINT"or go to menu 1; 6" Z$="123456":GOSUB 9800 ON ZX GOTO 7144,7152,7164,7176,7250,7064 INPUT"X, Y co-ordinates, and value. ego '1,4,100"';LX,LY,V UX=LX:UY=LY:GOTO 7184 PRINT"X co-ordinate of column, Y co-ord of bottom and" INPUT"top, and value. ego '3,2,5,100"';LX,LY,UY,V UX=LX:GOTO 7184 PRINT"Y co-ordinate of row, X co-ord. of left & right" INPUT"end, and value. ego '3,2,5,100"';LY,LX,UX,V UY=LY:GOTO 7184 PRINT"Left X, right X, bottom Y, top Yeo-ordinates," INPUT"and value. ego '1,4,2,6, 100"'iLX,UX,LY,UY,V IF LXNC THEN PRINT"Upper column No. to greatl":GOTO 7108 IF UY>NR THEN PRINT"Upper row No. to greatl":GOTO 7108 IF PA=2 THEN INPUT"Duration of discharge/recharge"iT1 ON PA GOTO 7216,7220 FOR J=LY TO UY:FOR I=LX TO UX:H(I,J)=V:NEXT I:NEXT J:GOTO 7108 FOR J=LY TO UY:FOR I=LX TO UX

Matrix model 7224 7228 7250 7260 7270 7280 7290 7300 7310 7320 7330 7340 7350 7360 7370 7380 7390 7400 7410 7420 7430 7440 7450 7460 7470 8500 8510 8520 8530 8540 8550 8560 9800 9810 9820 9830 9840 9850 9999

Q(I,J)=V:TI(I,J)=T1/TF NEXT I:NEXT J:GOTO 7108 REM === Video and printer output CLS ON PA GOTO 7280,7340 PRINT"Heads":IF LP=1 THEN LPRINT:LPRINT"Heads" FOR J=NR TO 0 STEP-1:FOR I=O TO NC PRINT USING G1jH(I,J)j:IF LP=1 THEN LPRINT USING G1jH(I,J)j NEXT I PRINT" R"jJ:IF LP=1 THEN LPRINT" R"jJ NEXT J:GOTO 7470 PRINT"Discharge, (metres cubed/day)":IF LP=1 THEN LPRINT: LPRINT"Discharge, (m cubed/day)" FOR J=NR TO 0 STEP-1:FOR I=O TO NC PRINT USING G1jQ(I,J)j:IF LP=1 THEN LPRINT USING G1jQ(I,J)j NEXT I PRINT" R"jJ:IF LP=1 THEN LPRINT" R"jJ NEXT J PRINT SCj:Z$="Cc":GOSUB 9800 CLS:PRINT"Duration of operation of discharge/recharge"j:IF LP=1 THEN LPRINT:LPRINT"Duration of operation of discharge/recharge."j IF TF=1 THEN PRINT" (Days)":IF LP=1 THEN LPRINT" (Days)" IF TF=1440 THEN PRINT" (Minutes)":IF LP=1 THEN LPRINT" (Minutes)" FOR J=NR TO 0 STEP-1:FOR I=O TO NC:PRINT USING G1jTI(I,J)/TFj :IF LP=1 THEN LPRINT USING G1jTI(I,J)/TFj NEXT I:PRINT" R"jJ:IF LP=1 THEN LPRINT" R"jJ NEXT J:GOTO 7470 PRINT SCj:Z$="Cc":GOSUB 9800:PRINT:PRINT:GOTO 7076 REM === Input constants === PT=4*3.141593*T:SY="YyNn" G1="ifIIHIIiI. 1111": SC="Press C to continue " CO=-.57721566H:C1=.99999193H:C2=-.24991055H:C3=5.51996 8E- 02 C4=-9.76004E-03:C5=1.07857E-03:C6=.250621:C7=2.334733 C8=1.681534:C9=3.330657 RETURN REM === Await valid key --ZZ$=INKEY$ ZZ$=INKEY$:IF ZZ$="" THEN 9820 FOR ZX=1 TO LEN(Z$):IF MID$(Z$,ZX,1)=ZZ$ THEN RETURN NEXT ZX GOTO 9820 END

8.9