EXPERIMENTAL AND NUMERICAL ANALYSIS OF GROUNDWATER ? SURFACEWATER INTERACTION IN HORIZONTAL AND SLOPING UNCONFINED AQUIFERS Except where reference is made to the work of others, the work described in this thesis is my own or was done in collaboration with my advisory committee. This thesis does not include proprietary or classified information. __________________________ Gopal Chandra Saha Certificate of Approval: ________________________ ________________________ Joel G. Melville T. Prabhakar Clement, Chair Profesor Profesor Civil Engineering Civil Engineering ________________________ ________________________ Xing Fang George T. Flowers Associate Professor Dean Civil Engineering Graduate School EXPERIMENTAL AND NUMERICAL ANALYSIS OF GROUNDWATER ? SURFACEWATER INTERACTION IN HORIZONTAL AND SLOPING UNCONFINED AQUIFERS Gopal Chandra Saha A Thesis submitted to the Graduate Faculty of Auburn University in Partial fulfillment of the requirements for the degree of Master of Science Auburn, Alabama August 10, 2009 iii EXPERIMENTAL AND NUMERICAL ANALYSIS OF GROUNDWATER ? SURFACEWATER INTERACTION IN HORIZONTAL AND SLOPING UNCONFINED AQUIFERS Gopal Chandra Saha Permission is granted to Auburn University to make copies of this thesis at its discretion, upon the request of individuals or institutions and at their expense. The author reserves all publication rights. Gopal Chandra Saha Signature of Author August 10, 2009 Date of Graduation iv VITA Gopal Chandra Saha, son of Santosh Kumar Saha and Niva Rani Saha, was born on July 1 st , 1979. In March 2003, he graduated from Bangladesh University of Engineering and Technology (BUET) with a Bachelor?s of Science degree in Civil Engineering. After graduation, he worked in a construction company for one year. Then he did his Master?s of Science degree in Water Resources Engineering and Management from Stuttgart University, Germany in December 2006. In January 2008, he started coursework to pursue his Master?s of Science degree in Civil Engineering in Auburn University under the guidance of Prof. T. Prabhakar Clement. He married Pinki Saha, daughter of Utkott Saha and Jyotti Saha, on December 6 th , 2008. v THESIS ABSTRACT EXPERIMENTAL AND NUMERICAL ANALYSIS OF GROUNDWATER ? SURFACEWATER INTERACTION IN HORIZONTAL AND SLOPING UNCONFINED AQUIFERS Gopal Chandra Saha Master of Science, August 10, 2009 (B.Sc., BUET, Bangladesh, 2003) 84 Typed Pages Directed by T. Prabhakar Clement We investigated the interaction between the water moving from an unconfined aquifer and a surfacewater body. This interaction can occur during or after a rainfall event, which can cause the groundwater table to change. We developed a laboratory-scale aquifer model, which is bounded by a stream and a water divide. The model simulates both horizontal and sloping aquifer conditions. We conducted both steady-state and transient experiments. First, we set the water table at a steady state using a constant recharge rate and then the aquifer was allowed to drain. The experimental data were vi recreated using a numerical model, which simulated the observed water table variations. The transient water table data indicated highly non-linear changes. It was found that under similar recharge conditions, water table in the horizontal unconfined aquifer is higher than that in the sloping unconfined aquifer (2.03 degree slope). Also, the horizontal unconfined aquifer required more time to drain. The numerical model was able to predict the observed groundwater table data. The proposed numerical model is a useful tool for managing groundwater-surfacewater management problems. vii ACKNOWLEDGEMENTS First of all, I would like to thank my advisor Prof. T. Prabhakar Clement for his enormous guidance, suggestion and inspiration in completing the thesis successfully. His constructive criticism has helped me to organize the work. I would also like to thank Professors Joel Melville and Xing Fang for serving as my committee members, and for their support. Special thanks to Jagadish, Shyam, Sunwoo and Eric for their help in doing lab work. I would also like to show my deepest gratitude to my parents and my family members for encouraging me during my whole degree program. This research was, in part, funded by the Auburn University WRC-Wolf Bay Basin Project # 101002-107311- 700-2055. viii Style manual or journal used Auburn University Graduate School Guide to Preparation of Master?s Thesis, Water Resources Research Journal. Computer software used Microsoft Office 2003 (Microsoft Word, Microsoft Excel), Visual Basic. ix TABLE OF CONTENTS LIST OF FIGURES ........................................................................................................... xi LIST OF TABLES........................................................................................................... xiv CHAPTER 1. INTRODUCTION .......................................................................................1 1.1 Background ......................................................................................................1 1.2 Objectives .........................................................................................................2 1.3 Organization of thesis .......................................................................................3 CHAPTER 2. LITERATURE REVIEW .............................................................................4 2.1 Unconfined hillslope flow .................................................................................4 2.2 Derivation of the governing flow equation......................................................10 2.3 Review of previous works related to groundwater-surfacewater Interactions.12 2.3.1 Review of analytical solutions to the governing flow equation ..........13 2.3.2 Review of numerical solutions.............................................................19 2.3.3 Review of experimental studies...........................................................21 CHAPTER 3. EXPERIMENTAL METHODS .................................................................23 3.1 Column experiments ........................................................................................23 3.2 Darcy experiments ...........................................................................................27 x CHAPTER 4. DEVELOPMENT AND TESTING OF A NUMERICAL MODEL .........30 4.1 Numerical solution strategy .............................................................................30 4.2 Initial and boundary conditions .......................................................................33 4.3 Code verification using published datasets......................................................34 4.3.1 Comparison against Kim et al (2001) horizontal aquifer solution......34 4.3.2 Comparison against Verhoest et al (2002) sloping aquifer solution...36 4.3.3 Comparison against Sanford et al (1993) solution..............................38 CHAPTER 5. LABORATORY DATA AND MODELING RESULTS...........................42 5.1 Steady state experiments...................................................................................42 5.2 Modeling steady-state water table profiles .......................................................48 5.3 Transient experiments.......................................................................................51 5.3.1 Modeling the transient head boundary at the left side of the tank.........51 5.3.2 Transient solution for horizontal unconfined aquifer ............................54 5.3.3 Transient solution for sloping unconfined aquifer.................................57 5.3.4 Accuracy of the numerical model..........................................................58 5.4 Correlation between the sloping unconfined aquifer with the horizontal unconfined aquifer ............................................................................................59 CHAPTER 6. SUMMARY, CONCLUSIONS AND RECOMMENDATIONS ..............61 REFERENCES ..................................................................................................................64 APPENDIX........................................................................................................................68 xi LIST OF FIGURES Figure 1: The schematic presentation of sloping bed unconfined aquifer...........................5 Figure 2: An arbitrary point in the porous medium of sloping bed unconfined aquifer.....6 Figure 3: Elevation of the water table in sloping bed unconfined aquifer...........................8 Figure 4: Horizontal bed unconfined aquifer under constant recharge............................10 Figure 5: Sloping bed unconfined aquifer under constant recharge ..................................10 Figure 6: 2-D control volume in an unconfined aquifer ...................................................11 Figure 7: Column experiment setup in the laboratory .......................................................24 Figure 8: Darcy experiment setup in the laboratory ..........................................................28 Figure 9: Comparison Kim et al. (2001) analytical model (linearized) results against our proposed model (non-linear) results for a horizontal unconfined aquifer example problem...............................................................................................35 Figure 10: Comparison between analytical model results of Kim et al. (2001) and the linearized numerical model for a horizontal unconfined aquifer example problem ............................................................................................................36 Figure 11: Comparison Verhoest et al. (2002) analytical results against linear and non- linear numerical models results for a sloping unconfined aquifer example problem ............................................................................................................37 xii Figure 12: Comparison Sanford et al. (1993) analytical model results against our proposed model (non-linear) results for a horizontal unconfined aquifer example problem...............................................................................................39 Figure 13: Comparison Sanford et al. (1993) analytical model results against our proposed model (non-linear) results for a sloping unconfined aquifer (theta = 0.089 degree)example problem ...........................................................40 Figure 14: Comparison Sanford et al. (1993) analytical model results against our proposed model (non-linear) results for a sloping unconfined aquifer (theta = 0.137 degree)example problem ...........................................................41 Figure 15: Laboratory tank with recharge generator used in this study for horizontal unconfined aquifer ..........................................................................................44 Figure 16: Laboratory tank with recharge generator used in this study for sloping unconfined aquifer (i = 2.03 degree)...............................................................44 Figure 17: A schematic laboratory set up for tank experiment..........................................45 Figure 18: Experimental Steady-state water table profiles in horizontal unconfined aquifer under different recharge rates...............................................................46 Figure 19: Experimental Steady-state water table profiles in sloping unconfined aquifer (i=2.03degree) under different recharge rates..................................................47 Figure 20: Modeling steady-state water table profiles using analytical equation in horizontal unconfined aquifer...........................................................................49 Figure 21: Modeling steady-state water table profiles using numerical model in sloping unconfined aquifer (i = 2.03 degree)...................................................50 xiii Figure 22: Modeling transient head boundary (left hand side of the tank) in horizontal unconfined aquifer ............................................................................................52 Figure 23: Modeling transient head boundary (left hand side of the tank) in sloping unconfined aquifer (i = 2.03 degree).................................................................53 Figure 24: Transient behavior of water table in horizontal unconfined aquifer (k = 90 cm/min (1296 m/day), y S = 0.35 and w = 5.4 cm/min at the beginning (steady state))................................................................................................................55 Figure 25: Transient behavior of water table in horizontal unconfined aquifer (k = 90 cm/min (1296 m/day), y S = 0.35 and w = 4.42 cm/min at the beginning (steady state))...................................................................................................56 Figure 26: Transient behavior of water table in sloping unconfined aquifer (i = 2.03 degree, k = 90 cm/min (1296 m/day), y S = 0.35 and w = 4.42 cm/min at the beginning (steady state))................................................................................57 Figure 27: Comparison of experimental steady state water table profile between horizontal and sloping unconfined aquifers under same condition .................60 Figure 28: Pressure at the bottom of a sloping bed unconfined aquifer ...........................68 xiv LIST OF TABLES Table 1: Results of Column experiment ............................................................................24 Table 2: Darcy velocities and Reynolds numbers in column experiments........................26 Table 3: Results from Darcy experiments .........................................................................28 Table 4: Darcy velocities and Reynolds numbers in Darcy experiments ..........................29 Table 5: Parameters used in Sanford et al. (1993) analytical model .................................38 Table 6: Recharge rates and constant head boundary conditions for both horizontal and sloping unconfined aquifers (i = 2.03 degree) .....................................................43 Table 7: Mass balance errors for the water table profiles generated by the numerical model for the horizontal unconfined aquifer ........................................................59 Table 8: Mass balance errors for the water table profiles generated by the numerical model for the sloping unconfined aquifer (i = 2.03 degree) .................................59 1 CHAPTER 1 INTRODUCTION 1.1 Background Interaction between groundwater and surfacewater systems is a common phenomenon observed in nature. This interaction can occur during or after a rainfall event. When rain water infiltrates into the ground it raises the water table which will then induce groundwater drainage to a surface water body. For sustainable groundwater management, it is important to understand the dynamics of groundwater table movement during and after recharge events. Groundwater discharge contribution to stream flow is known as groundwater runoff or base flow. Traditionally it is assumed that during precipitation, stream flow is generated primarily from surface runoff; while during dry period, stream flow is generated primarily from groundwater flow (Hall, 1968, Norum et al., 1968). However, recent studies have shown that groundwater flow (so called old water) can dominate stream flow even during a precipitation event (Shanley et al., 2002). Typically, base flow 2 is a good indicator of aquifer storage characteristics within a groundwater basin (Meyboom, 1961). In order to estimate base flow, a rating curve of base flow can be made by plotting mean groundwater table within a basin against stream flow during dry periods (Schneider, 1961). During or after a rainfall event in a small basin the groundwater table could rise rapidly and contribute to the base flow. The variations in base flow, with time during and after rainfall event within a basin is indicated by a recession curve (Chow et al., 1988). Hence, recession curve is a good measure of drainage rate from a basin (Werner et al., 1951). Groundwater table response to a recharge event may depend on the topography of the aquifer bottom. It is necessary to understand this response because sloping hillslopes are the basic landscape elements in many catchments (Torch et al., 2003). In addition, hillslope subsurface flow is one of the key elements that control the hydrology of upland watershed under both wet and dry conditions (Brutsaert, 1994). 1.2 Objectives We investigated the hydraulics of unconfined flow near a surface water body during and after a rainfall event under both horizontal and sloping aquifer conditions. We developed a laboratory-scale aquifer model to conduct multiple experiments. A numerical model, which is based on a fully implicit finite difference scheme with Picard iteration, was used to model the experimental results. 3 1.3 Organization of thesis This thesis is organized into six chapters including this introduction chapter. Chapter 2 provides the literature review and derivation of the governing equations for both horizontal and sloping-bed unconfined aquifers. It also summarizes previous works related to groundwater-surfacewater interactions. Chapter 3 describes all the experimental methods employed in this study. Chapter 4 details the development of a numerical model and the validation of the model results based on published analytical and numerical results. Chapter 5 discusses the modeling of laboratory data, and the transient model solution for both horizontal and sloping unconfined systems. Chapter 6 summarizes the conclusion of this study and provides some recommendation for future research in this area. An alternate analysis of depth-averaged Darcy flux in sloping bed unconfined aquifer is listed in the appendix. 4 CHAPTER 2 LITERATURE REVIEW 2.1 Unconfined hillslope flow Changes in the groundwater table in an unconfined aquifer after a recharge event has been studied by many researchers. Boussinesq (1877) first developed the general depth averaged flux equation for unconfined flow in a sloping aquifer based on depth averaged Darcy?s law. Figure 1 shows the schematic presentation of a sloping bed unconfined aquifer. In Figure 1, x is along the bed,  x is the Cartesian coordinate, and i is the slope. The relationship between x and  x is given by the following equation:  cos x i x = (1) 5 Figure 1: The schematic presentation of sloping bed unconfined aquifer. In Figure 2, an arbitrary point (x, Z) is considered in the porous medium of sloping bed unconfined aquifer. The elevation of point (x, Z) = z cos sinzZ ix i=+ (2) From Darcy?s law, H vK x ? =? ? (3) Where v is Darcy flux or Darcy velocity or specific discharge L T ? ? ? ? ? ? , K is the hydraulic conductivity L T ?? ?? ?? , and H x ? ? is hydraulic gradient. h i x=0 X=L 1 h 2 h D z  x x i 6 Figure 2: An arbitrary point in the porous medium of sloping bed unconfined aquifer. Total head is defined as: 2 2 pv H z g? = ++ (4) where p ? is pressure head, 2 2 v g is velocity head, and z is elevation head. Since the groundwater velocity (v) is typically a very small, we can neglect the velocity term. Hence, p Hz ? =+ (5) Darcy velocity in Z direction z H vK Z ? =? ? x Z z Z (x,Z) i 7 z P vK z Z ? ??? =? + ?? ? ?? (6) If the flow is parallel to the bed, 0 z v = . After substituting 0 z v = in (6) we get, 0 P Kz Z ? ??? ?+= ?? ? ?? (7) After differentiating we get 1 0 Pz ZZ? ?? += ?? (8) From Figure 2 (or differentiating (2)) we can obtain cos z i Z ? = ? ; and using this we get, 1 cos 0 P i Z? ? += ? (9) Darcy velocity in x direction, x P vK z x ? ??? =? + ?? ? ?? 1 x P z vK K xx? ?? =? ? ?? (10) From Figure 2,sin z i x ? = ? ; and using this we can get, 1 sin x P vK Ki x? ? =? ? ? (11) 8 x Z Z i h h = h(x) Figure 3: Elevation of the water table in sloping bed unconfined aquifer. The elevation of the water table function (h = h(x)) in sloping bed unconfined aquifer is shown in Figure 3. Using this figure we can develop the following boundary condition for the water table: at Z h= , 0P = From equation (8) cos P i Z ? ? =? ? (12) After integrating cos ( )P iZ f x?=? + (13) 9 Using the water table boundary condition we get, () cosf xih?= After substituting f(x) in (13) we get, cos cosP iZ ih? ?=? + (14) Substituting P in equation (11) we get, (cos cos) sin x K viZihK x ?? ? ? =? ? + ? ? cos sin x Kh viKi x ? ? ? =? ? ? cos sin x h vKi i x ? ?? =? + ?? ? ?? (15) This also shows that () xx vvZ? Depth-averaged Darcy flux can be defined as: x qvh ? = (16) where q ? is depth averaged flux along the x direction 2 L T ? ? ? ? ? ? . After substituting (15) in (16) we get the general depth averaged flux equation for unconfined flow in a sloping aquifer as: () cosqKhhiz x ? ? =? + ? cos hz qKh i x x ? ?? ?? =? + ?? ?? ?? cos sin h qKh ii x ? ? ?? =? + ?? ? ?? (17) 10 2.2 Derivation of the governing flow equation The Boussinesq equation for a sloping unconfined aquifer with a constant recharge rate can be formulated by combining equation (17) with the continuity equation. Schematic of the horizontal and sloping bed unconfined aquifers considered in this study are shown in Figures 4 and 5. h w i X=0 X=L 1 h 2 h h w 2 h 1 h X=0 X=L Figure 4: Horizontal bed unconfined aquifer under constant recharge Figure 5: Sloping bed unconfined aquifer under constant recharge 11 Now, let us consider a 2-D control volume for groundwater occurring in these systems, as illustrated in Figure 6. Figure 6: 2-D control volume in an unconfined aquifer Applying mass balance over this control volume for a system with a constant recharge rate (w), we get the following conservation equation: xxx yyy y qq ytqq xtwxytSxyh?? ?? ? ?? ?? +? +? ???????+??+??=?? ???? (18) where y S is the specific yield and w is the recharge rate L T ? ? ? ? ? ? . After dividing by x yt?? ?? and considering 0, 0, 0x y and t?? ?? ?? we get the following 2-D depth-averaged flow equation, x? x q? ? xx q? ? +? y q? ? yy q? ? +? y? 12 y hqq Sw txy ?? ??? =? ? + ??? (19) For one dimensional flow, the above equation reduces to, y hq Sw tx ? ?? =? + ?? (20) Inserting (17) in (20) and considering homogeneous and isotropic values of hydraulic conductivity and specific yield values, we get the following one-dimensional Boussinesq equation for a sloping unconfined aquifer: cos sin y hhh SKih w txxx ???? ?? ?? =++ ???? ?? ?? ?? (21) For horizontal (i = 0) unconfined aquifer, the above equation reduces to, y hh SKh w txx ????? ?? =+ ???? ??? ?? ?? (22) 2.3 Review of previous work related to groundwater-surfacewater interactions An approximate analysis of groundwater flow resting on an impermeable bed was initiated by Dupuit (1863). The common assumptions used in that study were: 1) the stream or drain boundaries completely cut through the entire depth of the aquifer and 2) groundwater flow occurs in a homogeneous, isotropic formation with hydraulic properties that remain constant. Furthermore, these analyses assume the following 13 assumptions known as the Dupuit-Forchheimer assumptions (as summarized in Reddi, 2003): a) For small changes in the slope of line of seepage, the hydraulic head is independent of depth. b) The hydraulic gradient causing flow is equal to the slope of the water table. 2.3.1 Review of analytical solutions to the governing flow equation Since equation (21) is a non-linear equation, it does not have a general analytical solution. Therefore, the equation is often simplified and solved by introducing additional approximations. The most common approximation used is linearization. Boussinesq (1877) suggested a simple linearization method that can be used when the water table variation is small compared to the aquifer depth. There are several analytical solutions for horizontal and sloping unconfined flow. Some of them assume constant recharge rate and others assume time-varying recharge rates. Since we used a constant recharge rate in our study, here we only review analytical solutions that considered a constant recharge rate. Marino (1974) derived the following analytical equation for transient water profiles in a horizontal aquifer. 22 2 2 0 00 2()(1) (,) 1 (1)4 (1)4 22 nn wDt l n x l n x hxt h i erfc i erfc S KDt KDt SS ?? == ?? ?? ?? ?? ++? ?= ? ? ? +? ? ?? ?? ?? ?? (23) 14 He used the following initial and boundary conditions to develop the above analytical equation: Initial condition: 22 0 0hh?= at t = 0 and 0x? Boundary conditions: 22 0 0hh?= at x= 0 and 0t > 22 0 ()0hh x ? ? = ? at 2 l x= and 0t > where 0 h is the initial depth of saturation of the aquifer, h is the height of the water table above the base of the aquifer after the incidence of recharge, l is the length of the aquifer which is bounded by two fixed head boundaries, 2 ()ierfcy is the second repeated integral of the error function of argument y, h(x, t) is the hydraulic head at a particular time and space, K is the hydraulic conductivity and D is the aquifer depth, w is the recharge rate, and S is the specific yield. Melville et al. (1987) derived steady-state dimensionless analytical equations for different scenarios of hillslope flow problem. They used the following boundary conditions to get those analytical equations: Boundary conditions: 0 dh dx = (no flow) at x= 0 and 0t > 1 hh= at x= L and 0t > When F = 0.25, 1 2 2 0 ln ln 1 2 2 X HXX H X HHH H ?? ?? ??? ? =??+ ???? ??? ? ??? ? ?? ??? (24) 15 where 2 R F KT = , R is vertical infiltration or recharge rate, K is hydraulic conductivity, tanT ?= , ? is bed slope, x X L = , x is the coordinate parallel to the bed, L is the length of aquifer, 2 hC H LS = , cosC ?= , sinS ?= and h is the vertical depth of the water table. When 0.25F < () () 1 2 2 0 2114 1 ln ln ln 1 21 4 2114 X F HXXH F XHHH F F H ???? ??? ?? ?? ? ? ?? ?? ?? ?? =??+ ? ??? ?? ?? ?? ??? ?? ?? ? ??? ? ? ?+? ?? ???? (25) When 0.25F > 1 2 2 1 0 2 1 ln tan ln 1 2 41 41 X HXXH F XHHH F F H ? ? ?? ???? ? ???? ?? ? ? ?? ?? ?? ?? ?? =? ??+ ? ??? ?? ?? ?? ?? ??? ?? ?? ?? ? ??? ??? ? ? ?? ?? ?? ?? (26) Sanford et al. (1993) derived analytical equation for transient water profiles in both horizontal and sloping aquifer. They considered a system with the following initial and boundary conditions. Initial condition: h = h (0, 0) at t = 0 and 0x? Boundary conditions: h = h (0, t) at x = 0 and 0t > 0 dh dx = (no flow) at x = L and 0t > 16 For horizontal aquifer, the solution is: () 2 2(0,) 2 hLxxht ? =?+ (27) where L is the length of the aquifer, h (0, t) is water table at the boundary (where x=0), and 3 3 (0, ) i VV htL Lsw ? ??? =? ?? ?? where i V is the initial volume of water contained in the rectangular region, V is the cumulative outflow volume at any time, s is the specific yield, and w is the width of the system. For sloping aquifer, the transient water profile can be predicted by the following analytical equation: () 2 2tan(0,) 2 hLxxx ht ? ?=??+ (28) where ? is slope of the impermeable layer, and 2 3 31 (0, ) tan 2 i VV htL L Lsw ? ? ??? =?+ ?? ?? Fetter (1994) derived the steady-state analytical equation for a horizontal bed unconfined aquifer using the following boundary conditions: Boundary conditions: 1 hh= at x= 0 and 0t > 2 hh= at x= L and 0t > 17 The steady-state analytical solution is: (){} (){ } 22 12 22 1 hhx wL xx hh LK ? ? =? + (29) where h is the hydraulic head in the unconfined aquifer, w is the recharge rate, L is the length of the aquifer, K is the hydraulic conductivity, 1 h and 2 h are left and right hand side boundary heads, respectively. The above equation was used in this study to analyze steady-state data. Kim et al. (2001) derived an analytical equation for transient water profiles in horizontal unconfined aquifer using the following initial and boundary conditions: Initial condition: 0 ()hhx= at t = 0 and 0x? Boundary conditions: 1 hh= at x= 0 and 0t > 2 hh= at x L= and 0t > The derived analytical solution becomes as: 2 2 () 21 1 1 (,) sin 22 n t L n n hhx Lnx hxt x h Ce LL ? ? ? ?? ?? ? ? = ? ?? =? + + + + ?? ?? ? (30) where h(x, t) is the hydraulic head at a particular time and space, L is the length of the aquifer, 1 h and 2 h is left and right hand side boundary heads, respectively. w S ? = where w is the recharge rate, and S is the specific yield. KD S ? = where K is the hydraulic conductivity and D is the aquifer depth. 18 2 21 01 0 2 sin 22 L n hhx Lnx Chh x d LLL ?? ? ?? ??? ?? =?+?+ ? ???? ?? ?? ? Verhoest et al. (2002) also developed an analytical equation for computing the transient state profiles of both horizontal and sloping aquifers. They used the following initial and boundary conditions to derive those analytical solutions. Initial condition: h = D at t = 0 and 0x? Boundary conditions: 1 D F y = at x = 0 and 0t > 2 D F y = at x = L and 0t > where D is the depth of the soil layer, F is the Laplace transform of the water table depth, y is the Laplace variable, 11 DDH= ? , and 22 DDH= ? , 1 H is the water table at x =0, 2 H is the water table at x =L, L is the length of the aquifer. For the sloping aquifer case, the analytical equation is: ( ) ( ) ()( ) () ()() 12 1 22 12 2 22 2 2 2 1 2 22 2 2 2 22 2 2 1 2 sinh( ) 2 sinh( ) 2 21 sin exp 21 sin exp 1 aL x n aL ax n n ax aL n n eLNafKHH ax Nx hH afK aL afK nDH DHe nx n eaKt aL n L L nLNe e nx n aKt LL fK a L n ? ?? ? ? ?? ? ?? ? ? = ? ? = +??? ?? =? + + ??? ? ? ? ?? ?+ ? ??? + ?? ? ? ?? ?? ? ?? ?? +?+ ?? ??+? ? ? ? ?? ?? (31) 19 where N is the recharge rate, f is the specific yield of the formation, 2 U a k =? , an sinki U f = , coskpD i K f = , k is the hydraulic conductivity, p is the linearization constant, and i is the slope of the aquifer bed. For horizontal unconfined aquifer i=0, U = 0, and a = 0 and equation (31) can be simplified as: ()() ()( ) () 12 21 1 1 2 22 22 233 2 1 21 2 211 sin exp sin exp n n n n DH DH xH H NxL x hH LfK LN nx n nx n KtKt LL fKn LL ? ?? ?? ? ? = ? = ?? ??? ? ?? ?? =+ + + ?? ?? ?? ?? ?? ?? + ? ?? ?? ?? ?? ? ? (32) 2.3.2 Review of numerical solutions Marino (1975) presented a digital computer model that simulated the response of an unconfined aquifer to changes in stream levels. The model is based on the one- dimensional Boussinesq equation for horizontal unconfined aquifer and was solved using a predictor-corrector scheme that used non uniform grid spacing. The numerical scheme was found to be unconditionally stable. 20 Beven (1981) solved non-dimensional form of one dimensional Boussinesq equation for a sloping unconfined aquifer with constant recharge. The solution method used was the implicit finite difference approximation. Sims (1986) presented another approach for solving hillslope flow problem. He developed a numerical model, which was based on explicit finite difference scheme to solve the nonlinear problem. The major disadvantage of that model was that discretization step must be small enough to generate a stable solution. Serrano (1998) presented a numerical model for modeling transient stream/aquifer interactions in an alluvial valley aquifer. The model is based on the one-dimensional Boussinesq equation for horizontal unconfined aquifer which was solved using a decomposition method. Parlange et al. (2001) developed a numerical model based on the one-dimensional Boussinesq equation for horizontal unconfined aquifer. The equation was solved using a finite element program written in PDE2D, a general purpose partial differential equation solver. PDE2D uses a second-order accurate backward Euler scheme and an adaptive time stepping scheme to generate transient solutions. Verhoest et al. (2002) presented a numerical model and compared the results against a transient analytical solution. The numerical scheme used the Crank-Nicholson approximation and a finite element mesh to solve the one-dimensional linearized form of Boussinesq equation. The finite element scheme employed a piecewise-linear Lagrange basis function and a piecewise-uniform weighting function. 21 Rocha et al. (2007) developed a numerical model that considered two- dimensional linearized Boussinesq equation for a horizontal aquifer receiving a constant recharge. The equation was solved using the MODFLOW employing the block-centered finite difference approach. 2.3.3 Review of experimental studies Vaculin et al. (1979) conducted laboratory experiment with recharge in a three- dimensional soil slab (6 m ? 5 cm ? 2 m). The slab was packed with fine rive sand between two perspex walls supported by a frame resting on impervious horizontal boundary. Both the sides of the slab were connected to constant head reservoirs. At the soil surface, a constant flux q = 0.148 m/hr was applied over a width of 1 m in the center. The remaining soil surface was covered to prevent evaporation losses. During the infiltration, water content was measured using gamma rays, and water pressure heads at different points were measured by 20 tensiometers, each connected to its own transducer mounted on one of the perspex walls. Another transducer measured outflow volume from the constant head reservoirs. Sanford et al. (1998) conducted a series of experiments in a three-dimensional glass-sided flume (245 cm ? 30 cm ? 30 cm). The flume was filled with a sand-gravel mixture to a depth of 23 cm. Several small monitoring wells were placed in the inside wall of the flume at various locations to measure the water table heights. One side of the flume was forced with a constant head boundary (x = 0), and other side was a no flow 22 boundary. Drainage holes were installed on the downslope end of the flume (x = 0) at various heights above the impermeable layer. At the beginning of the experiment, the flume was filled with water to a saturated depth, designated as h (0, 0), when the flume was in the horizontal position. Then the upslope end was raised to the desired height and water was added at the upslope end to maintain a constant water depth as close to h (0, 0) as long as possible until steady-state condition was reached. At steady state, water table heights were recorded in various wells. Then the drainage holes at h = 0 (x = 0) were opened and water addition was stopped to produce a sudden drawdown. The water table height in each well and the cumulative outflow were recorded periodically until the drainage stopped. The same procedure was applied for several slopes. Kim et al. (2001) performed laboratory experiments using a three-dimensional physical unconfined aquifer model (67 cm ? 15 cm ? 50 cm) that used acrylic plates to visualize the water table levels from outside. Seven manometers were placed at a regular spacing in the front wall to measure water table heights. The model consisted of an unconfined aquifer region bounded by two reservoirs. The unconfined aquifer region was packed with uniform fine sand materials. Two different flow experiments were completed using two different recharge conditions. In the rising head experiment, the water table was allowed to rise by applying a constant recharge. In the falling head experiment, the water table was allowed to drain until a steady-state condition was reached. Both the rising and falling head experiments were performed for two different types of boundaries, one with equal heads, and the other with unequal heads. Transient water table heights were recorded at each manometer at a regular time intervals for all cases. 23 CHAPTER 3 EXPERIMENTAL METHODS 3.1 Column experiments A schematic column experimental setup is shown in Figure 7. The column was 3-cm diameter and 60-cm long. The column was filled with uniform glass beads (mean diameter 1.1 mm) up to 36.9 cm from the bottom of the column. A constant head reservoir was used to provide constant flow rate for the column experiment. Flow rate and water levels at the column inlet and at exit boundaries were measured under steady- state conditions. Three column experiments with three different flow rates were conducted. Table 1 provides a summary of the experimental results. 24 Figure 7: Column experiment set up in the laboratory Experiment no. 1 h (cm) 2 h (cm) Q ( 3 /mincm ) Calculated K (cm/min) Average K (cm/min) 1 49.1 44.3 88 95.68 2 47.9 44.3 69 99.36 96.07 3 55.5 44.3 200 93.19 Table 1: Results from column experiments h1 h2 Q Holding stand Porous medium Screen 25 The hydraulic conductivity values shown in Table 1 are computed based on the following analysis: According to continuity equation we know, QAq= (33) where Q= flow rate 3 L T ?? ?? ?? , A = area of the medium 2 L? ? ? ? , and q = Darcy flux L T ?? ?? ?? . From Darcy?s law, we can write h qK x ? =? ? (34) Equation (34) can be written as 12 hh qK x ? = ? (35) where K is the hydraulic conductivity, x? is the thickness of porous medium layer, 1 h and 2 h are head levels in the column at the top and bottom boundaries, respectively. By substituting (35) in (33) we get, 2 12 4 hhd QK x ? ? = ? 2 12 4 x KQ dh h? ? = ? (36) where d is the diameter of the column [ ] L . Hydraulic conductivity of glass beads can be calculated using the Hazen method (Hazen 1911). This method uses the following empirical equation, 2 10 ()KCd= (37) 26 where K is hydraulic conductivity (cm/sec), 10 d is the effective grain size (cm), and C is a coefficient based on particle size and packing of medium. For well sorted clean coarse sand, C is 120 to 150. Considering this coefficient for glass beads, we get hydraulic conductivity 1.45 cm/sec (87.12 cm/min) to 1.815 cm/sec (109.9 cm/min). Darcy velocity and Reynolds number (R) for different flow rates in column experiments were calculated. Darcy velocity was calculated using equation (35). Reynolds number was calculated using the following equation (Bear 1978), qd R ? = (38) where q = Darcy flux L T ?? ?? ?? , d = diameter of the porous medium [ ] L , ? = kinematic viscosity of fluid 2 L T ?? ?? ?? . For room temperature, kinematic viscosity of water is 2 6 10 sec m ? . The calculated values of Darcy velocities and Reynolds numbers in column experiments are tabulated in Table 2. Flow Rate, Q ( 3 /mincm ) Darcy flux (cm/min) Reynolds Number 88 12.45 2.28 69 9.76 1.79 200 28.29 5.19 Table 2: Darcy velocities and Reynolds numbers in column experiments 27 3.2 Darcy experiments Several Darcy experiments were done in a three-dimensional flow tank (115 cm ? 2.5 cm ? 60 cm), which was filled by uniform glass beads (mean diameter 1.1 mm). A constant head reservoir was used to provide constant flow rate at the left hand side of the tank for the Darcy experiment. Water levels at the left and right hand side of the tank, and outflow from the right hand side of tank were measured under steady-state conditions. Three Darcy experiments with three different flow rates were conducted. A schematic Darcy experimental setup is shown in Figure 8. Table 3 provides a summary of the experimental results. Table 4 provides the calculated values of Darcy velocities and Reynolds numbers in Darcy experiments. 28 Figure 8: Darcy experiment set up in the laboratory Experiment no. 1 h (cm) 2 h (cm) Q ( 3 /mincm ) Calculated K (cm/min) Average K (cm/min) 1 23.7 23.2 107 85.6 2 24.4 23.5 200 88.9 90 3 25.6 23.8 422 93.8 Table 3: Results from Darcy experiments 115 cm 60 cm 2.5 cm Glass tube Glass beads Wire mesh 29 Flow Rate, Q ( 3 /mincm ) Darcy flux (cm/min) Reynolds Number 107 1.98 0.363 200 3.63 0.66 422 7.43 1.36 Table 4: Darcy velocities and Reynolds numbers in Darcy experiments From Tables 2 and 4, it is seen that Reynolds numbers are less than 10. Experimental results have found that Darcy?s law is valid when the Reynolds number is less than 1 to 10 (Lindquist 1933. Rose 1945, Bear 1978). 30 CHAPTER 4 DEVELOPMENT AND TESTING OF A NUMERICAL MODEL 4.1 Numerical solution strategy A fully implicit finite difference scheme with a Picard iteration scheme for the non-linear terms was used to solve the following sloping unconfined flow equation. cos sin y hhh SKih w txxx ???? ?? ?? =++ ???? ?? ?? ?? (39) The nonlinear spatial terms in the above equation (right hand side) can be written as: () () () 11 2 11 cos sin cos 22 sin 2 mm mm ii ii ii ii ii hh Kih i w xx x hh hh Ki hh hh x Kihh w x +? +? ??? ?? ?? ++= ?? ?? ?? ? ?? ?? ?? ?? ?? ?? ?? ?? ? ?? ?? ?? ?? ?? ?? ?? ? ?? ??? ++ ?? ? ?? (40) 31 In the above equation, m is the Picard iteration level, x? is grid size, i h and 1i h + are the water level at the current (m + 1 )-th Picard iteration level at the i-th node and (i +1)- th node, respectively. The first term in the equation is approximated using the forward difference approximation, and water levels (h) are averaged using the arithmetic mean of Picard updated head values of the adjacent nodes. The second term is approximated using the central difference approximation. By using () 1 1 2 mm ii avg hh h + + = ( ) 1 2 2 mm ii avg hh h ? + = We can write equation (40) as: ()(){} () 11 2 1 11 2 cos sin cos sin 2 avg i i avg i i ii hh Kih i w xx x Kihhhhhh Kihh w xx +? +? ??? ?? ?? ++= ?? ?? ?? ? ?? ?? ?? ?? ? ??? + +?? ?? ?? ?? ?? (41) The temporal term (left hand side term) in the equation (39) can be expressed using backward-Euler approximation as: () ii yy hph h SS tt ? ? = ?? (42) 32 where i ph is the water level at the previous time level, and t? is the time step. Using finite-difference expressions for the spatial and temporal terms of equation (39), we get the following equation: () ()(){ } () 11 2 1 2 11 cos sin 2 avg i i avg i i ii y ii Kihhhhhh hph S tx Kihh w x +? +? ?? ?? ? ? =? ? ?? ?? ?? ??? ++ ?? ? ?? (43) After rearranging all the terms, we get the following equation: ( ) ( ) ( ) 2112 11avg i avg avg i avg i i hhhhhhhph? ???? ?+ ?+ + ++ +?? = + (44) where (tan) 2 x i ? ? = , 2 () (cos) y Sx Kt i ? ? = ? 2 () (cos) wx Ki ? ? = Equation (44) can be written as: 11iii ah bh ch d ?+ ++ = (45) where, 2 , avg ah ?=? + 12avg avg bh h ?= ++ 1 , avg ch ?=? ? i dph? ?= + 33 In equation (45), all the unknowns are in the left hand side and all the known are in the right hand side. Using this equation, the approximate water level for the new time level at (m +1)-th Picard level is obtained after solving the linear matrix problem. The matrix equation (45) can be written in the following form: []{} {}A hd= (46) where [ ]A is the square matrix which contains all the coefficients (except d), { }h is the unknown water levels, and { }d is based on the known water levels at the previous time level. At every Picard step, a, b, c vales are updated. Once the convergence is established the computation is moved to the next time level by updating the { }d vector. 4.2 Initial and boundary conditions For initial conditions, water tables before applying to recharge or the steady-state water table condition after applying a constant recharge rate were used. For boundary conditions, constant head, known transient head, and no flow conditions were used. Equation (45) was appropriately modified to implement these boundary conditions. For no flow condition, cos sin 0 h Kh i i x ? ?? ?+= ?? ? ?? was used at the right hand side boundary in both horizontal and sloping bed unconfined aquifers, since cos sin h qKh ii x ? ? ?? =? + ?? ? ?? . 34 4.3 Code Verification using published datasets 4.3.1 Comparison against Kim et al. (2001) horizontal aquifer solution Before using our proposed numerical algorithm in our study, we verified the code against other numerical and analytical results that are published in journal articles. First, we compared our numerical results to the analytical results of Kim et al. (2001). The parameters used by Kim et al. (2001) model are 1 h = 14.5 cm, 2 h = 14.6 cm, theta or i= 0 degree (horizontal bed), aquifer depth (D) = 16 cm, hydraulic conductivity (K) = 6.41 cm/min, specific yield ( y S ) = 0.33, recharge (w) = 1.96 cm/min, ?x = 1 cm, ?t = 0.0001 min, length of aquifer (L) = 47 cm, and total simulation time to reach the steady state was 3.3 minutes. It is important to note that the Kim et al. model used a liner approximation to simplify the model. They assumed the aquifer depth was large enough compared to the water table variations and substituted the nonlinear water table term (h) by aquifer depth (D) in one-dimensional Boussinesq equation for horizontal unconfined aquifer. We found that the numerical results matched at the beginning, but with the passage of time some deviations were observed. At the steady state, the peak water table predicted by Kim et al. (2001) is about 0.29 cm higher than that our proposed non-linear model (Figure 9). But when we modified our nonlinear code by using aquifer depth (D) as the nonlinear water table term (h) in one-dimensional Boussinesq equation for horizontal aquifer to implement the Kim et al. (2001) linear approximation, the numerical results matched well with the analytical results with very little deviation (see Figure 10). 35 10 12 14 16 18 20 22 0 5 10 15 20 25 30 35 40 45 50 Distance (cm) W a t e r Ta bl e ( c m ) t=0.3 min(model) t=0.3 min (kim) t=0.8 min (model) t = 0.8 min (kim) t=1.8 min (model) t = 1.8 min (kim) t=3.3 min (model) t = 3.3 min (kim) Figure 9: Comparison Kim et al. (2001) analytical model (linearized) results against our proposed model (non-linear) results for a horizontal unconfined aquifer example problem 36 8 10 12 14 16 18 20 22 0 5 10 15 20 25 30 35 40 45 50 Distance (cm) W a t e r Tabl e ( c m ) t=0.3 min(mod) t=0.3 min(ana) t=0.8 min (mod) t= 0.8 (ana) t=1.8 min(mod) t=1.8 min(ana) t=3.3 min(mod) t=3.3 min(ana) Figure 10: Comparison between analytical model results of Kim et al. (2001) and the linearized numerical model for a horizontal unconfined aquifer example problem 4.3.2 Comparison against Verhoest et al. (2002) sloping aquifer solution We compared our numerical model results against the results of Verhoest et al. (2002). The parameters used in Verhoest et al. (2002) model are 1 h = 0.5 m, 2 h = 1.5 m, theta or i= 2 degree, aquifer depth (D) = 2 m, hydraulic conductivity (K) = 0.001 m/s, 37 specific yield ( y S ) = 0.34, recharge (w) = 3 mm/hr, and the total simulation time to reach the steady state was 200 days. Figure 11 compares Verhoest et al. (2002) analytical results, which is based on linearization, against the linearized and non-linearized numerical results. From the figure, it can be observed that Verhoest et al. (2002) model yielded higher peak (about 0.6 m more) than the non-linear numerical model under similar conditions. 0 0.5 1 1.5 2 2.5 3 0 10203040506070809010 Distance (m) Water Tabl e (m) t=200 days(Non-linear Num. model) t=200 days(Linear num model) t=200 days (Verhoest et al(2002) analytical) Figure 11: Comparison Verhoest et al. (2002) analytical results against linear and non-linear numerical models results for a sloping unconfined aquifer example problem 38 But, when we modified our numerical code to implement the linearization procedure, the numerical model results matched well with Verhoest et al. results (2002) (see Figure 11). 4.3.3 Comparison against Sanford et al. (1993) solution We also compared our numerical model results against the analytical results of Sanford et al. (1993) for three different cases. In Sanford et al. (1993) the right hand side was set to no-flow boundary condition and the left hand side was set to transient head boundary condition. The parameters used in Sanford et al. (1993) analytical model are shown in Table 5, and the comparison of Sanford et al. (1993) analytical results and our proposed model results are shown in Figures 12, 13, and 14. It is to be noted here that the analytical model of Sanford et al. (1993) was developed on equations 22 and 23, and the left side boundary of the experimental tank was suddenly dropped removing the cork from the bottom hole (h = 0) only. Case Slope (degree) h(0,0) cm K (cm/s) y S 1 0 22.3 0.76 0.238 2 0.089 22.2 0.78 0.238 3 0.137 21.8 0.64 0.238 Table 5: Parameters used in Sanford et al. (1993) analytical model 39 0 5 10 15 20 25 0 50 100 150 200 Distance (cm) H e i ght ( c m ) IC(O min) Sanford (1.6 min) proposed (1.6 min) sanford(5.5 min) proposed (5.5 min) sanford(12 min) proposed(12 min) sanford (33 min) proposed(33 min) Figure 12: Comparison Sanford et al. (1993) analytical model results against our proposed model (non-linear) results for a horizontal unconfined aquifer example problem 40 0 5 10 15 20 25 0 50 100 150 200 Distance (cm) H e i ght ( c m ) IC(O min) Sanford (0.67 min) proposed (0.67 min) sanford(1.5 min) proposed (1.5 min) sanford(3 min) proposed(3 min) sanford (5 min) proposed(5 min) Figure 13: Comparison Sanford et al. (1993) analytical model results against our proposed model (non-linear) results for a sloping unconfined aquifer (theta = 0.089) example problem 41 0 5 10 15 20 25 0 50 100 150 200 Distance (cm) Heig ht (cm) IC(O min) Sanford (0.33 min) proposed (0.33 min) sanford(2 min) proposed (2 min) sanford(3 min) proposed(3 min) sanford (5 min) proposed(5 min) Figure 14: Comparison Sanford et al. (1993) analytical model results against our proposed model (non-linear) results for a sloping unconfined aquifer (theta = 0.137) example problem From the Figure 12, we found Sanford et al. (1993) analytical results in horizontal aquifer matched well against our proposed numerical model results. For the sloping aquifer, the results matched quite well, except near the no-flow boundary (see Figures 13 and 14). 42 CHAPTER 5 LABORATORY DATA AND MODELING RESULTS 5.1 Steady state experiments A three-dimensional flow tank (115 cm ? 2.5 cm ? 60 cm) which was connected to a constant head reservoir and a constant recharge generator was used in this study. The flow tank was made of plexi glass plates hence the water level can be recorded from outside the tank. Uniform glass beads (mean diameter 1.1 mm) were used as the porous medium to construct an unconfined aquifer. Multiple 6 mm diameter glass tubes filled with red-colored water were used to delineate the location of the water table in the aquifer. Digital photos of the actual tank fitted with recharge generator in both horizontal and sloping bed unconfined aquifers are shown in Figures 15 and 16 respectively. A schematic of the laboratory setup is shown in Figure 17. Experimental works conducted with water flow from the constant head reservoir to laboratory tank through a recharge generator. In all the steady-state experiments, constant head boundary condition was used for the left side of the tank, and no flow 43 boundary condition was used for the right side of the tank. During infiltration, water table in unconfined aquifer was raised and when it became steady-state, the outflow from the right hand side of the tank was collected using a beaker for 2 minutes. The flow was measured by dividing the collected outflow volume by 2 minutes. The sloping aquifer was made by lifting the left hand side of the tank placing a wooden block of 3.8 cm height below the bottom. Observations from multiple steady-state experiments conducted with different recharge rates under horizontal and sloping aquifers are shown in Figures 18 and 19 respectively. Recharge rates were designed in such a way that there was no ponding. The recharge rates and constant head boundary conditions for both horizontal and sloping aquifers are tabulated in Table 6. Table 6: Recharge rates and constant head boundary conditions for both horizontal and sloping unconfined aquifers Slope (degree) Constant Boundary condition (cm) Recharge rate (cm/min) 0 22.0 2.9 0 22.75 4.2 0 25.0 5.4 2.03 23.8 4.42 2.03 24.3 4.83 2.03 26.0 5.55 44 Figure 15: Laboratory tank with recharge generator used in this study for horizontal unconfined aquifer Figure 16: Laboratory tank with recharge generator used in this study for sloping unconfined aquifer (i =2.03 degree) 45 Figure 17: A schematic laboratory set up for tank experiment 115 cm 60 cm 2.5 cm 107cm No flow boundary Glass tubes Glass beads Wire mesh Constant head boundary Recharge 46 0 5 10 15 20 25 30 35 40 2040608010120 Distance (cm) W a t e r Tabl e Hei ght , h ( c m ) exp(w=2.9 cm/min) expt(w=4.2 cm/min) expt(w=5.4 cm/min) No Flow Boundary Figure 18: Experimental steady-state water table profiles in horizontal unconfined aquifer under different recharge rates. 47 0 5 10 15 20 25 30 35 40 0 20 40 60 80 100 120 Distance (cm) Water Tabl e Height, h (cm) expt (w=5.55 cm/min) expt(w=4.83 cm/min) expt (w=4.42 cm/min) No flow Boundary Figure 19: Experimental steady-state water table profiles in sloping unconfined aquifer (i = 2.03 degree) under different recharge rates. 48 5.2 Modeling steady-state water table profiles Steady-state water table profiles observed were modeled using analytical and numerical models. For the horizontal aquifer data, we used the analytical equation (29), which was derived by Fetter (1994). The hydraulic conductivity estimated from Darcy experiments was used in analytical equation (29), and the resulting analytical model results matched well with the experimental steady-state results for different recharge rates (see Figure 20). For the sloping aquifer, steady-state water table profiles were simulated using the numerical model. The same numerical modeling results matched reasonably well with the experimental results (see Figure 21). 49 0 5 10 15 20 25 30 35 40 0 2040608010120 Distance (cm) W a t e r Ta bl e H e i ght , h ( c m ) expt(w=4.2 cm/min) ana(w=4.2 cm/min) expt(w=5.4 cm/min) ana (w=5.4 cm/min) exp(w=2.9 cm/min) ana(w=2.9 cm/min) No flow Boundary Figure 20: Modeling steady-state water table profiles using analytical equation in horizontal unconfined aquifer 50 0 5 10 15 20 25 30 35 40 0 2040608010120 Distance (cm) W a t e r Ta bl e H e i ght , h ( c m ) expt (w=5.55 cm/min) model(w=5.55 cm/min) expt(w=4.83 cm/min) model(w=4.83 cm/min) expt (w=4.42 cm/min) model (w=4.42 cm/min) No Flow Boundary Figure 21: Modeling steady-state water table profiles using numerical model in sloping unconfined aquifer (i = 2.03 degree) 51 5.3 Transient experiments Multiple transient experiments were conducted in both horizontal and sloping aquifers. For transient experiments, at first water tables was raised to a steady-state level by applying a constant recharge rate, and then the system was allowed to drain after removing the recharge. During transient experiments, water table at the left boundary changed with time. Therefore, we had to model the constant head boundary condition as a known transient head boundary. All the transient experiments were simulated using known transient head boundary condition for the left side of the tank, and no flow boundary condition for the right side of the tank. The specific yield ( y S ) was estimated by fitting the numerical model results to the transient experimental results. Water table readings were recorded at regular time intervals using a digital camera. All the experimental data sets and the modeling results are summarized in the following section. 5.3.1 Modeling the transient head boundary at the left side of the tank Transient head boundary for both horizontal and sloping aquifers was described using a best fit curve of the observed transient left hand boundary data (Figures 22 and 23). These data were collected by sampling the digital images at every half-minute from an initial steady-state condition to fully drained condition. Multiple steady-state experiments were done using different recharge rates in horizontal aquifer. Here only one case with a constant recharge of 5.4 cm/min is shown as an example. At first the observed 52 transient left hand boundary data were plotted against time, and then a best-fit curve was drawn for that profile of the experimental result. The best fit curve is given by equation (39), which was used to model the left hand side transient head boundary condition (i.e. x = 0 cm) in horizontal aquifer. Figure 22 shows the best fit model results against observed head data. 432 ( ) 0.0061 0.1356 1.1147 4.0671 25.348ht t t t t=?+?+ (47) h = 0.0061t 4 - 0.1365t 3 + 1.1147t 2 - 4.0671t + 25.348 R 2 = 0.99 18 19 20 21 22 23 24 25 26 0246810 Time (min) Left boundar y w a ter l e vel (cm ) h vs t (LHS) Poly. (h vs t (LHS)) Figure 22: Modeling transient head boundary (left hand side of the tank) in horizontal unconfined aquifer. 53 Similarly, the equation for transient head boundary condition at the left hand side in sloping aquifer was developed. Figure 23 also shows the best fit model results with observed head data. The equation of that best fit curve for recharge rate of 4.42 cm/min is: 432 ( ) 0.0713 0.6172 1.5002 0.3683 23.875ht t t t t=? + ? + + (48) h = -0.0713t 4 + 0.6172t 3 - 1.5002t 2 - 0.3634t + 23.875 R 2 = 0.9921 18 19 20 21 22 23 24 25 0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Time (min) L e f t b oun dar y w a t e r l evel ( c m ) h vs t (LHS) Poly. (h vs t (LHS)) Figure 23: Modeling transient head boundary (left hand side of the tank) in sloping unconfined aquifer (i = 2.03 degree) 54 5.3.2 Horizontal unconfined aquifer- transient experimental data and modeling results Figure 24 shows the transient experimental data and numerical results when the water table was allowed to drain from the steady-state condition (initial condition) with a constant recharge rate of 5.4 cm/min in the horizontal unconfined aquifer system. The specific yield ( y S ) was estimated as 0.35 by fitting these numerical model results to the transient experimental results. Then this estimated value was used in the numerical model to fit the numerical results with another experimental results (see Figure 25 where recharge rate= 4.42 cm/min at the steady-state). The numerical model results show good match with the experimental results at different time intervals. The results show that the water table drops fast at the beginning because of high hydraulic gradient, but drops slowly at later stage due to low hydraulic gradient. Kim et al. (2001) made similar observation. This indicates a non-linear drainage pattern for water-table movement in the unconfined aquifer. 55 10 15 20 25 30 35 40 0 2040608010120 Distance (cm) W a t e r Tabl e Hei g ht , h ( c m ) mod(st) exp(st) mod(0.25 min) exp(0.25 min) mod(1 min) exp(1 min) mod(2 min) exp(2 min) mod(4 min) exp(4 min) mod( 8 min) exp(8 min) No flow Boundary Figure 24: Transient behavior of water table in horizontal unconfined aquifer (k = 90 cm/min (1296 m/day), y S = 0.35 and w = 5.4 cm/min at the beginning (steady state)) 56 10 15 20 25 30 35 40 0 2040608010120 Distance (cm) W a t e r Ta bl e H e i g ht , h ( c m ) mod(st) exp(st) mod(0.25 min) exp(0.25 min) mod(1 min) exp(1 min) mod(2 min) exp(2 min) mod(5 min) exp(5 min) No flow Boundary Figure 25: Transient behavior of water table in horizontal unconfined aquifer (k = 90 cm/min (1296 m/day), y S = 0.35 and w = 4.42 cm/min at the beginning (steady state)) 57 5.3.3 Sloping unconfined aquifer- transient experimental data and modeling results Figure 26 also shows the transient experimental data and numerical results when the water table was allowed to drain from the steady-state condition (initial condition) with a constant recharge rate of 4.42 cm/min in the sloping unconfined aquifer (i = 2.03 degree). 0 5 10 15 20 25 30 35 2040608010120 Distance (cm) Water Table Height, h (cm) exp(st) mod(st) exp(0.5 min) mod(0.5 min) exp(1 min) mod(1 min) exp(3.5 min) mod(3.5 min) No flow Boundary Figure 26: Transient behavior of water table in sloping unconfined aquifer ((i = 2.03 degree, K = 90 cm/min (1296 m/day), y S = 0.35 and w = 4.42 cm/min at the beginning (steady state)) 58 The numerical model results also show good agreement with the experimental results. The experimental result shows that the water profile dropped very rapidly because of sloping effect. It is observed that the water table dropped fast at the beginning and slowly at the later stage, similar to the horizontal aquifer case. This indicates a non- linear pattern in the water-table movement in unconfined aquifer. 5.3.4 Mass balance accuracy of the Numerical results Our numerical model solved a set of non-linear equations using the Picard iterative scheme. The model employed 1x? = cm, and 0.00025t? = min to simulate the transient profiles. The mass balance errors were less than 1%, which are tabulated in Tables 7 and 8 for horizontal and sloping unconfined aquifers respectively. Elapsed Time (min) Mass balance error (%) 0.25 0.46 0.5 0.5 1 0.57 2 0.65 4 0.66 8 0.68 Table 7: Mass balance errors for the water table profiles generated by the numerical model for the horizontal unconfined aquifer 59 Elapsed Time (min) Mass Balance Error (%) 0.5 0.57 1 0.6 3.5 0.6 Table 8: Mass balance errors for the water table profiles generated by the numerical model for the sloping unconfined aquifer (i = 2.03 degree) 5.4 Correlation between the sloping unconfined aquifer with the horizontal unconfined aquifer Figure 27 compares the steady-state head profiles observed in the horizontal and sloping unconfined aquifers under similar conditions. The figure shows that the peak head values observed in horizontal and sloping (2.03 degree slope) unconfined aquifer occurred at the no Flow boundary. The peak in the horizontal unconfined aquifer is approximately 1.9 cm higher than that in the sloping unconfined aquifer under identical recharge rate of 4.42 cm/min. But the total draining time in sloping unconfined aquifer is approximately 1.5 minute less than that in horizontal unconfined aquifer because of the influence of the gravity effect in sloping aquifer. 60 10 15 20 25 30 35 40 0 2040608010120 Distance (cm) W a t e r Tabl e H e i ght , h ( c m ) Steady state (2.03 degree slope) Steady state (Horizontal) No flow Boundary Figure 27: Comparison of experimental steady state water table profile between horizontal and sloping unconfined aquifers (i = 2.03 degree) under same condition 61 CHAPTER 6 SUMMARY, CONCLUSIONS AND RECOMMENDATIONS We developed a laboratory-scale aquifer model to conduct steady-state and transient experiments in both horizontal and sloping unconfined aquifers. The experimental results were recreated using a numerical model. A three-dimensional flow tank (115 cm ? 2.5 cm ? 60 cm) connected to a constant head reservoir and a constant recharge generator was used to conduct steady- state and transient experiments in both horizontal and sloping unconfined aquifers. The hydraulics conductivity of the porous medium was measured using 1-D column and 2-D Darcy experiments. The average hydraulic conductivity was estimated to be 90 cm/min. At first, multiple steady-state experiments were conducted for different recharge rates. For steady-state experiments, constant head boundary at the left hand side and no flow boundary at the right hand side of the tank were used. Experimentally measured steady-state water table profiles were modeled using analytical and numerical models. Horizontal aquifer data was described using an analytical model. For sloping aquifer, a numerical model was developed to simulate the experimental results. During transient experiments, it was observed that the water table at the left hand boundary changed with time. Therefore, the left hand 62 boundary was modeled as a known transient head boundary. The transient head boundary variations observed for both horizontal and sloping aquifers were modeled by using a best fit curve of the observed head data. This data was collected by sampling the digital images at every half-minute from an initial steady-state condition to fully drained condition. All the transient experiments were simulated using transient head boundary condition for the left boundary, and no flow condition for the right boundary. The transient experiments and numerical simulations were done for known transient head boundary condition. The transient results show that water table drops faster at the beginning than that in the later stage because of high hydraulic gradient at the beginning. These results indicate the non-linear nature of water-table drainage patterns in unconfined aquifer system. It also shows that the peak water table in the horizontal unconfined aquifer is higher than that in the sloping unconfined aquifer (2.03 degree slope) under the same recharge condition, hence more time required for draining in horizontal unconfined aquifer. The novel experimental setup developed in our laboratory also allowed us to understand the interactions among an unconfined aquifer, a stream, and a water divide during or after rainfall event, when the water level in the near by stream varies. The numerical model predicted better groundwater table data than that in linearized analytical model. The proposed numerical model is a useful tool for predicting the variations in groundwater table during or after rainfall event. For future application, the model needs to be further validated using a field scale stream-aquifer interaction dataset. However, there are some challenges in the field for characterizing the hydraulic properties (e.g. hydraulic conductivity). In the field, 63 hydraulic conductivity variations will be non-homogenous and anisotropic. This would require a careful analysis of various averaging procedures used for interpretation the field data. 64 REFERENCES Bear J. 1978. Hydraulics of groundwater. McGraw-Hill , Inc. NY. Beven K. 1981. Kinematic subsurface strormflow. Water Resources Research. 17(5): 1419-1424. Boussinesq MJ. 1877. Essai sue la theorie des eaux courants. Memories de L? Academic des Sciences de L? Institute de France 23(1): 252-260. Brutsaert W. 1994. The unit response of groundwater outflow from a hillslope. Water Resources Research 30(10): 2759-2763 Chow VT, Maidment DR, and Mays LW. 1988. Applied Hydrology, McGraw-Hill, New York. Dupuit J. 1863. Etudes theoriques et pratiques sur le movement des eaux dans les canaux decouverts et a travers les terrains permeables, ed. 2, Dunod, Paris. Fetter CW. 1994. Applied Hydrogeology. Prentice Hall: Englewood Cliffs, NJ; 691. Hall FR. 1968. Base-flow recession- A review, Water Resources Research, v. 4, pp. 973- 983. Hazen A. 1911. Discussion: Dams on sand foundations. Transactions, American Society of Civil Engineers 73:199. 65 Kim DJ, and Ann MJ. 2001. Analytical solutions of water table variation in a horizontal unconfined aquifer: Constant recharges and bounded by parallel streams. Hydrological Processes 15: 2691-2699. Lindquist E. 1933. On the flow of water through porous soil. Premier Congres des grands barrages (Stockholm), 81-101. Marino MA. 1974. Rise and decline of the water table induced by vertical recharge. Journal of Hydrology. 23: 289-298. Marino MA. 1975. Digital simulation model of aquifer response to stream stage fluctuation. Journal of Hydrology. 25: 51- 58. Melville JG, and PN Sims. 1987. An analysis of hillside seepage. Ground Water. V. 25, Number 6. Meyboom P. 1961. Estimating ground-water recharge from stream hydrographs. Journal of Geophysical Research. V. 66, pp. 1203-1214. Norum DI, and Luthin JN. 1968. The effects of entrapped air and barometric fluctuations in the drainage of porous mediums. Water Resources Research. V. 4, pp. 417- 424. Parlange JY, Stagnitti F, Heilig A, Szilgyi J, Parlange MB, Steenhuis TS, Hogarth WL, Barry DA, and Li L. 2001. Sudden drawdown and drainage of a horizontal aquifer. Water Resources Research. 37(8):2097-2101. Rocha D, Feyen J, and Dassargues A. 2007. Comparative analysis between analytical approximations and numerical solutions describing recession flow in unconfined hillslope aquifers. Hydrogeology Journal 15:1077-1091. 66 Rose HE. 1945. An investigation into the laws of flow of fluids through beds of granular materials. Proceedings of the Institute of Mechanical Engineers 153:141-148. Reddi LN. 2003. Seepage in soils: principles and applications. John Wiley and Sons. New York. Sanford WE, Parlange JY, and Steenhuis TS. 1993. Hillslope drainage with sudden drawdown: closed form solution and laboratory experiements. Water Resources Research. 29(7):2313-2321. Schneider R. 1961. Correlations of Groundwater Levels and Air Temperatures in the Winter and Spring in Minnesota, U.S. Geological Survey Water-Supply Paper1539- D, pp.14. Serrano SE, and Workman SR. 1998. Modeling transient stream/aquifer interaction with the non-linear Boussinesq equation and its analytical solution. Journal of Hydrology. 206: 245- 255. Shanley JB, Kendall C, Smith TE, Wolock DM, and McDonnell JJ. 2002. Controls on old and new water contributions to stream flow at some nested catchments in Vermont, USA. Hydrological Processes 16: 589-609. Sims PN. 1986. Hillside seepage modeling. Master?s Thesis. Auburn University. Troch PA, Paniconi C, and Loon EV. 2003. Hillslope-stroage Boussinesq model for subsurface flow and variable source areas along complex hillslopes: 1. Formulation and characteristic response. Water Resources Research. (11), 1316, doi: 10.1029/2002WR001728. 67 Vaculin M, Khanji D, and Vachaud G. 1979. Experimental and numerical study of a transient, two-dimensional unsaturated-saturated water table recharge problem. Water Resources Research. 15(5), 1089-1101. Verhoest NEC, Pauwels VRN, Troch PA, and Troch FPD. 2002. Analytical solution for transient water table heights and outflows from inclined ditch-drained terrains. Journal of Irrigation and Drainage Engineering 128(6):358-36 Werner PW, and Sundquist KJ. 1951. On the ground-water recession curve for large watersheds, International Association of Scientific Hydrology Publication, 33 pp. 202-212. 68 APPENDIX An Alternate Analysis of Depth-averaged Darcy Flux in Sloping Bed Unconfined Aquifer: Figure 28: Pressure at the bottom of a sloping bed unconfined aquifer In Figure 28, the weight of fluid phase is f wnghxy?=?? (I)where ? is the density of water 3 M L ? ? ? ? ? ? , g is gravitational acceleration 2 L T ?? ?? ?? , h is i cos f wi f w x? h y? 69 height of water table perpendicular to the sloping bed, n is porosity, x? is the length [L] along the aquifer bottom, and y? is the width [L] of the aquifer. Component of this weight perpendicular to the sloping bed, cos f wi, is equal to cosnghxy i? ?? Pressure acting on perpendicular to the sloping bed, cos cos ff wiwi F P A Axy == = ? ? (II) where F is the force acting perpendicular to the sloping bed 2 ML T ? ? ? ? ? ? , and Anxy=?? is the area normal to the sloping bed 2 L? ? ? ? . After substituting equation (II) in equation (III), we get cos cos nghxy i P gh i nxy ? ? ?? == ?? (IV) Total head is defined as: 2 2 pv H z g? =+ + (V) where p ? is pressure head, 2 2 v g is velocity head, and z is elevation head. Since the groundwater velocity (v) is typically a very small, we can neglect the velocity term. Hence, p Hz ? =+ (VI) Substituting (IV) in (VI), we get 70 cosgh i Hz g ? ? =+ cosH hiz=+ (VII) From Darcy?s law, H vK x ? =? ? (VIII) Where v is Darcy flux or specific discharge L T ? ? ? ? ? ? , K is the hydraulic conductivity L T ? ? ? ? ? ? , and H x ? ? is hydraulic gradient. Depth-averaged Darcy flux can be defined as: qvh ? = (IX) where q ? is depth averaged flux along the x direction 2 L T ? ? ? ? ? ? . After substituting equations (VII) and (VIII) in equation (IX) we get the general depth averaged flux equation for unconfined flow in a sloping aquifer as: (cos )qKhhiz x ? ? =? + ? cos hz qKh i x x ? ?? ?? =? + ?? ?? ?? cos sin h qKh ii x ? ??? =? + ?? ? ?? (X)