PAGE 1
GROUNDWATER FLOW TO DRAINS: SEMIANALYTICAL MODEL FOR RESIDENCE TIME by Thomas A. Farkas A thesis submitted in partial fulfillment of the requirements for the degree of Master of Science in the Department of Geology in the University of South Florida May 1992 Major Professor: Len Vacher, Ph.D.
PAGE 2
Graduate Council University of South Florida Tampa, Florida CERIFICATE OF APPROVAL MASTER'S THESIS This is to certify that the Master's Thesis of Thomas Alex Farkas with a major in the Department of Geology has been approved by the Examining Committee on January 24, 1992 as satisfactory for the Thesis requirement for the Master of Science degree Thesis Committee: Major Professor: Len Vacher, Ph.D. Mark Stewart, Ph. D. Ph.D.
PAGE 3
ACKNOWLEDGEMENTS This project was made possible by my wife, Renee, whose confidence in me, gave me the strength to finish after all these years. Much appreciated are the encouragements of my parents, Erika and Alex Farkas, who instilled in me the benefits of a good education. Drs. Len Vacher, Mark Stewart, and Sam Upchurch are thanked for their guidance as well as their editorial and technical comments. The contributions of Ken Badgely deserve special mention as he lent invaluable CAD assistance on the figures. Finally, everyone at Schreuder & Davis, Inc is thanked for their encouragements and expertise which contributed to this project. ii
PAGE 4
TABLE OF CONTENTS LIST OF TABLES v LIST OF FIGURES vi ABSTRACT viii INTRODUCTION 1 THEORY: SOLUTIONS FOR V AND . . . . . 3 Potential Theory Solution for Flow to TileDrain . . . . 3 Potential Theory Solution for Flow to DitchDrain . . . 6 DupuitForschheimer Theory Solution for Flow to DitchDrain . . . . 10 Summary of the Three Problems . 13 THEORY: RESIDENCE TIME ...... Average Residence Time . . Streamline Residence Time . METHODS Tileand DitchDrain Problems . . DupuitForschheimer DitchDrain Problem RESULTS . . . . . . cauchyReimann Analysis . TileDrain Problem DitchDrain Problem DupuitForschheimer DitchDrain Problem DISCUSSION . . . . . CauchyReimann Analysis . TileDrain Problem . . . . DitchDrain Problem . . . . . DupuitForschheimer DitchDrain Problem . Comparison of Tileand DitchDrain Problems Comparison of Tileand DF DitchDrain Problems . . . . . . CONCLUSIONS iii 15 15 16 19 19 20 22 22 24 24 30 36 36 38 40 41 43 44 46
PAGE 5
REFERENCES 48 APPENDICES 50 Appendix A 51 Appendix B 5 5 Appendix c 63 iv
PAGE 6
LIST OF TABLES Table 1 Streamline residence times . v Page 37
PAGE 7
LIST OF FIGURES Figure 1 2 3 4 Flow to a tiledrain (after Kirkham & Powers, 1972) . . . . . . . . . . . Flow to a ditchdrain (after Kirkham, 1958) DupuitForschheimer flow to a ditchdrain. (A) Equal ditch water elevation (after Kirkham, 1967), (B) Halfsymmetrical section ...... Distribution of residence time within flow problem (after Vacher et al, 1990) .... 5 Graphs showing the distribution of oCR. (A) Standard tiledrain, (B) Nonpenetrating ditchdrain, (C) Fully penetrating ditchdrain, (D) Page 4 7 11 18 Tiledrain with depth equal to 200 m . . . 23 6 Tiledrain problems with D/L=0.04. (A) Standard tiledrain with D=20 m, L=500 m, R=O.l mjyr, n=0.1, and K=100 m jyr, (B) R=0.2 m jyr, (C) n=0.2, (D) K=200 m jyr . . . . . . . 25 7 Tiledrain problems with variable D and L. (A) D=40 m and L=500 m (D /L=0.08), (B) D=80 m and L=500 m (D/L=0.16), (C) D=200m and L=500m (D/L=0.4), (D) D=20 m and L=1000 m (D /L=0.02), (E) D=4 0 m and L=1000 m (D/L=O. 04) . 2627 8 Distribution of dimensionless residence time vs. D / L ratio . . . . . . . 28 9 Ditchdrain problems with standard parameters and D/L=0.04. (A) Nonpenetrating ditchdrain d=0.1 m, (B) Onequarter penetrating ditchdrain d=5 m, (C) Onehalf penetrating ditchdrain d=10 m, (D) Fully penetrating ditchdrain d=20 m . 29 vi
PAGE 8
10 DupuitForschheimer ditchdrain problems. (A) Problem with standard parameters (K=100 mjyr), (B) K=200 mjyr, (C) K=500 mjyr, (D) K=1000 mjyr, (E) K=2000 mjyr, (F) K=4000 mjyr, (G) K=SOOO m/yr . . . . . . . . . . . 3132 11 Comparison of the 0.80 V0 streamline residence time in the tiledrain and DupuitForschheimer ditchdrain problems; 1"1=DF, 1"2=Tile . 33 12 DupuitForschheimer ditchdrain problems with R=0.2 mjyr. (A) K=100 mjyr, (B) K=SOOO mjyr 34 13 DupuitForschheimer ditchdrain problems with L=1000 m. (A) K=100 mjyr, (B) K=SOOO mjyr 35 vii
PAGE 9
GROUNDWATER FLOW TO DRAINS: SEMIANALYTICAL MODEL FOR RESIDENCE TIME by Thomas A. Farkas An Abstract A thesis submitted in partial fulfillment of the requirements for the degree of Master of Science in the Department of Geology in the University of South Florida May 1992 Major Professor: Len Vacher, Ph.D. viii
PAGE 10
The residence time of groundwater, where groundwater flows from an upland to a drain, can be determined from three different analytical models: potentialtheory solution for flow to a tiledrain (Kirkham and Powers, 1972), potentialtheory solution for flow to a ditchdrain (Kirkham, 1958), and DupuitForschheimer flow to a ditchdrain (Kirkham, 1967). A CauchyRiemann analysis of the tileand ditchdrain models shows that neither model's analytical solutions exactly satisfy Laplace's equation. The errors in the and solutions, however, are larger in the ditchdrain model. These errors are largest near the drain and prevent the precise quantification of the relationship between drain penetration and streamline residence time. Comparison with the tiledrain model reveals that changes in streamline residence time with increasing drain penetration are relatively minor. The residence times for the entire suite of various penetrating ditchdrain problems may be calculated from the simpler but more accurate tiledrain potentialtheory solution. As long as the geometry of the tiledrain problem remains uniform, the residence time is linearly affected by the porosity, recharge and depth of the flow problem. If the ix
PAGE 11
geometry of the problem changes as a result of altering the depth/length ratio, both the streamline residence time and the position of the streamline are affected. The analysis of the DupuitForschheimer ditchdrain model reveals the model behaves differently depending upon the presence or absence of a watertable mound. In the presence of the watertable mound, which is controlled by the recharge/hydraulic conductivity ratio, the residence time is nonlinearly affected by recharge, hydraulic conductivity and length of the flow problem. If the watertable mound is absent, the residence time is linearly affected by porosity, recharge and the height of the problem. The residence times calculated under conditions where the watertable mound is negligible, very closely approximate those calculated in the tiledrain model. ABSTRACT APPROVED: Major Professor: Len Vacher, Ph.D Professor, Department of Geology Date Approved X
PAGE 12
INTRODUCTION Time nets showing streamlines and contours of groundwater age have been produced for idealized, steadystate flow systems such as the coastal wedge (Vacher et al., 1991), the strip island lens (Vacher et al., 1990) and the Toth regional flow system (Robinson, 1991). The intent of all of these studies has been to show the relationship between timerelated parameters (i.e., residence time, velocity, ground water ages) and underlying controls such as geometry and geology of the flow system, its dimensions, the distribution of recharge, and general hydrologic setting. The purpose of this research is twofold. First is to document relationships between timerelated parameters and underlying controls for each of three, twodimensional analytical models: (1) the potentialtheory solution for flow to a tiledrain (Kirkham and Powers, 1972), (2) the potentialtheory solution for flow to a ditchdrain (Kirkham, 1958) and (3) DupuitForschheimer flow to a ditchdrain (Kirkham, 1967). The second purpose is to determine the advantages andjor disadvantages of using any one of the three analytical models to calculate travel times for groundwater flow to draintype outlets. Each of these analytical models treat groundwater flow from an upland to
PAGE 13
2 an outlet such as a drain. Knowledge about groundwater residence times is important in the study of contaminated groundwater and plume migration. A better understanding of the hydrogeologic parameters affecting various flow systems and the analytical models representing these systems will help increase the ability to predict contaminant plume migration rates.
PAGE 14
THEORY: SOLUTIONS FOR V AND Potential Theory Solution for Flow to TileDrain The conceptual model and nomenclature for the tiledrain problem are shown in Figure 1 The origin is at the drain. The half width of the mound (distance from drain to the groundwater divide), the depth of the impermeable barrier below the drain, the width of the drain, and the recharge rate are denoted by L, D, A and R, respectively. The flow domain is the rectangle bounded by the x and y axes, the divide, and the impermeable barrier. The flow in the mound above the x axis is ignored to simplify the problem, so the solution is approximate. The solution becomes exact, however, as the height of the mound (HL) becomes small relative to D. The problem is solved by Kirkham and Powers (1972) in terms of the stream function (V, L3/ T /L). The governing equation is Laplace's flow.equation: 3 (1) The boundary conditions are: BC1 V=V 0 where x=O and O $ y $ D
PAGE 15
RAINFALL R l l l l l y Figure 1. Flow to a tiledrain (after Kirkham & Powers, 1972) 4 I' HL X
PAGE 16
BC2 lf=lf 0 where O:$x:$L and y=D BC3 lf=lf 0 where x=L and O:$y:$D BC4 1f={{xa)/(la)}lf0 where a:$x:$L and y=O BC5 lf={ (ax)ja}lf0 where O:$x:$a and y=O 1f'0 is the total flow into a unit length of the sink (drain) and is given by: 5 1f'0=R(La) {2) The solution for 1f is: 1f=lf0 21f'/7r L 1/m sin{m7rX/L) sinh[m7r(Dy)/L]/sinh(m7rD/L) (3) where limits of l: is between 1 and Kirkham and Powers {1972, pg. 115) make approximations in order to derive a solution of (3). These approximations include an impermeable barrier to prevent deep seepage, neglecting the loss of head in the watertable mound, and utilizing a drain as a strip sink of zero thickness. From equation (3) and the CauchyRiemann conditions, the solution for the specific discharge potential which is K, where K is hydraulic conductivity and is the hydraulic head) is derived (Kirkham and Powers, 1972). The hydraulic head is = 2RL/K7r { (1/2) ln[ert/L [cosh(7ry/L) cos(7rX/L) ]/ (4)
PAGE 17
6 + L 1/m ( eIIIWO/L) /sinh } where limits of r is between 1 and Potential Theory Solution for Flow to DitchDrain The conceptual model for the ditchdrain problem flow system is shown in Figure 2. The nomenclature is similar to the tiledrain system except d is introduced as the depth of penetration of the ditch, and h is the depth from the water level of the ditch to the impermeable barrier. The analytical equations needed to solve for V and for three cases of varying ditch penetration are developed by Kirkham (1958). The case of the fully penetrating ditch, where d=h, and the nonpenetrating ditch, where d approaches o, are the two end members of the problem. The solution for V and for these two cases is presented in Appendix A. The intermediate case in the ditchdrain problem is where the ditch (O
PAGE 18
7 RAINFAL L R l 1 1 l 1 r H L X 4 y Figure 2. Flow to a ditchdrain (after Kirkham, 1958)
PAGE 19
For this problem, V=RL. With these boundary conditions, the solution to Laplace's equation with some rearrangement is (Kirkham, 1958) where S1= tan[sin (7rX/L) 1 { cos (7rX/L) } ] S2= L ( 1/m) sin (mnx/L) emm/L sinh (mny /L) I sinh (mnh/L) S3= tancos(nyjh) }] S4= L ( 1/m) sin (mny /h) emwL/h sinh (mnxjh) I sinh (mnL/h) (5) S5= (hjnd) L (1/m"2) sin(mnd/h) sin(mnyjh) emwx/h S6= (h/7rd) L (1/m"2) sin(mnd/h) sin(mnyjh) e mwL/h sinh(mnxjh)/sinh(mnL/h) The limits of Lis between 1 From equation (5), the explicit solution for V is 8 V= 2Vo/7r {SlS2+S3S4S5+S6} (6) From equation (6), the CauchyRiemann conditions, and some rearrangement, Kirkham (1958) developed the following equation for = T1T2+T3+T4T5T6 (7)
PAGE 20
where: T1= (1/2) ln[2ewy/L (cosh('Jt'y/L) COS('Jt"X/L))] T2= }; (1/m) cos(m'Jt"X/L) eBh/L cosh(m'Jt"y/L)/sinh(m'Jt"h/L) T3= (1/2) ln[2en/h (cosh('Jt"X/h) cos('Jt"Y/h))] T4= }; (1/m) cos(m'Jt"yjh) eLih cosh(m'Jt"X/h)/sinh(m'Jt"L/h) TS= (h/'Jt"d) E(l/mA2) sin(m'Jt"d/h) cos(m'Jt"Y/h) T6= {h/,.d) r ( 1/mA 2) sin (m'Jt"d/h) cos (m'Jt"y /h) e L/h cosh(m,.xjh)/sinh(m,.L/h) 9 The limits of }; is between 1 and Because I= equation (7) can be rearranged to give an explicit solution for as: 2'ofK,. (T1T2+T3+T4TST6) +C/k (8) The solution for (equation 8) contains an unevaluated constant (C), chosen such that for the lowest point on the water table (Kirkham, 1958). To calculate C, therefore, a point (x,y)= (r,O) at the ditch drain is considered. The drain is assumed to be full to the level y=O and r is the semiwidth of the drain. Equation (8) may be rewritten to show C= (21lof"') V(r,O) {9) where V(r,o)= o and V(x,y) is equivalent to T1T2+T3+T4T5T6 (Kirkham, written communication 1989). With these stipulations, the solution for C in the partially
PAGE 21
10 penetrating ditch is: C/K = { r(1/m) (10) [ )1] ln(2e .. r/Zh + r (1/m) eL/h h/7rd r (1/m"2) sin er/h h/1fd r ( 1/m"2) sin eIIIII'L/h where limits of r is between 1 and Kirkham (1958, pg. 896897) rearranges the solutions to and in order to obtain more rapidly converging infinite series for values of x and y near the sinks. The result is the ability to calculate near the ditchdrain using a simpler solution. DupuitForschheimer Solution for Flow to DitchDrain The conceptual model for the DupuitForschheimer (DF) ditchdrain problem with discharge to fully penetrating ditches is shown in Figure 3a. The heights of water in the ditches above the impermeable barrier are h1 and h2 respectively. The distance between the two ditches is Recharge (R) is a constant rate. The elevation of the water table at the groundwater divide is h(D). For the case where h1=h2, the groundwater divide, which is a noflow boundary, occurs at the midpoint between the two rivers. The flow systemis symmetrical about the divide, so only one side of
PAGE 22
11 RAINFALL R (A) h ( d ) T 1 h 1 h2 l l&.1l! l RAINFALL R (B) l w h1 h (X) h (D) 1.()0+. Figure 3 DupuitForschheimer flow to a ditchdrain. (A) Equal ditch water elevation (after Kirkham, 1967), (B) Halfsymmetrical section.
PAGE 23
12 the flow system needs to be considered (Figure 3b) The origin of the x,w coordinate system is at the lower left corner of the ditch with x increasing to the right and w upward. Along the watertable surface, w is equal to the height of the water table, h(x). Fetter (1980) shows the derivation of the hydraulic head based on the DF assumptions and the boundary conditions that h=hl at x=O and h=h2 at x=2L. The DupitForschheimer assumptions are: 1) The hydraulic gradient is equal to the slope of the water table. 2) Flow is evenly distributed with depth (i.e., streamtubes at any one point are all the same thickness). The solution for the hydraulic head is: h2 = h12 {(hl2 h22)/2L}X + (R/K) (2Lx)x (11) For the case where the height of the water in both ditches is the same, the solution reduces to: h= {h12 + (R/K) (2LX)X}A0.5 (12)
PAGE 24
13 Kirkham (1967) uses the DF assumptions to develop an approximate analytical solution for the stream function. The solution for the stream function with some rearranqement is: != f0 (W/h) R(Lx) (13) where !0 is equal to RL. Summary of the Three Problems The three problems considered in this thesis represent idealized, twodimensional, steadystate flow systems which are drained by varyinqly penetrating drains. Recharge is uniform, and hydraulic conductivity is homoqeneous and isotropic. The drain water levels, regardless of drain penetration, are at an equal elevation above the impermeable barrier. Thus flow is symmetric about the groundwater divide. The tiledrainage and ditchdrainage potentialtheory problems assume that the volume of the watertable mound is negligible in comparison to the volume of the rectanglar area below the mound. This assumption, together with the boundary conditions following from the stipulated uniform distribution of recharge, enables solution of the Laplace equation in terms of the stream function. The solution for potential then follows by application of the CauchyRiemann conditions.
PAGE 25
14 The nonpenetrating ditchdrain and the tiledrain are geometrically identical. The solution for f and for the nonpenetrating ditch, however, is mathematically more complicated than for the tile drain. This results from the use of a horizontal stripdrain along the x axis in the tiledrain problem. In contrast, the solution for the ditchdrain problems use a vertical strip drain along the y axis. The solutions to and become increasingly more complicated, progressing from the nonpenetrating ditch to fully penetrating ditch to the partially penetrating ditch drain problems. The DupuitForschheimer ditchdrain flow system does not assume that the volume of the watertable mound is negligible. The volume of the mound increases as R/K increases. Thus as R/K decreases, the geometry of the flow system becomes increasingly similar to that of the ditchdrain problem (Kirkham, 1958). In spite of this geometric similarity, the solution to and for the DupuitForschheimer ditch drain is considerably less complicated than in Kirkham's potential theory solutions because of mathematical simplifications allowed by the DupuitForschheimer assumptions.
PAGE 26
15 THEORY: RESIDENCE TIME Average Residence Time In any steadystate reservoir flow system, total inflows to the reservoir equal total outflows. The average residence time of the ground water within the reservoir is a lumped parameter equal to the volume of ground water in storage divided by the recharge to the flow system (Vacher et al., 1990). For the tiledrainage problem (Figure 1), 0= nDL/RL = nD/R (14) where n is porosity and D, the depth to the impermeable barrier; represents the flow system thickness. This equation also applies to the ditchdrainage problem (Figure 2). The depth of ditch penetration doesnot alter the geometry of the flow system, and therefore, does not affect Hydraulic conductivity is from the solution because it does not alter the volume of groundwater in the flow system and waterparticle transport is dependent on recharge rate. The solution (Vacher, written. communication 1989) to in the DF ditchdrain problem is derived from the equation:
PAGE 27
16 T0= nh/RL (15) where his given by equation (12). Equation (15) is integrated between the limits of the problem (0 to L): T0 = {nf h1dx + nf (R/K) (Lxx2)0 5 dX}/(RL) The solution to T0 is: .,. 0= nh1/R + mrL/8 (RK) 0 5 (16) (17) The first and second terms of equation (17) represent the residence time calculated for the volume of groundwater in the rectangle and in the watertable mound respectively. The hydraulic conductivity is present in the second term of the solution because it affects the volume of the watertable mound. As the volume of the mound becomes negligible compared to the rectangle, the equation for the T0 reduces to the first term, which is equivalent to the T0 calculated for the tiledrain and ditchdrain problems. Streamline Residence Time Each individual water particle within a steadystate flow system has its own residence time, which is the time for the particle to travel the length of a streamline. Residence time, therefore, is also a distributed parameter
PAGE 28
that depends on the length and distribution of potential along the various streamlines (Vacher et al., 1990). As 17 shown in Figure 4, the residence time (ri) in each streamtube is the volume of ground water in the tube dlvided by the discharge through the tube per unit length of the ditch, or (18) where m is equal to the number of streamtubes. With an infinite number of streamtubes, the overall residence time (r0 ) becomes equal to the average of all the streamline residence times (Vacher et al., 1990); for example (Figure 4): 5 5 r f = nr Vi/ (RL/5) = 5nV /RL = 5r 0 (19) i=1 1 i=1 f = 0 r r;J5 (20)
PAGE 29
RL/5 RL/ 5 l l d=h RL/ 5 l vs RL/5 l 'Ti= nVJ(RL/m) 'T0= Tjm 18 Figure 4. Distribution of residence time within flow problem (after Vacher et al, 1990).
PAGE 30
19 METHODS Tileand DitchDrain Problems A tracking routine was developed to calculate streamlineresidence times. The tracking begins at the point where the particle first enters the flow system along the y axis. The Y of the streamline is calculated using equation BC 4 and BC 5 for the tileand ditchdrain problems respectively. Then with x incremented by toward the drain and y taking on successively larger values by the equation for the appropriate Y is solved for each trial combination of x and y until the calculated Y exceeds the Y for the streamline. This locates a point on the streamline. From the x and y coordinates of the newly determined point, the potential is calculated using the appropriate solution to The gradient is then calculated between the new point and the previous point on the streamline. With the solution to the gradient in hand, the flow velocity (v) between the two points is calculated using Darcy's Law : (21) where 0 and are the potentials at the new and previous pointson the streamline and d(x,y) is the Pythagorean
PAGE 31
distance between those two points. A travel time (rt) between the two points is then calculated using the equation: 20 (22) The travel times between successive points on the streamline are cumulated, in the tiledrain problem, to the point where the streamline discharges into the drain (x=0.1 m) and in the ditchdrain problem to the point x= 10 m from the drain. The sum of the travel times at these points is equal to the The and increments are 1 m and 0.1 m respectively. Computer runs with a increment of 0.1 m changed the solution for f; by less than 0.1 percent at the expense of increasing computer run time by a factor of 7. The computer programs for these calculations are listed in Appendix B. DupuitForschheimer DitchDrain Problem The calculation of r ; for the DF problem is similar to that for the two potentialtheory solutions. The tracking procedure is considerably easier because of the simpler mathematics of the analytical expressions for V and The water particle enters the flow system at the water table where V, using equation (13), is calculated for the
PAGE 32
21 streamline. The height of the water table at the entry point is given by equation (12). A new solution for the head (hn) is obtained by incrementing towards the ditch. The next point on the streamline is thus calculated by rearranging the stream function equation to solve for the exact position of coordinate w: w= (V0 V)hn/R(Lx) (23) where V is determined at the entry point for the streamline. With the new point on the streamline known, the procedure to calculate the hydraulic gradient, flow velocity and travel time between the new and previous points on the streamline is the same as in the tileand ditchdrain problems. The streamline is tracked to the point of discharge into the ditch (x=O m). The sum of the travel times gives the ri. The increment is one meter. An exact solution for the coordinate w is calculated.The computer program for these calculations is listed in Appendix c.
PAGE 33
22 RESULTS CauchyRiemann Analysis The mathematical imprecision of any analytical solution to Laplace's equation can be quantified in terms of the CauchyRiemann (CR) conditions, (24) These conditions hold if and only if V and both satisfy Laplace's equation. Therefore one can test the solution by calculating the parameter (25) at each point. If the analytical solution to and V exactly satisfies Laplace's equation, 6CR is 0. The magnitude of the deviation from CauchyRiemann conditions is represented by increasing positive values of 6CR. The graphs showing deviation from CauchyRiemann conditions are shown in Figures 5AD. Figure 5A is a graph showing the distribution of 6CR for the standard tiledrain problem. The graphs in Figures 5B and 5C show the 6CR distribution for the fully penetrating and non
PAGE 34
(A) Cll a: ... ... ... ::1 E ::r:: ... G.. ... 0 (B) "' a: ... ... ... ::I ,._ c :z: ... G.. ... c (C) lQ ... ... ... ::I E :z: ... G.. ... c (D) "' "' ... ... ... ::I ,._ c :z: ... a. ... 0 0 10 15 20 010 100 150 200 250 300 350 oo so 500 0 5 10 1 5 20 0 5 10 0 50 200 DISTANCE (X) 1otTERS s a 5 2 7.0 s o 2.1 .. J 010 50 100 150 200 250 300 350 400 so 500 DISTANCE (X) lotETERS S.l 5.2 7.0 3.0 2. 6 .3 01 o so 1 oo 150 200 2so 3oo 350 oo so soo DISTANCE (X) lotETERS ../ O.J O J o s o. o s O.J o.s 0.1 o.z O.J o s 010 50 100 150 200 250 300 350 400 450 500 DISTANCE {X) lotETERS Figure 5. Graphs showing the distribution of 6CR. (A) Standard tiledrain, (B) Nonpenetrating ditchdrain, (C) Fully penetrating ditchdrain, (D) Tiledrain with depth equal to 200 m. 23
PAGE 35
penetrating ditchdrain problems with standard parameters. Figure 50 shows the distribution of 6CR in the tile drain problem with 0=200 m. TileDrain Problem 24 The distribution of streamline residence times for a variety of tiledrain problems is shown in Figures 6 and 7. For the purpose of comparison, the system of Figure 6A is taken as a standard: L=500 m, 0=20 m, K=lOO mjyr, R=O.l mjyr, n=O.l The systems in Figures 6BD have the same depthlength (D/L) ratio, but the value of one of the controlling parameters has been doubled. The systems in Figures 7AE have different combinations of depth and length with the same n, K, and R. Figure 8 is a dimensionless graph which shows the streamline residence timejaverage residence time ratio for systems with varying D/L ratios. The Ti/T0 ratio is shown for the streamlines for which! is 0.20, 0.40, 0.60, 0.70, 0.80, and 0.90 of 'o DitchDrain Problem The distribution of streamline residence time in varyingly penetrating ditchdrain problems is shown in Figures 9AD. The D/L ratio is 0.04 and the controlling parameters in each figure are the same as the standard of comparison tiledrain problem. Figure 9A shows a non
PAGE 36
25 0 (A) C/1 5 "' .., :::1 E 1 0 "' .... 15 ... Q T c 20 0 20 0 100 200 300 soo DISTANCE {X) WtT R S 0 (B) C/1 5 "' .., ... :::1 E 10 = .... 1 5 .., Q T, =10.0 2 0 0 100 200 .300 5 0 0 DISTANCE {X) WETERS 0 (C) C/1 5 "' .., ... .., :::1 E 10 = .... 15 .., Q T, 20 0 100 200 300 SOD DISTANCE {X) WETERS 0 (D) C/1 5 "' .., .., :::1 E 10 = 0.. 15 .., Q T =20.0 20 0 100 200 .300 soo DISTANCE {X) WETERS Figure 6 T iledrain problem s with D/L=0.04. (A) Standard tiledrain with 0=20 m L=SOO m, R=O. l mfyr, n=O.l, and K=lOO mfyr, (B) R =0.2 mfyr, (C) n =0.2, (D) K =200 mfyr.
PAGE 37
(A) 0 Cll 10 a:: LoJ 1LoJ ::i ,..., 20 >:I: 1Q.. 30 ..... 0 To =40.0 40 I I 0 100 200 .300 400 500 (B) DISTANCE {X) METERS 0 Cll 20 a:: ..... .... ..... :::::E 40 >:I: .... Q.. 60 LoJ 0 To =80.0 80 0 100 200 300 400 500 DISTANCE {X} METERS (C) 0 Cll 50 a:: LoJ 1..... :::::E 100 >....... :I: .... Q.. 150 w 0 io=200.0 200 0 100 200 300 400 500 DISTANCE {X) METERS Figure 7. Tiledrain problems with variable 0 and L. (A) 0=40 m and L=SOO m (0/L=O.OS), (B) 0=80 m and L=SOO m (O/L=0.16), (C) 0=200m and L=SOOm (O/L=0.4), 26
PAGE 38
Figure 7 cont ... Tiledrain problems with variable D and L. (D) 0=20 m and L=lOOO m (D/L=O. 02), (E) 0=40 m and L=lOOO m (D/L=0.04). 27
PAGE 39
3.0 2.5 2 0 0.80't0 .._o 1.5 0.70't 0 1.0 0 .60't0 0.5 0 .40't0 .. 0.20'to :1 0.0 0 0 0 1 0 2 0.3 0.4 D/L Figure a Distribution of dimensionless residence time vs. D/L ratio. 28 0.5
PAGE 40
(A) (B) {C) (D) 0 5 ... ... ... :::1 10 !:. :z: ... ... ... c 20 010 0 "' 5 "' ... :a E :z: ... ... 15 ... c 20 0 "' 0:: 5 ... ... ... :::1 E 10 % ... ... 15 ... c 20 ::1 5 ... ... ... :::1 10 :z: ... ... ... c 0 10 010 T. 0 100 200 30 <100 500 DISTANCE (X) WTRS ,., :20 0 100 200 30 <100 500 OISTAHCC (X) WCTERS 100 200 30 400 500 DISTANC E (X) WCTEAS ,., .. 20. 0 010 100 200 <100 500 DISTANCE (X) WI:TEAS 29 Figure 9 Ditchdrain problems with standard parameters and D/L=0 .04. (A) Nonpenetrating ditchdrain d=O.l m, (B) onequarter penetrating ditchdrain d=S m, (C) onehalf penetrating ditchdrain d =lO m, {D) Fully penetrating ditchdrain d=20 m.
PAGE 41
30 penetrating ditchdrain (d=O.l m). Figure 9B shows onequarter ditchdrain penetration (d=5 m). Figures 9C and 90 show onehalf (d=lO m) and fullypenetrating (d=20 m) ditchdrains. The streamline residence times depicted in each of the four Figures correspond to the point at x=lO m along the streamline. This is in contrast to the streamline residence times depicted in the tiledrain problems which correspond to x=O.l m (discharge point at the drain). DupuitForschheimer DitchDrain Problem The distribution of streamline residence time for a variety of OF ditchdrain problems is shown in Figures 1013. The system in Figure lOA has controlling parameters equal to that of the standard of comparison tiledrain problem. The controlling parameter K is varied in Figures lOBG in order to show its effect on the system. Figure 11 shows the relationship of the residence time for water traveling the 0.80 V0 streamline. The comparison is between the standard tiledrain and the OF ditchdrain which has a changing watertable mound. Figures 12A and 12B show the effect of double recharge at K=100 and 8000 mjyr respectively. Figures 13A and 13B show the effect of double system length at K=100 and 8000 mjyr respectively.
PAGE 42
(A) 20 C/1 16 "' ... ,_ ... 12 :z ,_ 8 :z: ;:;; 4 :z: T.=26.2 0 0 100 200 30 400 500 DISTANCE (X) litTERS (B) 20 C/1 16 "' ... ,_ ... 12 :z .!. 8 ,_ :z: !2 ... 4 :z: r. =24.4 0 0 100 200 30 400 500 DISTANCE {X) METERS (C) 20 "' "' 16 w ,_ ... 12 :z 8 !r ;:;; 4 :z: r. =22.8 0 0 100 200 30 400 500 DISTANCE {X) METERS (D) 20 :Q 16 ... ,_ w 12 :z .!. 8 ,_ :z: ;:;; 4 :z: r. .0 0 10 100 200 30 400 500 DISTANCE (X) METERS Figure 10. DupuitForschheimer ditchdrain problems. (A) Problem with standard parameters (K=lOO mjyr), (B) K=200 mjyr, (C) K=500 mjyr, (D) K=lOOO mjyr, 31
PAGE 43
32 (E) 20 en 0:: 16 loJ ..... loJ 1 2 ::::E ........ ........ ..... :I: B C> w :I: 7;, =21." 0 0 100 200 30 "00 50 0 (F) DISTANCE {X) t.4ETERS 2 0 en 1 6 0:: loJ ..... loJ 12 ::::E ..... B :I: C> w :I: r. =21.0 0 0 10 0 2 0 0 30 "00 500 DISTANC E (X) t.4ETERS (G) 2 0 en 0:: 1 6 loJ ...... w 12 ::::E ........ !. ..... B :I: C> w :I: 7;,=20. 7 0 0 1 0 0 2 0 0 3 0 4 00 500 DISTANCE (X) t.4ETERS Figure 10 cont DupuitForschheimer ditchdrain problems. (E) K=2000 mfyr, (F) K=4000 mfyr, (G) K=SOOO mfyr.
PAGE 44
1 2 3 4 5 6 7 8 Figure 11. Comparison of the 0.80 streamline residence time in the tiledrain and DupuitForschheimer ditchdrain problems; T 1=DF, T 2=Tile. 33
PAGE 45
Figure 12. DupuitForschheimer ditchdrain problems with R=0.2 mjyr. (A) K=lOO mjyr, (B) K=8000 mjyr. 34
PAGE 46
Figure 13. DupuitForschheimer ditchdrain problems with L=lOOO m (A) K=lOO mjyr, (B) K=SOOO mjyr. 35
PAGE 47
DISCUSSION CauchyRiemann Analysis The occurrence of oCR greater than zero at all data points in Figures SAD, indicates that neither the tileor ditchdrain analytical solutions exactly satisfy Laplace's equation. One of the most striking features of Figures 6, 7 and 9 is the occurrence of waves in the streamlines of the tileand ditchdrain problems. Particle tracking at finer increments of x and y do not affect the waves. Comparison of the streamlines (Figures 6, 7 and 9) with the distribution of oCR (Figures SAC) shows that the waves result from the approximating assumption used in deriving the analytical solutions for and of these problems. Waves are absent from the streamlines of Figure 7C, in which the depth of the flow system is 200 m The CauchyRiemann graph (Figure 50) for this problem shows that oCR throughout the entire flow system is considerably less than the tiledrain problem with a depth of 20m (Figure SA) The waves are quite pronounced in the tiledrain problem which has a length of 1000 m (Figure 70). This analysis appears to indicate that the stability of the and solutions is related to the D/L ratio of the problem. The solution is more stable with an increasing D/L ratio. Similar results
PAGE 48
are found for the waves in the ditchdrain streamlines. The occurrence of large values of oCR near the drain indicates the analytical solutions to and are the most inaccurate near the drain. The parameter, oCR, is significantly larger near the drain in the ditchdrain 37 problems. Table 1 is a list of streamline residence times at the points x= 10 m and 0.1 m along a streamline for the standard tiledrain problem of Figure 6A and the nonpenetrating and fully penetrating ditchdrain problems of Figures 9A and 9D respectively. TABLE 1: Streamline Residence Times Tile Drain 0.20 '0 0.40 '0 0.60 '0 0 .80 '0 x= 10 m: 4.0 9.8 18. 0 32.0 x=0.1 m: 4.2 10.0 18.3 32.5 Non Penetr. 0.20 '0 0.40 '0 0.60 '0 0.80 '0 Ditch x= 10 m: 4.2 10.1 18.5 33.1 x=0.1 m: 5.0 11.0 19.5 34.6 Full Penetr.: 0. 20 '0 0.40 '0 0.60 '0 0.80 '0 Ditch x= 10 m: 4.7 10.3 18.4 32.9 x=0.1 m: 13.8 11.3 19.0 33.4 An analysis of Table 1 reveals that the streamline residence times of the ditchdrain problems are more severely affected in the region near the drain. The increase in residence time, just between x = 10 m and 0.1 m in each
PAGE 49
of the tiledrain streamlines, is less than five percent of the total streamline residence time. The increases in residence time, between the same x interval in the nonpenetrating and fully penetrating ditchdrain streamlines, are as much as 16 percent and 65 percent respectively of the total streamline residence time. These increases near the drain are disproportionate and must be identified whenever using the ditchdrain analytical solutions. TileDrain Problem An examination of Figures 6AC reveals that residence time (both 10 and 1i), varies proportionally with recharge and porosity. Residence time is halved if the recharge rate is doubled (Figure 6b) and doubled if the porosity is doubled (Figure 6C). Residence time is not affected by changes in hydraulic conductivity. In Figure 60, residence time remains the same while K is doubled. Hydraulic conductivity is not a controlling parameter because it affects neither the size of the system nor the total recharge. Groundwater velocity is not affected by hydraulic conductivity, because neither the groundwater age nor the length of the streamlines is changed by K. If the depth of the flow system (Figure 7A) is increased by a factor of two, the average residence time, 10 increases by the same factor. However, the 1 i's,
PAGE 50
39 residence time of individual streamlines, do not all increase by the same factor. For example, the residence time along streamlines 0.20 !0 0.40 !0 0.60 !0 and 0.80 !0 are larger by factors of 1.83, 1.95, 1.99 and 2.03 respectively. The reason for the nonuniform effect of depth or lenqth is that a change in the flowsystem depthjlenqth ratio causes alteration of that system's geometry. There is a redistribution of the streamlines. Because velocities are greatest in the upper part of the flow system and systematically diminish with depth, streamlines displaced into shallower portions of the flow system have residence times less than predicted and those streamlines displaced into deeper portions of the flow system have residence times greater than what is predicted. Fiqures 7B and 7C illustrate how the position of streamlines are displaced as a function of increasing flow system depth. Fiqure 70 shows how changing the length of the system by a factor of two has only a very small effect on the position of the streamlines and streamline residence times. The change in length, unlike depth, does not affect the. average residence time (T0). If the depth and length of the flow system are changed proportionally, such that the D/L ratio remains unchanged, then the T0 and ri's change directly in proportion to the change in depth. This can be seen_in Fiqure 7E where the depth and length both increase by a factor of two, the D/L ratio remains 0.04, and the r0 and r;'s are exactly doubled.
PAGE 51
40 The dimensionless graph in Figure a shows the redistribution of the streamlines with different D/L ratios. The r1/r0 ratio decreases in the shallow streamlines and increases in the deeper streamlines, owing to the velocity distribution within the flow system. One can use this graph to compute the streamline residence time for any tiledrain problem of reasonable shape (D/L 0.5), provided the average residence time is known from the overall dimension and recharge rate. DitchDrain Problem Varying ditch penetration affects streamline length and postion only in the region near the ditch. In Figures 9AD, all the changes between corresponding streamlines in each ditchdrain problem are limited to a maximum of 50 meters from the ditch. The residence time of the streamlines is cumulated only to the point x = lO m because of the existence of large errors in the and V solutions between the ditchdrain and x=lO m (Figures SB and SC). A comparison of streamline residence times of four varying ditch penetration problems reveals that the change in residence time is small due to varying ditch penetration. The maximum change occurs in the 0.20 V0 streamline, which experiences a 10 % increase in residence time as a result of g 'oing from no ditch penetration (d= O .1 m), to full ditch
PAGE 52
penetration (d=20 m). Changes in residence time in the longer streamlines is progressively less as a function of ditch penetration. DupuitForschheimer DitchDrain Problem 41 The OF ditchdrain problem in Figure lOA, with controlling parameters identical to the tiledrain problem in Figure 6A, has a f0 which is 24 percent larger than the tiledrain f0 Individual ri's are as much as 17 percent larger in the OF problem. The existence of a watertable mound in the OF problem results in a larger volume in the system. A longer period of time, therefore, is required for ground water to flow through this larger volume system, and this results in the larger average residence time, f0 Similarly, individual streamlines are longer so all the ri's are larger. Increasing K in the OF ditchdrain problem, as is done in Figures lOBG, results in shrinkage of the watertable mound and a reduction in the r0 and individual ri's. The relationship between r0 and K is empirically shown to be nonlinear in the OF figures and consistent with equation 17. The OF ditchdrain problem in Figure lOG is for the case where K is 8000 mjyr. The watertable mound has essentially disappeared, resulting in a system which is geometrically similar to the tiledrain problem of Figure 6a. The average residence time in Figure lOG is less than
PAGE 53
three percent larger than in Figure 6A. The 0.20 'o streamline residence time's differ by approximately 6 42 percent while the 0.80 'o streamline residence times are identical. The OF ditchdrain problem assumes a fully penetrating drain with flow evenly distributed along the entire y axis. The tile drain flow system assumes no penetration, resulting in flow that is diverted back towards the tile drain along the x axis. The tile drain exerts a larger relative influence on the paths of the shorter streamlines which is why their r;'s do not agree as closely to the OF r ;'s. Figure 11 shows the effect of the watertable mound on the 0.80 '0 streamline. As (proportion which shows height of the flow problem) in the OF problem approaches 1 and the water table mound disappears, the difference in the streamline residence time between the standard tiledrain (r2 ) and the OF ditchdrain (r1 ) approaches o. The watertable mound is also infuenced by recharge (R). In figure 12A, R (0.2 mjyr) has been increased by a factor of two while K (100 mjyr) is kept constant. The height of the watertable mound increases significantly causing the r and r .s to decrease by amounts which are 0 1 less than factor of a two. In Figure 12B, where R=0.2 mjyr and k=8000 mjyr, the watertable mound has essentially disappeared and r0 and r;'s exhibit a decrease of factor two to within one percent.
PAGE 54
43 The length in a DF ditchdrain problem, unlike in the tile drain, influences the residence time. At the R/K ratio of lx10"3 (R=O.l mjyr, K=lOO mjyr) increasing the length of the flow system by a factor of 2 results in a 20 percent increase in T0 (Figure 13A) For the same length and a R/K ratio of 1.3xl0"5 (R=O.l mjyr, K=SOOO mjyr) the increase in T0 is limited to approximately 3 percent and to less than 1 percent in the T 1's (Figure 13B). Analysis of the controlling parameters R, K and L shows the height of the water table mound is controlled by the R/K ratio. The mound height is large at high R/K ratios and the mound height decreases as the R/K ratio decreases. The length of the flow problem influences the residence time when the mound is present and is dropped from consideration when the mound is absent. Comparison of Tileand DitchDrain Problems Streamline residence times at x=lO m for the standard tiledrain problem do not equal those for the nonpenetrating ditchdrain problem at the same point. The differences are even greater at x=O.l m. The two drain problems, however, are geometrically and hydrogeologically identical. This suggests that the differences in residence time are caused by their unique analytical solutions to and V. In the CauchyRiemann analysis, the tiledrain analytical solutions are shown to be more accurate, especially in the region near
PAGE 55
the drain. This accounts for the larger disparity in streamline residence time at x=O.l m as opposed to x=lO m. Residence times of the longer streamlines show better agreement between the tiledrain and fully penetrating ditchdrain problems than between the tiledrain and nonpenetrating ditchdrain problems. The residence times of the shorter streamlines show the reverse with better agreement between the tiledrain and nonpenetrating ditchdrain problems. Attempts to attribute these results to effects associated with varying ditch penetration is subject to much uncertainty. It appears, therefore, that the errors in the ditchdrain and analytical solutions obscure the exact magnitude of the ditch penetration effects, which, as discussed above, are small (Figure 9AD). From these considerations, it can be concluded that the residence time is more easily and accurately calculated, for the entire suite of various penetrating ditchdrain problems, using the potentialtheory solution for flow to a tiledrain (Kirkham and Powers, 1972). Comparison of Tileand DF DitchDrain Problems The absence of a watertable mound in the DF ditchdrain problem of Figure lOG allows comparison with the tiledrain problem of Figure 6A which has a similar geometry. The paths of corresponding streamlines in the two figures, vary in the region near the drain. This effect is caused by the
PAGE 56
45 difference in drain penetration: full drain penetration in the case of the OF ditchdrain and no drain penetration in the tile drain. A comparison of ri's reveals that the corresponding 0.20 'o streamline differ by a maximum 6 percent. Increasingly longer streamlines differ by progressively smaller percentages. The 0.80 'o streamline residence time is 32.5 years in both problems. These results show that the OF ditchdrain problem can nearly duplicate the accuracy of the ri's from the much more complicated tiledrain solution. The observed differences in ri's is thought to be caused by drain penetration effects. These effects are relatively small and progressively decrease away from the drain. Numerous studies have shown that most natural uplandtodrain systems have very small watertable slopes. Thus the comparison of the results for the tiledrain and OF ditchdrain problems indicatethat watertable mounds are negligible and the simple OF technique can be effectively applied to solve for ground water residence time in most any ditchdrain problem.
PAGE 57
CONCLUSIONS 1. Provided the geometry of the tiledrain problem remains uniform, the residence time is linearly affected by the parameters n, R and D. Changes in D and L which alter the D/L ratio redistributes the streamlines and slightly alters the streamline residence times. 2. The controlling parameters affect the residence times in ditchdrain problems in the same manner as in the tiledrain problems. Drain penetration does not affect the average residence time but does affect the streamline residence times. Errors in the and V solutions near the ditch drain prevent precise quantification of the relationship between drain penetration and streamline residence time. However, the study does indicate that the change in streamline residence time with increasing drain penetration is relatively minor. 3. Residence times for DupuitForschheimer ditchdrain problems are linearly affected by the parameters n, R h0 in the absence of a watertable mound. In the presence of a watertable mound, which is controlled by the R/K ratio, the residence time is nonlinearly affected by R, K and L. 4. The residence times for the entire suite of various penetrating ditchdrain problems are more easily and
PAGE 58
47 accurately calculated using the tiledrain potentialtheory solution instead of the ditchdrain solution. 5. In the absence of a watertable mound, the residence times calculated from the DupuitForschheimer ditchdrain solution very closely approximate those calculated from the more complicated tiledrain solution.
PAGE 59
REFERENCES Bengtsson, T.O., 1987 Flow net analysis of strip island lenses and calculations of timerelated parameters; Masters Thesis, University of South Florida. Fetter, C.W.,Jr., 1980 Applied Hydrogeology. Columbus: Merrill Publ. Co., 588 p. Fogg, G.E., and Senger, R.K., 1985 Automatic generation of flow nets with conventional groundwater modeling algorithms: Ground Water, v. 23, p. 336345 Fukuda, H.,1957 Underdrainage into ditches in soil overlying an impervious substratum: Transactions, American Geophysical Union, v. 38, n. 5, p. 730739 Kirkham, D., 1950 Seepage into ditches of a plane water table and an impervious substratum: Transactions, American Geophysical Union, v. 31, n. 3, p. 425430 Kirkham, D., 1958 Seepage of steady rainfall through soil into drains: Transactions, American Geophysical Union, v. 39, n. 5, p. 892907 Kirkham, D., 1967 Explanation of paradoxes in DupuitForchheimer seepage theory: Water Resources Research, v. 3, p. 609622 Kirkham, D., and Powers, W.L., 1972 Advanced Soil Physics. New York: WileyInterscience, 534 p. Robinson, J.L., 1991 Toth regional flow residence time analysis; Masters Thesis, University of South Florida Vacher, H.L., 1988a, DupuitGhybenherzberg analysis of strip island lenses: Geological Society of America Bulletin, v. 100, n. 4, p. 580591. Vacher, H.L., Bengtsson, T.O. 1989 Effect of hydraulic conductivity on the residence time of meteoric groundwater in Bermudian and Bahamiantype islands. in Mylorie, J. E. (ed.), Proceedings of the 4th Symposium of the Geology of the Bahamas. San Salvador, Bahamas: Bahamian Field Station, p. 337351.
PAGE 60
Vacher, H.L., Bengtsson, T.O., and Plummer, L.N., 1990, Hydrology of meteoric diagenesis: Residence time of meteoric ground water in island freshwater lenses with application to aragonitecalcite stabilization rate in Bermuda: Geological Society of American Bulletin, v. 102, p. 659676. Vacher, H.L., Farkas, T.A., and Robinson, J.L., 1991 Time net for groundwater flow in idealized coastal wedge: Journal of Coastal Research, v. 7, n. 1, p. 3138 49
PAGE 61
50 APPENDICES
PAGE 62
51 APPENDIX A Potential Theory Solutions to DitchDrain Model
PAGE 63
Potential Theory Solution for Flow to Fully Penetrating DitchDrain (Kirkham, 1958). where S1= tan[sin(7rX/L)/{e7)'/L COS(7TX/L)}) S2= L ( 1/m) sin (m7rx/L) e mwh/L sinh (m7ry /L) ;sinh (m7rh/L) S3= tan[sin (1ry /h) I { exth cos (7ry /h) } ) S4= L ( 1/m) sin (m1ry /h) emwL/h sinh (m7rx/h) ;sinh (m7rL/h) 52 The limits of L is between 1 and oo. After some rearranging, the solution for the stream function is: lf= 21fc/7r {S1S2+S3S4} Using the Cauchy Reimann conditions, Kirkham (1958) developed the following equation which satisfies t. where: 7r(tc)/2lf0 = T1T2+T3+T4 T1= (1/2) ln[2eyY/L (cosh(7ry/L) COS(7TX/L))) T2= L ( 1/m) cos (m7rX/L) e mwh/L cosh (m1ry /L) I sinh (m7rh/L) T3= ( 1/2) ln [ 2exth (cosh (7rx/h) cos (1ry /h))] T4= L ( 1/m) cos (m1ry /h) emwL/h cosh (m7rxjh) /sinh (m7rL/h)
PAGE 64
APPENDIX A (Continued) C/K = 21fc/7rK { ln (2simrr/2L) L ( 1/m) cos (m7rr/L} [ (cosh(m7rL/h)/sinh(m7rL/h)) 1] ln(2eIT/2h sinh(1Tr/2h)) + L ( 1/m) eIIII'L/h cosh (m1rrjh) ;sinh (m7rL/h) The limits of L is between 1 and Because t= K, the solution may be rearranged to show: = 21fc/K1T (T1T2+T3+T4) +C/K 53 Potential Theory Solution for Flow to NonPenetrating DitchDrain (Kirkham, 1958) 7r ( 1f 0 ') I ( 2 '0) =S 1s 2 where S1= tan[sin (1TX/L) 1 { eyY/ L cos (1TX/L} } ] S2= L ( 1 /m) sin (rn1TX/L) ellll'h/ L sinh (m7ry /L) ;sinh (m7rh/ L) The limits of L is between 1 and oo. Afte r some rearranging, the solution for the stream function is: 1f= 21fc/7r {S1S2} Using the Cauchy R eimann conditions, Kirkham (1958)
PAGE 65
APPENDIX A (Continued) developed the following equation which satisfies t. where: = T1T2 T1= (1/2) ln[2e7)'/L (cosh(nyjL) T2= L(1/m) e mm/L C/K = { L(1/m) [ 1] ln(2enl2h 54 sinh + L ( 1/m) eIIII'L/h cosh /sinh The limits of L is between 1 and oo. Because t= the solution may be rearranged to show: 2 V ( T1T2) +C/K
PAGE 66
55 APPENDIX B Basic Computer Algorithm for Tileand DitchDrain Problems
PAGE 67
10 'Program modifed by Tom Farkas to calculate streamline residence times for tiledrain problem 11 start at input y=O and input x 12 increments x at 50 m 30 OPEN "VALlg.DAT" FOR OUTPUT AS #1 32 OPEN "VAL2g.DA':f" FOR OUTPUT AS #2 34 OPEN "VAL3g.DAT" FOR OUTPUT AS #3 36 OPEN "VAL4g.DAT" FOR OUTPUT AS #4 50 READ R, D, S, K, A, POR 51 DATA 0.1, 40, 1000, 100, .1, .1 60 STRMO=R*(SA): PI=3.14159 70 PRINT "Trace of streamline in tiledrain problem" 71 PRINT" R(mjyr), D(m), S(m), K(m/d), A(m),porosity are:" R;D;S;K;A;POR 90 RESTIME=S*D*POR/(R*(SA)) 'residence time 91 PRINT "residence time (days, yrs)="RESTIME/365; RESTIME 100 GOSUB 1000 find potl at x=s,y=O 200 'FIND STRM, POTL VS X,Y 210 INPUT "x (start at 100, then x will be incremented by lOO)";X 220 Y=O: XSTORE=X: TIMETOTL=O: DISTOTL=O: CNT=CNT+l 222 IF XSTORE>S.5 THEN GOTO 950 PRINT ">>>>> FOR x="X 230 290 Ul=PI*X/S: U2=PI*Y/S: U3=PI*D/S: U4=PI*(DY)/S: U5=PI*A/2/S: U7=PI*A/S GOSUB 1200 PRINT X PRINT X,Y,STRM,POTL STRMCOMP=STRM y XOLD=X: YOLD=Y: POTLOLD=POTL X=X1 IF X>S11 THEN Y=O: GOTO 420 STRM IF X<5 AND X>O THEN IF X<.Ol THEN Y=.5: Y=.5: GOTO 420 X=O: GOTO 425 Y=Y1 IF Y>>>>>> velocities and times between potentials 620 630 V=K*(POTLPOTLOLD)/DIST/POR 650 TIMEINC=DIST/V: PRINT," travel time (yrs)="TIMEINC 56
PAGE 68
APPENDIX B (Continued) 660 TIMETOTL=TIMETOTL+TIMEINC : PRINT ," residence time(yrs)="TIMETOTL 670 DISTOTL=DISTOTL+DIST: PRINT ," cumulated distance="DISTOTL 680 700 710 720 730 Y=Y*1 IF CNT=1 IF CNT=2 IF CNT=3 IF CNT=4 740 Y=Y*1 THEN WRITE THEN WRITE THEN WRITE THEN WRITE #1,X,Y,TIMETOTL,POTL #2,X,Y,TIMETOTL,POTL #3,X,Y,TIMETOTL,POTL #4,X,Y,TIMETOTL,POTL 800 IF X=O THEN X=XSTORE+200: PRINT X X X X X X " :print: 220 810 950 close #1 952 #2 954 #3 956 #4 970 end 1qqq i 1010 setup for nondimensioning: find pot'l at s,o 1q2q 1030 1040 u1="PI*X/S:" 1050 1060 potlo="POTL" 1070 1080 print,"strm, strmjstrmo, are:" 1090 print,strm, strm strmo, : 1095 return 12qq 1201 compute stream function 1202 1209 sum="O" 1210 m="1" to 50 1220 sin1="SIN(M*U1)" 1230 1240 sinh2="(EXP(M*U3)EXP(M*U3))/2" 1250 1251 'print m, 1260 next 1270 1285 1 strmo 13qq 1301 potential 1302 1303 1310 25 1320 cos1="COS(M*U7):" cos2="COS(M*U1)" 1330 cosh="(EXP(M*U2)+EXP(M*U2))/2" 57
PAGE 69
APPENDIX B (Continued) 1340 SINH=(EXP(M*U3)EXP(M*U3))/2 1350 SUM=SUM+(1/M)*(COS1COSH*COS2)*EXP(M*U3)/SINH 1351 'PRINT M,SUM 1360 NEXT M 1370 COSH=(EXP(U2)+EXP(U2))/2 1380 U6=(EXP(U2))*(COSHCOS(U1))/2/SIN(U5)/SIN(U5) 1390 POTL=2*R*S/K/PI*(.5*LOG(U6)+SUM) 1400 1 PRINT POTL 1500 RETURN 58
PAGE 70
APPENDIX B (Continued) 10 'Program written by Tom Farkas to calculate streamline 11 1 residence time for Kirkham's DitchDrain problem 12 1 increments x at 50 m 30 OPEN 11biTl.DAT11 FOR OUTPUT AS #1 32 OPEN 11biT2.DAT11 FOR OUTPUT AS #2 34 OPEN "biT3.DAT11 FOR OUTPUT AS #3 36 OPEN 11biT4.DAT11 FOR OUTPUT AS #4 50 READ R, S, K, H, A, POR 51 DATA 0.1, 500, 100, 20, .1, .1 55 D=20 60 STRMO=R*(S): PI=3.14159 65 POTL0=(2*R*S)/(PI*K) 70 PRINT "Trace of streamline in tiledrain problem" 71 PRINT 11 R(mjyr), S(m), K(m/d), H(m), A(m) ,porosity 59 are:" R;S;K;H;A;POR 90 RESTIME=S*H*POR/(R*S) 'residence time 91 PRINT "residence time (days, yrs)="RESTIME/365; RESTIME 100 GOSUB 1000 find potl at x=s,y=O 200 'FIND STRM, POTL VS X,Y 210 INPUT 11x (start at 100, then x will be incremented by 100)11;X 220 Y=O: XSTORE=X: TIMETOTL=O: DISTOTL=O: CNT=CNT+l 222 IF XSTORE>S.5 THEN GOTO 950 230 PRINT 11>>>>> FOR x="X 250 Fl=(PI*X/S) :F2=(PI*Y/S) :F3=(PI*H/S) :F4=(PI*X/H) 260 F5=(PI*Y/H) :F6=(PI*S/H) :F7=(PI*A)/(2*S) :F8=(PI*A)/(2*H) 270 F9=H/(PI*D): FlO=(PI*D)/H 300 GOSUB 1200 310 PRINT 11 X Y STRM POTL11 320 PRINT X,Y,STRM,POTL 350 IF Y=O THEN GOSUB 1800 400 STRMCOMP=STRM 410 XOLD=X: YOLD=Y: POTLOLD=POTL 412 X=X1 415 Y=Y1: IF Y>>>>>> velocities and times between potentials 620 DIST=SQR((XXOLD)A2+(YYOLD) A2)
PAGE 71
APPENDIX B (Continued) 630 V=ABS(K*(POTLPOTLOLD)/DIST/POR) 650 TIMEINC=DIST/V: PRINT," travel time (yrs)="TIMEINC 660 TIMETOTL=TIMETOTL+TIMEINC: PRINT ," residence time(yrs)="TIMETOTL 670 DISTOTL=DISTOTL+DIST: PRINT ," cumulated distance="DISTOTL 680 Y=Y*1 700 IF CNT=1 710 IF CNT=2 720 IF CNT=3 730 IF CNT=4 740 Y=Y*1 THEN WRITE THEN WRITE THEN WRITE THEN WRITE #1,X,Y,TIMETOTL,POTL,V #2,X,Y,TIMETOTL,POTL,V #3,X,Y,TIMETOTL,POTL,V #4,X,Y,TIMETOTL,POTL,V 60 800 IF X=.1 THEN X=XSTORE+100: PRINT 11 X X X X X X " :print: 220 810 410 950 close #1 952 #2 954 #3 956 #4 970 end 1000 i 1010 1 setup for nondimensioning: find pot11 at s,o 1020 1030 1035 f1="(PI*X/S)" 1045 1047 f9="H/(PI*D):" f10="(PI*D)/H" 1050 1070 1080 print,"strm, strmjstrmo, potl are:" 1090 print,strm, strmo, : 1095 return compute stream function 1201 1202 s2sum="O:" s4sum="O:" s5sum="O:" s6sum="O" 1204 t1sum="O:" t2sum="O:" t4sum="O:" t6sum="O:" t8sum="O" 1205 z1sum="O:" z2sum="O:" z3sum="O:" z4sum="O" 1210 s1="ATN(SIN(F1)/(EXP(F2)COS(F1)))" 1220 s3="ATN(SIN(F5)/(EXP(F4)COS(F5)))" 1230 m="1" to 10 1240 tsinh1="(EXP(M*F3)EXP(M*F3))/2" 1249 tsinh2>1E+09 THEN GOTO 1260 1250 TSINH2=(EXP(M*F6)EXP(M*F6))/2 1260 TSINH4=(EXP(M*F2)EXP(M*F2))/2 1270 S2=(1/M)*SIN(M*F1)*EXP(M*F3)*TSINH4/TSINH1 1279 IF TSINH5>1E+09 THEN GOTO 1290 1280 TSINH5=(EXP(M*F4)EXP(M*F4))/2 1290 S4=(1/M)*SIN(M*F5)*EXP(M*F6)*TSINH5/TSINH2 1292 S5=((1/(M*M))*SIN(M*F10)*SIN(M*F5)*EXP(M*F4))
PAGE 72
APPENDIX B (Continued) 1294 S6=((1/(M*M))*SIN(M*F10)*SIN(M*F5)*EXP(M*F6)* (TSINH5/TSINH2)) 1300 S2SUM=S2+S2SUM 1310 S4SUM=S4+S4SUM 1312 S5SUM=S5+S5SUM 1314 S6SUM=S6+S6SUM 1320 NEXT M 1322 S5S=F9*S5SUM 1324 S6S=F9*S6SUM 61 1330 STRM=((2*STRMO/PI)*(S1S2SUM+S3S4SUMS5S+S6S))+STRMO 1350 1 COMPUTE POTENTIAL 1360 FOR M=1 TO 10 1370 T1=(1/M)*COS(M*F1)*EXP(M*F2) 1380 TCOSH1=(EXP(M*F2)+EXP(M*F2))/2 1390 TSINH1=(EXP(M*F3)EXP(M*F3))/2 1400 T2=(1/M)*COS(M*F1)*EXP(M*F3)*(TCOSH1/TSINH1) 1410 TCOSH2=(EXP(F4)+EXP(F4))/2 1420 T3=.5*LOG(2*EXP(F4)*(TCOSH2COS(F5))) 1429 IF TSINH2>1E+09 THEN GOTO 1440 1430 TSINH2=(EXP(M*F6)EXP(M*F6))/2 1440 T4=(1/M)*COS(M*F5)*EXP(M*F6)*(TCOSH2/TSINH2) 1450 TCOSH3=(EXP(M*F3)+EXP{M*F3))/2 1460 T6=(1/M)*COS(2*M*F7)*((TCOSH3/TSINH1)1) 1470 TCOSH4=(EXP(2*M*F8)+EXP(2*M*F8))/2 1480 T8=(1/M)*EXP(.;..M*F6)*(TCOSH4/TSINH2) 1482 Z1=((1/{M*M))*SIN(M*F10)*COS(M*F5)*EXP(M*F4)) 1484 Z2=((1/(M*M))*SIN(M*F10)*COS(M*F5)*EXP(M*F6)* {TCOSH2/TSINH2) ) 1486 Z3={(1/(M*M))*SIN(M*F10)*EXP{M*2*F8)) 1488 Z4=((1/(M*M))*SIN(M*F10)*EXP(M*F6)*{TCOSH4/TSINH2)) 1490 1 SUMMATION 1500 T1SUM=T1+T1SUM 1510 T2SUM=T2+T2SUM 1520 T4SUM=T4+T4SUM 1530 T6SUM=T6+T6SUM 1540 T8SUM=T8+T8SUM 1542 Z1SUM=Z1+Z1SUM 1544 Z2SUM=Z2+Z2SUM 1546 Z3SUM=Z3+Z3SUM 1548 Z4SUM=Z4+Z4SUM 1550 NEXT M 1560 T5=LOG(2*SIN(F7)) 1570 TSINH3=(EXP(F8)EXP(F8))/2 1580 T7=LOG{2*EXP(F8)*TSINH3) 1582 Z1S=Z1SUM*F9 1584 Z2S=Z2SUM*F9 1586 Z3S=Z3SUM*F9 1588 Z4S=Z4SUM*F9
PAGE 73
APPENDIX B (Continued) 1590 POTL=(POTLO*(T1SUMT2SUMT3+T4SUMZ1SZ2S))(POTLO*(T5T6SUMT7+T8SUMZ3SZ4S)) 1600 RETURN 1650 END 1710 RETURN 1800 IF CNT=1 THEN WRITE #1,X,Y,TIMETOTL,POTL 1810 IF CNT=2 THEN WRITE #2,X,Y,TIMETOTL,POTL 1820 IF CNT=3 THEN WRITE #J,X,Y,TIMETOTL,POTL 1830 IF CNT=4 THEN WRITE #4,X,Y,TIMETOTL,POTL,STRM 1840 RETURN 62
PAGE 74
63 APPENDIX C Basic Computer Algorithm for DupuitForschheimer DitchDrain Problem
PAGE 75
64 10 'THIS PROGRAM (DUPUIT) WILL CALCULATE RESIDENCE TIMES IN 20 'A FLOW SYSTEM OF TWO FULLY PENETRATING RIVERS WITH EQUAL 30 'WATER ELEVATION. 40 OPEN 11REN1.DAT11 FOR OUTPUT AS #1 42 OPEN 11REN2.DAT11 FOR OUTPUT AS #2 44 OPEN 11REN3.DAT11 FOR OUTPUT AS #3 46 OPEN 11REN4.DAT11 FOR OUTPUT AS #4 50 READ R,L,DX,K,H1,H2,POR 55 DATA 0.1,1000,1,8000,20,20,0.1 200 1 CALCULATION OF TOTAL FLOW ON ONE SIDE OF DIVIDE 205 X=L/2 210 TOTFL=(R*L/2) 215 PRINT:PRINT 11X AT DIVIDE IS11;X 220 PRINT "TOTAL FLOW="; TOTFL 225 STRM=TOTFL 240 H=SQR(H1A2(((H1A2H2A2)/L)*X)+((R/K)*(LX)*X)) 245 PRINT "HEAD AT DIVIDE=";H 250 FOR J=1 TO 4 252 TIMEINC=O: TIMETOTL=O :DIST=O 255 X=X(L*.1) 260 STRM=STRM(R*L/10) 270 H=SQR(H1A2(((H1A2H2A2)/L)*X)+((R/K)*(LX)*X)) 275 PRINT 277 PRINT 11X IS11X 280 PRINT "STREAMLINE VALUE=11STRM 290 PRINT "HYDRAULIC HEAD=11H 300 GOSUB 500 310 NEXT J 312 CLOSE #1 313 CLOSE #2 314 CLOSE #3 315 CLOSE #4 350 END 500 1 CALCULATE X,W COORDINATES 502 PRINT: PRINT 11 X W DIST TIMEINC TIMETOTL 505 X1=X 510 FOR X1=X1 TO 0 STEP OX 520 H=SQR(H1 A2(((H1A2H2A2)/L)*Xl)+((R/K)*(LX1)*X1)) 530 W=((TOTFLSTRM)/(R*(L/2X1)))*H 590 IF H=W THEN GOTO 650 600 1 COMPUTE DISTANCE, VELOCITY AND TRAVELTIME 610 DIST=SQR((DX)A2+(WWOLD) A2) 620 V=K*(HHOLD)/DIST/POR 630 TIMEINC=DIST/V 635 IF TIMEINC
PAGE 76
APPENDIX C (Continued) 700 HOLD=H: WOLD=W 705 OIST=O 710 NEXT X1 720 RETURN 800 WRITE #1,X1,W,TIMETOTL,H :GOTO 700 810 WRITE #2,X1,W,TIMETOTL,H :GOTO 700 820 WRITE #3,X1,W,TIMETOTL,H :GOTO 700 830 WRITE #4,X1,W,TIMETOTL,H :GOTO 700 65
