Citation |

- Permanent Link:
- http://ufdc.ufl.edu/AA00008997/00001
## Material Information- Title:
- Two-dimensional analytical and three-dimensional finite-element method modeling of the interactions between wetlands and groundwater
- Creator:
- Zhong, Jinhua
- Publication Date:
- 2004
- Language:
- English
- Physical Description:
- xvi, 126 leaves : ill. ; 29 cm.
## Subjects- Subjects / Keywords:
- Aquifers ( jstor )
Groundwater ( jstor ) Groundwater flow ( jstor ) Modeling ( jstor ) Parametric models ( jstor ) Peat ( jstor ) Saturated flow ( jstor ) Three dimensional modeling ( jstor ) Water tables ( jstor ) Wetlands ( jstor ) Civil and Coastal Engineering thesis, Ph. D ( lcsh ) Dissertations, Academic -- Civil and Coastal Engineering -- UF ( lcsh ) - Genre:
- bibliography ( marcgt )
theses ( marcgt ) non-fiction ( marcgt )
## Notes- Thesis:
- Thesis (Ph. D.)--University of Florida, 2004.
- Bibliography:
- Includes bibliographical references.
- General Note:
- Printout.
- General Note:
- Vita.
- Statement of Responsibility:
- by Jinhua Zhong.
## Record Information- Source Institution:
- University of Florida
- Holding Location:
- University of Florida
- Rights Management:
- Copyright [name of dissertation author]. Permission granted to the University of Florida to digitize, archive and distribute this item for non-profit research and educational purposes. Any reuse of this item in excess of fair use or other copyright exemptions requires permission of the copyright holder.
- Resource Identifier:
- 003163679 ( ALEPH )
779971760 ( OCLC )
## UFDC Membership |

Downloads |

## This item has the following downloads: |

Full Text |

TWO-DIMENSIONAL ANALYTICAL AND THREE-DIMENSIONAL FINITE-ELEMENT METHOD MODELING OF THE INTERACTIONS BETWEEN WETLANDS AND GROUNDWATER By JINHUA ZHONG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2004 TABLE OF CONTENTS ACKNOWLEDGMENT ................................................. ii LIST OF TABLES ....................................................... v LIST OF FIGURES ..................................................... vi LIST OF SYM BOLS .................................................... ix A B STR A CT .......................................................... xiii CHAPTER ........................................................... . 1 INTRODUCTION ...................................................... 1 1.1 W etland ...................................................... 1 1.2 Wetland Water Budget ........................................... 2 1.3 Goals and Objectives ............................................ 4 2 ANALYTICAL GROUNDWATER MODEL RESPONSE TO THE PERIODIC OR NONPERIODIC WATER STAGE FLUCTUATION IN A CIRCULAR LAKE OR WETLAND ............................................. 7 2.1 Introduction ................................................... 7 2.2 Problem Formulation: Groundwater Flow Around a Circular Wetland ..... 8 2.3 Analytical Solution for the Periodic Boundary Condition ................ 9 2.4 Analytical Solution for the Non-periodic Boundary Condition ........... 12 2.5 Application .................................................. 14 2.6 Comparison of Analytical Solution and Numerical Solutions ............ 17 2.7 Sum m ary .................................................... 18 3 TRANSIENT THREE-DIMENSIONAL DEFORMABLE FINITE-ELEMENT SATURATED GROUNDWATER MODEL .......................... 19 3.1 Introduction .................................................. 19 3.2 Groundwater Flow Models Based on Finite-Difference and Finite-Element Approximation .................... .. ..... 21 3.3 Three-Dimensional Deformable Finite-Element Saturated Groundwater Model ................................. 24 3.3.1 Governing Equation and Boundary Conditions ............... 27 3.3.2 Finite Element Model of Saturated Groundwater Flow ......... 28 3.3.3 Adaptation of the Saturated Finite-Element Model to Describe the Free Surface .............................. 36 3.4 Validating the Code of the Three-Dimensional Deformable Finite-Element Groundwater Flow Model ....................... 37 3.4.1 Theis Solution ........................................ 37 3.4.2 Groundwater Flow around a Circular Wetland ............... 41 3.5 A Numerical Model Application to Simulate the Interaction Between a W etland and an Aquifer .................................... 42 3.6 Application of Three-Dimensional Deformable Saturated Groundwater Flow Model to SV5 Wetland ...................... 50 3.6.1 Site Description-SV5 ................................... 50 3.6.2 Field Methods and Wetland Aquifer Interaction Test ........... 52 3.6.3 Stage-Volume Relationships .............................. 54 3.6.4 Creation of Finite-Element Model Mesh .................... 56 3.6.5 Solution Procedure ..................................... 57 3.6.6 Model Calibration and Results ..................... .... 58 3.7 Comparison with Results of Wise and others ........................ 69 3.8 Sum m ary .................................................... 71 4 INVERSE MODEL OF THREE-DIMENSIONAL FINITE-ELEMENT METHOD SATURATED GROUNDWATER MODEL IN SEARCHING FOR THE HYDRAULIC CONDUCTIVITY OF THE PEAT LAYER ....... 69 4.1 Introduction .................................................. 69 4.2 Adjoint Problem for Three-Dimensional Finite-Element Saturated Groundwater Model ............................... 73 4.3 Searching Objective Function .................................... 80 4.4 Adjoint State Controlling Equations ............................... 80 4.5 Searching Method and Stop Criterion .............................. 82 4.6 Field Application on SV5 ....................................... 82 4.8 Summary ................................ .................... 85 5 THREE-DIMENSIONAL FINITE-ELEMENT VARIABLY-SATURATED GROUNDWATER MODEL ....................................... 87 5.1 Introduction ................................ .................. 87 5.2 Three-Dimensional Finite-Element Variably-Saturated Groundwater Model ......................................... 88 5.2.1 Governing Equation .................................... 88 5.2.2 Linearizing the Governing Equations ....................... 89 5.2.3 Boundary Conditions for the Variably-Saturated Model ........ 91 5.3 Finite-Element Method for Modeling Three-Dimensional Variably-Saturated Flow ..................................... 94 5.4 Transient, Variably-Saturated Water-Table Recharge Example Test ...... 98 5.5 Application on SV5 W etland .................................... 101 5.5.1 Finite-Element Mesh for Three-Dimensional Variably-Saturated Flow Model ......................... 101 5.5.2 Modeling Results and Summary .......................... 103 6 SUMMARY AND CONCLUSIONS ..................................... 111 REFERENCES ....................................................... 113 BIOGRAPHICAL SKETCH ............................................. 126 LIST OF TABLES Figure page 2-1 Parameters and their values ............................................ 17 3-1 Selected early references for the finite-difference method .................... 22 3-2 Selected early references for the finite-element application to model groundwater 25 3-3 Parameters and their values .................................. .......... 42 3-4 Data for calculating the discharge to the wetland from the aquifer at 10 hours .... 63 3-5 Data for calculating the discharge to the wetland from the aquifer at 20 hours .... 64 3-6 Data for calculating the discharge to the wetland from the aquifer at 30 hours .... 64 3-7 Data for calculating the discharge to the wetland from the aquifer at 40 hours .... 65 3-8 Data for calculating the discharge to the wetland from the aquifer at 50 hours .... 65 3-9 Sensitivities calculated by saturated groundwater flow model ................. 66 4-1 Error and sensitivities based on the adjoint method model .................... 84 5-1 Data for calculating the discharge to the wetland from the aquifer at 10 hours ... 107 5-2 Data for calculating the discharge to the wetland from the aquifer at 20 hours ... 108 5-3 Data for calculating the discharge to the wetland from the aquifer at 30 hours ... 108 5-4 Data for calculating the discharge to the wetland from the aquifer at 40 hours ... 109 5-5 Data for calculating the discharge to the wetland from the aquifer at 50 hours ... 109 v 5-6 Sensitivities calculated by variably saturated groundwater flow model ......... 110 LIST OF FIGURES Figure page 2-1 Setup of the groundwater problem around a circular wetland ................... 9 2-2 Transient groundwater level changes at radial distances 1,2,3 and 4 meters away from the wetland shoreline ..................................... 15 2-3 Groundwater distribution over the distance to the wetland shoreline at 6,7,8,9,10 h .................................................... 15 2-4 Three-dimensional view of the groundwater table around a circular wetland at t=3hr .......................................... 16 2-5 Comparison of the analytical solution and numerical solution at 2 and 5 meters from the wetland shore line ........................... 18 3-1 Conceptual profile model of a wetland/groundwater system .................. 26 3-2 Comparison of finite element mesh and finite difference mesh ................ 26 3-3 Setup of the saturated groundwater problem around a wetland ................. 27 3-4 Flowchart illustrating the iteration procedure of finding the coupling boundary condition ..................................... 29 3-5 Unsteady flow to a well in phreatic aquifer at time t ......................... 38 3-6 Finite-Element Method mesh of the wedge used to simulate of a well pumping in phreatic aquifer ................................. 40 3-7 Comparison of the numerical and analytical head for the example problem of a transient pumping test in a phreatic aquifer ......................... 40 3-8 Comparison of the numerical and analytical simulated groundwater table ........ 41 vii 3-9 Three-dimensional finite-element mesh for the simulation of wetland and aquifer 43 3-10 Relationship between the area of wetland water surface ..................... 44 3-11 Relationship between the deepest water depth and the volume of wetland ....... 45 3-12 Illustration of the initial condition where it was assumed a wetland surface-water stage and the phreatic surface elevation were the same ........ 45 3-13 Plan view of the numerical mesh at time zero. ............................ 46 3-14 Three-Dimensional Finite-Element mesh and simulated head contour at hour 5 .. 46 3-15 Simulated head contour at hour 5 ...................................... 47 3-16 Simulated iso-surface plot of 10.6 meter head at hour 5 ..................... 47 3-17 Three-dimensional Finite-Element mesh and head contour on a cross section on a cross section through wetland ..................... 48 3-18 Simulated head contour on the vertical cross section through wetland .......... 48 3-19 Modeled iso-surface head plot of 11.1 m,l1 m and 10.5 m ................... 49 3-20 Mesh and modeled head contour with vertical cross section through the pumping well .......................................... 49 3-21 Head contour (m) on the vertical cross section through the pumping well ....... 50 3-22 W ell location map (Switt et al. 1998) ................................... 51 3-23 Top surface of peat layer ............................................. 54 3-24 Bottom surface of peat layer .......................................... 55 3-25 Relationship between wetland volume and water stage ..................... 55 3-26 Finite-element mesh of wedge of wetland and aquifer for the saturated flow model ........................................ 57 3-27 Vertical cross-section view of three-dimensional finite-element mesh and head from the saturated flow model ............................... 61 3-28 Piezometric head at positions beneath wetland with time .................... 61 3-29 Variation in the inundated area of wetland during simulation ................ 62 3-30 Discharge between wetland and aquifer with respect to time ................. 62 3-31 Comparison between calculated head and observed head .................... 63 3-32 Comparison results using Walser's parameter values ....................... 67 3-33 3-D numerical model simulated results over specific storage ................. 67 4-1 Sketch showing setup of the saturated groundwater inverse problem around a wetland .................................... 74 4-2 Relation between error and peat layer hydraulic conductivity .................. 84 4-3 Flowchart illustrating the iteration procedure of inverse problem .............. 85 4-4 Vertical cross-section of the adjoint state distribution at T = 10 hours ........... 86 5-1 Sketch showing setup of the variably saturated groundwater flow problem around a wetland ...................................... 92 5-2 Schematic of relative evapotranspiration (ET/PET), as affected by soil water potential,j .................................. 93 5-3 Schematic diagram of the flow domain ................................... 99 5-4 Water moisture distribution (x:y=2:1) at 4 hours .......................... 100 5-5 Water moisture distribution at 8 hours (x:y=2:1) .......................... 101 5-6 Finite-element mesh of the wedge part for variably saturated flow model (Unite is m) ............................ 102 5-7 Calculated piezometric head variation beneath wetland using the three-dimensional variably saturated groundwater flow model .......... 104 5-8 Model predicted discharge from the aquifer to the overlaying wetland ......... 104 5-9 Inundated area of wetland during simulation ............................. 105 5-10 Simulated and observed heads in wetland and in aquifer ................... 105 5-11 Simulated head from the saturated flow model and the variably saturated flow model ............................................. 106 5-12 Horizontal view of the three-dimensional finite-element mesh and piezometric head distribution from the variably saturated flow model at time of 9 hours .. 106 5-12 Enlarged top and middle part of the mesh and head contour at time of 9 hours .. 107 LIST OF SYMBOLS C = the conductance of the wetland confining unit Co = coefficient di(t) = the unknown values of hydraulic head and at time t ET = evapotranspiration [ULT] F(t) = the Fourier coefficient G = ground-water level [L] Gi = groundwater inflow [L3/T] Go = groundwater outflow [L3 /T] GWNs = the net flow or seepage of groundwater to the wetland surface [L3fT] h = piezometric head [L] h(e) = the approximate solution for hydraulic head within element e [L] Ho = initial saturated aquifer thickness [L] H' = Sobolev space Ho(t) = the Hankel function order 0 hoo = the regional average water table elevation, [L] i = reference integer spatial interval Jo = the Bessel function of the first kind order 0 k = a constant parameter K = the saturated hydraulic conductivity tensor [ IT] Kh/KI = anisotropic ratio Ks = hydraulic conductivity [LIT] Kxx = horizontal hydraulic conductivity in x direction [L/T] Kyy = horizontal hydraulic conductivity in y direction [L/T] K= vertical hydraulic conductivity [L/T] Ki. = the saturated hydraulic conductivity tensor [L/T] K = the relative permeability while the degree of water saturation is Sw L = average peat thickness [L] n = current integer time level Nie) = the interpolation functions for each node within element e n = a unit outward normal vector to the boundary No = the number of observation wells N, = the total number of nodes used in the numerical solution p = the number of nodes within element e Pn = net precipitation ( total precipitation minus intercepted precipitation[ L3/T] Pc = the capillary pressure ( Pc = 1) [MILT2] PET = the daily potential evapotranspiration [L/T] Q = pumping discharge [L3T] r = the radial coordinate [L] radial distance from the well [L] ro = the radius of the circular lake [L] s = drawdown [L] S = the specific storage Si = surface inflows, including flooding stream [ L/T] So = surface outflows [L3T] Ss = the specific storativity [L-'] SA, = wetland wetted surface area [L2] SAv = wetland dry surface area [L2] SP = exchanging flux [L3/T] S = the degree of saturation, equal to 06/6 S = the residual saturation Se = effective saturation t = time [T] tf = the time duration of the simulation [T] t = the observation times [T] T = the transmissivity in the aquifer [L2IT] T, = tidal inflow (+) or outflow (-) [L3/T] u = space function Vw = the wetland volume [L3] Vp = the volume above the bottom of the peat layer [L3] W = wetland water level [L] Wi( = the weighting function x, y, z = the axes in the physical domain X = the function of time Yo = the Bessel function of the second kind order 0 z = the vertical coordinate or elevation head ( L) V = the gradient operator AH = the average difference in hydraulic head [L] Ar = the radial step [L] At = the time step [T] AV = change in the volume of water storage in wetland [L3] Ax = the mesh space step in x-direction [L] Az = the mesh space step in y-direction [L] 60 = the angular frequency S,r] = the axes in the computational domain r = the boundary of volume 0 F, = flux boundary F2 = head boundary y = specific weight of a unit volume of water = pg [M/L2T2] i = the pressure [M/LT2] Q, = the exclusive subdomain of node i 8K = variation of K [L/T] 8h = variation of head h [L] '(x,y,tf- = the adjoint state of h T) S= tf- t [T] a and m = parameters obtained from a fit of above equations to experimental results 6 = the water content, [L/L3] 6s = the saturated water content, [L3/L3] Or = residual moisture content, [L3/L3] = pressure head; [M/LT2] l* = the value of i when ET = 0.5PET, [M/LT2] wi = the i corresponds to PE/PET=0.05, [M/LT2] ,W, = the value of i when PE=0.95 PET [M/LT2] Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy TWO-DIMENSIONAL ANALYTICAL AND THREE-DIMENSIONAL FINITE ELEMENT METHOD MODELING OF THE INTERACTIONS BETWEEN WETLANDS AND GROUNDWATER By Jinhua Zhong December, 2004 Chairman: Louis H. Motz Cochair: Kirk Hatfield Major Department: Civil and Coastal Engineering Simulating saturated/unsaturated flow between wetlands and the contiguous subsurface is complex, because it involves the modeling of wetland storage changes that in turn cause changes in pertinent hydrologic boundary conditions. Existing groundwater models are incapable of simulating the flow between a wetland and an aquifer, if small changes in wetland storage cause significant vertical and lateral movements of the wetland boundary. In addition to the simulation complexities introduced by a change in wetland geometry and storage, periodic and nonperiodic fluctuations of the wetland stage induce corresponding elevation changes in the phreatic surface of the contiguous aquifer. As these induced waves propagate through an aquifer, friction causes a loss of energy, which is manifested as spatially dampened potentiometric head perturbations along the inland direction. An analytical model was developed that describes subsurface flows around a wetland induced by any nonperiodic fluctuations in surface-water stage, such as a flood wave. Not considered, however, are any changes in wetland storage, or wetland geometry as a function of surface water stage. The analytical model was validated using a finite-difference numerical groundwater model that is based on the governing equation expressed in radial coordinates. Comparisons of analytical and numerical results show excellent agreement. Two numerical models were specifically designed to simulate wetland-aquifer interactions. The first model is a saturated groundwater flow model that incorporates adaptive simulation technologies to permit real-time simulation of moving wetland boundaries and their associated local influence on the phreatic surface. This model is numerically efficient and uses deformable hexahedral finite elements to trace phreatic surface changes in real-time. The second model also uses deformable hexahedral finite elements; however, this model simulates variably saturated groundwater flow. For both models, the numerical elements deform as required to characterize the horizontal and vertical extent of the moving wetland boundary (which serves to define a location within the porous system, where the pressure is known). Simulation results from both models were applied to a field site where water was rapidly withdrawn from an isolated wetland to observe system response and evaluate the wetland aquifer hydraulic connection. CHAPTER 1 INTRODUCTION 1.1 Wetland Wetlands are transitional lands between terrestrial and aquatic systems where the water table is usually at or near the surface, or the land is covered by shallow water (Cowardin et al. 1979). Wetlands are also called marsh, fen, peatland or water and can be natural or artificial and permanent or temporary. The water can be static or flowing and fresh or brackish, or it can include areas of marine water the depth of which at low tide does not exceed 6 meters (Ramsar Convention 1971, Middleton 1998). Wetlands perform various environmental functions that include increasing flood mitigation, enhancing storm abatement (i.e., in the case of coastal wetlands), adding groundwater recharge, providing water quality treatment, serving as wildlife habitats, and performing aesthetic functions (Mitsch and Gosselink 1993). These functional attributes can be classified into three basic types: hydrologic, water quality, and life support. Measures of the life-support function must integrate several wetland attributes, including water quality, hydroperiod and hydrodynamics, habitat structure, food chain support, biogeographic setting, nutrient status, corridors for migration, and primary production (Preston and Bedford 1988). In recent years, increasing interest has focused on the study and protection of wetlands. Development, drainage, population pressure, and pollution have destroyed extensive areas of these sensitive ecosystems. Another increasingly significance to wetlands is municipal well fields supplying water to large urban areas. These well fields affect 1 2 wetlands by reducing surficial aquifer levels over large areas that encompass isolated and interconnected wetlands (Sonenshein and Hofstetter 1990). As a result of depressing the water table, historically wet ecosystems have been drained (Switt et al. 1998). 1.2 Wetland Water Budget Hydrologic conditions are extremely important for the maintenance of a wetland's structure and function. A hydrologic budget is often used to quantify the volume and rate of water flowing into and out of a wetland. The primary hydrologic inflows include precipitation, streams, groundwater, and, in coastal wetlands, tides (Mitsch and Gosselink 1986, Rosenberry 2000). The wetland water budget equation as defined by Mitsch and Gosselink (1993) is: AV/At = Pn +S, +G,-ET-So-G,T (1-2) where AV = change in the volume of water storage in wetland [L3] At = change in time [T] Pn = net precipitation (total precipitation minus intercepted precipitation (I)) [ L3/T] Si = surface inflows, including flooding stream [ L3T] Gi = groundwater inflow [L3/T] ET = evapotranspiration [L3/T] So = surface outflows [L3 /T] Go = groundwater outflows [L3/T] T = tidal inflow (+) or outflow (-) [L3T] 3 Terms in the wetland water budget vary in importance depending on the wetland type, and not all terms are required to characterize a specific type of wetland. In most wetland systems, precipitation and evapotranspiration are major components of the water budget (Winter and Woo 1990). Groundwater can also have a significant influence on some wetland systems, but not always. Three of the most common relationships that characterize a wetland are groundwater discharge, groundwater inflow and outflow, and groundwater recharge. Groundwater inflow occurs when the contiguous wetland surface water (or groundwater) stage is lower than the phreatic surface of the surrounding aquifer. When the surface water stage or groundwater levels inside a wetland are higher than the phreatic surface in the surrounding contiguous aquifer, groundwater and surface water will flow out of wetland. When a wetland is situated well above the groundwater, the wetland is identified as being perched. This type of hydrologic system, referred to as a surface water depression wetland by Novitzki (1979), loses water only through infiltration into the ground and through evapotranspiration. The rate of water exchange between a wetland and local aquifer is a function of the hydraulic gradient between the two hydrologic subsystems. The discharge, in accordance with Darcy's law (Truesdell 1995), is proportional to the hydraulic gradient, the hydraulic conductivity, and the cross-sectional area through which the flow occurs. Periodic or nonperiodic fluctuations of the wetland stage in a flow-through wetland will induce fluctuations in the elevation of the phreatic surface in the contiguous aquifer. As these induced waves propagate through an aquifer, friction causes a loss of energy, which is manifested as spatially dampened head perturbations. This damping effect progressively diminishes the amplitude of piezometric head along the inland direction (Carr and Van Der 4 Kamp 1969). Many studies have examined the problem through one-dimensional analytical solutions (Cooper 1963, Gregg 1966, Nielsen 1990, Serfes 1992, Li, Barry et al. 2000, 2001, Boufadel and Peridier 2002, Seerano 2003, Swamee and Singh 2003). A two-dimensional analytical solution was obtained by considering a shift in phase and a change in amplitude along a straight coastline (Sun 1997). Three-dimensional flow regimes near a shallow circular lake was studied using numerical and analytical mathematical models of aquifer flow (Townley and Trefry 2000). To further the understanding of the interactions between lakes or wetlands and groundwater, without the restrictive assumptions typically associated with analytical solutions, numerical models are usually used. For example, it is possible to obtain numerical solutions for the case of anisotropic and nonhomogeneous conditions. The accuracy of solutions obtained by numerical methods can be excellent; however, it depends on several factors, including the type of numerical method used, the complexity of boundary and initial conditions, and the computational precision of the computer used to implement the method. Regardless of these advantages, few investigators have chosen to focus on the numerical simulation of interactions between wetlands and groundwater ( McDonald and Harbaugh 1988, Schaffranek 1996, Lee 2000). 1.3 Goals and Objectives Analytical models to simulate the complex groundwater surface water interactions that occur with isolated and flow-through wetlands would be valuable tools for developing hydrologic budgets for wetlands. Given that most surface water fluctuations are nonperiodic, there is a need to develop analytical groundwater models for simulating groundwater table fluctuations induced by both periodic and nonperiodic surface water flood waves. However, 5 with the exception of the one-dimensional study by Cooper and Rorabaugh (1963), most investigations addressed surface-water fluctuations that were periodic. Thus, the first goal of this dissertation was to develop an analytical method for simulating groundwater table fluctuations induced by both periodic and nonperiodic surface-water flood waves. This goal was achieved by meeting the single objective of developing an analytical model using Fourier transforms. A limited number of three-dimensional numerical groundwater models have been developed to simulate groundwater wetland interactions; however, these models are generally inadequate tools for simulating the complex interactions between a wetland and a contiguous unconfined aquifer. For example, Walser, Annable and Wise (1998) used FEMWATER (Yeh, 1992), a variably saturated flow model, to simulate the exchange of water between a wetland and an aquifer. McDonald and Harbaugh (1988) used MODFLOW and developed a three-dimensional finite difference groundwater model and studied the dynamic coupling between wetlands and the subsurface. Lee (2000) used the HYDRUS-2D flow model and simulated vertical two-dimensional groundwater flow through Lake Barco in Florida. Neither FEMWATER (1996) nor MODFLOW can simulate or incorporate a moving real-time wetland boundary caused by a change in surface water stage. Furthermore, as a saturated three-dimensional finite-difference model, MODFLOW can not simulate fluctuations of the phreatic surface as a moving transient boundary as required to define fully the three- dimensional flow. Finally, the water stage of a wetland itself is a function of the transient exchange of water between a surface and subsurface systems; which is not fully modeled in any existing groundwater flow codes. Thus, to simulate groundwater surface water interaction without the above identified shortcomings of existing models, a new model is needed. The second goal of the dissertation was to develop two numerical models capable of simulating 1) a transient and moving wetland boundary, 2) three-dimensional flow, 3) the movement of the phreatic surface, and 4) the coupling that occurs between the transient surface-water stage and the phreatic surface. To achieve this goal, two objectives were pursued. The first objective was to develop two numerical groundwater models using finite elements. The second objective was to apply both numerical models to simulate groundwater and surface water fluctuations over an isolated wetland in southern Florida. Of the two models developed, the first is a three-dimensional saturated groundwater flow model that employs deformable hexahedral elements. The finite-element mesh deforms in space and time to trace the location of the phreatic surface of an unconfined aquifer and the stage of the contiguous surface water body. Based on the above model and the adjoint state method, an inverse model was developed to estimate the hydraulic conductivity of peat layers which typically exist as a wetland bottom that serves to isolate surface wetlands from subsurface groundwater systems. The second numerical model, a three-dimensional variable saturated flow model, was also developed using deformable mesh. This model was used to investigate prediction differences between the saturated flow code and the variable saturated model. From the groundwater flow pattern given by both models, an improved understanding was obtained for the interactions occurring between surface waters and groundwaters. In addition, an improved understanding was obtained of how short-term pumping from a phreatic aquifer affects wetlands and lake water stages. CHAPTER 2 ANALYTICAL GROUNDWATER MODEL RESPONSE TO THE PERIODIC OR NON-PERIODIC WATER STAGE FLUCTUATIONS IN A CIRCULAR WETLAND 2.1 Introduction Periodic or nonperiodic fluctuations of the surface water stage in wetlands or lakes will induce fluctuations in the elevation of the phreatic surface of contiguous unconfined aquifers. As these induced waves propagate through the aquifer, friction causes a loss of energy, which is manifested as spatially dampened head perturbations. This damping effect progressively diminishes the amplitude of fluctuations in the piezometric head along the inland direction (Carr and Van Kamp 1969). Many studies have examined the problem through one-dimensional analytical solutions (e.g.,Cooper and Rorabaugh 1963, Gregg 1966, Nielsen 1990, and Serfes 1992, Li, Barry et al. 2000, 2001, Boufadel and Peridier 2002, Seerano 2003, Swamee and Singh 2003). Cooper and Rorabaugh (1963) obtained the analytical solution for groundwater movement and bank storage due to flood stages in surface streams. Gregg (1966) found a one-dimensional relation between tidal fluctuations and groundwater table perturbations. Nielsen (1990), Serfes (1992) and Carr (1969) also published papers investigating the similar problems. A two-dimensional analytical solution was obtained by considering a phase shift in water-table perturbations and the associated change of amplitude along a straight coastline (Sun 1997). An analytical solution representing vertical two-dimensional groundwater-surface water interaction was obtained by Anderson (2003). Three-dimensional flow regimes near a shallow circular lake was 7 8 studied using numerical and analytical mathematical models of aquifer flow (Townley and Trefry 2000, Smith and Townley 2002). With the exception of Cooper and Rorabaugh (1963), all of the above studies addressed periodic surface-water fluctuations, in spite of the fact that most stage variations are not periodic i.e., consider, for example, the pervasive nonperiodic flooding of rivers, lakes, and wetlands. Cooper's (1963) nonperiodic analytical model for a river/groundwater system was a one-dimensional solution along an axis perpendicular to the river. However, this analytical model is not suitable for simulating the nonperiodic interactions between a wetland and an aquifer. In this chapter, two analytical groundwater models are developed to describe transient groundwater levels affected by periodic and nonperiodic fluctuations of the surface- water stage in a wetland or a lake. 2.2 Problem Formulation: Groundwater Flow Around a Circular Wetland Considering a circular wetland that is hydraulically contiguous to a surfacial aquifer and ignoring recharge from rainfall, the governing equation for two-dimensional groundwater flow in a homogeneous isotropic aquifer can be expressed in radial coordinates as: 2h 1 Dh S Oh 2 h (2-1) ar 2 r r Tat where h is the piezometric head, r is the radial coordinate, t is time, S is the specific storage and Tis the transmissivity in the aquifer. When a thick unconfined aquifer is considered, the gradient and small fluctuations relative to saturated thickness are assumed to be small, T is approximately equal to Khoo (Figure 2-1) in which K is the hydraulic conductivity and ho is 9 the average unconfined aquifer thickness equal to the regional average water-table elevation above an underlying impermeable confining unit. The boundary condition along the shore line of the wetland is Eq. 2-2, where r o is the radius of the circular wetland. It is also assumed in Eq. 2-3 that as r--o, the groundwater table elevation approaches the regional average water table elevation, hoo. B.C: h(r,t) ,=r, = hoo + f(r ,t) @ r= r. (2-2) h(r,t) ,r= = hoo @ r= o (2-3) wetland or lake _ h_ //777/7- //7777//77///77/7 77777/ Figure 2-1. Setup of the groundwater problem around a circular wetland 2.3 Analytical Solution for the Periodic Boundary Condition Since the boundary condition is a periodic function of time, substituting forfl(rt), Equation 2-4 is obtained: B.C: h(r,t) r=ro = hoo + h ei't @ r= r, (2-4) h(r,t)| ,- = hoo @ r= co (2-5) where ho is the half height of the periodic constituent and to is the angular frequency. No initial condition is necessary, since here the steady-state oscillatory solution is sought as t -,. Because of the oscillatory boundary condition, a solution is sought in the form h(r,t) =ho +hog(r) e (2-6) (2-6) Substitution of Eq. 2-6 into Eq. 2-1, Eq. 2-2, and Eq. 2-3 yields d2g +l.O g S.i-g=0 (2-7) Br2 r ar T where the boundary conditions are: g(r)=] @ r=ro (2-8) g(r)=O @ r=- (2-9) Equation 2-7 can now be rewritten to obtain: v2 +vg +v2g= av2 Ov (2-10) v 2= T (2-11) The solution to Equation 2-10 subject to the boundary conditions in Equation 2-8 and 2-9 depends on the value of w. Considering the boundary condition along the shoreline, groundwater around a circular water body will have the following forms. When (o is greater than 0: g(r)=Co( +iY)= CoH v)- H() (2-12) H' (vC) and when wo is less than 0: g(r) = C -iY) = CH (u)= (2-13) where Jo is the Bessel function of the first kind order 0; Yo is the Bessel function of the second kind order 0; Ho(t) is the Hankel function order 0; and CO is a coefficient (Lebedev 1965). Substitution of Eq. 2-12 and Eq. 2-13 into Eq. 2-6 gives the following solution for the groundwater flow around a circular lake affected by periodic function: h,(rt)=hoo+h0o )e H.(.vo) (2-14) for 0 > 0 and: h.(r,t)=h0 +h0 -e Nt( S (o2)(o) (2-15) for c < 0. The asymptotic approximations of the functions Ho')(z) and Ho(2)(z) are: H (z)= 2 i*) [ 1-32**(2k- 1)2(2 iz)-k+O(z --1) (2-16) S7 z k= 22k-k! I i 2 /) 1/2 (-i e( )~ 132..-(2k 1)2 2)(z)= 2 / 4 ( iz)-k+O(|zi -n-1) VZ k=o 22k-M! uo= '*(-1-O'i)-ro V0 = -*(i-1)ro (2-17) u= .*(-1-i)-r 2v= T F2T 2.4 Analytical Solution for the Non-Periodic Boundary Condition A flood wave will typically induce a nonperiodic variation in the wetland stage. Here, f(ro,t), which is a nonperiodic function of time, is assumed to describe transient variations in the wetland stage. Using the Fourier transformation, f(ro ,t) can be represented as: J(ro,t)= fF()e s"'d" (2-18) where F(wo) is the Fourier coefficient, and F(o)dw is the amplitude of oscillation for the component with a frequency in the range of (o,(+do). For this F(x), the corresponding fluctuations in the groundwater table are described by: (W >0) HO 2)(v)) h.(r,t)=hoo + F(w) *e w (2-19) H2)(u) h,(r,t)=hoo + F(o) "*eu ,, k) ( o < 0 ) (2-20) as given by Eq. 2-14 and Eq. 2-15. The inverse Fourier transformation of Equation 2-19 and Eq. 2-20 leads to: h(r,t)=h oo+ + 2n =h0 2x -h' H'(v) e e Hf ~(v0) J H (u0) (2-21) wl.dj From the above analytical model, it is clear that the unsteady groundwater table around a circular wetland depends not only on the instantaneous values of the wetland water stage but also on the accumulative effect represented by the integration. However, physically it is not possible for the wetland water stage at a later time to affect the groundwater table at an earlier time. Hence, the phreatic surface depends on its current stage and historical stage. The _na =ooh +- 14 final analytical solution to the nonperiodic boundary condition 2-2 is then h(r,t)=h o+- In 2x (2-22) (1),.a f (v) H()(o) 0 2.5 Application If it is assumed that wetland stage fluctuations can be described by a nonperiodic Gaussian function of time, then: 1 202 (ro, t) =- e 2a V2-n G (2-23) The analytical solution obtained when Eq. 2-23 is combined with Eq. 2-21 is: 1 r 1 ( t H(2) h(r,t)=h o+ e 202 2 e) e(t-)i*d In 2 f H2)(Uo) (2-24) The transient changes of groundwater table described by Eq. 2-24 at four different distances from the wetland are shown in Figure 2-2. The groundwater distributions along the radius at four different times are shown in Figure 2-3. From both figures, it is evident that groundwater fluctuations are greatest near the wetland. From Figure 2-3, it can be seen that the variation of the groundwater table over time is no longer symmetrical although the " '(-)'-d*o dr variation of the wetland water stage is a Gaussian distribution which is symmetrical over time. 101 I 100.8 | I 100.6 100.4 100.2 100 99.8 0 5 10 15 20 25 Time (hr) lake d=1m d=2m d=3m d=4m Figure 2-2. Transient groundwater level changes at radial distances 1,2,3, and 4 meters away from the wetland shoreline. 100.3 - H00.25 :100.2 3- 00.15 100.1 900.05 (o 100 t 99.95 50 t=6hr t=7hr t=8hr t=9hr t=10hr 55 60 65 Distance to lake shoreline (m) Figure 2-3. Groundwater distribution over the distance to the wetland/lake shoreline at 6,7,8,9, and 10 h. 16 Figure 2-4 shows the groundwater table elevations around a circular wetland after 3 hours. The radius of the wetland is 50 meters. All other parameters are listed in Table 2-1. Theoretically, when o is very large, it means that the frequency of the corresponding boundary condition is very large and the period is very short. When the absolute value of a is very large, its effect on the groundwater around the wetland is very small. This is because the period of the wave component is too small to have any impact on the groundwater around the wetland. That is to say, when the integration of a is taken, a certain negative number could be used to replace the negative infinite number, and a certain positive number could also be used to replace the positive infinite number. Accurate model results are obtained when w takes the range of [-5, 5]. Thus, it may be assumed that the boundary condition has no effect on the ground water table if the frequency of the boundary water level is very large. H (M) Figure 2-4. Three-dimensional view of the groundwater table around a circular wetland at t=3 h. 17 2.6 Comparison of Analytical Solution and Numerical Solutions A comparison was made between the analytical and the numerical solutions to the radial flow problem specified through Eq. 2-1-Eq. 2-3. An implicit finite-difference scheme was used to represent the radial form of the groundwater flow equation Eq. 2-1 as: h,- -2h" i+h,' I h11 -hi-I _S h, -h(-5) (Ar)2 r, 2Ar T At Here i is the reference integer spatial interval; n is the current integer time level; At is the time step; and Ar is the radial step. A circular wetland which has a radius of 50 meters was assumed. All other system parameters are given in the Table 2-1. By considering the different boundary conditions (Eq. 2-3 and Eq. 2-4), the corresponding numerical results are found [Figure 2-5]. Analytical and numerical results illustrate time variations in groundwater levels. Figure 2-5 shows that the analytical solution and numerical solution compare very well. Table 2-1. Parameters and their values Parameter Value Parameter Value Hoo (m) 100.0 T (m2/h) 0.635 Ar (m) 0.25 S 0.2 At (hr) 0.02 to (h) 5.0 k (m/hr) 0.00635 o 0.4 0.4 S0.3 analy(d=2m) 0.2 0.2 ............. analy(d=5m) 0.1 -o 0num (d=2m) 2. num (d=5m) -0.1 7 __ 0 5 10 15 20 25 Time (hr) Figure 2-5. Comparison of analytical solution and numerical solution at 2 and 5 meters from the wetland shoreline. 2.7 Summary The analytical model presented in this chapter can be used to simulate aquifer surface fluctuations caused by the transient wetland stages regardless of the temporal function assumed to describe variations of water stage. In other words, the model can deal with any general transient of the boundary function. The analytical model was validated using a finite- difference numerical groundwater model for a homogenous aquifer. This first numerical method is not recommended for simulations involving heterogenous aquifers. Rather, a new method is presented for simulating wetland/aquifer interactions in such aquifers, which requires a simple transformation of the equations. The excellent agreement between two models shows that the analytical solution presented is valid. This model is also useful in a situation where groundwater and the surface water components are being considered in a system-wide water balance. CHAPTER 3 TRANSIENT THREE-DIMENSIONAL DEFORMABLE FINITE-ELEMENT SATURATED GROUNDWATER MODEL 3.1 Introduction Recent research has focused on the details of various wetland functions such as water-quality treatment, flood-wave suppression, and groundwater/surface water interactions (Wise et al. 2000, Krabbenhoft and Webster 1995, Anderson and Cheng 1993, Ruddy and Williams 1991, and Mills and Zwarich 1986). For example, Ruddy and Williams (1991) investigated the hydrologic relation between the South Fork Williams Fork stream in the Colorado River basin and four adjacent subalpine wetlands, and they also evaluated the potential effects on these wetlands when stream flow is diverted. They concluded that stream flow was the primary source of water for two of the four wetlands, whereas precipitation, snowmelt, valley side-slope flow (including overland and small-channel flow), and groundwater inflow served to define the transient hydrological behavior of all wetlands. The above studies showed that wetland functions (i.e., water quality treatment, floodwave suppression, etc.) are correlated with recharge and discharge areas, ground-water- flow paths and rates, and soil-profile development. The location of recharge and discharge areas can change in response to seasonal or weather-related factors and can affect the groundwater flow direction. As a result, groundwater and surface water quality and the presence of various wetland-plant communities may all vary in response to these changing conditions. 20 Analytical models such as the one presented in Chapter 2 can be used to simulate deep wetlands or lakes. However, the primary limitation of analytical methods is that solutions can only be obtained by imposing severely restrictive assumptions on aquifer properties and boundary and initial conditions. For example, an assumption commonly made to obtain analytical solutions is that the aquifer is isotropic and homogeneous for hydraulic conductivity. These assumptions, however, are not valid and often not adequate to describe groundwater flow or solute transport in most field situations. Numerical methods have the primary advantage over analytical approaches in that they do not require such restrictive assumptions. For example, it is possible to obtain numerical solutions for the case of anisotropic and nonhomogeneous conditions. The accuracy of solutions obtained by numerical methods can be excellent. However, this accuracy depends on several factors, including the type of numerical method used, the complexity of boundary and initial conditions, and the computational precision of the computer used to implement the method. Regardless of these advantage, only a few investigators have chosen to focus on the numerical simulation of interactions between wetlands and groundwater (e.g., Schaffranek 1996 and McDonald and Harbaugh 1988, Lee 2000, Rosenberry 2000). In this chapter, a general discussion of the advantages/disadvantages of two numerical methods commonly used to model hydrologic systems is briefly discussed. Next, the mathematical development of a three-dimensional deformable finite-element saturated groundwater model is presented. This model and the associated code are then validated by comparing model predictions with analytical modeling results. A more complicated application of the three-dimensional model is presented to illustrate the type of information provided through the numerical approach taken. Finally, the 21 model is used to simulate an actual wetland using field data collected when the hydrologic system was responding to induced stress. 3.2 Ground-Water Flow Models Based on Finite-Difference and Finite-Element Approximation Several types of numerical methods have been used to solve groundwater flow problems, the two principal ones being the finite-difference method and the finite-element method. The finite-difference method was initially applied to the flow of fluids in petroleum reservoirs, and it wasn't until the mid-1960's that problems of groundwater flow and solute transport were examined. The method has a number of advantages that contribute to its continued widespread use and popularity: (1) for simple problems such as one-dimensional, steady-state groundwater flow in an isotropic and homogeneous aquifer, mathematical formulation and computer implementation are easily understood; and (2) the accuracy of solutions to steady-state and transient groundwater flow problems is generally quite good. The finite-difference method also has disadvantages: (1) the method works best for rectangular or prismatic aquifers of uniform composition (for example, it is difficult to incorporate irregular or curved aquifer boundaries, anisotropic and heterogeneous aquifer properties, or sloping soil and rock layers into the numerical model without introducing numerous mathematical and computer programming complexities); and (2) the accuracy of solutions to solute transport problems is less than can be obtained through application of the finite-element method. Table 3.1 lists selected early references where finite-difference methods were applied to model groundwater flow. More recent applications of hydrologic modeling using finite-difference for simulating of surface water /groundwater interactions include studies by Mc Donald and Harbaugh (1998), Winter (1983, 1986, and 1989), Sacks 22 (1992) and Krabbenhoft and Anderson (1990). Table 3-1. Selected early references for the finite-difference method. Topic Reference Early Developments in Petroleum Bruce et al. (1953), Peace and Reservoir Modeling Rachford (1962) Saturated Groundwater Flow Remson et al. (1965), Freeze and Whitherspoon (1966), Pinder and Bredehoeft (1968). Unsaturated Groundwater Flow Philip (1957), Ashroft et al. (1962), Freeze (1971), Brutsaert (1973), Healy (1990), Kirkland et al. (1992). Solute Transport Stone and Brian (1963), Oster et al. (1970), Tanji et al. (1967), Wierenga (1977) Application to Field Problems Orlob and Woods (1967), Gambolati et al. (1973), Fleck and McDonald (1978), Jain 1992 Comprehensive Reference Trescott and Larson (1977), Ames (1977), Mitchell and Griffiths (1980), Lapidus and Pinder (1982) Computer Program Trescott et al. (1976), Konikow and Bredehoft (1978), Healy (1990) Winter (1983, 1986, 1989), Sacks (1992) and Krabbenhoft and Anderson (1990) conducted numerical investigations on the interactions between groundwaters and lakes. Both Winter and Krabbenhoft and Anderson used fixed finite-difference meshes. For lakes, a finite-difference mesh may be suitable; however, because the bed of most wetlands exhibits significant slope, a poorly designed fixed finite-difference mesh may not adequately characterize water fluxes across a bed situated between an overlying body of water and an underlying aquifer. The issue of properly modeling fluxes across a sloping bed may not be 23 as important as the fact that accurate simulations of wetland-flow systems require measurements of hydrologic properties at finer scales and frequencies than are commonly made (Hunt 1996, Hunt and others 1996 1997). Assuming the hydrological model is correct, the accuracy of any computer simulated hydrologic system, such as a wetland, depends on the quality of the data used to calibrate model. Through a well instrumented study-site in the Kankakee River Valley in northwest Indiana, the USGS and the USEPA conducted detailed measurements of many various variables required to characterize a wetland hydrologic system. McDonald and Harbaugh (1998) developed a groundwater model using the available data and conducted the simulations of observed hydrologic process using the computer code commonly known as MODFLOW. The focus of the research effort was to study the effects of various wetland restoration activities on local groundwater levels. Computer simulations were used to predict the effects of wetland restoration on the existing hydrology of the wetlands and adjacent farmlands and to suggest residence times and path ways for surface water and groundwaters. The finite-element method was first used to solve groundwater flow and solute transport problems in the early 1970s. The method has several advantages: (1) irregular or curved aquifer boundaries, anisotropic and heterogeneous aquifer properties, and sloping soil and rock layers can be incorporated easily into the numerical model; and (2) the accuracy of solution to groundwater flow and solute transport problems is quite good. The main disadvantages of the finite-element method include: (1) for simple problems, the finite element method requires a greater amount of mathematical and computer programming sophistication than does the finite-difference method; and (2) there are fewer well- documented computer programs. 24 The finite-element method provides approximations of much higher order than does finite-difference methodology. The whole modeling domain can be subdivided into irregular elements as dictated by the physical geometry of the problem. The size of each element can vary. Small elements can be used in areas where rapid change occurs with model parameters and predicted state variables, whereas the large elements can be used where such variations are less severe. Inhomogeneities and anisotropy are easily accommodated with the finite- element method. The advantage of irregular elements and higher-order approximations is certainly important when the "space" domain is complex. However, this advantage is less obvious in problems involving time (evolution) (Gupta 1975). Selected references of early finite-element applications to model groundwater flow are listed in Table 3-2. 3.3 Three-Dimensional Deformable Finite-Element Saturated Groundwater Model To simulate three-dimensional groundwater flow in phreatic aquifer systems and the exchange of water between a wetland and the contiguous aquifer, a numerical model is needed that is capable of modeling the free groundwater surface. This requires first that the model assume a tentative location of the free surface and then solve for the distribution of the dependent variable, *ii, over a fixed domain. All this requires that the flux boundary condition across the surface be satisfied and that the pressure be atmospheric. To model the exchange of water between a wetland and an aquifer, it is necessary to recognize that significant vertical and lateral movements of the wetland boundary typically occur with changes in wetland storage (see Figure 3-1). This transient boundary suggests that a deformable finite-element mesh may be needed in the model to characterize adequately the moving groundwater/surface-water interface. 25 Table 3-2. Selected references of early finite-element applications to model groundwater. Topic Reference Early Developments in Petroleum Price et al. (1968) Reservoir Modeling Saturated Groundwater Flow Zienkiewicz et al. (1966), Javandel and Witherspoon (1968), Zienkiewicz and Parekh (1970), Pinder and Frind (1972). Unsaturated Groundwater Flow Neuman (1973), Gureghian et al. (1979), Pickens and Gillham (1980) Yeh and Ward (1981),Huyakom (1984), Yeh (1992). Solute Transport Price et al. (1968), Guymon et al. (1970), Neuman (1973), Van Genuchten et al. (1977), Kirkner et al. (1984). Application to Field Problems Pinder (1973), Gupta and Tinji (1976), Senger and Fogg (1987). Comprehensive References Ziekiewicz (1971), Pinder and Gray (1977), Lapidus and Pinder (1982), Huyakorn and Pinder (1983). Computer Programs Neuman and Witherspoon (1970), Reeves and Duguid (1975), Segol et al. (1975), Pickens et al. (1979) Huyakom (1984), Yeh (1992). Numerical groundwater modeling that is focused on simulating the interaction between surface-water and groundwater systems, where the groundwater/surface-water interface moves, must use well-defined linking conditions to describe the exchange of water across the interface. These linking conditions are unknown in the beginning of the modeling time step and found through a flux type interaction that then forecasts the linking boundary conditions at the end of each time step. Thus, the model has transient boundary conditions. Wetland SWater Surface at Time t, SWater Surface at Tim -- Phreatic Surface at Time t, Phreatic Surface at Time t2 Groundwater Aquifer Impermeable Aquacludc Figure 3-1. Conceptual profile model of a wetland/groundwater system Finally, typical finite-difference modeling is not able to trace the curvature of the bottom of a wetland (see Figure 3-2). Finite elements, however, are easily configured to fit non-linear spatial geometries. In conclusion, to cope with the above described three technical problems of modeling, i.e., the free surface location, incorporating dynamic changes in the vertical and horizontal extent of the ground/surface-water interface, and tracing the curvature of a wetland bottom, a special three-dimensional finite-element model is needed. Finite-Element Method Finite-Difference Method Figure 3-2. Comparison of finite-element mesh and finite-difference mesh 27 3.3.1 Governing Equation and Boundary Conditions Eq. 3-1 below can be used to describe incompressible, saturated groundwater flow beneath a wetland and in a phreatic aquifer: V.(KVh) Q = S, ah ij= 1,2,3 (3-1) at where K is the saturated hydraulic conductivity tensor ( LT -'); h is equal to z+* (L); i is pressure head (L); z is elevation head (L); and S, is the specific storativity (L'1). From Figure 3-3, the boundary of volume 0 is characterized through a combination of specified flux boundaries r, and specified head boundaries P2. Pertinent specific no-flux boundaries include the groundwater phreatic surface (EC and DF) and the bottom impermeable boundary (GH). Relevant specified head boundaries r2 include the wetland head boundary (AB and CD) and lateral boundaries (EG and FH). E F Peat Layer Groundwater Aquifer G HI GH, EC, and DF are no-flow boundaries F EG, FH, ACDB, and CD are specified head boundaries rF Figure 3-3. Setup of the saturated groundwater problem around a wetland 28 As to the boundary condition controlled by the wetland water table, the initial water table of the wetland and the pumping discharge were taken as the only inputs of the coupling model. During pumping, the pumping discharge lowered the wetland water table, and it also triggered an inflow from the aquifer to the wetland through the peat layer. In each simulation time step, the pumping discharge, estimated inflow from the aquifer to the wetland, initial water table of the wetland, and relation between the volume of the wetland and wetland water table were used as inputs in a mass balance equation that was solved to determine the water table at the end of the time step. The estimated water table at the end of the time step was put back into the saturated groundwater model to adjust the inflow from the aquifer to the wetland estimated initially. Therefore, several iterations in each time step were necessary to get an accurate estimated inflow from the aquifer to the wetland and a convergent solution of the water table in the wetland at the end of the time step. When pumping was stopped, the pumping discharge was taken as zero. Figure 3-4 shows the procedure of finding the inflow from the aquifer to the wetland and the coupling boundary. 3.3.2 Finite Element Model of Saturated Groundwater Flow To obtain a finite element approximation of Eq. 3-1, two spaces, S and V, are defined (Hughes 1987): S= {uIuEH1, u(xyz)2=qr2 } (3-2) V= uuEHuH, u(x,y,z)r2=O} (3-3) where H' is Sobolev space comprised of functions which are square-integrable generalized Predict a wetland water table hl(t+At) Simulate groundwater flow based on the predicted wetland water table Find the groundwater inflow to the wetland Wetland pumping Based on the water mass balance equation of the discharge Q wetland to find the wetland water table h2(t+At) Yes Next time step Figure 3-4. Flowchart illustrating the iteration procedure of finding the coupling boundary condition derivatives through order 1; u is a space function that is member of Sobolev space; qr` is the vertical flux along P2; and x, y, z are three spatial coordinates. Any function from S space is equal to the flux boundary value along specified head boundary (F2). Also, any function u from V space is equal to zero along head boundary (F2). Assuming a basis function w is chosen from V space, then the following integration of Eq. 3-1 over space domain Q is always true: w[V-(KVh)]dQ wQdO= wS, -'dQ, ij= 1,2,3 (3.4) J fa f at where Eq. 3-4 expresses the governing groundwater flow equation integrated over the modeling domain. Manipulation of the first term on the left hand side of Eq. 3-4 begins with the application of Green's theorem, which is an expanded form of the divergence theorem. This converts the integral into two terms, one of which is evaluated only on the surface of the region to be simulated. Green's theorem is: (v.-P)AdQ= j ( -.)AdT- (-VA)d (3-5) (3-5) o r a where A is a scalar, and P is a vector quantity. The boundary of volume 0 is defined by F, and n is a unit outward normal vector to the boundary. Application of the Galerkin criterion (Hughes 1986) and Green's theorem to Eq. 3-4 leads to: f w[(KVh)-n]dP+ f w[(KVh)-n]dd-f (Vw)-(KVh)dQ j w-QdQ = wS- dQ a 11 adt (3-6) 31 Because w is a member of V space, it is equal to zero at all specified head boundaries ("2); thus, the second term in Eq. 3-6 equals zero. If it is further assumed that Q equals zero, then Eq. 3-6 reduces to: w[(KVh)-n]dr- (Vw)-(KVh)dQ= wS, -dQ (3-7) For an approximate solution of h(x, y, z, t), it is assumed that: h(e)(x,y,z,t)= Nie)(x,y,z) d,(t) N,e V heS (3-8) i=1 where h(e) is the approximate solution for hydraulic head within element e, Ni(e) are interpolation functions for each node within element e, p is the number of nodes within element e, and d,(t) are the unknown values of hydraulic head for each node within element e and at time t. At the specified head boundaries, F2, di(t) = h(x,y,z)I2. wi(e) is the weighting function for node i, and the limits of the integration are chosen to represent the volume of the elements. In Galerkin's method, the weighting function for each node in the element is set equal to the interpolation function for the node, wi(e) = Ne). When the approximate solutions (Eq. 3-8 and Eq. 3-9) are substituted into Eq. 3-7, the resulting equation is not satisfied exactly; consequently, an error or residual occurs at every point in the problem domain. However, if it is also assumed that the aquifer is 32 isotropic, that is to say, values of saturated hydraulic conductivity in the three coordinate directions are constant within an element (but can vary from one element to the next), the contribution of any element e to the residual at a node i associated with the element can be described through the following system of equations: R (e)(t) dite)It W =[K(e) [~) Re)(t) d(e)(t) adi (e)(t) at ad(e)(t) at (3-10) -[Fe] where p = the number of nodes in each element (i.e., triangle element: p = 3), and Ri(e) = the residual error at node i in element e: (e ) =fff V(9 aN(ie) ax aN(e) ax aN(e) ay aN(e) ay aN (e) az a N(e) az K= 0 0 0 K 0 0 0 K and: [F(e) ffrace 1 p J aN (e) ax aN (e) ay aNle) az (3-11) aN.(e) ax aN (e) ay aN(e) az (3-12) and: N (e) C (e). e) ... e)] (3-13) V(e) NP where V(e) is the volume of element e. Eq. 3-10 represents a system of equations pertinent to a single element, and there is an equation for every node in the element. When equations for all elements in the mesh are combined, the global equations are obtained for the system: ad,(t) RI(t) dd(t) t =[K] + [C] -[F] (3-14) at By setting the residuals equal to zero, we have: ad,(t) [K] + + [C] I =[F] (3-15) at Eq. 3-15 is a system of ordinary differential equations, whose solution provides values of head d and the derivative with respect to time, ad/at, at each node in the finite-element mesh. 34 Although several methods are available for solving this system of equations, it has become standard practice to use the finite-difference method: ah d(t+At)-d(t) at At (3-16) If a variable, co is defined, such that: A-t At Then: (3-18) which can be extended to the vector of [d]: (3-17) [d(e)] = (1 -)[d(t)]+ [d(+At)] (3-19) If Eq. 3-19 and Eq. 3-18 are substituted into Eq. 3-15, the finite-difference formulation for the transient saturated flow equation is as follows: [C] + w At[K] )[d(t+At)] =( [C] (1 -w) At [K]) [d(t)]+At((1-o)[F(t)]+([F(t+At)] (3-20) d(e) = (1-o)d(t) +a d(t+AOt) 35 Eq. 3-20 represents a system of linear equations with unknown variables of d,(t), and where N,(x,y,z) are the basis functions in a tri-linear hexahedral element domain. Ni(x,y,z) does not change with time. The most common basis functions used are hat functions, which are linear functions of space variables. The solution of Eq. 3-20 begins by specifying the initial values of [d] (i.e., the values of head at time to= 0). Next, the system Eq. 3-20 is solved to obtain values of {d(t+At)} at the end of the first time step 1. This solution process is repeated for the next time step, and so on until the specified time duration of simulation is satisfied. Depending on the choice of w, several different subsets of the finite-difference formulations are defined. When o=0, it is the forward-difference method is obtained: ( [C] )[d(t+At)] = ( [C] At[K])[d(t)] +At[F(t)] (3-21) When w=1/2, the central-difference, or Crank-Nicholson, method is obtained: ([C] + 1 At[K] [d(t+At)] = [C] 1 -1 At [K] [d(t)] + (F(t)]+[F(t)]) (3-22) When (=1, it is the backward-difference method is obatined: ( [C] + At[K] )[d(t+At)] = [C] [d(t)] +At[F(t)] (3-23) 36 3.3.3 Adaptation of the Saturated Finite-Element Model to Describe the Free Surface Several papers have been written on the groundwater free-surface problem (e.g.,Taylor 1967, Neuman 1970, Finn 1976, France 1976, and Huyakom 1986). Most of these papers discussed the steady-state free-surface problem except for France (1976). This thesis examines the transient problem starting from the basic ideas that Neuman (1970) presented for steady flow. Basically, both France (1976) and Neuman (1970) tried to satisfy the two boundary conditions and are consistent. This chapter is also based on the two conditions and adapts the method they used to locate the free surface. The primary dependent variable to be used in modeling the free surface presented here is the hydraulic head h defined as: h = z + (3-24) Y where z is the vertical coordinate; y is specific weight of a unit volume of water pg and ip is the pressure. Along the free surface, there are two conditions which must be satisfied. First, qr=0 (pressure atmospheric), and hence for a given position, h is prescribed on it. That is: h =z (3-25) For the steady-state scenario, the second condition is simply that there is a specified flux across the free surface boundary (i.e., a zero flux or a flow imposed equal to net recharge). Thus, in addition to Eq. 3-24 at the free surface, there is: ah h h - K-i + K-j + K- k = q. (3-26) ax 'y zz If steady state has not been reached and the boundary is moving with a velocity normal to its instantaneous configuration of V., then the quantity of liquid entering its unit area is modified, and: Ah- h- 9h- K.,hax i+ K qj + K- k = q ,+V =q (3-27) where i, j, k are unit orthogonal vectors; S is the specific yield coefficient (or the effective porosity) relating the total volume of material to the quantity of liquid which can be drained from it, and q, equals the net flux across the free surface. The transient problem is solved by considering the solution to be a number of steady- state problems at small intervals of time At. For a given time level, the governing equation is solved to obtain values of nodal heads and velocities of free surface. The components of actual velocities are determined by dividing by the specific yield coefficient. By repeating this process for all modal points on the free surface, a new configuration and a new set of boundary conditions are defined. The actual nodal coordinates of the free surface are found using the Newton-Raphson iterative method. Thus, a new finite-element mesh is generated by modifying the coordinates of the nodes in the new flow domain. 38 3.4 Validating the Code of Three-Dimensional Deformable Finite-Element Groundwater Flow Model 3.4.1 Theis Solution To validate the code of the three-dimensional deformable finite-element groundwater flow model, a comparison was made of results generated by the numerical model and those created with an analytical model describing unsteady flow to a well in a phreatic aquifer. Figure 3-5 shows a schematic of unsteady flow to a well in a phreatic aquifer. Ho 7/ ////777 //// Figure 3-5. Unsteady flow to a well in a phreatic aquifer at time t. where s = drawdown(L), Ho= initial saturated aquifer thickness (L), Q = pumping discharge (L3/T), and r = radial distance from the well (L). As shown in Figure 3-2, the boundary condition along the phreatic surface is nonlinear and the shape and position of the phreatic surface are unknown. As a result, an exact analytical solution that describes groundwater flow is not known. However, by the employing the Dupuit assumptions and assuming recharge is zero, the unsteady flow to a fully penetrating well in a phreatic aquifer can be described in radial coordinates by: 2h 2 1 9h2 2S 9h ar 2 r ar K at (3-28) Initial condition: h(to) = Ho Boundary condition: h(t) (p) = Ho Here, S is the phreatic aquifer specific yield (or effective porosity, ne). If drawdowns are sufficiently small compared to the average depth of flow, (say 0.02Ho), then the approximate solution to this problem is the Theis solution (Bear 1976): s W(u) 41fT (3-29) where u = r2S/(4T t); W(u) is the well function, and T=KHo= initial transmissivity. The conditions that define the test problem used to validate the numerical groundwater flow model include: a pumping discharge, Q, of 1.5 m3/hr; an initial saturated phreatic aquifer thickness, Ho, of 20 m; a specific yield S (or ne), of 0.1; and a hydraulic conductivity, K, of 15 m/day. For the three-dimensional finite-element model, the specific storativity, Ss, is 0.0001 m-1; where S,=pg(a+n3); n equals porosity; a equals solid matrix compressibility [LT2M-']; g is acceleration of gravity [LT2]; and 3 is fluid compressibility [LT2M-1]. Because of problem symetry, modeling was limited to simulating hydraulic head variations over a wedge section of the aquifer using a 21x11x3 mesh (Figure 3-6). The boundary condition controlled by well water table, the initial water table of well and the pumping discharge were taken as the only input. The pumping pulled water out of the well, on the other hand, the groundwater supplied well. In each simulation time step, based on the known water table in the beginning of time step, and an estimated water table of well, a new discharge supplied from groundwater to well was calculated from the groundwater modeling. Based this new calculated supplying discharge to the well, the pumping discharge, and the circular area of the well, a more accurate well water table could 40 be estimated again. The groundwater flow model was solved again. After many iterations, a convergent and accurate well water table could be obtained before the modeling of next time step. z 20 1 35 5 Figure 3-6. Finite-element mesh of the wedge used to simulate of a well pumping in phreatic aquifer Because the mesh is deformable, the known pressure at the free surface and flux were used to locate the free surface. During calculation, a slice successive over-relaxation (SSOR) method was used as an iterative approach to locate the free surface. The results illustrated in Figure 3-7 validate the coding of the numerical model. 3.4.2 Groundwater Flow around a Circular Wetland To validate the three-dimensional finite-element saturated groundwater flow model again, the numerical model was used to simulate the groundwater flow around a circular wetland and induced by a flood wave in the wetland. It is assumed that the wetland stage fluctuations can be described by a nonperiodic Gaussian fucntion of time as Eq. 3-30: 1 2o2 (3-30) V~lO 0.1 i T....... 0.08 .. .. .0 .- FEM (t=5.0hr) 0.06 - 3 .ana (t=5.0hr) 0.04 . 0.04 FEM (t=2.5hr) 0.02 ana (t=2.5hr) 0 0 5 10 15 20 25 30 35 The radius to well Figure 3-7. Comparison of the numerical and analytically simulated head for the example problem of a transient pump test in a phreatic aquifer Unlike the coupling boundary condition that was used in the last case, a flood wave boundary condition used in this numerical model here. The solid line and the dash line in Figure 3-8 are the three-dimensional finite-element model simulated results. The triangle points and the diamond points are the groundwater analytical solutions from Eq. 2-24 obtained in Chapter 2. The circular wetland which has a radius of 50 meter was assumed. All other system parameters are given in the Table 3-3. The numerical model and analytical solution matched very well. 3.5 A Numerical Model Application to Simulate the Interaction Between a Wetland and an Aquifer The finite-element model was applied to simulate the interaction between a wetland and a contiguous phreatic aquifer stressed by an active productive well. A horizontal, irregular quadrilateral mesh was used to define the irregular boundaries of the wetland and the area beyond the wetland. Figure 3-9 illustrates the numerical mesh corresponding to the 100.9 100.8 -- 3D num d-2m 100.7 ....... 3D num d=5m 3 100.6 A analytical d=2m 100.5 analytical d=5m 100.4 100.3 -- J 100.2 100.1 -s 100 77 7 99.9 2 7 12 17 Time (hr) Figure 3-8. Comparison of the numerically and analytically simulated groundwater table. Table 3-3. Parameters and their values Parameter Value Parameter Value Hoo (m) 100.0 T (m2/h) 0.635 Ar (m) 0.25 S 0.01 At (hr) 0.02 to (h) 5.0 k (m/hr) 0.00635 a 0.4 horizontal and vertical extent of the phreatic aquifer. Beyond the wetland, the vertical extent of the mesh defines the thickness of the phreatic aquifer. Over the inundated portions of wetland, the thickness of the mesh defines the vertical distance between the wetland bottom and the underlying aquiclude. The wetland is shown in the middle of this area. The model domain is 104x104x10 meters. The four lateral boundaries and bottom boundary are assumed to be impermeable. An active production well is located within the model domain. It was assumed that the initial groundwater table and the water stage of wetland were the same. The deepest water depth in the wetland is 1.0 meter. Model runs were conducted to 43 simulate the exchange of water between wetland and aquifer. Such interactions between the wetland and the aquifer were examined under transient pumping conditions and precipitation. Since the real-time stage of the wetland is used as a boundary condition over the wetted part of wetland, the horizontal and vertical extent of the mesh must be adjusted to trace the vertical and horizontal movement of the transient wetland water surface. z Upper-extent of the phreatic aquifer defined x vY by the bottom elevation inundated wetland Top entent of the phreatic aquifer P Mg well 3.9999E+01 .9999E+01 S 9999E0i1 5 9999E-+0-1 7.99 999E01 "E+0I aquiclude 9 9999E+01 9.9999E+01 Figure 3-9. Three-dimensional finite-element mesh for the simulation of wetland and aquifer The modeled area was divided into a 21x21x10 hexahedron mesh. In the middle, a 9x9 quadrilateral mesh was used to represent the wetland area (Figure 3-9). Hydrologic conditions that define the test problem include: a net precipitation rate of 0.02 m/hour; a specific yield, S, of 0.075. a hydraulic conductivity, K, of 10 m/day; and a pumping discharge, Q, of 1.554 m3/hr. Finally, data from Figures 3-10 and Figure 3-11 were used to define the relationship between the wetland water stage, the surface water area of wetland, and the corresponding wetland volume. Figure 3-12 illustrates the initial condition where it was assumed that the wetland surface-water stage and phreatic surface elevation were the same. Figure 3-13 illustrates a 44 plan view of the numerical mesh at time zero. Figure 3-14 through Figure 3-21 represent various depictions of head contours and iso-surface plots of head predicted at hour 5 of the simulation. These figures show that water moves laterally into the wetland and pumping well because the same net rainfall forces the water table to increase faster than in the water stage in the wetland. The directions of flow also depend on the active production well. If production is large, more water flows to the well instead of flowing to the wetland. Otherwise, some water flows towards and into the wetland as groundwater inflow. Beneath the wetland, the hydraulic head increases with time and with increasing depth (Figure 3-17 and Figure 3-18), which shows that water flows from the aquifer into the wetland. As to the pumping well, the contours are almost vertical and parallel to each other near the well (Figure 3-20 and Figure 3-21). Groundwater moves towards the well in a horizontal direction because the well fully penetrates the aquifer. ,700 n 600 ... ........ . CU 500 .... 400 .. .......... CU 300 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 The water depth of wetland D (m) Figure 3-10. Relationship between the area of wetland water surface. 500 --- ,----------j 450 ........ .. .. . ?400 -F - E 3 5 0 -- !- --" > 300 - 1- 250 __......_ _- 200 400 450 500 550 600 650 The surface area of wetland S (s-m) Figure 3-11. Relationship between the deepest water depth and the volume of wetland. N Figure 3-12. Illustration of the initial condition where it was assumed the wetland surface water stage and the phreatic surface elevation were the same 5 9999E+-01 -9.9999E+01 9.9999E+01 4.9999E+01 -5.2005E-04 Figure 3-13. Plan view of the numerical mesh at time zero. x , Figure 3-14. Three-dimensional finite-element mesh and simulated head contours at hour 5. 2 4999EOI I 2005E-04 49E 9999E 01 7 4999E,01 5 9999E+01 7 9999E+01 9 9999E+0-1 9 9999E+01 Figure 3-15. simulated head contour at hour 5. zi Figure 3-16. Simulated iso-surface plot of 10.6 meter head at hour 5. 48 -52005E-04 -2005E-04 19999E4 0 1 9999E+01 3.9999E+01 39999ME 1 5.9999E+01 5.9999E+01 7.9999E+01 7.9999E+01 9.9999E+01 9.9999E+01 Figure 3-17. Three-dimensional finite-element mesh and head contour on a cross section through wetland 0.-046 04 .5 200SE.04 95+2005E.04 Figure 3-18. Simulated head contour on the vertical cross section through the wetland x Y -5 2005E-04 --5 2005E-04 9999E+01 L9999E.o 3.9999E+01 3,9999E+01 5.9999E+01 '.9999EE1I 7.9999E+01 7 9999E+01 9.9999E+01 p9f9999E401 Figure 3-19. Modeled iso-surface plot of head equal to 11.1 m, 11 m and 10.5 m Figure 3-20. Mesh and modeled head contour with vertical cross section through the pumping well A computer with a Pentium 1.4GHz processor took 4 hours to complete the 5-hour 1 901< E.01 59999 oo 9 9S99 99 9EI 9 (Si + Figure 3-21. Head contour (m) on the vertical cross section through the pumping well A computer with a Pentium 1.4GHz processor took 4 hours to complete the 5-hour imulation usimulng a time step of 0.1 hours and 3 iterations per time step. To use this model over a real system, a finer mesh would be desirable. Unlike most rectangular finite-difference grids, the finite-element grid used here precludes the use of more efficient solution algorithms that divide simulation domain into three individual directions. In order to use this model to investigate wetland and aquifer interaction on a large scale, it would be desirable to introduce a simplified mesh that would in turn be more amenable to efficient solvers. With computer hardware developing, the simulation speed can be accelerated consequently. To improve the simulation speed of the model from fundamental, a kind of hybrid method which combines finite-element method and finite-difference method could be an effective choice for a large scale practical problem. 3.6 Application of Three-Dimensional Deformable Saturated Groundwater Flow Model to SV5 Wetland 3.6.1 Site Description-SV5 SV5 is an isolated wetland found in pine flatwoods located in southeast Florida. It is part of the Savannas wetland system located south of the city of Jensen Beach in Martin 51 County. SV5 is located approximately two km west of the Atlantic Ocean and is adjacent to a large municipal well field. The South Florida Water Management District (SFWMD) monitors this wetland system. There are 13 surficial aquifer municipal water-supply wells located from a few hundred meters to four kilometers from SV5 (Figure 3-22). The wetland is roughly circular in shape, with a diameter of approximately 60 meters (Wise et al. 2000, Walser 1998). N A SV5 A A A A AA 0* 0 0 AA SFWMD Well \ / A Monitor Well o Interior Well Figure 3-22. Well Locations (Switt et al. 1998) The surficial aquifer that underlies the wetland is composed of three layers. The first layer, which is located from land surface to 12-18m below the surface, consists of white, gray and brown, largely fine- to coarse-grained quartz sand intermixed with shell beds. It receives recharge directly from precipitation. The second layer, which is 3-6 m thick, is comprised of fine to very fine sand with small amounts of shell and clay. The third layer, which extends to a depth of between 37-43 m, consists of limestone intermixed with sand and shell. The upper confining layer lies beneath the surficial aquifer. This layer consists of low- 52 permeability rocks that are primarily clastic. The Floridian aquifer system lies under the upper confining layer. The system is quite thick and consists of limestone of variable permeability (Fetter 1994). The climate in the area is considered humid sub-tropical, and the area receives an average of 1.5 m of rainfall per year (Walser 1998). 3.6.2 Field Methods and Wetland Aquifer Interaction Test The SV5 wetland aquifer interaction test was conducted (Wise et al. 2000). Eighteen monitoring wells were installed. Wells were placed along transects running from east to west and north to south. Six wells were placed outside the wetland on the western transect and two on each of the other three transects. The wells were identified by the transect designation followed by the distance location from the center of the wetland, i.e., well 37 is the well located 37 m from the center of the wetland along the west transect. SFWMD installed two wells and a staff gage at the site. The wells measure groundwater. The groundwater well is screened at a depth of six meters and grouted from the surface to a depth of five meters. The water and peat depths in the wetland were measured along four transects through the center of the wetland at 450 angles from one another. Peat and water depth measurements were taken at three-meter increments. The peat depth was determined by pushing a piece of 1 cm rebar into the peat at a steady rate until significant resistance was felt. This method is most easily applied to sites such as SV5 that are free of significant tree roots. The depth of the water table was determined by measuring the water level in the exterior monitoring wells. Elevation data for the sites were measured utilizing a Topcon Self-Leveling Rotating laser level. The South Florida Water Management District staff gage on the site was used as a relative reference for the elevation data, and all elevation data were expressed in terms of NGVD of 1929. 53 The water and peat depth information was entered into Surfer, a three-dimensional graphical interpolation program, to develop three-dimensional contours of the wetland. A stage -volume curve to describe the wetland was developed by utilizing Surfer to determine the wetland volume at varying water levels. A regression curve was fitted to this data to determine the equation that describes the stage-volume curve (Switt et al. 1998, Wise et al. 2000). This curve describes the top of the peat layer. The process was repeated to develop a curve that describes the bottom of the peat layer. The goal of the test was to stress hydraulically a wetland by pumping the standing water to lower the water level a predetermined amount, which varies from wetland to wetland, and then monitoring the response in wells and subsequent recovery of the water in the wetland. For the wetland test to be effective, the wetland must be a discharge wetland at the end of the pumping phase. The wetland test was conducted utilizing a 10.2 cm centrifugal trash pump to pump the water from the wetland until the water level in the wetland dropped 30.5 cm. The water was pumped into a drainage ditch located approximately 100 m from the wetland. This distance was chosen to ensure that the water pumped from the wetland did not reenter the system. The discharge rate of the pump was approximately 70 m3/hr. The discharge rate was measured by determining the amount of time it took to fill a five-gallon plastic bucket at the end of the discharge hose. It took nine hours to lower the wetland water level 30.5 cm. During the testing, water levels in all of wells were manually measured every hour, with the exception of the surface water well, which was measured every 30 minutes, utilizing a Slope Indicator water tape. Once the pump was turned off, water levels in all of the wells were manually measured every hour for 22 hours. Exterior and interior well water levels 54 were recorded for the wells in the wetland area. 3.6.3 Stage-Volume Relationships Water and peat depth data were used to produce a three-dimensional contour plot of SV5 (Figure 3-23 and Figure 3-24). The regression equation that describes the volume of water above the peat layer is V =3454S2-148.45s (3-30) where Vw is the wetland volume (m3) at s, the stage level (m) (Switt et al. 1998, Wise et al. 2000). The bottom of the peat layer is described by the follow regression equation (Figure 3.25): V =-116s4 +607s3 +938s2+108s (3,31) where Vp is the volume above the bottom of the peat layer (m3); and Vp-V, volumee of saturated peat. The coefficient of determination equals unity for both of these regressions. Figure 3-23. Top surface of peat layer of wetland Figure 3-24. Bottom surface of peat layer Stage-volume curve 1600 1400 -- - 1200 - 1000 800 600 400 200 0 0 0.1 0.2 0.3 0.4 Stage Vw~S Vp-~S 0.5 0.6 0.7 Figure 3-25. Relationship between wetland volume and water stage 56 Switt et al. (1998) used a similar formula (Eq. 3-32) to calculate the interaction between a wetland and the underlying groundwater: SP-K G-W) SA (3-32) where SP= exchanging flux (m3/day), K,= hydraulic conductivity (m/day), G = groundwater level (m), W = wetland water level (m), L = average peat thickness (m), and SA, = wetland wetted surface area (m2). Over the de-watered area of the wetland, seepage decreases from the above value at the edge of the water to zero at the edge of the wetland. Switt et al. (1998) assumed a standing water head difference equal to 1/3 of the value use in Equation 3-8. Thus, seepage for the de-watered area of the wetland was SP=Ki( W SA, (3-33) where SArw= wetland dry surface area (m2). The hydraulic conductivity of the organic soil in the peat layer is critical in wetland systems because it controls the amount of groundwater that flows from the aquifer into the wetland or from the wetland into the aquifer (Truesdell 1995). Determining a value of hydraulic conductivity for the soils in wetland system is particularly challenging. 3.6.4 Creation of Finite-Element Mesh The model domain was approximately three times the size of the wetland. Because the wetland is approximately radically symmetrical it was possible to develop a highly refined 4.00 numerical wedge to represent the wetland aquifer system. The aquifer was 57 divided into 10 layers to obtain a 11x26x3 mesh (Figure 3-26). The peat layer was divided into three layers producing a 3x11x3 mesh. The mesh was divided and fitted along the shape of peat layer. z -Wtland Pel Layer 15 5 40 20 0 0 Figure 3-26. Finite-element mesh of wedge of wetland and aquifer for the saturated flow model 3.6.5 Solution Procedure To simulate the wetland-aquifer interaction observed at SV5, the model must correctly simulate the movement of the free surface as changes occur with the wetland water stage. Two boundary conditions must be satisfied on free surfaces; first, there can be the flux across the surface must be represented, and, second, the pressure must be atmospheric. The initial fixed domain solution usually will not satisfy both boundary conditions simultaneously. Therefore, the trial free surface is adjusted and a solution obtained for the revised fixed domain. This process is repeated until the required adjustment in the free surface is considered negligible. Thus, the solution of the non-linear free surface problem is achieved by a series of solutions to linear problems with prescribed boundaries. Finally, because the wetted part of wetland may change significantly with variations of wetland water stage, the finite-element mesh must be adjusted horizontally and vertically between time steps. Before the pumping test began, the SV5 wetland water surface was about 60 meters in diameter. However, after pumping stopped it was less than 30 meters. Thus, the horizontal inundated area changed from 2750 m2 to 750 m2. To model the influence of a changing horizontal free surface required that the finite-element mesh be adjusted horizontally. The zone of mesh adjustment would be within the peat layer and include an upper boundary representing the free surface or the surface water stages and a lower boundary corresponding to the vertical limit of the peat layer. Adjusting the mesh horizontally keeps the same number of nodes on the boundary. Thus, the physical domain of the mesh may change horizontally and vertically. However, the calculation domain does not. After each simulation, information defining the location of every node is updated. Since the wetland boundary was assumed to vary over time, the total head was taken as a boundary condition. Thus, at any time, the wetland boundary can be viewed as an equipotential line. Water moves vertically to and from the equipotential line represented by the free surface boundary. The vertical hydraulic gradient at the wetland boundary times the product of the wetland area and hydraulic conductivity equals the volume of water exchanged between the wetland and the aquifer. 3.6.6 Model Calibration and Results Because the aquifer underlying SV5 is composed of sand, a specific storativity of S = 10-5 m-1 was used (Bear 1979). The saturated water content and residual water content in the peat layer were 0s = 0.7 and Or = 0.6. In the aquifer, 6, = 0.34 and Or = 0.18, based on Walser (1998). Water stage level observation data and groundwater monitoring data collect during the pump test and many hours after pumping stopped were used to calibrate the numerical model. In Krabbenhoft's (1996) interaction model between groundwater and lakes, an anisotropy ratio of 100 was obtained after several trial and error simulations. Also an anisotropy ratio of horizontal to vertical hydraulic conductivity of 10 was obtained by Winter (1978) in his numerical simulation of steady state three-dimensional groundwater flow model near lakes. In this finite-element model, it was assumed that the hydraulic conductivity of the peat layer was different from the hydraulic conductivity of the underlying aquifer, and the aquifer was also assumed anisotropic. Furthermore, it was assumed that the hydraulic conductivities and anisotropies of both layers were spatially uniform within each layer throughout the aquifer zone. By adjusting the anisotropy ratio (Kh/Kv) and checking for the best agreement between modeled and observed data, it was determined that the best results were obtained when Kh/Kv= 1.0, Kxx= Kyy = 8.7 m/day, and Kzz = 8.7 m/day in the aquifer, while Kp = 0.2 m/day in the peat layer. This number is different from the value of 0.03m/day obtained by Switt (1998) and 0.087m/day obtained by Wise (2000). The possible difference may be caused by several reasons. First, the peat layer depth used in the numerical model is not a constant. When the water stage in the wetland is lowered, the peat layer depth in the inundated area is much greater than the average value of 0.25 m used in the analytical model. The variation of groundwater head immediately beneath the peat layer may not be the same value estimated using a groundwater head variation at a depth of 6 meters beneath of the wetland. Thus, the second reason may be that the head distribution beneath the wetland obtained from the numerical model is different from the values used by the analytical model. Since the aquifer beneath the SV5 wetland is composed of sand, the specific 60 storativity is extremely small. The heads along the depth are very different. Figure 3-27 shows the head distribution beneath the wetland. Figure 3-28 shows the head variation with time at different depths beneath the wetland. The piezometric head changes relatively quickly in the deep zone, but the head does not change much near the wetland. Figure 3.29 shows the variation of the inundated area of the wetland. The area changes from 2,600 m2 to 700 m2just after pumping. Figure 3-30 shows the variation of the interaction discharge. The discharge supplied from the aquifer to the wetland increases quickly and then decreases slowly after pumping was stopped. Figure 3-31 shows excellent agreement between the simulated results given by the saturated flow model and observed groundwater levels. Table 3-4 to Table 3-8 show the modeled data used to calculate the interaction discharges at 5 different times. Table 3-8 shows the sensitivities of the head relative to changes in the hydraulic conductivities in the pear layer and the aquifer and to changes in the specific storativity. The hydraulic conductivities in both the peat layer and the aquifer play the most important part in the change of head distribution and interaction discharge. The sensitivity of the specific storativity is relatively small. 61 Lx -50 0 50 Figure 3-27. Vertical cross-section view of three-dimensional finite-element mesh and head from the saturated flow model at a time of 87 hours. ._ 15.9 E . 15.8 15.7 0 20 40 60 Time (hour) 80 100 Figure 3-28. Piezometric head at positions beneath the wetland with time 2.5 2.5 ---------------------- S1.5. 0 S0.5 .......-- .. ........ .. 4.-^-.......-- 0 20 40 60 80 Time (hour) Figure 3-29. Variation in the inundated area of wetland during simulation 3000 < .... .. ...... ............... i .. .. 2500 --- S2000 - S1500 1000 - 500 I- 0 20 40 Time (hour) 60 80 Figure 3-30. Discharge between wetland and aquifer with respect to time 0 20 40 60 Time (hour) WETL OBS WETL CAL G-W OBS G-W CAL 80 100 Figure 3-31. Comparison between calculated head and observed head Table 3-4. Data for calculating the discharge to the wetland from the aquifer at 10 hours Elements Area Head in Peat Head on Peat Distance of Two of Peat Layer (m2) Layer (m) Surface (m) Nodes (m) 1 6.95828 15.7810 15.7223 0.301929 2 20.8748 15.7813 15.7223 0.305786 3 34.7914 15.7819 15.7223 0.306822 4 48.7075 15.7827 15.7223 0.303675 5 62.6240 15.7843 15.7223 0.297326 6 76.5414 15.7868 15.7223 0.289907 7 90.4590 15.7898 15.7223 0.281774 8 104.377 15.7936 15.7223 0.271725 9 118.296 15.7979 15.7223 0.260863 10 132.215 15.8028 15.7223 0.248279 Q to wetland at 10 hr 1.48466 mn/hr Table 3-5. Data for calculating the discharge to the wetland from the aquifer at 20 hours Elements Area Head in Peat Head on Peat Distance of Two of Peat Layer (m2) Layer (m) Surface (m) Nodes (m) 1 8.36449 15.7953 15.7427 0.302115 2 25.0935 15.7956 15.7427 0.306344 3 41.8223 15.7960 15.7427 0.306027 4 58.5508 15.7964 15.7427 0.302183 5 75.2802 15.7982 15.7427 0.294043 6 92.0102 15.8003 15.7427 0.287151 7 108.741 15.8036 15.7427 0.275314 8 125.472 15.8070 15.7427 0.264839 9 142.204 15.8114 15.7427 0.250437 10 158.937 15.8162 15.7427 0.236900 Q to wetland at 20hr 1.65871 n3/hr Table 3-6. Data for calculating the discharge to the wetland from the aquifer at 40 hours Elements Area Head in Peat Head on Peat Distance of Two of Peat Layer (m2) Layer (m) Surface (m) Nodes (m) 1 10.9288 15.8222 15.7797 0.302417 2 32.7864 15.8225 15.7797 0.307252 3 54.6436 15.8225 15.7797 0.304707 4 76.5010 15.8229 15.7797 0.298378 5 98.3597 15.8243 15.7797 0.289097 6 120.219 15.8264 15.7797 0.277635 7 142.080 15.8288 15.7797 0.266032 8 163.941 15.8321 15.7797 0.249581 9 185.804 15.8359 15.7797 0.233460 10 207.669 15.8406 15.7797 0.216548 Q to wetland at 40hr 1.86649 m'/hr Table 3-7. Data for calculating the discharge to the wetland from the aquifer at 60 hours Elements Area Head in Peat Head on Peat Distance of Two of Peat Layer (m2) Layer (m) Surface (m) Nodes (m) 1 13.1219 15.8463 15.8115 0.302649 2 39.3657 15.8466 15.8115 0.307800 3 65.6089 15.8464 15.8115 0.303651 4 91.8530 15.8469 15.8115 0.295196 5 118.098 15.8480 15.8115 0.285844 6 144.345 15.8498 15.8115 0.271268 7 170.593 15.8519 15.8115 0.255422 8 196.843 15.8544 15.8115 0.239015 9 223.095 15.8579 15.8115 0.219572 10 249.349 15.8626 15.8115 0.199788 Q to wetland at 60hr 1.94458 m'/hr Table 3-8. Data for calculating the discharge to the wetland from the aquifer at 80 hours Elements Area Head in Peat Head on Peat Distance of Two of Peat Layer (m2) Layer (m) Surface (m) Nodes (m) 1 15.0334 15.8681 15.8391 0.302835 2 45.1003 15.8683 15.8391 0.307522 3 75.1664 15.8680 15.8391 0.302700 4 105.234 15.8685 15.8391 0.292596 5 135.303 15.8695 15.8391 0.280887 6 165.374 15.8709 15.8391 0.266576 7 195.446 15.8727 15.8391 0.247805 8 225.521 15.8749 15.8391 0.228349 9 255.599 15.8777 15.8391 0.208494 10 285.679 15.8819 15.8391 0.186883 Q to wetland at 80hr 1.93622 m3/hr 66 Table 3-9. Sensitivities calculated by saturated groundwater flow model Parameter (PR) Percentage varied APR I hi-h211 (m) II hi-h211/APR Ka(m/day) 0.10 -0.85 0.9924E-02 -1.1675E-02 Ka(m/day) +0.10 +0.85 0.9055E-02 1.0653E-02 Kp(m/day) -0.10 -0.02 0.6410E-02 -0.3205 Kp(m/day) +0.10 +0.02 0.5564E-02 0.2782 So -0.50 0.000005 0.2068E-04 -4.136 So +0.50 +0.000005 0.1265E-04 2.53 3.7 Comparison with Results of Wise's (2000) and others Parameter values obtained by Wise (2000), kpeat= 0.087 (m/day), by Walser (1998), kand = 0.75 (m/day), kpeat = 0.03 (m/day) for SV5 were put in the three-dimension finite- element model in this dissertation, and the simulated results were plotted on the Figure 3-32. When the specific storage was set on 0.01 (1/m), both the simulated wetland water table and groundwater are pretty close to the observed conditions. The simulated wetland water table and the groundwater table corresponding to the specific storage of 0.001 (1/m) are much different from the observed results. There are two main differences between Walser's numerical simulation and the models in this dissertation. First, unlike a coupling boundary condition used in this dissertation, the fixed observed water table process over time was used in Walser's numerical simulation. The other difference is that a fixed three-dimension mesh was used in Walser's numerical simulation. When kaqu at 8m/day and kpeacat 0.2m/day were used, and specific storage was varied, the simulated results are plotted on the Figure 3-33. The simulated wetland water tables for both specific storage values are very close to the observed wetland water table results. The 67 groundwater table for the specific storage of 0.01(1/m), however, is above the observed results. 4.3 - 4.25 S 4.2 --. -- gwobs S4.15 .- ----gw sim(Ss=0.01) .) 4.1 ..- ''' gw sim(Ss=0.001) S4.05 -- wetlobs ...-----wetilsim(Ss=0.01) 4 S- -... wet sim(Ss=0.001) 0 3.95 ,,-_ _ 3.9 __............. ... - 3.85 - 0 10 20 30 40 50 60 70 80 90 Time (hr) Figure 3-32. Comparison results using Walser's parameter values. Figure 3-33. 3-D numerical model simulated results over specific storage Table 3-10 shows the hydraulic conductivities in sand layer and peat layer. The variations are rather large. 68 Table 3-10. Hydraulic Conductivity values for sand and peat layers k in sand layer k in peat layer (m/day) (m/day) Switt et al.(1998) 7-12 0.03 Slug test 2 Rycroft et al. (1975) 10-7-102 Wise et al. (000) 0.083 Walser Numerical Model (1998) 0.75 0.03/0.3 Models in this dissertation 8 0.2 3.8 Summary A three-dimensional finite-element numerical saturated groundwater flow model using deformable mesh and real-time wetland-aquifer coupling boundary is developed in this chapter. Considering the area of the water surface changes quickly because most wetlands have a shallow bank slope, the model has the capability of tracking the interaction the interaction boundary between groundwater and a wetland. The model is verified using two analytical solutions, one is Theis solution, the other is the analytical solution found in Chapter 2 of this dissertation, which describes groundwater flow around a circular wetland. The model is further used to simulate the interaction problem of groundwater and the SV5 wetland. The pumping from the SV5 wetland and the recovery of the SV5 wetland surface afterwards are closely simulated by this model, the interaction boundaries are effectively traced by the deformable mesh. CHAPTER 4 INVERSE MODEL OF THREE-DIMENSIONAL FINITE-ELEMENT METHOD SATURATED GROUNDWATER MODEL IN SEARCHING FOR THE HYDRAULIC CONDUCTIVITY OF THE PEAT LAYER 4.1 Introduction 4.1.1 Inverse Problem During the past twenty years, mathematical modeling has been extensively used in the study of groundwater problems. However, mathematical models require a number of parameters to characterize the aquifer. Usually, prior knowledge of the values of these parameters is limited, and they must be quantified by calibration which requires observations of hydraulic heads or solute concentrations. Given a parametric model of the physical system and values of the model parameters, the prediction of system output (response) to any input is referred to as the direct problem, or the forward problem (simulation). The forward problem predicts unknown system states by solving appropriate governing equations. Conversely, the estimation of model parameters given the parameteric system model and an input-output relationship is known as the inverse problem (calibration). The typical inverse problem involves estimating the values of model parameters, solving the direct problem with this set of parameters to calculate the predicted system response, and then adjusting model parameters until deviations between observed and model predicted responses are in some specified sense minimized (Sun 1994). 70 Trial and error and graphical matching techniques are the most basic and traditional methods to solve inverse problem. They are widely used because they are simple to implement. Because the number and combinations of parameter adjustments are not bounded, "trial and error" is rather flexible but also time-consuming. Additionally, the solution is strongly dependent on the skill of the practitioner (Keidser 1991). 4.1.2 Parameterization Most groundwater inverse algorithms adopt either a blocked or a geostatistical (random field) description of spatial variability. The first approach divides the region of interest into a number of discrete blocks or zones that are believed to correspond to distinct geological units (Yeh 1986, Carrera and Neuman 1986a,b,c, and Cooley 1977,1979,1982,1983). Each block is characterized by a set of spatially uniform hydrogeologic properties which are treated as parameters in an appropriate inverse problem. Thus, the unknown distributed paramters can be represented by a number of constants. The weakness of this method is obvious, i.e., although we have some of the geological information is available, it is still difficult to divide the flow region into a limited number of zones. If the zonation pattern is incorrect, the identified paramter values also will be incorrect (Sun and Yeh 1985). The geostatistical alternative views the properties of interest as stationary random fields that vary relatively smoothly over space (Hoeksma and Kitanidis 1984, Dagan 1985,Yeh 1996, and Zhang 1997, Michalak 2004). In the geostatistical interpolation method of kriging, the unknown parameter is regarded as a stochastic field described by some statistical parameters, such as the mean, variance, and correlation distance. Then the identified parameter vector consists of only a few statistical parameters if the aquifer is 71 basically homogeneous. The problem of overparameterization may be avoided, and the inverse solution obtained by the maximum likelihood estimate and cokriging is generally stable. The correctness of the inverse solution, however, depends on the correctness of some statistical assumptions and the structure of covariance functions. Although the two approaches are based on different parameterizations, they both treat hydrogeologic properties as spatial functions. This suggests that it should be possible to formulate a general inverse theory which encompasses both the blocked and geostatistical alternatives, as well as hybird methods which lie between these extremes. 4.1.3 Parameter Identification Methods Many methods are available for the solution of the common groundwater modeling inverse problem whereby the conductivity or transmissivity of a porous medium are estimated using available observations of head or pressure and other information (Distefano 1975, Chang 1976, Loaiciga 1987, Mishra 1989, Knopman 1989, Keidser 1991, Sun 1995, Yeh 1996, Zhang 1997, Lamorey 1998, Anderman 1999, Michalak 2004). Neuman (1973) classified inverse modeling solution techniques as either "direct" or "indirect". The direct method or equation error method is used when head variations and derivatives (usually estimated) are known over the entire flow region. If the measurement and model errors are negligible, the original governing equation becomes a linear first-order partial differential equation of the hyperbolic type in terms of the unknown parameters. With the aid of boundary conditions and flow data, a direct solution for the unknown parameters may be possible. In practice, observation wells are sparsely distributed in the flow region and only a limited number of observation wells are available. To formulate the inverse problem by the equation error criterion, missing data have to be estimated by interpolation. The 72 interpolation data contain errors resulting from the interpolation, and, thus, cause additional errors in the results of parameter identification. The solution approach of the indirect method to solve inverse problems is to formulate the inverse problem in term of an optimization framework, based on the most common output error criterion, i.e., output error least squares criterion. Using this approach, a set of parameters is identified which best describes observed system states. With inverse modeling, the structure and parameters of a model are adjusted simultaneously or sequentially so as to obtain a parameter input system and output relations that characterize observed excitation-response functions of the real system. Since the model structure for the wetland/groundwater system is assumed to follow Eq.5-1-5-3, the problem of determining the hydraulic conductivity of the peat layer from the observed water stage and total head in the aquifer may be taken as a typical inverse problem. The investigation of the inverse problem based on the three-dimensional saturated groundwater flow model and the adjoint state method is presented in this chapter. 4.1.4 Computation of Sensitivity Coefficients Sensitivity coefficients, i.e., the partial derivatives of head with respect to each of the parameters, play an important role in the solution of the inverse problem. In the Gauss- Newton algorithm, elements of the Jacobian matrix are represented by the sensitivity coefficients, ahi/aTI, i=1,....,M, 1=1,...,L. If h is the head vector, the sensitivity coefficients are 8h,1/8T, 1=1 ,...,L. There are three methods used in calculation of sensitivity coefficients. There are the influence coefficient method, sensitivity equation method, and variational method. The influence coefficient method determines the sensitivity coefficients by perturbing each parameter once at a time. 73 In the sensitivity equation method, a set of sensitivity equations is obtained by taking the partial derivatives with respect to each parameter in the governing equation and in the initial and boundary conditions. The sensitivity coefficients are obtained by solving the new governing equations. The variational method was first used for solving the inverse problem of parameter identification by Jacquard and Jain (1965) and then by Carter et al. (1974, 1982) associated with finite-difference schemes. Sun and Yeh (1985) extended the method to the case of a finite-element scheme. The sensitivity coefficients can be computed using the following equation: S q(xt-) Vh(xy,) ddyd (4-1) 0 0, j=1,2,..., No i=1,2,...,N, where Q is the exclusive subdomain of node i; V is the gradient operator; h(x,y,t) is the solution of the governing equation; No is the number of observation wells; Nn is the total number of nodes used in the numerical solution; and q(x,y,t) is the time derivative of q(x,y,t), which is the solution of the set of adjoint equations. 4.2 Adjoint Problem for Three-Dimensional Finite-Element Saturated Groundwater Model The focus of this section is the formulation of the adjoint problem which defines sensitivities to be used in determining the hydraulic characteristics of the peat layer and the underlying aquifer. It is assumed that the thickness of the peat and aquifer are known. As 74 discussed in Chapter 3, three dimensional groundwater flow in an unconfined aquifer is governed by: V (KVh)+ Q = S,- (xyz)eO, Ot: tf (4-2) subject to the initial and boundary conditions: h(x,y,z,0) = fo (xy,z)EF (4-3) hir, =f, KVh-nlr =f2 (4-4) fo is the initial head; f, is the specific head condition; f2 is the specified flux condition; and t f is the time duration. The boundary of volume Q shown in Figure 4-1 is denoted by r including wetland head boundary (r2), no flow boundary (rF), and lateral head boundary (r2). E A B F -- I Peat Layer Groundwater Aquifer G I HI AB, EG, FH or CD are head boundaries r2 GH or EC, DF are no flow boundaries r Figure 4-1. Sketch showing setup of the saturated groundwater inverse problem around a wetland 75 The hydraulic conductivities of the peat layer and the underlying aquifer are the unknown parameters. Any variation 6K of K where K=hydraulic conductivity must cause a variation 6h of head h, because h is dependent on K through the governing equation and subsidiary conditions. Taking the variation of Eq. 4-2 and subsidiary conditions, the following is obtained: v.(6 KVh) + V *( KV6h ) = S, (xvz)Q, Ot5tf (4-5) subject to the conditions: 6h(x,y,z)l o = 0 (xy,z)eQ (4-6) 6hIr2 = 0, [6KVh+KV8h]-nr = 0 (4-7) Eq. 4-5-Eq. 4-7 represent the variational problem of the primary problem in Eq. 4-2. The variational problem relates the variation 6k and 6h. Assuming h is an arbitrary function having continuous second-order space derivatives on 0 and first-order derivatives in [0, tf], then multiplying Eq. 4-5 by and integrating the result over time-space domains yields: |(ST68h)I d f [S_6h + TV-(8kVh) + TV'(kV6h) dOdt = 0 (4-8a) 0 00 If we define: V(x,y,z,tf) = 0 (x0y,z)eQ (4-8b) then the first term of Eq. 4-8a vanishes because of condition 4-7. Thus, Eq. 4-8b reduces to: h a + TV-(6kVh) + V.-(kV6h) ]ddt = 0 (4-9) 0 Q where: YV-(8kVh)dQ = T(8kVh)-nds - 8kVh-VWdV 'V-(kV8h)dD = Y(kVbh)nd - - kV8h-Vfd = I a a 8h[V-(kVT)]ad - 6hkVWTn ds Using identities of Eq. 4-10, Eq. 4-11 and Eq. 4-12, Eq. 4-9 can be written as: and: (4-10) (4-11) (4-12) ffS36h + h[V-(kVT)]-65kVh-VIY d dt 0 0 + f If6kVh*nds+ fIkV6h-nds f 8hkVW-nds dt= 0 (4-13) 0 s s s Using conditions 4-7 and defining: Yr2 = 0, kV'-nlr, = 0 then the second term of Eq. 4-13 is equal to zero. Eq. 4-13 produces: [s +V-(kVT) dhd} dt (V .-Vh)dkdOdt = 0 (4-14) 00 0 0 In order to identify the hydraulic conductivity k(x,y,z) that would lead to the best prediction of observed heads, a performance criterion must be defined. A general form of the performance function may be written as: E(h,k)= fh,k;x,y,z,t)dDdt = 0 (4-15) 00 a Taking the variation of E(h, k) yields: I E(h, k)= 0 f 0 a The summation of Eq. 4-16 and Eq. 4-14 gives: - s5 -V(kVW) dhdQdt + (4-17) S ++VT-Vh dkdQdt 0 a If W(x,y,z,t) is selected to satisfy the following PDE: S -t +V-(kVT)- -a=0 a T h subject to conditions: (4-18) t ((x,y,z,t ) = 0, r = 0, T Ir, " kVT'nlr =0 2* then, the first integral on the right hand side of Eq. 4-17 related to unknown 5h will vanish. Thus: f68k+ f 6h/d dt = 0 9k 9h ) (4-16) 8E = 0 0 t 8E= ff( +VY-Vh dkdQdt 00 From Eq. 4-19, the following partial derivative can be defined: 6E Ak I f f-+VT-Vh d(dt (8k ) 00 By introducing the following time transformation: T = tf- t the adjoint problem can be expressed as i -_ V-(kVY) + at subject to: I lr=o=0 kV'nr, =0 where T = J(x,y,z,tf -T). Solving each equation once, all components of the gradient VE(k) can be obtained. (4-19) (4-20) (4-21) (4-22) (4-22a) S1|=o=0, (xy)EQ 80 4.3 Searching Objective Function To find the best parameter value, the model structure and model parameters are solved simultaneously or sequentially so as to make the input-output relation of the model fit any observed excitation-response relationship of the real system. Since the model structure of the problem in this study is determined, the problem of determining the hydraulic conductivity of the peat layer and the underlying aquifer from the observed water stage and total head in the aquifer may be taken as a typical inverse problem: E(Kest) =min E(K) (4-23) where: L E(k) = [ha-hobs (4-24) 1=1 The traditional method to solve this kind of problems is using the trial and error method based on the output least squares criterion (Eq. 4-23). 4.4 Adjoint State Controlling Equations In this study the adjoint state method was used to solve the inverse problem. From the governing equation (Eq. 4-2) and the associated boundary conditions, the adjoint formulation of the primary problem may be found using Green's first and second theorems and other transforms: S -V(k VT)+ -f =0 s rT h (4-25) subject to: TPL=0=0, (x,y)eO 1r2=0=0 kV'nIr, =0 where h is the total head; T = T (x,y,tf-T) is the adjoint state of h; tf is the time duration of the simulation criterion and f(h,k;x,y,t) expresses the 'goodness' of fit between the predicted head, h, and the observed head, hobs expressed through the ordinary least squares (OLS): f(h,k; x,y,z,t) = 11 h(k) -hob o6(X-Xobs)8(y-y obs) (Z-Zobsb)6(t-tj) i=1 (4-28) where ( xobs, Yobs, Zobs) is the location of the observation well; and t, (j =1,2,---,J) are the observation times. Integrating, Eq. 4-28 over time space domains produces: E(h,k)= f(h,k;x,y,z,t)dQdt (4-29) By solving the governing Equation 4-2 and the adjoint Equation 4-25 once, we may obtain the following functional partial derivative required to solve Equation 4-23 is obatined: and (4-26) (4.27) E f f 1Vhddt (4-30) 4.5 Searching Method and Stopping Criterion Since the gradient VE(k) is known, many searching methods including the gradient method, quasi-Newton method, and the Fletcher-Reeves method could be used to solve the inverse problem 4-23 (Sun 1985). Two stopping criteria used here are: IIK,-K-1il and: IIE(kn)-E(k,1) 1 Using the adjoint method and the Fletcher-Reeves algorithm to update the search sequence, a highly efficient numerical scheme for solving the inverse problem can be formulated. 4.6 Field Application on SV5 The inverse model was used to estimate the hydraulic conductivity of the peat layer underlying the SV5 wetland (Wise et al. 2000). Since the piezometric head at a depth of 5 meters beneath the wetland was observed, f is defined as: J f(h,k;x,y,z,t) = I h (k) hobI 8(x-xobs)(y-yobs)6(z-zbs)8(t-tj) (4-33) Water-level data were measured every hour for a total duration, of 87 hours. An example Water-level data were measured every hour for a total duration, tf, of 87 hours. An example 83 of how the hydraulic conductivity, K, is updated between iterations is as follows. Assuming an initial estimate of the vertical hydraulic conductivity K of the peat layer was K,=1.00m/day, then aE/9K = -0.0125 by Eq.4-29 and E(K,)= 0.48998. Because it is not possible to match the modeling result with the observed result exactly, the final E(k) is not zero. The sensitivity was used to decide the direction of searching. Assuming E(K) is supposed to be the smallest value, the updated value of K=K2 is assumed to be 0.6 m/day as generated from the formula: E(K) = E(K) + K, (K2 K, ) (4-34) Based on the new K = K2 = 0.6 m/day, the updated gradient 9E/dK = 0.0188, and E(K2)= 0.41467. Then a new K is assumed until E(K) is minimized. The final K is 0.2 m/day. Figure 4-2 and Table 4-1 are the calculated values from the adjoint method model. In this model, the observed groundwater data were used in the output least squares criterion. The observed water stage in the wetland was not used in the model. This is the major reason that the final result from the inverse model is a little bit different from the result obtained in Chapter 3. The saturated flow model in Chapter 3 uses both the observation of groundwater head beneath the wetland and the surface water stage of wetland to calibrate the model. In the inverse model, however, it is hard to use the surface-water observation result, since it is a boundary of the model. Because the wetted area of surface water associated with the wetland changes with time, the space location of each finite-element point tied to the location of the water surface must be saved before the adjoint state equations are solved. Thus, for each assumed hydraulic 84 conductivity in the peat layer, the forward problem, and adjoint state problem must be solved first, and then aE/ak for this assumed value of K could be found. Figure 4-3 shows the procedure for solving an inverse problem using adjoint method. 0 0.2 0.4 0.6 K in peat layer (m/day) 0.8 1 Figure 4-2. Relation between error and peat hydraulic conductivity Table 4-1. Error and sensitivity based on the adjoint method model K (m/day) E(k) aE/aK 1 1.0 10.5698 0.26965 2 0.6 8.94509 0.40555 3 0.4 7.17941 1.03330 4 0.2 3.47742 0.31495 5 0.1 4.82392 3.72980 6 0.08 5.59932 3.14520 |

Full Text |

PAGE 1 TWO-DIMENSIONAL ANALYTICAL AND THREE-DIMENSI ONAL FINITE-ELEMENT METHOD MODELING OF THE INTERACTIONS BETWEEN WETLANDS AND GROUNDWATER By JINHUA ZHONG A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF DOCTOR OF PHILOSOPHY UNIVERSITY OF FLORIDA 2004 PAGE 2 TABLE OF CONTENTS ACKNOWLEDGMENT . . . . . . . . . . . . . . . . . . . . . . . . 11 LIST OFT ABLES ...................... ................................ V LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . vi LIST OF SYMBOLS .. ...... .... ............ . .... . ................... i x ABSTRACT ........................... .............................. xii i CHAPTER .............................................. ............. 1 INTRODUCTION ...... . .......... ................................... 1 1.1 Wetland .... . ............................................... 1 1.2 Wetland Water Budget .................. ............... .......... 2 1 3 Goals and Objectives .......... ....... .......... ............... 4 2 ANALYTICAL GROUNDWATER MODEL RESPONSE TO THE PERIODIC OR NONPERIODIC WATER ST AGE FLUCT U ATION IN A CIRCULAR LAKE OR WETLAND . ............ . ............ . ............. 7 2 1 Introduction .. ................................... ....... .... ... 7 2.2 Problem Formulation: Groundwater Flow Around a C i rcular Wetland ..... 8 2.3 Analytical Solution for the Periodic Boundary Condi t ion ............... 9 2.4 Analytical Solution for the Non-periodic Boundary Condition ........... 1 2 2.5 Application ............................. ..................... 1 4 2 6 Comparison of Analytical Solution and Num e rical Solutions ....... .... 1 7 2 7 Summary ........ ............... ............................. 18 3 TRANSIENT THREE DIMENSIONAL DEFORMABLE FINITE ELEMENT SATURATED GROUNDWATER MODEL ........................... 19 3.1 Introduction . ...... ...... . . ............ ................ . 19 11 PAGE 3 3.2 Groundwater Flow Models Based on Finite-Difference and Finite-Element Approximation ............................. 21 3 3 Three-Dimensional Deformable Finite-Element Saturated Groundwater Model ............. ......... ......... 24 3.3.1 Governing Equation and Boundary Conditions ........ . . . 27 3 3 2 Finite Element Model of Saturated Groundwater Flow ......... 28 3.3.3 Adaptation of the Saturated Finite-Element Model to Describe the Free Surface ............... .............. 36 3.4 Validating the Code of the Three-Dimensional Deformable Finite-Element Groundwater Flow Model ........ ............. 3 7 3.4.1 Theis Solution ........................................ 37 3.4.2 Groundwater Flow around a Circular Wetland ............... 41 3.5 A Numerical Model Application to Simulate the Interaction Between a Wetland and an Aquifer ........................ ....... .... 42 3.6 Application of Three-Dimensional Deformable Saturated Groundwater Flow Model to SV5 Wetland ...................... 50 3.6.1 Site Description-SYS .............. .................... 50 3 6.2 Field Methods and Wetland Aquifer Interaction Test ........... 5 2 3.6.3 Stage-Volume Relationships .................... .......... 54 3.6.4 Creation of Finite-Element Model Mesh ....... ............. 56 3.6.5 Solution Procedure .... ............ ..... .... .......... 5 7 3.6.6 Model Calibration and Results ................ . . . .... 58 3.7 Comparison with Results of Wise and others . .................. . 69 3.8 Summary ....................... ............................. 7 1 4 INVERSE MODEL OF THREE-DIMENSIONAL FINITE-EIEMENT METHOD SATURATED GROUNDWATER MODEL IN SEARCHING FOR THE HYDRAULIC CONDUCTIVITY OF THE PEAT LA YER ....... 69 4 1 Introduction .................................................. 69 4 2 Adjoint Problem for Three-Dimensional Finite Elem e nt Saturated Groundwater Model . . . . . . . . . . . . . . . 7 3 4.3 Searching Objective Function ....................... . ..... ... .. 80 lll PAGE 4 4.4 Adjoint State Controlling Equations ........ . ..................... 80 4 5 Searching Method and Stop Criterion ....................... ....... 82 4.6 Field Application on SV5 . . . ................................. 82 4.8 Summary ............................. . ...... . .......... 85 5 THREE-DIMENSIONAL FINITE-ELEMENT VARIABLY-SATURATED GROUNDWATER MODEL .......................... ...... ...... 87 5 1 Introduction ................................................. 8 7 5 2 Three-Dimensional Finite-Element Variably-Saturated Groundwater Model ................. ...... ........... ..... 88 5.2.1 Governing Equation ..... ... . .......................... 88 5.2.2 Linearizing the Governing Equations ............. ....... . 89 5.2.3 Boundary Conditions for the Variably-Saturated Model ........ 91 5.3 Finite-Element Method for Modeling Three-Dimensional Variably-Saturated Flow ............ .......... ............... 94 5.4 Transient, Variably Saturated Water Table Recharge Example Test ...... 98 5 5 Application on SV5 Wetland .... .... ................ .... ...... 101 5.5.1 Finite-Element Mesh for Three-Dimensional Variably-Saturated Flow Model ......................... 10 1 5.5.2 Modeling Results and Summary ......................... 103 6 SUMMARY AND CONCLUSIONS ......... . ........ ...... ..... ..... 111 REFERENCES .............................................. ......... 113 BIOGRAPIDCAL SKETCH ... ........................................ 126 IV PAGE 5 LIST OF TABLES Figure 2-1 Parameters and their values .................................. .......... 17 3-1 Selected early references for the finite-difference method .................... 22 3-2 Selected early references for the finite-element application to model groundwater 25 3-3 Parameters and their values .................... ..... ..... ............. 42 3-4 Data for calculating the discharge to the wetland from the aquifer at 10 hours .... 63 3-5 Data for calculating the discharge to the wetland from the aquifer at 20 hours ... 64 3-6 Data for calculating the discharge to the wetland from the aquifer at 30 hours .... 64 37 Data for calculating the discharge to the wetland from the aquifer at 40 hours .... 65 3 8 Data for calculating the discharge to the wetland from the aquifer at 50 hours .... 65 3 9 Sensitivities calculated by saturated groundwater flow model .... ............. 66 4-1 Error and sensitivities based on the adjoint method model .................... 84 5-1 Data for calculating the discharge to the wetland from the aquifer at 10 hours ... 10 7 5-2 Data for calculating the discharge to the wetland from the aquifer at 20 hours ... 108 5-3 Data for calculating the discharge to the wetland from the aquifer at 30 hours . 108 5 4 Data for calculating the discharge to the wetland from the aquifer at 40 hours ... 109 5-5 Data for calculating the discharge to the wetland from the aquifer at 50 hou r s ... 109 V PAGE 6 5-6 Sensitivities calculated by variably saturated groundwater flow model ....... . 110 VI PAGE 7 LIST OF FIGURES Figure 2-1 Setup of the groundwater problem around a circular wetland ................... 9 2-2 Transient groundwater level changes at radial distances 1,2,3 and 4 meters away from the wetland shoreline ........ ............. ..... .......... 15 2-3 Groundwater distribution over the distance to the wetland shoreline at 6 7 8,9 10 h ................................ .... ........... ... 15 2-4 Three-dimensional view of the groundwater table around a circular wetland at t=3hr ........................................ . 16 2-5 Comparison of the analytical sol ution and numerical solution at 2 and 5 meters from the wetland shore line .......................... 18 3-1 Conceptual profile model of a wet l and/gro und water system .... .............. 26 3-2 Comparison of finite element mesh and finite difference mesh ............... 26 3-3 Setup of the saturated groundwater problem around a wetland ................. 27 3-4 Flowchart illustrating the iteration procedure of finding the coupling boundary cond ition ................. ...... ... .. ........ 29 3-5 Unsteady flow to a well in phreatic aquifer at time t .... ........... ........ .. 38 3-6 Finite-Element Method mesh of the wedge u sed to simulate of a well pumping in phreatic aquifer .......... ....... ............ ... 40 3 7 Comparison of the numerical and analytical head for the example problem of a transient pumping test in a phreatic aquifer ......................... 40 3 8 Comparison of the numerical and ana lytical simulated groundwater table .... .... 41 Vll PAGE 8 3-9 Three-dimensional finite-element mesh for the simulation of wetland and aquifer 43 3-10 Relationship between the area of wetland water surface ..................... 44 3-11 Relationship between the deepest water depth and the volume of wetland ....... 45 3-12 Illustration of the initial condition where it was assumed a wetland surface-water stage and the phreatic surface elevation were the same ....... 45 3 -13 Plan view of the numerical mesh at time zero ............................. 46 3-14 Three-Dimensional Finite-Element mesh and simulated head contour at hour 5 . 46 3-15 Simulated head contour at hour 5 ...................... . ............. 47 3-16 Simulated iso-surface plot of 10.6 meter head at hour 5 ............. ....... 47 3-17 Three-dimensional Finite-Element mesh and head contour on a cross section on a cross section through wetland .................... 48 3-18 Simulated head contour on the vertical cross section through wetland ......... 48 3-19 Modeled iso-surface head plot of 11.1 m,llm and 10 5 m ................ . 49 3-20 Mesh and modeled head contour with vertical cross section through the pumping well .......................... . ...... ...... . 49 3-21 Head contour (m) on the vertical cross section through the pumping well ....... 50 3-22 Well location map (Switt et al. 1998) ........... ....................... 51 3 23 Top surface of peat layer ............................................. 54 3 24 Bottom surface of peat layer .......................................... 55 3-25 Relationship between wetland volume and water stage .......... ... ....... 55 3-26 Finite-element mesh of wedge of wetland and aquifer for the saturated flow model ...... ...... ............................ 57 3 27 Vertical cross section view of three dimensional finite element m e sh and head from the saturated flow model ................... .......... . 61 Vlll PAGE 9 3-28 Piezometric head at positions beneath wetland with time ................... 61 3-29 Variation in the inundated area of wetland during simulation ................ 62 3-30 Discharge between wetland and aquifer with respect to time ................. 62 3-31 Comparison between calculated head and observed head .................... 63 3-32 Comparison results using Walser s parameter values . . ............. ..... 67 3-33 3 D numerical model simulated results over specific storage ... . ............ 67 4 1 Sketch showing setup of the saturated groundwater inverse problem around a wetland ...................... .............. 74 4-2 Relation between error and peat layer hydraulic conductivity . ................ 84 4-3 Flowchart illustrating the iteration procedure of inverse problem ... .......... 85 4-4 Vertical cros s -section of the adjoint state distribution at -r = 10 hour s ........... 86 5-1 Sketch showing setup of the variably saturated groundwater flow problem around a wetland .................... . .............. . 92 5 2 Schematic of relative evapotranspiration (ET/PET) as affected by soil water potential ,\jJ ................ .................. 93 5-3 Schematic diagram of the flow domain ............ . . ... ... . ....... . 99 5-4 Water moisture distribution (x : y=2 : 1) at 4 hours ....................... ... 100 5 5 W a ter moi s ture di s tribution at 8 hour s ( x : y =2: 1) .......................... 101 5-6 Finite-element mesh of the wedge part for variably saturated flow model (Unite is m) ........... ... ... ......... 102 5 7 Calculated piezometri c head variation beneath wetland u s in g the three dimensional variably saturated groundwater flow model ....... . 104 5 8 Model predicted dischar g e from th e aquifer to the ov e rlayin g w e tland ......... 104 5 9 Inundated area of wetland durin g simulation ............................. 105 l X PAGE 10 5-10 Simulated and observed heads in wetland and in aquifer ................... 105 5-11 Simulated head from the saturated flow model and the variably saturated flow model ..................................... ....... 106 5-12 Horizontal view of the three-dimensional finite-element mesh and piezometric head distribution from the variably saturated flow model at time of 9 hours .. 106 5-12 Enlarged top and middle part of the mesh and head contour at time of 9 hours .. 107 X PAGE 11 LIST OF SYMBOLS C the conductance of the wetland confining unit c o coefficient dlt) the unknown values of hydraulic head and at time t ET evapotranspiration [UT] F(w) the Fourier coefficient G ground-water level [L] G i groundwater inflow [L3ff] G o groundwater outflow [L3ff] GWNS the net flow or seepage of groundwater to the wetland surface [L3ff] h = piezometric head [L] h(e) the approximate solution for hydraulic head within element e [L] H o initial saturated aquifer thickness [L] H 1 Sobolev space H0(t) the Hankel function order 0 hoo the regional average water table elevation, [L] reference integer spatial interval J o the Bes sel function of the first kind order 0 k a constant parameter K the sa turated hydraulic conductivity tensor [ UT] K/Kv anisotropic rati o K s hydraulic conductivity [UT] Kxx horizontal hydraulic conductivity in x direction [UT] XI PAGE 12 -horizontal hydraulic conductivity in y direction [UT] -vertical hydraulic conductivity [UT] K = the saturated h ydraulic conductivity tensor [UT] IJ = the relative permeability while the degree of water saturation is S w L = average peat thickness [L ] n -current integer time level N.< e > -the interpolation functions for each node within element e I n -a unit outward normal vector to the boundary N o -the number of observation wells Nn -the total number of nodes used in the numerical solution p = the number of nodes within element e pn = net precipitation ( total precipitation minus intercepted precipitation[ L3fT] p c = the capillary pressure (Pc = lfT) [M/LT2 ] PET -the daily potential evapotranspiration [UT] Q -pumping discharge [L3fT] r -the radial coordinate [L] radial distance from the well [L] r o -the radius of the circular lake [L] s -drawdown [L] s -the specific storage S i -surface inflows, including flooding stream [ L3rr] s o -surface outflows [L3fT] s s -the specific storativity [L -1 ] SAW -wetland wetted surface area [L2 ] SANW -wetland dry surface area (L2 ] SP -exchanging flux [L3fT] S W -the degree of sa turation equal to 8/85 Xll PAGE 13 S wr -the residual saturation S e -effective saturation t -time [T] tr -the time duration of the simulation [T] t j -the observation times [T] T -the transrnissivity in the aquifer [L2/T] T i -tidal inflow ( +) or outflow (-) [L3/T] u -space function v w -the wetland volume [L3 ] Vp -the volume above the bottom of the peat layer [L3 ] w -wetland water level [L] w . PAGE 14 r1 -flux boundary r2 -head boundary y -specific weight of a unit volume of water= pg [MIL 2T2 ] 'V -the pressure [MIL T 2 ] Qi -the exclusive subdomain of node i oK = variation of K [IJT] oh -variation of head h [L] 'P(x,y, tr-the adjoint state of h i-) 't t r t [T] a andm = parameters obtained from a fit of above equations to experimental results 0 -the water content, [L3/I}] 0 s -the saturated water content, [L3/L3 ] e r -residual moisture content [L3/L3 ] 'V -pressure head; [MILT 2 ] w -the value of \JI when ET= 0.5PET, [MILT 2 ] 'Vwmm -the \JI corresponds to PE/PET=0.05, [MILT 2 ] l.Jrwmax -the value of \JI when PE=0. 95 PET [MILT 2 ] XIV PAGE 15 Abstract of Dissertation Presented to the Graduate School of the University of Florida in Partial Fulfillment of the Requirements for the Degree of Doctor of Philosophy TWO-DIMENSIONAL ANALYTICAL AND THREE-DIMENSIONAL FINITE ELEMENT METHOD MODELING OF THE INTERACTIONS BETWEEN WETLANDS AND GROUNDWATER By Jinhua Zhong December 2004 Chairman: Louis H. Motz Cochair : Kirk Hatfield Major Department: Civil and Coastal Engineering Simulating saturated/unsaturated flow between wetlands and the contiguous subsurface is complex because it involves the modeling of wetland storage changes that in turn cause changes in pertinent hydrologic boundary conditions. Existing groundwater models are incapable of simulating the flow between a wetland and an aquifer, if small changes in wetland storage cause significant vertical and lateral movements of the wetland boundary. In addition to the simulation complexities introduced by a change in wetland geometry and storage, periodic and nonperiodic fluctuations of the wetland stage induce corresponding elevation changes in the phreatic surface of the contiguous aquifer. As these induced waves propa gate throu g h an aquifer, friction causes a lo ss of e ner gy, which is manife ste d as spatially dampened potentiometric head perturbations along the inland direction. An analytical model was developed that describes subsurface flows around a wetland induced by any nonperiodic fluctuations in surface-water stage, such as a flood wave Not considered, xv PAGE 16 however are any changes in wetland storage or wetland geometry as a function of surface water stage. The analytical model was validated using a finite difference numerical groundwater model that is based on the governing equation expressed in radial coordinates. Comparisons of analytical and numerical results show excellent agreement. Two numerical models were specifically designed to simulate wetland-aquifer interactions. The first model is a saturated groundwater flow model that incorporates adaptive simulation techno l ogies to permit real-time simulation of moving wetland boundaries and their associated local influence on the phreatic surface. This model is numerically efficient and uses deformable hexahedra1 finite elements to trace phreatic surface changes in real time. The second model also uses deformable hexahedral finite elements ; however this model simulates variably saturated groundwater flow For both models the numerical element s deform as required to characterize the horizontal and vertical extent o f the moving wetland boundary (which serves to define a location within the porous system where the pressure is known). Simulation results from both models were applied to a field site where water was rapidly withdrawn from an isolated wetland to observe system response and evaluate the wetland aquifer hydraulic connection. XVI PAGE 17 CHAPTER 1 INTRODUCTION 1.1 Wetland Wetlands are transitional lands between terrest1ial and aquatic systems where the water table is usually at or near the surface, or the land is covered by shallow water (Cowardin et al. 1979) Wetlands are also called marsh fen peatland or water and can be natural or artificial and permanent or temporary The water can be static or flowing and fresh or brackish or it can include areas of marine water the depth of which at low tide does not exceed 6 meters (Ramsar Convention 1971, Middleton 1998). Wetlands perform various environmental functions that include increasing flood mitigation enhancin g storm abatement (i.e ., in the case of coastal wetlands) adding groundwater recharge providing water quality treatment serving a s wildlife habitats, and performing aesthetic functions (Mitsch and Gosselink 1993) These functional attributes can be classified into three basic types: hydrologic water quality and life support. Measures of the life-support function must integrate several wetland attributes including water quality hydroperiod and hydrodynamics habitat structure food chain support bio g eographic setting nutrient status, corridors for migration and primary production (Preston and Bedford 1988). In recent years increasing interest has focused on the study and protection of wetlands. Developm e nt drainage population pre s sure and pollution have destroyed extensive areas of these sensitive ecosystems. Another increasing! y significance to wetlands is municipal well fields supplying water to large urban areas. These well fields affect 1 PAGE 18 2 wetlands by reducing surficial aquifer levels over large areas that encompass isolated and interconnected wetlands (Sonenshein and Hofstetter 1990). As a result of depressing the water table, historically wet ecosystems have been drained (Switt et al. 1998). 1.2 Wetland Water Budget Hydrologic conditions are extremely important for the maintenance of a wetland's structure and function. A hydrologic budget is often used to quantify the volume and rate of water flowing into and out of a wetland. The primary hydrologic inflows include precipitation, streams, groundwater, and, in coastal wetlands, tides (Mitsch and Gosselink 1986, Rosenberry 2000). The wetland water budget equation as defined by Mitsch and Gosselink (1993) is: where VI ~t = P + S + G ET-S G T n I I O 0 V = change in the volume of water storage in wetland [L3 ] ~t = change in time [T] (1-2) P0 = net precipitation ( total precipitation minus intercepted precipitation en) [ L3/T] Si = surface inflows including flooding stream [ L3/T] Gi = groundwater inflow [L3/T] ET= evapotranspiration [L3/T] SO= surface outflows [L3/T] G0 = groundwater outflows [L3/T] T = tidal inflow ( +) or outflow (-) [L3/T] PAGE 19 3 Terms in the wetland water budget vary in importance depending on the wetland type and not all terms are required to characterize a specific type of wetland. In most wetland systems, precipitation and evapotranspiration are major components of the water budget (Winter and Woo 1990). Groundwater can also have a significant influence on some wetland systems, but not always. Three of the most common relationships that characterize a wetland are groundwater discharge groundwater inflow and outflow, and groundwater recharge. Groundwater inflow occurs when the contiguous wetland surface water (or groundwater) stage is lower than the phreatic surface of the surrounding aquifer. When the surface water stage or groundwater levels inside a wetland are higher than the phreatic surface in the surrounding contiguous aquifer, groundwater and surface water will flow out of wetland When a wetland is situated well above the groundwater, the wetland is identified as being perched. This type of hydrologic system, referred to as a surface water depression wetland by Novitzki (1979), loses water only through infiltration into the ground and through evapotranspiration. The rate of water exchange between a wetland and local aquifer is a function of the hydraulic gradient between the two hydrologic subsystems. The discharge, in accordance with Darcy's law (Truesdell 1995) is proportional to the hydraulic gradient, the hydraulic conductivity, and the cross-sectional area through which the flow occurs Periodic or nonperiodic fluctuations of the wetland stage in a flow through wetland will induce fluctuations in the elevation of the phreatic surface in the contiguous aquifer. As these induced waves propagate through an aq uife r, friction cau ses a loss of e n ergy, which is manifested as spatially dampened head perturbations This damping effect progressively diminishes the amplitude of piezometric head along the inland direction (Carr and Van Der PAGE 20 4 Kamp 1969). Many studies have examined the problem through one-dimensional analytical solutions (Cooper 1963, Gregg 1966, Nielsen 1990 Serfes 1992, Li Barry et al. 2000, 2001, Boufadel and Peridier 2002 Seerano 2003, Swamee and Singh 2003). A two-dimensional analytical solution was obtained by considering a shift in phase and a change in amplitude along a straight coastline (Sun 1997) Three-dimensional flow regimes near a shallow circular lake was studied using numerical and analytical mathematical models of aquifer flow (Townley and Trefry 2000) To further the understanding of the interactions between lakes or wetlands and groundwater, without the restrictive assumptions typically associated with analytical solutions numerical models are usually used. For example, it is possible to obtain numerical solutions for the case of anisotropic and nonhomogeneous conditions. The accuracy of solutions obtained by numerical methods can be excellent; however it depends on several factors including the type of numerical method used the complexity of boundary and initial conditions and the computational precision of the computer used to implement the method. Regardless of these advantages few investigators have chosen to focus on the numerical simulation of interactions between wetlands and groundwater ( McDonald and Harbaugh 1988 Schaffranek 1996 Lee 2000). 1.3 Goals and Objectives Analytical models to simulate the complex groundwater surface water interactions that occur with isol a ted and flow throu g h wetlands would be valuable tools for developin g hydrologic budgets for wetlands. Given that most surface water fluctuations are nonperiodic there is a need to develop analytical groundwater models for simulating groundwater table fluctuations induced by both periodic and nonperiodic surface water flood waves. However PAGE 21 5 with the exception of the one-dimensional study by Cooper and Rorabau gh (1963), most investigations addressed surface-water fluctuations that were periodic. Thus, the first goal of this dissertation was to develop an analytical method for simulating groundwater table fluctuations induced by both periodic and non periodic surface-water flood waves. This goal was achieved by meeting the single objective of developing an analytical model using Fourier transforms A limited number of three dimensional numerical groundwater models have been developed to simulate groundwater wetland interactions; however, these models are generally inadequate tools for simulating the complex interactions between a wetland and a contiguous unconfined aquifer. For example Walser, Annable and Wise (1998) used FEMW ATER (Yeh, 1992), a variably saturated flow model, to simulate the exchange of water between a wetland an d an aquifer. McDonald and Harbaugh (1988) used MODFLOW and developed a three-dimensional finite difference groundwater model and studied the dynamic coupling between wetlands and the subsurface. Lee (2000) used the HYDRUS-2D flow model and simulated vertical two-dimensional groundwater flow through Lake Barco in Florida Neither FEMW ATER (1996) nor MODFLOW can simulate or incorporate a moving real-time wetland boundary caused by a change in surface water stage Furthermore as a saturated three-dimensional finite-difference model MOD FLOW can not simulate fluctuations of the phreatic surface as a moving transient boundary as required to define fully the three dimensional flow. Finally the water stage of a wetland itself is a function of the transient exchange of water between a surface and subsurface sys tems; which is not fully modeled in any existing groundwater flow codes. Thus, to simulate groundwater surface water interaction without the above identified s hortcomin gs of existing models a new model is PAGE 22 6 needed. The second goal of the dissertation was to develop two numerical models capable of simulating 1) a transient and moving wetland boundary, 2) three-dimensional flow 3) the movement of the phreatic surface, and 4) the coupling that occurs between the transient surface-water stage and the phreatic surface. To achjeve trus goal, two objectives were pursued. The first objective was to develop two numerical groundwater models using firute elements. The second objective was to apply both numerical models to simulate groundwater and surface water fluctuations over an isolated wetland in southern Florida. Of the two models developed the first is a three dimensional saturated groundwater flow model that employs deformable hexahedral elements. The finite-element mesh deforms in space and time to trace the location of the phreatic surface of an unconfined aquifer and the stage of the contiguous surface water body. Based on the above model and the adjoint state method an inverse model was developed to estimate the hydraulic conductivity of peat layers wruch typically exist as a wetland bottom that serves to isolate surface wetlands from subsurface groundwater systems. The second numerical model a three dimensional variable saturated flow model was also developed using deformable mesh. This model was used to investigate prediction differences between the saturated flow code and the variable saturated model. From the groundwater flow pattern given by both models, an improved understanding was obtained for the interactions occurring between surface waters and groundwaters. In addition an improved understandin g was obtained of how short term pumping from a phreatic aquifer affects wetlands and lake water stages PAGE 23 CHAPTER 2 ANALYTICAL GROUNDWATER MODEL RESPONSE TO THE PERIODIC OR NON-PERIODIC WATER STAGE FLUCTUATIONS IN A CIRCULAR WETLAND 2.1 Introduction Periodic or nonperiodic fluctuations of the surface water stage in wetlands or lakes will induce fluctuations in the elevation of the phreatic surface of contiguous unconfined aquifers As these induced waves propagate through the aquifer, friction causes a loss of energy, which is manifested as spatially dampened head perturbations. This damping effect progressively diminishes the amplitude of fluctuations in the piezometric head along the inland direction (Carr and Van Kamp 1969). Many studies have examined the problem through one-dimensional analytical solutions (e.g. Cooper and Rorabaugh 1963, Gregg 1966 Nielsen 1990 and Serfes 1992 Li, Barry et al 2000 2001 Boufadel and Peridier 2002, Seerano 2003 Swamee and Singh 2003). Cooper and Rorabaugh (1963) obtained the analytical solution for groundwater movement and bank storage due to flood stages in surface streams Gregg (1966) found a one-dimensional relation between tidal fluctuations and groundwater table perturbations. Nielsen (1990), Serfes (1992) and Carr (1969) also published papers investigating the similar problems. A two-dimensional analytical solution was obtained by considerin g a phase shift in water table perturbations and the associated change of amplitude along a straight coastline (Sun 1997) An analytical solution representing vertical two-dimensional groundwater-surface water interaction was obtained by Anderson (2003). Three-dimensional flow regimes near a shallow circular lake was 7 PAGE 24 8 studied using numerical and analytical mathematical models of aquifer flow (Townley and Trefry 2000, Smith and Townley 2002). With the exception of Cooper and Rorabaugh (1963), all of the above studies addressed periodic surface-water fluctuations, in spite of the fact that most stage variations are not periodic i.e., consider, for example, the pervasive nonperiodic flooding of rivers, lakes, and wetlands. Cooper's (1963) nonperiodic analytical model for a river/groundwater system was a one-dimensional solution along an axis perpendicular to the river. However, this analytical model is not suitable for simulating the nonperiodic interactions between a wetland and an aquifer. In this chapter, two analytical groundwater models are developed to describe transient groundwater levels affected by periodic and nonperiodic fluctuations of the surface water stage in a wetland or a lake. 2.2 Problem Formulation: Groundwater Flow Around a Circular Wetland Considering a circular wetland that is hydraulically contiguous to a surfacial aquifer and ignoring recharge from rainfall, the governing equation for two dimensional groundwater flow in a homogeneous isotropic aquifer can be expressed in radial coordinates as: a2h 1 ah s ah -+--= ar 2 r ar Tat (2-1) where his the piezometric head r is the radial coordinate tis time, Sis the specific storage and Tis the transmissivity in the aquifer. When a thick unconfined aquifer is considered, the gradient and small fluctuations relative to saturated thickness are assumed to be small Tis approximately equal to Kh00 (Figure 2-1) in which K is the hydraulic conductivity and h00 is PAGE 25 9 the average unconfined aquifer thickness equal to the regional average water-table elevation above an underlying impermeable confining unit. The boundary condition along the shore line of the wetland is Eq. 2-2, where r O is the radius of the circular wetland. It is also assumed in Eq 2-3 that as r----.00, the groundwater table elevation approaches the regional average water table elevation h00 B.C: h(r,t)j r=r o = hoo + f(ro ,t) h(r,tJI r=.., = hoo w e tland or lake @ r= r 0 @ r= oo .......... -2_ ........ Figure 2-1. Setup of the groundwater problem around a circular wetland 2.3 Analytical Solution for the Periodic Boundary Condition (2-2) (2-3) Since the boundary condition is a periodic function of time, substituting for f(r 01t) Equation 2 4 is obtained: B.C: h(r,t)I r=r o = hoo + h o e iw t h ( r ,tJI r=.., = hoo @ r = oo (2 4) (2 5) where h0 is the half height of the periodic constituent and w is the angular frequency. No initial condition is necessary since here the steady-state oscillatory solution is sought as t--00 Because of the oscillatory boundary condition a solution is sought in the form PAGE 26 10 Substitution of Eq. 2-6 into Eq. 2-1, Eq. 2-2, and Eq. 2-3 yields a2g 1 ag Sw -+----zg=O ar2 r ar T where the boundary conditions are: g(r)=l g(r)=O @ r=r 0 @ r=oo Equation 2-7 can now be rewritten to obtain: v = ~ 8"'-(i-l)r 2T (2-6) (27) (28) (2-9) (2-10) (2-11) The solution to Equation 2-10 subject to the boundary conditions in Equation 2-8 and 2-9 depends on the value of w. Considering the boundary condition along the shoreline, groundwater around a circular water body will have the following forms. When w i s greater than 0: H(1)(v) g(r)= C (J +iY ) = C H (1)(v) = -0 -o o o o o ( I ) Ho (vo) (2-12) PAGE 27 11 and when w is less than 0: H(2)(u) g(r) = C (J -iY ) = C H(2)(u) = -0 -O O O O O (2) H0 (u) (2-13) where J0 is the Bessel function of the first kind order 0; Y0 is the Bessel function of the second kind order 0; H0(t) is the Hankel function order 0; and C0 is a coefficient (Lebedev 1965). Substitution ofEq. 2-12 and Eq. 2-13 into Eq. 2-6 gives the following solution for the groundwater flow around a circular lake affected by periodic function: H~1)(v) h (r t)= h +h ---e wiz W 00 0 ( ) H0 (v) (2-14) for w > 0 and: H;1)(u) h (r t) = h +h ---e wu w 00 0 (2) H0 (u) (2 15) for w < 0. The asymptotic approximations of the functions H0<1l (z) and ~ <2l (z) are: (2 16) PAGE 28 12 Hi2)(z) = e 4 L (-Ii -(2izfk+O(lzl n l) ( 2 l 1/2 ( l)i(z-~)[ n 1 3 2 (2k 1)2 l 1tZ k = O 22kk! u =~ S l w l (1-1)-r 2T v = ~ Sw (i l)r 2T u =~ S l w l (-1-1)-r o 2T o v = ~ Sw (i l)r o 2T o 2.4 Analytical Solution for the Non-Periodic Boundary Condition (2-17) A flood wave will typically induce a non periodi c variation in the wetland stage. Here, f(r0,t) whic h is a n o np erio dic function of time, is assumed to describe transient variations in the wetland stage. Usi n g the Fou rier transformation f(r0 t) can be represented as: 1 J f{_r 0,t) = 2 TT F( w )e zwt dw (2-18) where F(w) is the Fourier coefficient, and F(w ) dw is the amplitude of oscillation for the compon e nt with a freq u ency in the range of ( w ,w+dw ) For this F(w), the corresponding fluctuations in the groundwater table are described by: H ~1)(v) h (r t) = h + F(w) ---e wu w oo (I) Ho (vo) ( w > 0 ) ( 2 19) PAGE 29 13 or H;1)(u) h (r t)=h + ~(w) ---e wti (i) 00 (2) H0 (u) ( w < 0) (2-20) as given by Eq. 2-14 and Eq. 2-15. The inverse Fourier transformation of Equation 2-19 and Eq. 2-20 leads to: 1 h(rt)=h +-, 00 21t = h + _l J Jf(r 1:)e -w,id,: oo 27t o' O -oo =h + -1 Jf(r 1:) oo 27t o' J H;1)(v) J o H;2)(u) ---e w(t-,)i dw + ---e w(M)i dw d,: H;1)(v) H;1)(u) 0 (2-21) From the above analytical model, it is clear that the unsteady groundwater table around a circular wetland depends not only on the instantaneous values of the wetland water stage but also on the accumulative effect represented by the integration However physically it is not possible for the wetland water stage at a later time to affect the groundwater table at an earlier time. Hence, the phreatic surface depends on its current stage and historical stage. The PAGE 30 14 final analytical solution to the nonperiodic boundary condition 2-2 is then I h(r,t)=h +-1 fj(r ;t) 00 21t 0 (2-22) 2.5 Application If it is assumed that wetland stage fluctuations can be described by a nonperiodic Gaussian function of time, then: (11 / 1 202 f(r0 t)=--e {i:rr.a (2-23) The analytical solution obtained when Eq. 2-23 is combined with Eq. 2-21 is: JI (tt ,)2 1 1 h(r,t) =hoo + ---e 202 2n {i:rr.a (2-24) The transient c h anges of groundwater table described by Eq. 2-24 at four different distances from the wetland are shown in Figure 2 2. The groundwater distributions along the radius at four different times are shown in Figure 2 3 From both figures it is evident that groundwater fluctuations are greatest near the wetland. From Figure 2 3 it can be seen that the variation of the groundwater table over time is no longer symmetrical although the PAGE 31 15 variation of the wetland water stage is a Gaussian distribution which is symmetrical over time 101 ]: 100 8 Q) I----:0 100.6 CCI .... Q) 100.4 -CCI 3: -a 100 2 C: ::::, 0 100 .... C} 99 8 0 -5 10 15 Time (hr) ---------~20 25 lake d=1m d=2m d = 3m d=4m Figure 2-2. Transient groundwater level changes a t radial distances 1 2,3 and 4 meters away from the wetland shoreline. 100.3 :[oo. 25 ~100. 2 <1l -I I I I ~00. 15 iii 3:100.1 "O I ., / 05 0 I // / (5 100 ~ 99. 95 50 t =6hr I\ \ \ \ t =7hr \ t =8hr -'-"""-' ~---..:: t =9hr t = 1 Ohr 55 60 65 7 0 Distance to lake shoreline (m) Figure 2-3. Groundwater distribution over the distance to the wetland/lake s horeline at 6 7 8 ,9, and 10 h. PAGE 32 16 Figure 2-4 shows the groundwater table elevations around a circular wetland after 3 hours The radius of the wetland is 50 meters. All other parameters are listed in Table 2-1. Theoretically when w is very large it means that the frequency of the corresponding boundary condition is very large and the period is very short. When the absolute value of w is very large, its effect on the groundwater around the wetland is very small. This is because the period of the wave component is too small to have any impact on the groundwater around the wetland. That is to say when the integration of w is taken, a certain negative number could be used to replace the negative infinite number and a certain positive number could also be used to replace the positive infinite number. Accurate model results are obtained when w takes the range of [-5 5]. Thus, it may be assumed that the boundary condition has no effect on the ground water table if the frequency of the boundary water level is very large. Figure 2-4 Three-dimensional view of the groundwater table around a circular wetland at t = 3 h. PAGE 33 17 2 6 Comparison of Analytical Solution and Numerical Solutions A comparison was made between the analytical and the numerical solutions to the radial flow problem specified through Eq. 2-1-Eq. 2-3. An implicit finite-difference scheme was used to represent the radial form of the groundwater flow equation Eq. 2-1 as: hn+l_2hn+ l hn+ l hn+l_hn+ l hn+l_hn i-1 I + I + I 1 I+ 1 1 1 S I I -------+-----(/:l.r)2 T l:l.t (2-25) Here i is the reference integer spatial interval; n is the current integer time level ; D.t is the time step; and D.r is the radial step. A circular wetland which has a radius of 50 meters was assumed. AJI other system parameters are given in the Table 2-1. By considering the different boundar y conditions (Eq. 2-3 and Eq. 2-4), the corresponding numerical results are found [Figure 2-5]. Analytical and numerical results illustrate time variations in groundwater levels Figure 2-5 shows that the analytical solution and numerical solution compare very well. Table 2-1. Parameters and their values I Parameter I Value II Parameter I Value I H00 (m) 100.0 T (m 2 /h) 0.635 D.r (m) 0.25 s 0.2 D.t (hr) 0.02 t o (h) 5.0 k (m/hr) 0.00635 a 0.4 PAGE 34 18 0.4 A --0 3 t----<--01 0.2 t----+---t--.-1---t---+--,l---t---+--,f----1 l -0.1 0 5 -------.... 10 15 Time (hr) 20 -25 analy(d=2m) analy( d=5m) num (d=2m) num (d=5m) Figure 2-5. Comparison of analytical solution and numerical solution at 2 and 5 meters from the wetland shoreline. 2.7 Summary The analytical model presented in this chapter can be used to simulate aquifer surface fluctuations caused by the transient wetland stages regardless of the temporal function assumed to describe variations of water stage. In other words the model can deal with any general transient of the boundary function. The analytica l model was validated using a finite difference numerical groundwater model for a homo genous aquifer. This first numerical method is not recommende d for simulations involving hetero ge nou s aquifers. Rather a new method is presented for simulating wetland/aquifer interactions in such aquifers, which requires a simple transformation of the equations. The excellent agreement between two models shows that the analytica l solution presented is valid. This model is also useful in a situation where groundwater and the surface water components are being considered in a system-wide water balance. PAGE 35 CHAPTER 3 TRANSIENT THREE-DIMENSIONAL DEFORMABLE FINITE-ELEMENT SATURATED GROUNDWATER MODEL 3.1 Introduction Recent research has focused on the details of various wetland functions such as water-quality treatment, flood-wave suppression, and groundwater/surface water interactions (Wise et al. 2000, Krabbenhoft and Webster 1995 Anderson and Cheng 1993 Ruddy and Williams 1991 and Mills and Zwarich 1986). For example, Ruddy and Williams (1991) investigated the hydrologic relation between the South Fork Williams Fork stream in the Colorado River basin and four adjacent subalpine wetlands and they also evaluated the potential effects on these wetlands when stream flow is diverted They concluded that stream flow was the primary source of water for two of the four wetlands whereas precipitation snowmelt valley side-slope flow (including overland and small-channel flow) and groundwater inflow served to define the transient hydrological behavior of all wetlands. The above studies showed that wetland functions (i.e. water quality treatment floodwave suppression, etc.) are correlated with rechar g e and discharge areas ground water flow paths and rates and soil profile development. The location of recharge and discharge areas can change in response to seasonal or weather-related factors and can affect the groundwater flow direction As a result groundwater and surface water quality and the presence of various wetland plant communities may all vary in response to these changing conditions. 19 PAGE 36 20 Analytical models such as the one presented in Chapter 2 can be used to simulate deep wetlands or lakes. However, the primary limitation of analytical methods is that solutions can only be obtained by imposing severely restrictive assumptions on aquifer properties and boundary and initial conditions. For example, an assumption commonly made to obtain analytical solutions is that the aquifer is isotropic and homogeneous for hydraulic conductivity. These assumptions, however, are not valid and often not adequate to describe groundwater flow or solute transport in most field situations. Numerical methods have the primary advantage over analytical approaches in that they do not require such restrictive assumptions For example it is possible to obtain numerical solutions for the case of anisotropic and nonhomogeneous conditions The accuracy of solutions obtained by numerical methods can be excellent. However, this accuracy depends on several factors, including the type of numerical method used the complexity of boundary and initial conditions and the computational precision of the computer used to implement the method Regardless of these advantage, only a few investigators have chosen to focus on the numerical simulation of interactions between wetlands and groundwater (e.g., Schaffranek 1996 and McDonald and Harbaugh 1988 Lee 2000, Rosenberry 2000). In this chapter, a general discussion of the advantages/disadvantages of two numerical methods commonly used to model hydrologic systems is briefly discussed. Next the mathematical development of a three dimensional deformable finite -e lement saturated groundwater model is presented This model and the associated code are then validated by comparing model predictions with analytical modeling results. A more complicated application of the three-dimensional model is presented to illustrate the type of information provided through th e numerical approach taken. Finally, th e PAGE 37 21 model is used to simulate an actual wetland using field data collected when the hydrologic system was responding to induced stress. 3.2 Ground-Water Flow Models Based on Finite-Difference and Finite-Element Approximation Several types of numerical methods have been used to solve groundwater flow problems, the two principal ones being the finite-difference method and the finite-element method. The finite-difference method was initially applied to the flow of fluids in petroleum reservoirs, and it wasn t until the rnid-1960's that problems of groundwater flow and solute transport were examined. The method has a number of advantages that contribute to its continued widespread use and popularity: (1) for simple problems such as one-dimensional, steady-state groundwater flow in an isotropic and homogeneous aquifer mathematical formulation and computer implementation are easily understood; and (2) the accuracy of solutions to steady-state and transient groundwater flow problems is generally quite good. The finite-difference method also has disadvantages: (1) the method works best for rectangular or prismatic aquifers of uniform composition (for example it is difficult to incorporate irregular or curved aquifer boundaries, anisotropic and heterogeneous aquifer properties or sloping soil and rock layers into the numerical model without introducing numerous mathematical and computer programming complexities); and (2) the accuracy of solutions to solute transport problems is less than can be obtained through application of the finite-element method Table 3.1 lists selected early references where finite-difference methods were applied to model groundwater flow. More recent applications of hydrologic modeling using finite difference for simulating of surface water /groundwater interactions include studies by Mc Donald and Harbaugh (1998), Winter (1983 1986 and 1989), Sacks PAGE 38 22 (1992) and Krabbenhoft and An d erson (1990). Table 3-1. Selected early references for the finite-difference method. Topic Reference Early Developments in Petroleum Bruce et al. (1953), Peace and Reservoir Modeling Rachford (1962) Saturated Groundwater Flow Remson et al. (1965), Freeze and Whitherspoon (1966), Pinder and B redehoeft (1968). Unsaturated Groundwater Flow Philip (1957) Ashroft et al. (1962), Freeze (1971), Brutsaert (1973), Healy (1990), Kirkland et al. (1992) Solute Transport Stone and Brian (1963), Oster et al. (1970), Tanji et al. (1967), Wierenga (1977) Application to Field Problems Orl ob and Woods (1967), Gambolati et al. (1973), Fleck and McDonald (1978), Jain 1992 Comprehensive Reference Trescott and Larson (1977), Ames (1977), Mitchell and Griffiths (1980), Lapidus and Pinder (1982) Computer Program Trescott et al. (1976), Konikow and Bredehoft (1978), Healy (1990) Winter (1983 1986, 1989), Sacks (1992) and Krabbenhoft and Anderson (1990) conducted numerical investigations on the interactions between groundwaters and lakes Both Winter and Krabbenhoft and Anderson used fixed finite-difference meshes. For lakes, a finite-difference mesh may be suitable; however, because the bed of most wetlands exhibits significant slope a poorly designed fixed finite difference mesh may not adequately characterize water fluxes across a bed situated between an overlying body of water and an underlying aquifer. The issue of properly modeling fluxes accross a sloping bed may not be PAGE 39 23 as important as the fact that accurate simulations of wetland-flow systems reqmre measurements of hydrologic properties at finer scales and frequencies than are commonly made (Hunt 1996, Hunt and others 1996 1997). Assuming the hydrological model is correct, the accuracy of any computer simulated hydrologic system, such as a wetland, depends on the quality of the data used to calibrate model. Through a well instrumented study-site in the Kankakee River Valley in northwest Indiana, the USGS and the USEPA conducted detailed measurements of many various variables required to characterize a wetland hydrologic system. McDonald and Harbaugh (1998) developed a groundwater model using the available data and conducted the simulations of observed hydrologic process using the computer code commonly known as MODFLOW. The focus of the research effort was to study the effects of various wetland restoration activities on local groundwater levels. Computer simulations were used to predict the effects of wetland restoration on the existing hydrology of the wetlands and adjacent farmlands and to suggest residence times and path ways for surface water and groundwaters. The finite-element method was first used to solve groundwater flow and solute transport problems in the early 1970s. The method has several advantages: (1) irregular or curved aquifer boundaries, anisotropic and heterogeneous aquifer properties and sloping soil and rock layers can be incorporated easily into the numerical model; and (2) the accuracy of solution to groundwater flow and solute transport problems is quite good. The main disadvantages of the finite element method include : (1) for simple problems the finite element method requires a greater amount of mathematical and computer programming sophistication than does the finite-difference method ; and (2) there are fewer well documented computer programs. PAGE 40 24 The finite-element method provides approximations of much higher order than does finite-difference methodology. The whole modeling domain can be subdivided into irregular elements as dictated by the physical geometry of the problem. The size of each element can vary. Small elements can be used in areas where rapid change occurs with model parameters and predicted state variables, whereas the large elements can be used where such variations are less severe. Inhomogeneities and anisotropy are easily accommodated with the finite element method. The advantage of irregular elements and higher-order approximations is certainly important when the space" domain is complex. However this advantage is less obvious in problems involving time (evolution) (Gupta 1975). Selected references of early finite-element applications to model groundwater flow are listed in Table 3-2. 3.3 Three-Dimensional Deformable Finite-Element Saturated Groundwater Model To simulate three-dimensional groundwater flow in phreatic aquifer systems and the exchange of water between a wetland and the contiguous aquifer, a numerical model is needed that is capable of modeling the free groundwater surface This requires first that the model assume a tentative location of the free surface and then solve for the distribution of the dependent variable tt,, over a fixed domain. All this requires that the flux boundary condition across the surface be satisfied and that the pressure be atmospheric. To model the exchange of water between a wetland and an aquifer, it is necessary to recognize that significant vertical and lateral movements of the wetland boundary typically occur with changes in wetland storage (see Figure 3-1) This transient boundary suggests that a deformable finite-element mesh may be needed in the model to characterize adequately the moving groundwater/surface water interface PAGE 41 25 Table 3-2. Selected references of early finite element applications to model groundwater. Topic Reference Early Developments in Petroleum Price et al. (1968) Reservoir Modeling Saturated Groundwater Flow Zienkiewicz et al. ( 1966), Javandel and Witherspoon (1968), Zienkiewicz and Parekh (1970), Pinder and Frind (1972). Unsaturated Groundwater Flow Neuman (1973), Gureghian et al. (1979), Pickens and Gillham (1980) Yeh and Ward (1981),Huyakom (1984), Yeh (1992) Solute Transport Price et al. (1968), Guymon et al. (1970) Neuman (1973), Van Genuchten et al. (1977), Kirkner et al (1984). Application to Field Problems Pinder (1973) Gupta and Tinji (1976) Senger and Fogg (1987). Comprehensive References Ziekiewicz (1971), Pinder and Gray (1977) Lapidus and Pinder (1982), Huyakom and Pinder (1983) Computer Programs Neuman and Witherspoon (1970) Reeves and Duguid (1975) Segal et al. (1975) Pickens et al. (1979) Huyakom (1984) Yeh (1992). Numerical groundwater modeling that is focused on simulating the interaction between surface water and groundwater systems where the groundwater/surface water interface moves must use well-defined linking conditions to describe the exchange of water across the interface These linking conditions are unknown in the beginning of the modeling time step and found through a flux type interaction that then forecasts the linking boundary conditions at the end of each time step. Thus the model has transient boundary condit i ons. PAGE 42 Phreatic Surface at Time t1 26 Phreatic Surface at Time ti Groundwater Aquifer Impermeable Aquacludc ////////////////////////////////// Figure 3 -1. Conceptual profile model of a wetland/groundwater system Finally typical finite-difference modeling is not able to trace the curvature of the bottom of a wetland (see Figure 3-2). Finite elements, however are easily configured to fit non-linear spatial geometries. In conclusion to cope with the above described three technical problems of modeling i e., the free surface location incorporating dynamic changes in the vertical and horizontal extent of the ground/surface water interface and tracing the curvature of a wetland bottom a special three-dimensional finite-element model is needed. :,,.b, /1 -___.,cC';l. -v Finit e -Elem ent M e thod Finite-Difference Method Figure 3 -2. Comparison of finite-element mesh and finite difference mesh PAGE 43 27 3 3 1 Governing Equation and Boundary Conditions Eq 3-1 below can be used to describe incompressible, saturated groundwater flow beneath a wetland and in a phreatic aquifer: V(KVh) Q=S ah i,j= l,2,3 s at (3-1) where K is the saturated hydraulic conductivity tensor ( LT -1); his equal to z+$ (L); $ is pressure head (L); z is elevation head ( L); and S s is the specific storativity (L-1). From Figure 3-3, the boundary of volume Q is characterized through a combination of specified flux boundaries r1 and specified head boundaries r2 Pertinent specific no-flux boundaries include the groundwater phrea t ic surface (EC and DF) and the bottom impermeable boundary (GH). Relevant specified head boundaries r2 include the wetland head boundary (AB and CD) and lateral boundaries (EG and FH). /\ E F Peat Layer Groundwater Aquifer GH, EC, and DF are no-flow boundaries r 1 EG, FH, ACDB, and CD are specified head boundaries r2 Figure 3 3. Setup of the saturated groundwater problem around a wetland PAGE 44 28 As to the boundary condition controlled by the wetland water table, the initial water table of the wetland and the pumping discharge were taken as the only inputs of the coupling model. During pumping, the pumping discharge lowered the wetland water table, and it also triggered an inflow from the aquifer to the wetland through the peat layer. In each simulation time step, the pumping discharge, estimated inflow from the aquifer to the wetland, initial water table of the wetland, and relation between the volume of the wetland and wetland water table were used as inputs in a mass balance equation that was solved to determine the water table at the end of the time step The estimated water table at the end of the time step was put back into the saturated groundwater model to adjust the inflow from the aquifer to the wetland estimated initially. Therefore, several iterations in each time step were necessary to get an accurate estimated inflow from the aquifer to the wetland and a convergent solution of the water table in the wetland at the end of the time step When pumping was stopped the pumping discharge was taken as zero. Figure 3-4 shows the procedure of finding the inflow from the aquifer to the wetland and the coupling boundary. 3.3.2 Finite Element Model of Saturated Groundwater Flow To obtain a finite element approximation of Eq. 3-1, two spaces, Sand V, are defined (Hughes 1987): (3-2) V = { ulu E H1 u(x,y,z)/r2 =0} ( 3 3) where H1 is Sobolev space comprised of functions which are square-integrable generalized PAGE 45 Wetland pumping discharge Q 29 Predict a wetland water table hi(t+dt) Simulate groundwater flow based on the predicted wetland water table Find the groundwater inflow to the wetland Based on the water mass balance equation of th f---------,3>-i wetland to find the wetland water table hi(t+dt) No Next time step Figure 3-4. Flowchart illustrating the iteration procedure of finding the coupling boundary condition derivatives through order 1 ; u is a space function that is member of Sobolev space ; ClI'2 is the vertical flux along r2 ; and x, y z are three spatial coordinates. Any function from S space is equal to the flux boundary value along specified head boundary (I'2). Also any function u from V space is equal to zero along head boundary (I'2). Assuming a basis function w is chosen from V space then the following integration of Eq. 3 1 over space domain Q is always true : PAGE 46 30 i w[\7(K\7h)]dQ i w QdQ "i wS, :~, dQ, ija 1,2,3 (3.4) where Eq. 3-4 expresses the governing groundwater flow equation integrated over the modeling domain Manipulation of the first term on the left hand side of Eq. 3-4 begins with the application of Green's theorem, which is an expanded form of the divergence theorem. This converts the integral into two terms, one of which is evaluated only on the surface of the region to be simulated Green s theorem is: I ( \7 )A dQ a I ( u )A ar-I( \7A)dQ (3-5) n r n where A is a scalar and J: is a vector quantity The boundary of volume Q is defined by r and n is a unit outward normal vector to the boundary. Application of the Galerkin criterion (Hughes 1986) and Green's theorem to Eq. 3-4 leads to: i w[(K\7h)u]ar.f. w [(K\7h)u]arf (vw)(KVh)dO r1 r2 n f 1 ah, wQdO= wSdQ s a t n n (36) PAGE 47 31 Because w is a member of V space, it is equal to zero at all specified head boundaries (I'2); thus, the second term in Eq. 3-6 equals zero. If it is further assumed that Q equals zero, then Eq. 3-6 reduces to: (3-7) For an approximate solution of h(x, y, z, t), it is assumed that: h(e)(x,y,z,t)=tN?)(x,y,z)d/t), NiEV, hES (38) i = I w(e) = t N;(e)(x,y,z) ci NiE V (3-9) i = 1 where h PAGE 48 32 isotropic, that is to say, values of saturated hydraulic conductivity in the three coordinate directions are constant wit hjn an element (but can vary from one element to the next), the contribution of any element e to the residual at a node i associated with the element can be d escribed through the following system of eq u ations: adl(e)(t) R?) (t) dl(e)(t) at =[K PAGE 49 33 and: (3-13) where y < e > is the volume of element e. Eq. 3-10 represents a system of equations pertinent to a single element, and there is an equation for every node in the element. When equations for all elements in the mesh are combined, the global equations are obtained for the system: ad,(t) R,(t) d,(t) at = [K] + [c ] [F] (3-14) Rit) dit) adit) at By setting the residuals equal to zero, we have : ad,(t) d,(t) at [K] + [c] = [F] (3-15) dit) Bdit) at Eq. 3 15 is a system of ordinary differential equations, whose solution provides values of head d and the derivative with respect to time Bd/Bt, at each node in the finite-element mesh. PAGE 50 34 Although several methods are available for solving this system of equations, it has become standard practice to use the finite-difference method: a\e) = d(t+/1t)-d(t) at 11t If a variable, w is defined such that: Then: e-t W = -/1t d(e) = (I-w)d(t)+wd(t+l1t) which can be extended to the vector of [d]: [d (e)] = (1w) [ d(t)]+w [d(t+/1t)] (3-16) (3-17) (3-18) (3-19) If Eq. 3-19 and Eq. 3-18 are substituted into Eq. 3-15, the finite-difference formulation for the transient saturated flow equation is as follows: : [ c] + w 11 t [K] )[d(t+l1t)] = ( [ c] ( 1 w) 11t [K] ) [d(t)] +/1~(1w ) [ F(t) ] +w[F(t+l1t)] (3-20) PAGE 51 35 Eq. 3-20 rep r esents a system of linear equations with unknown variables of ~(t), and where N / x,y,z) are the basis functions in a tri-linear hexahedral element domain. Nlx, y,z) does not change with time The most common basis functions used are hat funct i ons which are linear functions of space variables. The solution of Eq 3-20 begins by specifying the initial values of [d] (i.e., the values of head at time t0= 0). Next the system Eq 3-20 is solved to obtain values of { d(t+~t)} at the end of the first time step 1 This solution process is repeated for the next time step and so on until the specified time duration of simulation is satisfied. Depending on the choice o f w several different subsets of the finite-difference formulations are defined When w=O it is the forward-difference method is obtained : ( [ c ] ) [d(t+At)] = ( [ c ] At[KJ)[d(t)] +~t[F(t)] (3-21 ) When w=l/2, the central-difference or Crank-Nicholson method is obtained: ( [c]+ ~At [ K ] ] [d(t+~t)]=( [ c ]-( 1-~]At[KJ][d(t)]+ ~t 1FU)]+[F(t)]) (3-22 ) When w=l, it is the backward-difference method is obatined: ( [ c ] + dt[ K ] ) [d(t+dt)] = [ c ] [ d(t)] +At[F(t)] ( 3 23 ) PAGE 52 36 3.3.3 Adaptation of the Saturated Finite-Element Model to Describe the Free Surface Several papers have been written on the groundwater free-surface problem (e.g.,Taylor 1967, Neuman 1970, Finn 1976, France 1976, and Huyakom 1986). Most of these papers discussed the steady-state free-surface problem except for France (1976). This thesis examines the transient problem starting from the basic ideas that Neuman (1970) presented for steady flow. Basically both France (1976) and Neuman (1970) tried to satisfy the two boundary conditions and are consistent. This chapter is also based on the two conditions and adapts the method they used to locate the free surface The primary dependent variable to be used in modeling the free surface presented here is the hydraulic head h defined as: h = z + "1 y (3-24) where z is the vertical coordinate; y is specific weight of a unit volume of water pg and \JI is the pressure. Along the free surface, there are two co nditions which must be satisfied. First, \Jf=0 ( pressure atmospheric), and hence for a given position, his prescribed on it. That IS: h =z (3-25) For the steady-state scenario, the second condition is simply that there is a specified flux across the free surface boundary (i.e. a zero flux or a flow imposed equal to net recharge). Thus, in addition to Eq. 3-24 at the free surface, there is: PAGE 53 37 ah-; ah-; ah K-z +K-J +K-k= q X ax y ay Z az O (3-26) If steady state has not been reached and the boundary is moving with a velocity normal to its instantaneous configuration of Vn, then the quantity of liquid entering its unit area is modified and: ah-; ah-; ah -K z + K -J + K -k = q + VS= q X ax y ay Z az O n n (3-27) where i, j k are unit orthogonal vectors; S is the specific yield coefficient (or the effective porosity) relating the total volume of material to the quantity of liquid which can be drained from it, and PAGE 54 38 3.4 Validating the Code of Three-Dimensional Deformable Finite-Element Groundwater Flow Model 3.4.1 Theis Solution To validate the code of the three-dimensional deformable finite-element groundwater flow model a comparison was made of results generated by the numerical model and those created with an analytical model describing unsteady flow to a well in a phreatic aquifer. Figure 3-5 shows a schematic of unsteady flow to a well in a phreatic aquifer. i Q ---------! f-~.i.------------------------S? ---....... s -----Figure 3-5. Unsteady flow to a well in a phreatic aquifer at time t. wheres= drawdown(L), H0= initial saturated aquifer thickness (L), Q = pumping discharge (L3/T), and r = radial distance from the well (L). As shown in Figure 3-2, the boundary condition along the phreatic surface is nonlinear and the shape and position of the phreatic surface are unknown. As a result an exact analytical solution that describes groundwater flow is not known. However, by the employing the Dupuit assumptions and assuming recharge is zero the unsteady flow to a fully penetrating well in a phreatic aquifer can be desc1ibed in radial coordinates by: a2h2 1 ah2 2s ah --+---=-ar2 r ar K at (3-28) PAGE 55 39 Initial condition: Boundary condition: h(t) I (r=<) = H0 Here, S is the phreatic aquifer specific yield (or effective porosity, n e). If drawdowns are sufficiently small compared to the average depth of flow, (say 0.02H0), then the approximate solution to this problem is the Theis solution (Bear 1976): s=_Q_W(u) 41tT where u = r2S/(4T t); W(u) is the well function, and T=KH0= initial transmissivity. (3-29) The conditions that define the test problem used to validate the numerical groundwater flow model include: a pumping discharge, Q of 1.5 m3/hr; an initial saturated phreatic aquifer thickness H0 of 20 m; a specific yield S ( or n e ), of 0.1; and a hydraulic conductivity, K of 15 m/day. For the three-dimensional finite-element model the specific storativity, S5 is 0.0001 m -1 ; where S5=pg(a+nP); n equals porosity ; a equals solid matrix compressibility [LT2M -1]; g is acceleration of gravity [LT2 ] ; and pis fluid compressibility [LT2M -1]. Because of problem symetry, modeling was limited to simulating hydraulic head variations over a wedge section of the aquifer using a 21xl lx3 mesh (Figure 3-6). The boundary condition controlled by well water table, the initial water table of well and the pumping discharge were taken as the only input. The pumping pulled water out of the well, on the other hand the groundwater supplied well. In each simulation time step, based on the known water table in the beginnin g of time step, and an estimated water t a ble of well, a new discharge supplied from groundwater to well was calculated from the groundwater modeling. Based this new calculated supplying discharge to the well, the pumping discharge and the circular area of the well, a more acc urate well water table could PAGE 56 40 be estimated again. The groundwater flow model was solved again. After many iterations, a convergent and accurate well water table could be obtained before the modeling of next time step z x--f\ 20 1 5 N 1 0 5 0 Figure 3-6. Finite-element mesh of the wedge used to simulate of a well pumping in phreatic aquifer Because the mesh is deformable, the known pressure at the free surface and flux were used to locate the free surface. During calculation a slice successive over-relaxation (SSOR) method was used as an iterative approach to locate the free surface. The results illustrated in Figure 37 validate the coding of the numerical model. 3.4.2 Groundwater Flow around a Circular Wetland To validate the three dimensional finite element saturated groundwater flow model again the numerical model was used to simulate the groundwater flow around a circular wetland and induced by a flood wave in the wetland. It is assumed that the wetland stage fluctuations can be described by a nonperiodic Gaussian fucntion of time as Eq. 3 30 : ( t t,l /{ 1 2o2 ro, t)=--e lino (3 30) PAGE 57 0 1 0 0.08 :[ 0 .06 C 3: 0 t \\ 'O 3: 0.04 0 0 .02 0 0 41 ....... t'.. b..,. lb--'--;--.._ -. ---5 10 15 20 25 30 35 The radius to well FEM ( 1 = 5 0hr ) ana (1=5. 0 hr) FEM (1=2. 5hr) ana (1=2. 5hr ) Figure 37. Comparison of the numerical and analytically simulated head for the example problem of a transient pump test in a phreatic aquifer Unlike the coupling boundary condition that was used in the last case, a flood wave boundary condition used in this numerical model here The solid line and the dash line in Figure 3-8 are the three-dimensional finite element model simulated results. The triangle points and the diamond points are the groundwater analytical solutions from Eq. 2-24 obtained in Chapter 2. The circular wetland which has a radius of 50 meter was assumed. All other system parameters are given in the Table 3-3 The numerical model and analytical solution matched very well. 3.5 A Numerical Model Application to Simulate the Interaction Between a Wetland and an Aquifer The finite-element model was applied to simulate the interaction between a wetland and a contiguous phreatic aquifer stressed by an active productive well. A horizontal irregular quadrilateral mesh was used to define the irregular boundaries of the wetland and the area beyond the wetland. Figure 3-9 illustrates the numerical mesh corresponding to the PAGE 58 42 100 9 100 8 100 7 --3D num d=2m ---3D num d=5m QI :c 100 6 Cll I100 5 ai ... analytical d=2m ~ analytical d=5m iii 100 4 :r: "C 100.3 C: :, 100.2 e C, 100.1 100 r... .J. 99.9 2 7 12 17 Time (hr) Figure 3-8 Comparison of the numerically and analytically simulated groundwater table Table 3-3. Parameters and their values Parameter Value Parameter Value H00 (m) 100.0 T (m2/h) 0.635 L\r (m) 0.25 s 0.01 L\t (hr) 0.02 t o (h) 5.0 k (m/hr) 0.00635 (J 0.4 horizontal and vertical extent of the phreatic aquifer. Beyond the wetland, the vertical extent of the mesh defines the thickness of the phreatic aquifer. Over the inundated portions of wetland, the thickness of the mesh defines the vertical distance between the wetland bottom and the underlying aquiclude The wetland is shown in the middle of this area The model domain is 104x104x10 meters. The four lateral boundaries and bottom boundary are assumed to be impermeable. An active production well is located within the model domain. It was assumed that the initial groundwater table and the water stage of wetland were the same. The deepest water depth in the wetland is 1.0 meter. Model runs were conducted to PAGE 59 43 simulate the exchange of water between wetland and aquifer. Such interactions between the wetland and the aquifer were examined under transient pumping conditions and precipitation Since the real-time stage of the wetland is used as a boundary condition over the wetted part of wetland, the horizontal and vertical extent of the mesh must be adjusted to trace the vertical and horizontal movement of the transient wetland water surface. -5 2005E-04 Upper-ex t e nt of the phreatic aquife r defin e d by the bottom elevation in u ndated wetland 10 ~5. 2005E-04 Figure 3-9. Three-dimensional finite-element mesh for the simulation of wetland and aquifer The modeled area was divided into a 21x2lxl0 hexahedron mesh. In the middle, a 9x9 quadrilateral mesh was used to represent the wetland area (Figure 3-9). Hydrologic conditions that define the test problem include: a net precipitation rate of 0.02 m/hour; a specific yield S of 0.075. a hydraulic conductivity K of 10 m/da y ; and a pumpin g discharge Q of 1.554 m3/hr. Finally data from Figures 3-10 and Figure 3-11 were used to define the relationship between the wetland water stage, the surface w a ter area of wetl a nd and the corresponding wetland volume. Figure 3 12 illustrates the initial condition wh e re it was assumed that the wetland surface-water stage and phreatic surface elevation were the same Figure 3-13 illustrates a PAGE 60 44 plan view of the numerical mesh at time zero. Figure 3-14 through Figure 3-21 represent various depictions of head contours and iso-surface plots of head predicted at hour 5 of the simulation These figures show that water moves laterally into the wetland and pumping well because the same net rainfall forces the water table to increase faster than in the water stage in the wetland. The directions of flow also depend on the active production well. If production is large, more water flows to the well instead of flowing to the wetland. Otherwise, some water flows towards and into the wetland as groundwater inflow. Beneath the wetland, the hydraulic head increases with time and with increasing depth (Figure 3-17 and Figure 3-18) which shows that water flows from the aquifer into the wetland. As to the pumping well the contours are almost vertical and parallel to each other near the well (Figure 3-20 and Figure 3-21). Groundwater moves towards the well in a horizontal direction because the well fully penetrates the aquifer. 700 -,---,---,---,---,---,---,---,---,---,---,---,---,---,---,----.----, .!!!., //_ 500 / Q) L,_V-+--+--+--+--; 3: 400 v -,____,___,,___,___,---, 0 .. ,../"--+-+-+---+-+---t--f--f C"d 300 / / : 200 f--t---f--it---+70 +v--+---+---+---+---+--+--+--+--+--+--I u :r--~-+----+---+--+--+--t---t---t---f--it----1--t--f {g 1 oo Vt--+--+--+--+--+--+---+--+-+--+---+--___, Q t<-/-+/-+--+--+--+--+--+--+--+--+--+--+--+--+---+---; .s=: f0 0 2 0.4 0 6 0 8 1 1 2 1.4 1 6 The water depth of wetland D (m) Figure 3 10. Relationship between the area o f wetland water surface PAGE 61 500 E45o I u -400 > 350 ::J 300 Q) .c I-250 200 >--V V / 45 V / /,-/ / / ,( v-r400 450 500 550 600 650 The surface area of wetland S (s-m) Figure 3-11. Relationship between the deepest water depth and the volume of wetland. Figure 3-12. Illustration of the initial condition where it was assumed the wetland surface water stage and the phreatic surface elevation were the same PAGE 62 46 9 9999E+0 I 4 9999E+0I 5 .2 005E-04 Figure 3-13. Plan view of the numerical mesh at time zero z ,A, 5 2005E-04 1 0 0 5 .2005~ Figure 3-14 Three-dimensional finite element mesh and simulated head contours at hour 5. PAGE 63 47 Figure 3-15. simulated head contour at hour 5 .\ .200'6- PAGE 64 48 Figure 3-17. Three-dimensional finite element mesh and head contour on a cross section through wetland Figure 3-18. Simulated head contour on the vertical cross section through the wetland PAGE 65 49 Figure 3-19. Modeled iso-surface plot of head equal to 11.1 m, 11 m and 10.5 m Figure 3-20. Mesh and modeled head contour with vertical cross section through the pumping well PAGE 66 50 Figure 3-21. Head contour (m) on the vertical cross section through the pumping well A computer with a Pentium 1.4GHz processor took 4 hours to complete the 5-hour simulation using a time step of 0.1 hours and 3 iterations per time step. To use this model over a real system, a finer mesh would be desirable. Unlike most rectangular finite-difference grids, the finite-element grid used here precludes the use of more efficient solution algorithms that divide simulation domain into three individual directions. In order to use this model to investigate wetland and aquifer interaction on a large scale, it would be desirable to introduce a simplified me s h that would in tum be more amenable to efficient solvers. With computer hardware developing, the simulation speed can be accelerated consequently. To improve the simulation speed of the model from fundamental a kind of hybrid method which combines finite-element method and finite difference method could be an effective choice for a large scale practical problem. 3 6 Application of Three-Dimensional Deformable Saturated Groundwater F l ow Model to SVS Wetland 3.6.1 Site Description-SVS SV5 is an i so lat e d wetland found in pine flatwoods located in southeast Florida. It is part of the Savannas wetland system located south of the city of Jensen Beach in Martin PAGE 67 51 County SV5 is located approximately two km west of the Atlantic Ocean and is adjacent to a large municipal well field. The South Florida Water Management District (SFWMD) monitors this wetland system. There are 13 surficial aquifer municipal water-supply wells located from a few hundred meters to four kilometers from SV5 (Figure 3-22). The wetland is roughly circular in shape, with a diameter of approximately 60 meters (Wise et al. 2000, Walser 1998) N 6 r E 6 SV5 0 6 6 6 6 6 6 o;ro o 0 66 <> M SFWMDWell 6 Monitor Well 0 Interior Well 6 6 Figure 3-22. Well Locations (Switt et al. 1998) The surficial aquifer that underlies the wetland is composed of three layers. The first layer, which is located from land surface to 12-18m below the surface, consists of white, gray and brown, largely fineto coarse-grained quartz sand intermixed with shell beds. It receives recharge directly from precipitation. The second l ayer, which is 3-6 m thick, is comprised of fine to very fine sand with small amounts of she ll and clay. The third layer, which extends to a depth of between 37 43 m, consists of limestone intermixed with sand and shell. The upper confining layer lies beneath the surficial aquifer. This layer consists of low - PAGE 68 52 permeability rocks that are primarily elastic. The Floridian aquifer system lies under the upper confining layer. The system is quite thick and consists of limestone of variable permeability (Fetter 1994 ). The climate in the area is considered humid sub-tropical, and the area receives an average of 1.5 m of rainfall per year (Walser 1998). 3.6.2 Field Methods and Wetland Aquifer Interaction Test The SV5 wetland aquifer interaction test was conducted (Wise et al. 2000). Eighteen monitoring wells were installed. Wells were placed along transects running from east to west and north to south Six wells were placed outside the wetland on the western transect and two on each of the other three transects. The wells were identified by the transect designation followed by the distance location from the center of the wetland i.e., well 37 is the well located 37 m from the center of the wetland along the west transect. SFWMD installed two wells and a staff gage at the site. The wells measure groundwater. The groundwater well is screened at a depth of six meters and grouted from the surface to a depth of five meters. The water and peat depths in the wetland were measured along four transects through the center of the wetland at 45 angles from one another. Peat and water depth measurements were taken at three-meter increments. The peat depth was determined by pushing a piece of 1 cm rebar into the peat at a steady rate until significant resistance was felt. This method is most easily applied to sites such as SV5 that are free of significant tree roots. The depth of the water table was determined by measuring the water level in the exterior monitoring wells. Elevation data for the sites were measured utilizing a Topcon Self-Leveling Rotating laser level. The South Florida Water Management District staff gage on the site was used as a relative reference for the elevation data and all elevation data were expressed in terms of NGVD of 1929 PAGE 69 53 The water and peat depth information was entered into Surfer a three-dimensional graphical interpolation program, to develop three-dimensional contours of the wetland A stage -volume curve to describe the wetland was developed by utilizing Surfer to determine the wetland volume at varying water levels. A regression curve was fitted to this data to determine the equation that describes the stage-volume curve ( Switt et al 1998, Wise et al. 2000). This curve describes the top of the peat layer. The process was repeated to develop a curve that describes the bottom of the peat layer. The goal of the test was to stress hydraulically a wetland by pumping the standing water to lower the water level a predetermined amount which varies from wetland to wetland, and then monitoring the response in wells and subsequent recovery of the water in the wetland For the wetland test to be effective, the wetland must be a discharge wetland at the end of the pumping phase. The wetland test was conducted utilizing a 10.2 cm centrifugal trash pump to pump the water from the wetland until the water level in the wetland dropped 30.5 cm. The water was pumped into a drainage ditch located approximate ly 100 m from the wetland. This distance was chosen to ensure that th e water pumped from the wetland did not reenter the system The discharge rate of the pump was approximately 70 m3/hr. The discharge rate was measured by determining the amount of time it took to fill a five-gallon plastic bucket at the end of the discharge hose. It took nine hours to lower the wetland water level 30.5 cm. During the testing water levels in all of wells were manually measured every hour with the exception of the surface water well which was measured every 30 minutes, utilizing a Slope Indicator water tap e. Onc e the pump was turn e d off, water levels in all of the wells were manually m eas ured eve ry hour for 22 hours Exterior and interior well wat e r levels PAGE 70 54 were recorded for the wells in the wetland area. 3.6.3 StageVolume Relationships Water and peat depth data were used to produce a three-dimensional contour plot of SV5 (Figure 3-23 and Figure 3-24). The regression equation that describes the volume of water above the peat layer is Vw = 3454S2 -148.45s (3-30) where V w is the wetland volume (m3 ) at s, the stage level (m) (Switt et al. 1998 Wise et al. 2000). The bottom of the peat layer is described by the follow regression equation (Figure 3.25): V =-116s4+607s3+938s2+108s p (3,31) where Vp is the volume above the bottom of the peat layer (m3 ) ; and V P-V w ::=volume of saturated peat. The coefficient of determination equals unity for both of these regressions. Figure 3-23. Top surface of peat layer of wetland PAGE 71 Q) E ::::, 0 > 55 F i g u re 3-2 4. Bott om s u rfac e o f peat layer Stage-volume curve 1600 --~~~~~~~~~~~-,--,--,---. / 1400 -----+--+--+-+-+-+--+--+--+--+--+--+----Y--/----i 1200 +-+-+--+--+--+--+--+--+--+---1--t--,i,.<1~--r_ -__ 1 000 1----f-f-+--+--+-+-+-+--+--+-----J.,lit'-/ ... ----------800 1----+-+-+--+--+-+-+-+--+-_,-.,/r--t----t--l--l 600 ~ =t==t== t =t= t = t = t =~./2'.~l("==::==t::;;~2""'::'.ll,/'.::::~ 400 +-+-+-+--+---+--+-__..,,-.L--'/~f----J.--1--.C:::~'eb --++---+-l----+l----l-l 200 t--t---t----t---t-==lr"" ~ -= -0 t=::t;~~---ctp=t=t=!=t=:j 0 0 1 0 2 0 3 0 4 0 5 0.6 0.7 St a ge VwS VpS F i g ure 3 -2 5 R elationship be t ween we tl an d vol u me and water sta g e PAGE 72 56 Switt et al. (1998) used a similar formula (Eq 3-32) to calculate the interaction between a wetland and the underlying groundwater: (3-32) where SP= exchanging flux (m3/day), K5= hydraulic conductivity (m/day) G = groundwater level (m) W = wetland water level (m), L = average peat thickness (m) and SAw = wetland wetted surface area (m2). Over the de-watered area of the wetland seepage decreases from the above value at the edge of the water to zero at the edge of the wetland. Switt et al. (1998) assumed a standing water head difference equal to 1/3 of the value use in Equation 3-8 Thus seepage for the de-watered area of the wetland was (3-33) where SANW= wetland dry surface area (m2 ) The hydraulic conductivity of the organic soil in the peat layer is critical in wetland systems because it controls the amount of groundwater that flows from the aquifer into the wetland or from the wetland into the aquifer (Truesdell 1995) Determining a value of hydraulic conductivity for the soils in wetland system is particularly challenging 3.6.4 Creation of Finite-Element Mesh The model domain was approximately three times the si z e of the wetland. Because the wetland is appro x imately radically symmetrical it was possibl e to develop a highly refined 4.0 numerical wedge to represent the wetland aquifer system. The aquifer was PAGE 73 57 divided into 10 layers to obtain a 1 lx26x3 mesh (Figure 3-26) The peat layer was divided into three layers producing a 3xl lx3 mesh. The mesh was divided and fitted along the shape of peat layer Wetland Peat Layer ,1 15 .I 10 ii ii 5 V 0 tT 4 0 2 0 X Figure 3 26 Finite-element mesh of wedge of wetland and aquifer for the saturated flow model 3 .6. S Sol u tion Procedure 5 0 0 .... To simulate the wetland-aquifer interaction observed at SV5 the model must correctly simulate the movement of the free surface as changes occur with the wetland water stage. Two boundary conditions must be satisfied on free surfaces; first there can be the flux across the surface must be represented and second the pressure must be atmospheric. The initial fixed domain solution usually will not satisfy both boundary conditions simultaneously. Therefore, the trial free surface is adjusted and a solution obtained for the revised fixed domain This process is repeated until the required adjustment in the free surface is considered negligible. Thus the solution of the non-linear free surface problem is achieved by a series of solutions to linear problems with prescribed boundaries. Finally because the wetted part of w e tland may chan g e significantly with variations of wetland water stage the finite element mesh must be adjusted horizontally and vertically between time PAGE 74 58 steps Before the pumping test began, the SV5 wetland water surface was about 60 meters in diameter. However, after pumping stopped it was less than 30 meters. Thus, the horizontal inundated area changed from 2750 m2 to 750 m2 To model the influence of a changing horizontal free surface required that the finite-element mesh be adjusted horizontally The zone of mesh adjustment would be within the peat layer and include an upper boundary representing the free surface or the surface water stages and a lower boundary corresponding to the vertical limit of the peat layer. Adjusting the mesh horizontally keeps the same number of nodes on the boundary. Thus the physical domain of the mesh may change horizontally and vertically. However, the calculation domain does not. After each simulation information defining the location of every node is updated. Since the wetland boundary was assumed to vary over time the total head was taken as a boundary condition Thus at any time, the wetland boundary can be viewed as an equipotential line. Water moves vertically to and from the equipotential line represented by the free surface boundary The vertical hydraulic gradient at the wetland boundary times the product of the wetland area and hydraulic conductivity equals the volume of water exchanged between the wetland and the aquifer. 3.6.6 Model Calibration and Results Because the aquifer underlying SV5 is composed of sand a specific storativity of S s = 10 -5 m -1 was used (Bear 1979) The saturated water content and residual water content in the peat layer were e s = 0.7 and e r = 0.6 In the aquifer, e s = 0.34 and er= 0.18, based on Walser (1998). Water stage level observation data and groundwater monitoring data collect during the pump test and many hours after pumping stopped were used to calibrate the PAGE 75 59 numerical model. In Krabbenhoft's (1996) interaction model between groundwater and lakes an anisotropy ratio of 100 was obtained after several trial and error simulations. Also an anisotropy ratio of horizontal to vertical hydraulic conductivity of 10 was obtained by Winter (1978) in his numerical simulation of steady state three-dimensional groundwater flow model near lakes. In this finite-element model it was assumed that the hydraulic conductivity of the peat layer was different from the hydraulic conductivity of the underlying aquifer, and the aquifer was also assumed anisotropic Furthermore it was assumed that the hydraulic conductivities and anisotropies of both layers were spatially uniform within each layer throughout the aquifer zone. By adjusting the anisotropy ratio (K/KJ and checking for the best agreement between modeled and observed data it was determined that the best results were obtained when K/Kv = 1.0 K x x = = 8.7 m/day, and = 8 7 m/day in the aquifer while~= 0 2 m/day in the peat layer. This number is different from the value of 0 03m/day obtained by Switt (1998) and 0.087m/day obtained by Wise (2000). The possible difference may be caused by several reasons. First the peat layer depth used in the numerical model is not a constant. When the water stage in the wetland is lowered the peat layer depth in the inundated area is much greater than the average value of 0 25 mused in the analytical model. The variation of groundwater head immediately beneath the peat layer may not be the same value estimated usin g a groundwater head variation at a depth of 6 meters beneath of the wetland Thus the second r e ason may be th a t the he a d distribution beneath th e we tland obtained from the num e rical model is differ ent from the valu e s used by th e a nalytical model. Since the aquifer beneath the SV5 wetland is composed of sand the specific PAGE 76 60 storativity is extremely small The heads along the depth are very different. Figure 3-27 shows the head distribution beneath the wetland Figure 3-28 shows the head variation with time at different depths beneath the wetland. The piezometric head changes relatively quickly in the deep zone but the head does not change much near the wetland. Figure 3.29 shows the variation of the inundated area of the wetland The area changes from 2,600 m2 to 700 m2just after pumping. Figure 3-30 shows the variation of the interaction discharge. The discharge supplied from the aquifer to the wetland increases quickly and then decreases slowly after pumping was stopped. Figure 3-31 shows excellent agreement between the simulated results given by the saturated flow model and observed groundwater levels Table 3-4 to Table 3-8 show the modeled data used to calculate the interaction discharges at 5 different times. Table 3 8 shows the sensitivities of the head relative to changes in the hydraulic conductivities in the pear layer and the aquifer and to changes in the specific storativity. The hydraulic conductivities in both the peat layer and the aquifer play the most important part in the change of head distribution and interaction discharge. The sensitivity of the specific storativity is relatively small. PAGE 77 61 z h_x IJ. W -10 0 1 "' --1n I \ 7 : -5 l:'i} l -0 -50 0 50 Figure 3-27 Vertical cross-section view of three-dimensional finite-element mesh and head from the saturated flow model at a time of 87 hours. -~ 16H \>-+---+--+--+--+---t--t-; ~ ; t .,ETD;m-:[ -g D1pth=, .5m I _ Q)!,!Q) 15.9 -11-:Dijr::J::d~;;t::1::J~::t=~r::1~:J D1pth 11.Sm E L---.--2 ----15. 8--------:~7V 15 7 1---+--+---,f---+--+------+--+--+---+----, 0 20 40 60 80 100 Time ( h our) Figure 3 -28. Piezometric head at positions beneath the wetl a nd with time PAGE 78 2 5 -0 1 5 i B 1 a -0 Q) ai 0 5 0 {\ ~ / \ I I 0 62 l--~ .----...-20 40 60 80 Time (hour) Figure 3 -2 9. Variation in the inundated area of wetland during simulation 3000 2500 \ r \ E \ ';;;-2000 \ <{ -0 --\ 1500 -0 C -----i----::, .!: 1000 ___ l..---"' t----t---~=---t----+--+--t --------500 I----+--+--+---+--+----,,----+--+---< 0 20 40 60 80 Time (hour) Figure 3-30 Discharge between wetland and aquifer with respect to time PAGE 79 63 \\ 4 2 \I' -tt-----+--+----+---+--t----+---+---------- PAGE 80 64 Table 3-5. Dat a for calculating the discharge to the wetland from the aq ui fer at 20 hours Elements Area Head in Peat Head on Peat Distance o f Two of Peat Layer ( m 2 ) Layer (m) Surface (m) Nodes ( m ) 1 8 36449 15 7953 15. 7427 0.302115 2 25. 09 35 15. 7956 15. 7427 0.306344 3 41.8223 15 7960 15. 7427 0 306027 4 58 5508 15. 7964 15 7427 0 302183 5 75 .2802 15.7982 15. 7427 0 294043 6 92. 0102 15. 8003 15. 7427 0 287151 7 108 741 15. 8036 15. 7427 0 275314 8 125.472 15.8070 1 5 7427 0 264839 9 142.204 15.811 4 15.7427 0 250437 10 158 937 1 5 8162 15. 7427 0 236900 Q to wetland at 20hr 1.65871 m1/hr Table 3-6. Data for calculating the discharge to the wetland from the aquifer at 40 hours Elements Area Head in Peat Head on P e at Distance o f Two of Peat Layer (m2) Layer ( m ) Surfac e ( m ) Nodes ( m ) 1 10. 9288 15. 8222 15.7797 0.302417 2 32.7864 15. 8225 1 5 7797 0 307252 3 54. 6436 15. 8225 1 5.7797 0 304707 4 76.5010 15 8229 1 5 779 7 0 298378 5 98 3597 15 8243 15 7797 0 28909 7 6 120 219 15. 8264 15. 7 797 0 277635 7 142 .080 15. 8288 15. 7797 0 266032 8 163 941 15. 8321 15.7797 0.249581 9 185 804 15. 8359 15 779 7 0 233460 10 207 669 15. 8406 15. 7797 0 .2 16 5 48 Q to w e tland at 40hr 1 86649 m1/hr PAGE 81 65 Table 37. Data for calculating the discharge to the wetland from the aquifer at 60 hours Elements Area Head i n Peat Head on Peat Di stance of Two of Peat Layer ( m2) Layer ( m) Surface ( m ) odes ( m ) I 13 1 2 19 15 .8463 1 5 8115 0 .302649 2 39.3657 15 .8466 15.8115 0 .307800 3 65 .6089 15 .8464 15.8115 0 303651 4 91.8530 1 5 8469 15.8115 0 .295196 5 118 .098 15 .848 0 15. 8115 0 .285844 6 144 .345 15.8498 1 5 8115 0 .27 126 8 7 170.593 15.8519 15 8115 0.255422 8 196 843 15 .8544 15. 8115 0 .239015 9 223 .095 15.8579 1 5.8115 0.219572 10 249.349 15 .8626 1 5.8115 0 199788 Q to wetland at 60hr 1.94458 m3/hr Table 3-8. D ata for calculating the dischar ge to the wetland from the aquifer at 80 hours Elements Area Head in Peat Head on Peat Distance of Two of Peat Layer (m2 ) Layer (m) Surface ( m) Nodes ( m ) I 15.0334 15 8681 1 5 8391 0 .302835 2 45.1003 15 .8683 1 5.83 91 0.307522 3 75.1664 1 5 .8680 1 5 839 1 0.302700 4 105 .234 15 .8685 1 5 8391 0 .292596 5 1 35 .303 15 8695 1 5 839 1 0 .280887 6 165 .374 15 .8709 1 5 8391 0 .266576 7 19 5 .446 15 8727 15.839 1 0 .247805 8 225 .521 15 .8749 15.839 1 0 .228349 9 255.599 1 5.8777 1 5 8391 0 .208494 10 285 .679 15 8819 1 5 839 1 0 18688 3 Q to w e tland at 80hr 1.93622 m3/hr PAGE 82 66 Table 3-9. Sensitivities calculated by saturated groundwater flow model Parameter (PR) Percentage varied ~PR II h1-h2II (m) II h1-h2ll/~PR K .( mlday) -0.10 -0.85 0.9924E-02 -1.1675E-02 K0( m/day) +0 10 +0. 85 0 9055E-02 l.0653E-02 Kp(m/day) -0.10 -0.02 0.6410E-02 -0.3205 ~(m/day) +0. 10 +0 .02 0.5564E-02 0.2782 s o -0.50 -0. 000005 0.2068E-04 -4.136 s o +0.50 +0.000005 0 1265E-04 2.53 3.7 Comparison with Results of Wise's (2000) and others Parameter values obtained by Wise (2000), ~t= 0.087 (m/day), by Walser (1998) ~and = 0 75 (m/day), kpeat = 0.03 (m/day) for SV5 were put in the three-dimension finite element model in this dissertation, and the simulated results were plotted on the Figure 3-32 When the specific storage was set on 0.01 ( 1/m), both the simulated wetland water table and groundwater are pretty close to the observed conditions. The simulated wetland water table and the groundwater table corresponding to the specific storage of 0.001 (1/m) are much different from the observed results. There are two main differences between Walser's numerical simulation and the models in this dissertation. First unlike a coupling boundary condition used in this dissertation, the fixed observed water table process over time was used in Walser' s numerical simulation. The other difference is that a fixed three-dimension mesh was used in Walser s numerical simulation. When ~qu at 8m/day and ~at at 0.2m/day were used and specific storage was varied the simulated results are plotted on the Figure 3 -33 The simu l ated wetland water tables for both specific storage values are very close to the observed wetland water table results. The PAGE 83 67 groundwater table for the specific storage of 0.01(1/m), however is above the observed results. 4.3 4.25 Q) 4 2 :c m I4.15 ... Q) 4.1 j 4 05 'C C: :::, 4 e CJ 3 95 3 9 3 85 -",;,J \ .------.. .. 0 1 0 20 30 40 50 60 70 80 90 Time (hr) --gwobs gw sim(Ss=0.01 ) gw sim(Ss=0 .001) --wetlobs wetl sim(Ss=0 .01) wetsim(Ss=0 .001) Figure 3-32. Comparison results using Walser s parameter values. 4 3 4.25 I 4 2 i Q) ~ -:c --gw obs m 4 .15 \'-; -gw n(Ss = 0 0001) If.=x ... .~ --~ .,_;/.? gw n ( Ss=0 01) Q) 4.1 -1 --well o bs ; 4 .05 w l-n (Ss=0.000 1) "C y ;::..,-C: .-;,...:. w l n ( Ss = 0 .01) :::, . -...:;:; 0 4 ... -..;...,---CJ -~ 3 .95 ::: ......... ,. 3 9 0 10 20 30 40 50 60 70 80 90 Time (hr) Figure 3-33. 3-D numerical model simulated results over specific storage Table 3-10 shows the hydraulic conductivities in sand layer and peat layer. The variations are rath e r lar ge PAGE 84 68 Table 3-10. Hydraulic Conductivity values for sand and peat layers k in sand layer k in peat layer (m/day) (m/day) Switt et al.(1998) 7-12 0.03 Slug test 2 Rycroft et al. (1975) 10-1-10 2 Wise et al. (000) 0.083 Walser Numerical Model (1998) 0.75 0.03/0.3 Models in this dissertation 8 0.2 3.8 Summary A three-dimensional finite-element numerical saturated groundwater flow model using deformable mesh and real-time wetland-aquifer coupling boundary is developed in this chapter. Considering the area of the water surface changes quickly because most wetlands have a shallow bank slope, the model has the capability of tracking the interaction the interaction boundary between groundwater and a wetland. The model is verified using two analytical solutions one is Theis solution the other is the analytical solution found in Chapter 2 of this dissertation which describes groundwater flow around a circular wetland The model is further used to simulate the interaction problem of groundwater and the SV5 wetland. The pumping from the SV5 wetland and the recovery of the SV5 wetland surface afterwards are closely simulated by this model, the interaction boundaries are effectively traced by the deformable mesh. PAGE 85 CHAPTER4 INVERSE MODEL OF THREE-DIMENSIONAL FINITE-ELEMENT METHOD SATURATED GROUNDWATER MODEL IN SEARCHING FOR THE HYDRAULIC CONDUCTIVITY OF THE PEAT LA YER 4.1 Introduction 4.1.1 Inverse Problem During the past twenty years, mathematical modeling has been extensively used in the study of groundwater problems However, mathematical models require a number of parameters to characterize the aquifer Usually, prior knowledge of the values of these parameters is limited, and they must be quantified by calibration which requires observat i ons of hydraulic heads or solute concentrations. Given a parametric model of the physical system and values of the model parameters the prediction of system output (response) to any input is referred to as the direct problem or the forward problem (simulation). The forward problem predicts unknown system states by solving appropriate governing equations. Conversely the estimation of model parameters given the parameteric system model and an input-output relationship is known as the inverse problem (calibration). The typical inverse problem involves estimating the values of model parameters solving the direct problem with this set of parameters to calculate the predicted system response and then adjusting model parameters until deviations between observed and model predicted responses are in some specified sense minimized (Sun 1994). 69 PAGE 86 70 Trial and error and graphical matching techniques are the most basic and traditional methods to solve inverse problem. They are widely used because they are simple to implement. Because the number and combinations of parameter adjustments are not bounded trial and error is rather flexible but also time-consuming. Additionally the solution is strongly dependent on the skill of the practitioner (Keidser 1991). 4.1.2 Parameterization Most groundwater inverse algorithms adopt either a blocked or a geostatistical (random field) description of spatial variability The first approach divides the region of interest into a number of discrete blocks or zones that are believed to correspond to distinct geological units (Yeh 1986 Carrera and Neuman 1986a b c and Cooley 1977 1979 1982 1983). Each block is characterized by a set of spatially uniform hydrogeologic properties which are treated as parameters in an appropriate inverse problem Thus the unknown distributed paramters can be represented by a number of constants. The weakness of this method is obvious i.e., although we have some of the geological information is available it is still difficult to divide the flow region into a limited number of zones. If the zonation pattern is incorrect the identified paramter values also will be incorrect (Sun and Yeh 1985) The geostatistical alternative views the properties of interest as stationary random fields that vary relatively smoothly over space (Hoeksma and Kitanidis 1984 D a gan 1985 Yeh 1996 and Zhan g 1997 Michalak 2 004). In the g eostati stica l interpolation m e thod of kriging the unknown param e t e r is regarded as a stochastic field described by some statistical parameters, such as the mean variance, and correlation distance. Then the identified parameter vector consists of only a few statistical parameters if the aquifer is PAGE 87 71 basically homogeneous. The problem of overparameterization may be avoided, and the inverse solution obtained by the maximum likelihood estimate and cokriging is generally stable. The correctness of the inverse solution however, depends on the correctness of some statistical assumptions and the structure of covariance functions. Although the two approaches are based on different parameterizations they both treat hydrogeologic properties as spatial functions. This suggests that it should be possible to formulate a general inverse theory which encompasses both the blocked and geostatistical alternatives as well as hybird methods which lie between these extremes. 4.1.3 Parameter Identification Methods Many methods are available for the solution of the common groundwater modeling inverse problem whereby the conductivity or transmissivity of a porous medium are estimated using available observations of head or pressure and other information (Distefano 1975, Chang 1976 Loaiciga 1987 Mishra 1989 Knopman 1989 Keidser 1991 Sun 1995 Yeh 1996, Zhang 1997 Lamorey 1998 Anderman 1999 Michalak 2004). Neuman (1973) classified inverse modeling solution techniques as either "direct or "indirect The direct method or equation error method is used when head variations and derivatives (usually estimated) are known over the entire flow region If the measurement and model errors are negligible the original governing equation becomes a linear first order partial differential equation of the hyperbolic type in terms of the unknown parameters. With the aid of boundary conditions and flow data, a direct solution for the unknown parameters may be possible. In practice, observation wells are sparsely distributed in the flow region and only a limited number of observation wells are available. To formulate the inverse problem by the equation error criterion missing data have to be estimated by interpolation. The PAGE 88 72 interpolation data contain errors resulting from the interpolation, and thus, cause additional errors in the results of parameter identification. The solution approach of the indirect method to solve inverse problems is to formulate the inverse problem in term of an optimization framework, based on the most common output error criterion i.e. output error least squares criterion. Using this approach a set of parameters is identified which best describes observed system states. With inverse modeling the structure and parameters of a model are adjusted simultaneously or sequentially so as to obtain a parameter input system and output relations that characterize observed excitation-response functions of the real system. Since the model structure for the wetland/groundwater system is assumed to follow Eq.5-1-5-3 the problem of determining the hydraulic conductivity of the peat layer from the observed water stage and total head in the aquifer may be taken as a typical inverse problem. The investigation of the inverse problem based on the three dimensional saturated groundwater flow model and the adjoint state method is presented in this chapter. 4.1.4 Computation of Sensitivity Coefficients Sensitivity coefficients i e. the partial derivatives of head with respect to each of the parameters, play an important role in the solution of the inverse problem. In the Gauss Newton algorithm elements of the Jacobian matrix are represented by the sensitivity coefficients, ah/of1 i=l, ... M, l=l, ... ,L. If his the head vector the sensitivity coefficients are ah/of1 1=1, ... ,L. Th e re are thre e m e thods used in calculation o f sensitivity coefficients. There ar e the influence coefficient method, sensitivity equation method and variational method The influence coefficient method determines the sensitivity coefficients by perturbing each parameter once at a time. PAGE 89 73 In the sensitivity equation method, a set of sensitivity equations is obtained by taking the partial derivatives with respect to each parameter in the governing equation and in the initial and boundary conditions. The sensitivity coefficients are obtained by solving the new governing equations. The variational method was first used for solving the inverse problem of parameter identification by Jacquard and Jain (1965) and then by Carter et al. (1974, 1982) associated with finite-difference schemes. Sun and Yeh (1985) extended the method to the case of a finite-element scheme. The sensitivity coefficients can be computed using the following equation: (4-1) j=l,2, ... N0 i=l,2, ... ,N0 where Q ; is the exclusive subdomain of node i; Vis the gradient operator; h(x,y t) is the solution of the governing equation; N0 is the number of observation wells; N0 is the total number of nodes used in the numerical solution; and q(x,y t) is the time derivative of q(x,y,t) which is the solution of the set of adjoint equations 4.2 Adjoint Problem for Three-Dimensional Finite-Element Saturated Groundwater Model The focus of this section is the formulation of the adjoint problem which defines sensitivities to be used in determining the hydraulic characteristics of the peat layer and the underlying aquifer. It is ass umed that the thickness of the peat and aquifer are known. As PAGE 90 74 discussed in Chapter 3, three dimensional groundwater flow in an unconfined aquifer is governed by: V (KV h ) Q = S ah s at (XJl,Z)E O subject to the initial and boundary conditions: h(x,y,z,0) = fo (XJl,Z)E Q h J r = J;, KVhnJr = /2, 2 I (4 -2) (4-3) (4-4) f0 is the initial head ; f1 is the specific head condition; f2 is the specified flux condition ; and t r is the time duration. The boundary of volume n shown in Figure 4-1 is denoted by r including wetland head boundary (I'2), no flow boundary (I'1), and lateral head boundary /\ E Peat Layer Groundwater Aquifer /// /// F I I H I AB EG, FH or CD are head boundaries r2 GH or EC, OF are no flow boundaries r 1 Figure 4-1. Sketch showing setup of the saturated groundwater inverse problem around a wetland PAGE 91 75 The hydraulic conductivities of the peat layer and the underlying aquifer are the unknown parameters. Any variation oK of K where K=hydraulic conductivity must cause a variation oh of head h, because h is dependent on K through the governing equation and subsidiary conditions. Taking the variation of Eq. 4-2 and subsidiary conditions, the following is obtained: v-(oKvh) +V ( K'voh) = S a(oh) s at subject to the conditions,;_ oh(x,y,z)l,= o = 0 (XJ',Z)EO oh I r = 0 [ oK 'i1 h + K 'i1 oh]-n I r = 0 2 1 (4-5) (4-6) (4-7) Eq. 4 5 -Eq. 4-7 represent the variational problem of the primary problem in Eq. 4-2. The variational problem relates the variation ok and oh. Assuming h is an arbitrary function having continuous second-order space derivatives on Q and first-order derivatives in [O, trL then multiplying Eq. 4-5 by 'P and integrating the result over time-space domains yields: 't Jes, 'P &h) I 'tiKi J J [s,&h aa~ + 'PV(&Wh) + 'PV(kV&h) ]iKia1 = o (4 8a) a o n PAGE 92 76 If we define_;_ (4-8b) then the first term of Eq. 4-8a vanishes because of condition 4-7. Thus, Eq. 4-8b reduces to: 't J f [s,oh aa~ + 'l'V( 6kVh) + 'l'V(kV6h) ]a'Oat = o (4-9) 0 0 where: I 'l'V(6kVh)a'Q = I 'l'(6kVh)rids -I 6kVhV'l'a'Q (4-10) 0 s 0 and : I 'l'V(kVoh) dQ = I 'I' (kV6h)rids -I kV6hV'l' dQ (4-11) 0 s 0 I kV6hV'P dQ = I 6h[V(kV'l')]a'Q -I 6hkV'l' rids (4-12) 0 0 s Using identities of Eq. 4-10, Eq. 4-11 and Eq. 4-12 Eq. 4-9 can be written as: PAGE 93 77 '1, If [s,Oh aa~ + oh[V(kv'')]-Ok'vhV'P r dt o n '1, + J f Okv'hnds + f kv'Ohrids -f OhkV'P rids dt = 0 (4-13) 0 s s s Using conditions 4 7 and defining: 'Pl r = O, 2 then the second term of Eq 4-13 is equal to z ero. Eq. 4-13 produces : (4 14) In order to identify the hydraulic conductivity k(x,y z) that would lead to the best prediction of observed heads a performance criterion must be defined. A general form of the performance function may be written as : '1, E(h,k) = J f j(h,k; x y,z,t)a'Qdt = O o n ( 4 15) PAGE 94 78 Taking the variation of E(h, k) yields: (4-16) The summation of Eq. 4-16 and Eq. 4-14 gives: If 'P(x y ,z,t) is selected to satisfy the following PDE: s a'P + 'v(k 'v'P) -af = o s ar ah (4-18) subject to conditions: 'P(x,y,z,t1 ) = 0, 'Plr = 0 I then the first integral on the right hand side of Eq. 4-17 related to unknown oh will vanish. Thus: PAGE 95 79 From Eq. 4-19, the following partial derivative can be defined: By introducing the following time transformation : r:=t-t f the adjoint problem can be expressed as s a'P -'v(k v''I') + af = o s ar: ah subject to: 'I' li:=0= 0 (x,y)d1 (4-19) (4-20) (4-21) (4 22) (4 22a) where 'I'= 'l'(x,y,z ,1t-r:). Solving each equation once, all components of the gradient v'E(k) can be obtained. PAGE 96 80 4.3 Searching Objective Function To find the best parameter value, the model structure and model parameters are solved simultaneously or sequentially so as to make the input-output relation of the model fit any observed excitation-response relationsh i p of the real system. Since the model structure of the problem in this study is determined, the problem of determining the hydraulic conductivity of the peat layer a nd the underlying aquifer from the observed water stage and total head in the aquifer may be taken as a typical inverse problem: E(Kest) =minE(K) (4-23) where: (4-24) The traditional method to solve this kind of problems is using the trial and error method based on the output least squares criterion (Eq. 4-23). 4.4 Adjoint State Controlling Equations In this study the adjoint state method was used to solve the inverse problem From the governing equation (Eq. 4-2) and the associated boundary conditions the adjoint formulation of the primary problem may be found using Green's first and second theorems and other transforms : s a'P v'(kv''P) + af = o s a1: ah (4-25) PAGE 97 81 subject to: (4-26) and (4.27) where his the total head; = (x ,y,tr-1") is the adjoint state of h; tr is the time duration of the simulation criterion and f(h,k;x,y,t) expresses the 'goodness' of fit between the predicted head h, and the observed head, hobs expressed throu gh the ordinary least squares (OLS) : J t(h,k;x,y,z,t) = I: 11 h(k)-h0b~I o(x-x0b)o(y-y0b)o(z-z0bs)ou-9 i = I (4-28) where ( x o b s Yobs zobs) is the location of the observation well; and t j (j =1,2;,J) are the observation times. Integrating, Eq. 4-28 over time space domains produces: l 'Jl J E(h,k) = ~f(h,k;x,y,z,t)dD.dt 0 (0) J I (4 29) By solving the governing Equation 4 2 and the adjoint Equation 4 25 once, we may obtain the following functional partial derivative required to solve Equ ation 4-23 is obatined: PAGE 98 82 (4-30) 4.5 Searching Method and Stopping Criterion Since the gradient VE(k) is known, many searching methods including the gradient method quasi-Newton method, and the Fletcher-Reeves method could be used to solve the inverse problem 4-23 (Sun 1985) Two stopping criteria used here are: (4-31) and: (4-32) Using the adjoint method and the Fletcher-Reeves algorithm to update the search sequence, a highly efficient numerical scheme for solving the inverse problem can be formulated 4.6 Field Application on SVS The inverse model was used to estimate the hydraulic conductivity of the peat layer underlying the SV5 wetland (Wise et al. 2000). Since the piezometric head at a depth of 5 meters beneath the wetland was observed, f is defined as: J t(h,k;x,y,z,t) = :E 11 h(k)-h0bJI o(x x0bs)o(y y0b)o(z-z0b)o(t-t) ; : I (4-33) Water level data were measured every hour for a total duration tr, of 87 hours An example PAGE 99 83 of how the hydraulic conductivity K, is updated between iterations is as follows. Assuming an initial estimate of the vertical hydraulic conductivity K of the peat layer was K1=1.00m/day, then aEJa K = -0.0125 by Eq.4-29 and E(K1) = 0.48998. Because it is not possible to match the modeling result with the observed result exactly the final E(k) is not zero The sensitivity was used to decide the direction of searching Assuming E(K) is supposed to be the smallest value, the updated value of K=K2 is assumed to be 0 6 m/day as generated from the formula: (4-34) Based on the new K = K2 = 0 6 m/day the updated gradient aEJa K = 0.0188 and E(K2)= 0.41467. Then a new K is assumed until E(K) is minimized. The final K is 0 2 m/day Figure 4-2 and Table 4-1 are the calculated values from the adjoint method model. In this model the observed groundwater data were used in the output least squares criterion. The observed water stage in the wetland was not used in the model. This is the major reason that the final result from the inverse model is a little bit different from the result obtained in Chapter 3. The saturated flow model in Chapter 3 uses both the observation of groundwater head beneath the wetland and the surface water stage of wetland to calibrate the model. In the inverse model, however it is hard to use the surface-water observation result, since it is a boundary of the model. Because the wetted area of surface water associated with the wetland changes with time the space location of each finite-element point tied to the location of the water surface must be saved before the adjoint state equations are solved. Thus for each assumed hydraulic PAGE 100 84 conductivity in the peat layer the forward problem, and adjoint state problem must be solved first, and then BE/Bk for this assumed value of K could be found. Figure 4-3 shows the procedure for solving an inverse problem using adjoint method. 12 10 -.:8 ..c: 1 w 6 :z: 4 2 0 0 2 OA OB 0~ K in peat layer (m/day) Figure 4-2. Relation between error and peat hydraulic conductivity Table 4-1. Error and sensitivity based on the adjoint method model I I K (m/da~'.) I E(k) I aEJa K I 1 1.0 10.5698 0 .2 6965 2 0 6 8 94509 0 40555 3 0.4 7 17941 1.03330 4 0.2 3.47742 0.31495 5 0 1 4 82392 3.72980 6 0.08 5.59932 3.14520 PAGE 101 85 Model groundwater flow based on the assumed parameters Solve the adjoint state governing equations Calculate aE/aK based on the distributions of groundwater head and adjoint state Calculate E and new ~ + 1 END No Figure 4 3. Flowchart illustrating the iteration procedure of inverse problem 4.8 Summary The inverse problem for three-dimensional transient groundwater flow in a phreatic aquifer is more complicated because the free surface changes with time such that the simulation zone is not stationary. In order to solve the adjoint state governing equations the position of each node must be saved at all time steps In other words for a deformable mesh it is not suitable to use the adjoint method to solve the inverse problem. The second problem is that the observed information related with interaction discharge, such as the surface water stage in the wetland in the SV5 case is difficult to be used. PAGE 102 86 As shown in Figure 4-4, to finish each iteration o f searching paramet ers four ste p s including solving the forward problem and the adjoint state equations must be executed. For one parameter in the case here the above method seems excessive. The "trial and error" method is clearer and faster. I I Iv f-I I I\' I II I I I/ I \ 1~ '\ I -[: ) I I I 1-'r-i I -50 \ I\ I z b_x talcl 0 X I I/ I I II I / I I \ r... I I -\\ I I -2 0 II I -II I I I I I I -10 N B p -I -~\ 0 50 Figure 4-4. Vertical cross -se ction of the adjoint stat e distribution at -r=lO hours PAGE 103 CHAPTER 5 THREE-DIMENSIONAL V ARIABL Y-SATURA TED GROUNDWATER MODEL 5.1 Introduction The limitations associated with sat urated subsurface flow models have inspired hydrologists to develop variably saturated groundwater models that consider both saturated and unsaturated flow (e.g. Voss 1984, van Genuchten 1980, Kirkland et al. 1992, Cox et al. 1994 and Boufadel et al. 1996). Flow in variably-saturated media has been investigated by a number of investigators over the past several decades. Rubin (1968) used the finite difference and alternating-direction-implicit (ADI) methods to solve the two-dimensional Richards equation. Freeze (1971) used the finite-difference and line successive over relaxation methods to simulate three-dimensional variably saturated flow. Neuman (1973) initiated the use of the finite-element method to model two-dimensional variably saturated flow in 1973. Huyakom el al.(1984, 1986) developed the SATURN two-dimensional and three-dimensional variably-saturated flow models using finite elements. Yeh (1992) used the finite-element method and Gauss elimination (or point iteration) methods to solve the three dimensional Richards' equation. His model is the well known three-dimensional FEMW ATER model. Other significant contributions to the area of variably-saturated flow can be found in Brutsaert (1971) Narasimhan ( 1978 ), Cooley (1983) Kirkland et al.(1992) Clement et al.(1994) and Bates et al. (2000), Zhang (2000), Thorenz 2002, and Russo (2003). 87 PAGE 104 88 The general modeling approach has consisted of extending Darcy's law to unsaturated media and combining it with the continuity equation to obtain Richard s equation. Empirical relations are then used to relate unsaturated hydraulic conductivity to the soil moisture content and pressure head [Brooks and Corey 1966; and Van Genuchten 1980]. Although there are many papers on variably-saturated flow modeling, few have investigated variably saturated flow using deformable finite elements. This chapter describes an efficient three-dimensional finite-element variably saturated groundwater model that is suitable for characterizing wetland/aquifer interactions using a deformable mesh and real-time boundary conditions. The model was tested using experimental data from Vauclin (1979) and then applied to simulate surface stage and groundwater level variations measured in an isolated wetland in South Florida A comparison was made between the simulation results predicted by the variably saturated flow model and that of the saturated flow model (described in Chapter 3) for the isolated wetland field study. 5.2 Three-Dimensional Finite-Element Variably-Saturated Groundwater Model 5.2.1 Governing Equation Eq. 5-1 describes variably-saturated groundwater flow assuming that the fluid is incompressible the system is isothermal, and the air phase is infinitely mobile: ~(K.. K B(lf,+z)l Q = ~ (0 S) i = 12 3 a 1J rw a s w' ,J '' x i t (5-1) where S w is the degree of saturation equal to 0/8s; 0 is the water content ; 0 s is the saturated water content; K;j is the saturated hydraulic conductivity tensor [UT]; ~=~(Sw ) represents the relative permeability of the medium which is a function of the degree of water PAGE 105 89 saturation Sw; lj1 is pressure head [L]; and z is elevation head [L]. In order to solve this nonlinear flow equation, constitutive relations must be established that relate the primary unknowns lj1 and S w Such relations are often given by a prescribed functional relationship that is determined experimentally. The following functional relation given by Van Genuchten (1980) was used in this development: O PAGE 106 90 (5-5) Combining Eq. 5-5 and Eq 5-2 produces : a(0Jw) =(0 -0 )m __ l_1m1 (_-_l)_ annPn -1_aP_c +S S-aw at S r n n 2 C a t e S at 1 +an p c (1 +an p c ) ( 5-6) where: aw=acw +z) at at (5-7) Thus Eq. 5-5 Eq. 5-6 and Eq. 5-7 yield: (5-8) By substituting Eq. 5-6 into the left hand side of Eq. 5-1, the following linearized equation was obtained : i,j = 1,2,3 (5-9) PAGE 107 91 The governing Eq. 5-9 was solved using the finite-element method. 5.2.3 Boundary Conditions for Variably-Saturated Model Figure 5-1 shows a hypothetical wetland consisting of a surface-water body lined with a peat layer, and beneath the wetland is an aquifer. The lower boundary of the aquifer is assumed to be an impermeable unit. As to the top boundary condition (IA BJ) constant moisture or a moisture flux can be specified. In the case of rainfall the net rainfall would be a specified flux as long as Ppreci p itatio n < K5 For the case of no precipitation, the moisture flux is equal to actual evapotranspiration. As to the coupling boundary condition controlled by wetland water level, the same way used in saturated flow model was used in this model. The initial water table of the wetland and pumping discharge were taken as the only inputs. During pumping, the pumping discharge lowered the wetland water level, and it also triggered the discharge from aquifer supplying to wetland through the peat layer. In each simulation time step, the pumping discharge, the estimated discharge from aquifer to wetland the initial water table of wetland and the relation between the volume of wetland and wetland water level were put together and the water table at the end of the time step was estimated. The estimated water table at the end of the time step was again put back to the model of the variably-saturated groundwater in the aquifer and peat layer; it corrected the discharge from aquifer to wetland estimated before. Therefore several iterations in each time step were necessary to get the accurate estimated discharge from aquifer to wetland and the convergent solution of water table in wetland at the end of the time step. While the pumping was stopped the pumping discharge was taken as zero. Thus the boundary condition at any time unknown in the beginning depended on the modeling results of the coupling problem itself before this time. PAGE 108 E I I G I /\ 92 Peat Layer Groundwater Aquifer // 1// F I I H I AB, EG, FH or CD are head boundaries r2 GH or EC, DF are no flow boundaries r 1 EA, BF are head, flow boundaries r 1 Figure 5-1. Sketch showing setup of the variably saturated groundwater flow problem around a wetland The determination of evapotranspiration (ET) is a two-step process in the model based on Fraisse and Campbel] (1997). First the daily potential evapotranspiration (PET) is calculated in terms of atmospheric data. Maximum evaporation depends on climatological factors which include net radiation, temperature humidity and wind velocity. The PET represents the maximum amount of water that will leave the s oil system by evaporation when there is a sufficient supply of soil water. After PET is calculated, checks are made to determine if ET is limited by soil water conditions If soil water conditions are not limiting ET is set equal to PET. When PET is higher than the amount of water that can be supplied from the soil system ET is set equ a l to the smaller amount. Maximum evaporation is equal to potential e v aporation. The followin g rel a tionship b e tw e en ET and PET de v eloped by Norero (1969) was us e d in the groundwat e r flow model: PAGE 109 93 ET=-P_E_T_ 1+( :.rl (5-10) and (5-11) where k is a constant ljJ is the soil water potential or metric pressure in the root zone lj, is the value of ljJ when ET= 0.5PET $wmm is the value of ljJ that corresponds to PE/PET=0.05, and $wmax is the value of ljJ that corresponds to PE/PET=0 95. Norero (1969) used $ wmm =1350 joules/kg and $ wmax=-200 joules/kg Figure 5-2 illustrates the relative evapotranspirat i on as a function of soil water potential. 10 8 6 4 2 ET/PET 1 0 -.75 50 .25 0 SOIL WATERPOTENTTIAL, 1/f, BARS Figure 5 -2. Schematic of relative evapotranspiration (ET/PET) as affected by soil water potential ljJ [Norero 1969]. PAGE 110 94 5.3 Finite-Element Method for Modeling Three-Dimensional Variably-Saturated Flow Given the complexity of the unsaturated flow Eq. 5-9, it is not readily seen that it has the same form as the governing equation for satuarted flow (Eq. 3-1). However, if f(05 0r, m, Cl, P c n) from Eq. 5-8 is used in lieu of specific storativity S s in the saturated groundwater equation the form of both the saturated and unsaturated flow equations becomes the same. As a result, the strategy taken to develop an finite-element model for saturated flow (Chapter 3) was applied again to develop an appropriate equation for an finite element model for variably saturated flow. Similar to what was done in Chapter 3, assuming that the basis function w is chosen from V space, the follow equation is always true: i ah, .. = wJ;(0 ,0 ,m,a,P ,n) dO., l,J = 1,2,3 S r C at n (5-12) Manipulation of the first term in the left hand side of Eq. 5-12 begins with the application of Green's theorem, which is an expanded form of the divergence theorem. This converts the integral into two terms one of which is evaluated only at the surface of the region to be simulated. Application of the Galerkin criterion and Green's theorem to Eq. 5 12 leads to PAGE 111 95 (5-13) The first term whic h depends on the moisture flux on the surface is no longer zero as it is in the saturated model. For head boundaries that are used along r2 w is equal to zero, because it belongs to V space Thus with the second term set equal to zero, Eq. 5-13 can be rewritten as: 1 ah, = wf,1 (0 ,0 ,m,a,Pc,n)-dQ s r a t n (5-14) An approximate solution for h(x y, z, t) can be assumed to have the following form: w PAGE 112 96 h(e)(x,y,z,t)=fN?)(x,y,z)d;(t), N/=V, hES (5-16) i = I where h PAGE 113 97 where: aNCe) aN(e) aNCe) aNCe) aNCe) l p -1 l l ax ax [K''>K~]= ff f ---k k 0 0 ax ay az xx rw aN(e) aNCe) 0 kylrw 0 l p --aN(e) aN(e) aN(e) 0 0 kzirw ay ay v PAGE 114 98 ad1(t) dl(t) at Ri(t) [KKrw] + [c] -[F] = (5-21) d/t) ad/t) R/t) at By setting the residuals equal to zero, the final system of ordinary differential equations is obtained: ad1(t) dl(t) at 0 [KKrw] + [c] [F] = 0 (5-22) d/t) ad/t) 0 at The system of Eq. 5-22 representing the variably saturated flow system was solved using finite differences. 5.4 Transient, Variably-Saturated Water-Table Recharge Example Test An example problem was formulated using laboratory data to verify the performance of the model. Details concerning the following laboratory data and the experiment were presented by Vauclin et al. (1979) The flow domain consisted of a rectangular soil slab 6.00 m by 2.00 m (Figure 5 3), with an initial horizontal water table located at a height of 0.65 m At the soil surface a constant flux of q = 3.55 m/day was applied over an infinite zone 1 m wide. The remaining soil surface was covered to prevent evaporation losses. Because of symmetry, it was necessary to model only one side (for this test the right side) of the flow domain. No-flow boundaries were set on the bottom and along the axis of symmetry PAGE 115 99 ( Figure 5-3). At the soil surface, the constant flux of 3.55 m/day was applied over the left 0 5 0 m of the top of the modeled domain. The remaining soil surface represented a no-flow boundary condition. A no-flow boundary was used above the water table on the right side of the domain (Figure 5-3). Infiltration Zone Soil Sufa A B C I I Impervious Substratum DE is head boundary AB is flux boundary BC, CD AF,FE are no flux boundary 1/ Figure 5 3. Schematic diagram of the flow domain From Vauclin et al. (1979) the characteristics of the experimental soil properties incl u ded a saturated hydraulic conductivity of 8.40 m/day, a porosity of T] = 8 s = 0 30, and a residual saturation of 8r= 0.01. When Eq. 5 2 was fitted to the water retention and the relative hydraulic conductivity data given by Vauclin et al (1979) estimated values for model parameters a and n were a = 3.3 m -1 and n = 4.1. Specific storage was neglected for this problem because changes in storage are facilitated by the filling of pores which overshadows the effects of compressibility for the vertical extent of this flow domain Hence, the specific storage coefficient was set to zero. The mesh spacing used in the test problem m the horizontal xand vertical z-directions was ~x = 0.10 m and ~ z = 0 05m respectively The soil system was assumed PAGE 116 100 to be initially at hydrostatic equilibrium with respect to the water table throughout the flow domain. Transient simulations were performed using time steps that varied from ~t = 180 secondes to ~t = 360 secondes. Transient positions of the water table are compared with the experimental results presented by Vauclin et al. (1979) in Figure 5-4 and Figure 5-5. These figures show excellent agreement between the transient water-table positions predicted by the variably saturated model and those obtained from Vauclin et al. (1979). Z ( m ) 0 .30 0 28 0 .24 0 .20 0 1 6 0 .12 0 .00 0 .50 1.00 1 .50 2 .00 2 .50 3 .00 X(2m) Moisture Content ... Observed point Figure 5-4. Water moisture distribution (x : y = 2 : 1) at 4 hours (Observed points are from Vauclin et al.1979) PAGE 117 Z(m) 101 0 .30 0 2 8 0.2 4 0 .20 0 1 6 0 1 2 0 .00 0 50 1 .00 1 .50 2 .00 2 .50 3 .00 X ( 2 m) Moisture Content .a. Observed point Figure 5-5. Water moisture distribution at 8 hours (x: y = 2 :1) (Observed points are from Vauclin et al.1979) 5.5 Application on SV5 Wetland 5.5.1 Finite-Element Mesh for Three-Dimensional Variably-Saturated Flow Model Unlike the saturated groundwater flow model (Chapter 3), matric pressure and water content may vary rapidly in time and space in the variably-saturated model, which produces rapid spatial and temporal variation in the variably-saturated flow hydraulic conductivity. In order to simulate the moisture or pressure variation accurately a fine numerical mesh must be used over zones of rapid change to determine flow parameter variations accurately Like the saturated flow model, the SV5 (Wise et al. 2000) simulation domain was divided into two main areas. One area was the aquifer; the other area was the wetland peat layer. Two parameters, ~ q u i f e r and eat> were used to represent their hydraulic conductivities. The numerical mesh was 26x21x3. The top 10 layers of mesh represent about top 1.55 m of the PAGE 118 102 aquifer. Among the toplO layers, the upper 4-5 layers were given a finer descretization, and the remaining layers of the top 10 were given a coarse descretization. The bottom 10 layers represented about 15 meters of saturated aquifer. The mesh in the saturated portion of the flow field was considerably coarser. z '~ 15 IO 5 0 0 Figure 5-6. Finite-element mesh of the wedge part for variably saturated flow model (Unit ism) Unlike the model for saturated flow, certain of the boundaries surrounding the simulation area are fixed. The top boundaries are the aquifer surface or the wetland bottom. The mesh deforms during a simulation because the elevation of the surface water changes The model seeks to locate the same elevation as the surface-water stage and at a location where a horizontal line defines where the surface-water surface intersects with the bottom of the wetland. As the mesh deforms during a simulation, a subroutine program is used to determine whether a cell is in the peat layer or the underlying aquifer and then the corresponding hydraulic conductivity and other parameters of the element are assigned. As the mesh deforms some elements are enlarged and some of them are compressed, but the PAGE 119 103 total area of the model domain does not change. 5.5.2 Modeling Results and Summary For the SV5 wetland simulation, the same hydraulic conductivity values used in the saturated model were used in the variably saturated flow model. Thus, Kxx=Kyy=8.5 m/day; =8.5 m/day; the anisotropy ratio K/Kv = 1.0 in aquifer; and = 0.2 m/day in the peat layer were used throughout the model domain. The specific storativity of S s =0 00001 (1/m) and the saturated water content and residual water content in the peat layer e s = 0 .7; er= 0.6 ; n=5.0; m=0.8; ex= 0.001 (11cm) = 0 1 (1/m) ; and in the aquifer e s = 0.34; e r = 0.18; n=9.l; m=0.89; a= 0.01(1/cm) = 1 (1/m) were also used, based on Walser et al. (1998) and Wise et al. (2000) Figure 57 shows the simulated piezometric head variation at selected depths beneath the wetland. When surface water was pumped from the wetland, the water stage decreased, and the associated piezometric head changes were greatest immediately beneath the wetland. The time needed for the wetland to recover from the stress was much longer. Figure 5-8 shows that the aquifer discharge to the wetland increased as the surface-water component of the wetland was drained and that the aquifer discharge decreased a little after the pumping stopped Figure 5-9 shows the transient variation of the inundated wetland area. During pumping, the inundated area decreased from 2,800 m2 to 700 m2 During the recovery period the inundated area increased s lowly Figure 5-10 shows that the model was ab le to fit the observed aquifer well which is the SFWMD well screened at a depth of six meters and grouted from the surface to a depth of five meters. In a comparison prediction, results from both the saturated model and the variably saturated flow model are consistent (Figure 5-11) The piezometric head distribution beneath the wetland is shown in Figure 5-12 and Figure PAGE 120 104 5-13. Tables 5-1 through 5-5 show the simulated heads, the corresponding areas, and the distances which were used to calculate the exchange discharge. 16 .05 16 \ --, :[ 15 .95 'O l PAGE 121 3000 fi' 2500 E ';; 2000 Q) < -0 1500 -0 C :, C 1000 500 \ \ ---\ \ \ -----0 20 105 .. ------v--i--40 60 80 Time (hour) Figure 5-9. Inundated area of wetland during simulation 4 2 \\'. i4.+--+---+---t--t E 4 15 -o-1, ~\\~------------+----t ---~,\l:5:::7"~""' ----" 4 1 --"""-~I/ 0 ;;: Q) 4 05 1--1-1---+---l--+--+---+--t---+-=+---1 E 0 N a: -----~ 4 +--H--t--+-~-+~ ~-~A--+--+--+---l 3.95 -I--H--------+,..,..,..,,..."-i ______________ -+---+--+--+--+------1 c---_/ 3 9 -+---+--+--+--+--+-->---+---+---+--< 0 20 40 60 80 100 Time (hour) WETLOBS WETLCAL G -WOBS G W CAL Figure 5 10. Simulated and observed heads in the wetland and in aquifer PAGE 122 106 4 .25 ~--------------4 2 ~ -+----+--+ - -1---t--+---I 4 1 s .. -+---+--+---+--+---+--+-~ ai rJ' ~~5~~l~~f=f=+~ I\ 4 1 r u ..... J .... .. -+-----I .: -, Qi 4 .05 !--+-+--+----+--!--;---,--t----t-:::::---1 E ---2 Q) a: -------4 -t---Hr--t----t----t---:::.l -------c...+------t----t-----i ---t;:;,-'i::--------+---3.95 -l---1+--b---=-----+--l---+---+---+--1----l ___..,..----3 9 +----,f----+--t--+----,--+--t--+----, 0 20 40 60 80 Time (h our) SM (WETL ) V M (WETL) S-M ( GW) V-S ( GW) Figure 5-11. Simulated heads from the saturated flow model and the variably saturated flow model z 6-. I 1 5 79 Ip. 89 11, 0 6 ( 3 ---1~.o s:: $ P9 1,;9B I). ( 5.' 7 5 tl .9\ -IC.0 K -I Y : 5 7 .50 0 50 1 5 1 0 0 Figure 5-12 Horizontal view of the three-dimensional finite element mesh and piezometric head distribution from the variably saturated flow model at a time of 9 hours PAGE 123 I .,J 107 ~ \ I\ "''-----V I IT. r, \ 1 ~ 1 lt / I\ I / \ V "l \ f'\ LI I \ I'\ V 7 Figure 5 -13. Enlarged top and middle part of the mesh and head contour at a time of 9 hours. I Ol Table 5-1. Data for calculating the discharge to the wetland from the aquifer at 10 hours Elements Area Head in Peat Head on Peat Distanc e of Two of Peat Layer (m2) L aye r (m) Surface ( m ) Nodes ( m ) l 6 82502 15.7300 15.72 03 0 04715 2 20.4751 15. 7301 1 5.7203 0 04718 3 34 1 25 1 15 7301 15. 7203 0 04721 4 47 7747 15.7291 15. 7203 0 0472 3 5 61.4246 15 7285 15.72 0 3 0 04726 6 75.0755 15. 7289 15. 7203 0 04729 7 88 7265 15. 7292 1 5.72 0 3 0 04732 8 1 0 2 378 15.7296 1 5.72 0 3 0 04736 9 116. 0 3 0 15.7301 1 5.72 0 3 0 04741 10 1 2 9.68 3 15. 7308 15. 7203 0 0474 7 Q to wetland at 10 hr 1 .6 954 m3/hr PAGE 124 108 Table 5-2. D ata for calculating the discharge to the wetland from the aquifer at 20 hours Elements Area Head in Peat Head on Peat Distance of Two of Peat Layer (m2 ) Layer ( m) Surface ( m) odes ( m) 1 8 39429 15.7513 15 7431 0 04715 2 25. 1 829 15.75 1 4 1 5.7431 0 047 1 8 3 41.97 1 3 15. 7515 15 7431 0 04721 4 58 .7594 15.7506 15. 7431 0 .04724 5 75.5484 15. 7501 15. 7431 0 04727 6 92.3 379 1 5.7505 15. 7431 0 04731 7 109.128 15. 7508 15. 7431 0 04735 8 125.919 15. 75 1 2 15. 7431 0 04740 9 1 42 711 15. 75 1 7 15. 7431 0 04746 10 1 59 .504 15.7524 15 7431 0 04753 Q to wetland at 20hr 1.8041 m1/hr Tab l e 5-3 D ata for calculati n g t h e d ischarge t o t h e wetland from th e a qu ifer at 4 0 hours Elements Area Head in P eat H ea d on P eat Di stance of Two of Peat Layer ( m 2 ) Layer ( m ) Surface ( m ) Nodes ( m ) 1 I 1.036 0 1 5.7876 15.7 8 1 3 0 04715 2 33 108 I 15.7876 1 5 7813 0 04719 3 55 1798 15.7870 1 5 7813 0 04722 4 77.2516 15 .7866 15. 7813 0 .04726 5 99 .3248 15. 7869 15 7813 0 04729 6 121.399 15.7871 15. 7813 0 .04734 7 1 43 .474 1 5 .7874 1 5 781 3 0 .04740 8 1 65 .55 0 1 5.7877 1 5.78 1 3 0.04747 9 187.6 28 15.7 882 15. 7813 0 .04754 10 209 707 15: 7889 15. 7813 0 04763 Q to w e tland at 40hr 1 8693 m1/hr PAGE 125 109 Table 5-4 Data for calculating the discharge to the wetland from the aquifer at 60 hours Elements Area Head in Peat Head on Peat Distance of Two of Peat Layer (m2) Layer ( m ) Surface ( m ) Nodes ( m) 1 13. 1576 1 5.8171 15.8120 0 04715 2 39 4727 1 5.8171 15. 8120 0 04719 3 65.7872 15.8166 15.8120 0 04723 4 92. 1052 15.8163 15. 8120 0 04723 5 118 .419 15. 8165 15.8120 0 04727 6 144.737 15.8167 15. 8120 0 04731 7 171.057 15. 8170 15. 8120 0 04737 8 197.378 15. 8173 15. 8120 0 04744 9 223.701 15.8177 15. 8120 0.04752 10 250.027 15. 8184 15. 8120 0 04772 Q to wetland at 60hr 1.8396 m3/hr Table 5-5. D ata for calc u lating t h e d ischarge to the wetland from the aq u ifer at 8 0 hours Elements Area Head in Peat Head on Peat Distance of Two of Peat Layer (m2 ) Layer (m) Surface ( m) Nodes ( m) 1 14 9321 15.8419 15. 8377 0 04715 2 44. 7963 15.8419 15. 8377 0 04720 3 74. 6597 15.8415 15.8377 0 04724 4 104 525 15. 8412 15. 8377 0 04728 5 134 391 15. 8415 15.8377 0 04733 6 164.259 15.8417 15.8377 0 04739 7 194 129 15.8419 15. 8377 0 04747 8 224 .001 15. 8422 15. 8377 0 04756 9 253 875 15. 8426 15. 8377 0.04767 10 283 753 15. 8432 15. 8377 0 047 7 9 Q to wetland at 80hr 1.7 653 m3/hr PAGE 126 110 The saturated model and the variably saturated model are different in several aspects. First of all, the variably saturated flow model requires that the moisture flux be specified on the soil surface as the boundary condition. This moisture flux on the soil surface induces the head distribution in the vadose zone that is not described in the saturated flow model. Second, the empirical formula used to characterize the matric pressure as a function of the degree of saturation may not be accurate, and as a result it displaces the phreatic surface (where pressure head is zero) to a location not consistent with observations. Table 5-6 shows the sensitivities of the heads related to charges in its hydraulic conductivities in the peat layer and aquifer and to changes i n the specific storativity. The hydraulic conductivities in both the peat layer and in the aquifer are important in the change of head distribution. Changes in the specific storativity do not affect the head very much relatively Table 5-6 Sensitivity calculated by variably saturated groundwater flow model Parameter (PR) Percentage varied dPR II h1-h2II (m) II h1-h2ll/dPR K (mlday) 0.1 -0.85 0 5700E-02 -0.6706E-02 K.(mlday) +0 1 +0.85 0.5090E-02 0.5988E 02 ~(m/day) -0.1 -0.02 0 .5 155E-02 -0 .25 78 ~(m/day) +0.1 +0.02 0.4665E-02 0.2333 s o -0.5 0 000005 0 .2 990E-04 -5.980 s o +0 .5 +0.000005 0 1998E 04 3 .996 PAGE 127 CHAPTER6 SUMMARY AND CONCLUSIONS This dissertation describes an analytical formula for groundwater variation affected by a flood wave in a wetland and also describes the development of two three-dimensional finjte-element method numerical models, one of which is a saturated groundwater flow model and the other a variably saturated groundwater flow model. Both models use a deformable mesh and real-time coupling for the wetland aquifer boundary along the interface between the surface water and groundwater. In the saturated groundwater flow model a predictor -corrector method is used to find the real-time boundary condition that satisfies both the he ad condition and flux condition at the same time The saturated groundwater flow model was verified using two analytical solutions; one is the Theis solution, and the other is the analytical solution found in Chapter2 of the dissertation, which describes groundwater flow around a circular wetland. The saturated groundwater flow model was further used to simulate the interaction problem of gro undw ater and the SV5 wetland in south Florida. The three-dimensional finite-element variably-saturated flow model was developed based on the linearized ruchard equation. This model is different from previous variablysaturated flow models in the deformable three-dimension firute-element mesh and dynarrucally interacting boundary conditio ns applied in this model. The predictor-corrector method is also used in the variably saturated grou ndw ater flow model to find the real time boundary conditions that satisfies both the h ead condition and flux condition at the same 111 PAGE 128 112 time. The variably-saturated groundwater flow model was verified based on the field and laboratory work described by Vauclin et aJ. (1979). The following specific conclusions can be drawn from this dissertation : 1. The analytical model describing the interaction between a flood wave in a lake and surrounding groundwater is an original development that can be used to simulate the groundwater fluctuations caused by a transient wetland surface water stage regardless of the function used to characterize the transient variation of the wetland stage. 2. Considering that the area of the water surface changes quickly because most wetlands have a shallow bank slope, the saturated groundwater flow model has the capability of successfully tracking the interaction boundary between groundwater and a wetland. The pumping from the SV5 wetland and the recovery of the wetland water surface afterwards were closely simulated by this model and thus the interacting boundaries were effectively traced by the deformable mesh. 3. The variably saturated groundwater flow model also successfully simulated the groundwater flow affected by the pumping from the SV5 wetland and the water level recovery in the wetland after pumping. The rapidly moving boundary i.e., the special phenomenon that occurs in the groundwater and wetland interaction problem is effectively dealt with by this model also using a deformable mesh PAGE 129 REFERENCES Allen, M B. And Murphy, C. L., 1986. A finite element collocation method for variably saturated flow in two space dimensions. Water Resour. Res. 22(11): 1537-1547. Ames W. F. 1977 Numerical methods for partial differential equations, 2nd ed. Academic Press. New York. Anderman Evan R. and Hill Mary C. 1999. A new multistage groundwater transport inverse method: Presentation evaluation and implications, Water Resour. Res. 35(4): 1053-1063. Anderson, E. I., 2004. An analytical solution representing groundwater-surface water interaction. Water Resour Res. 39(3), 1071 doi: 10.1029/2002WR001536 Anderson M P. and Cheng Xiangxue. 1993. Long-term and short-term transience in a groundwater/lake system in Wisconsin USA. Journal of Hydrology 145: 1-18. Asherson M P. 1979. Using models to simulate the movement of contaminants through groundwater flow system. Critical Reviews in Environmental Control, 9 (2):97-1 56. Barry, D A. Li L., and Jeng D.-S ., 2001. Tidal fluctuations in a leaky confined aquifer: Dynamic effects of an overlying phreatic aquifer. Water Resour. Res ., 37(4): 10951098 Barry D. A. Li L., Stagnitti F., and Parlange, J.-Y. 2000. Groundwater waves in a coastal aquifer: A new governing equation including vertical effects and capillarity. Water Resour. Res. 36(2): 411-420 Bates P D. Stewart M. D. Desitter A. Anderson, M. G ., Ronaud, J.-P. and Smith J. A., 2000. Numerical simulation floodplain hydrology WaterResour. Res. 36(9): 25172529. Bear, Jacob. 1979 Hydraulics of groundwater, McGraw -Hi ll Inc ., New York Boufadel Michel C ., and Peridier Vallorie 2 002. Exact a nalytical expressions for the piezometric profile and water exchange between stream and groundwater during and after a uniform rise of the stream level. Water Resour Res ., 38(7). 113 PAGE 130 114 Boufadel, Michel C. and Suidan Mak.ram 1998 2-D variably-saturated flow: physical scaling and Bayesian estimation. Journal of Hydrologic Engineering, Vol.3 No. 4 Boufadel, Michel C., Suidan Mak.ram T., and Venosa, Albert D 1997a Density-dependent flow in one-dimensional variably-saturated media. Journal of Hydrology, 202: 280301. Boufadel, Michel C., Suidan Makram T. and Venosa Albert D 1999 A numerical model for density-and-viscosity-dependent flows in two-dimensional variably saturated porous media Journal of Contaminant Hydrology 37:1-20 Boufadel Michel C., Suidan Mak.ram T. Venosa Albert D. and Bowers Mark T.1997b Steady seepage in trenches and dams : effect of capillary flow Journal of Hydrologic Engineering 125 (3):302-323 Bruce, G H., D W Peaceman and H H. Rachford Jr. 1953 Calculation of unsteady-state gas flow through porous media. Petrol. Trans. AIME 198. Pp74-92. Bruch, John C. 1975. Finite element solutions for unsteady and unsaturated flow in porous media California Water Resources Center Davis Brutsaert, W. F. 1971. A functional iteration technique for solving the Richards' equation applied to two-dimensional infiltration problems Water Resour. Res., 7(6): 15381596. Brutsaert W 1973 Numerical solution of multiphase well flow. Journal of the Hydraulic Division, American Society of Civil Engineers, 99(HY1): 1981 2001. Carr, P. A., and Van Der Kamp, G S. 1969. Determining aquifer characteristics by tidal method Water Resour Res., 5 1023-1031. Carrera J and Neuman S. P. 1986a. Estimation of aquifer parameters under transient and steady state conditions, 2 Uniqueness stability and solution algorithms Water Resour. Res 22(2) : 211 227. Carrers J. and Neuman S. P 1986b. Estimation of aquifer parameters under transient and steady state conditions, 3 Application to synthetic and field data Water Resour. Res. 22(2):228-242. Chang Shin and Yeh William W-G 1976 A proposed algorithm for the solution of the large scale inverse problem in groundwater Water Resour. Res 1 2 (3) : 368 371. Clement T. P ., Wise William R. and Molz, Fred J. 1994. A physical based, two dimensional finite -difference algorithm for modeling variably saturated flow Journal of Hydrology 161: 7190 PAGE 131 115 Cooley R. L. 1971. A finite difference method for unsteady flow in variably saturated porous media: application to a single pumping well. Water Resour. Res., 7: 1607-1625. Cooley, R. L., 1977. A method to estimating parameters and assuming reliability for models of steady state groundwater flow, 1, Theory and numerical properties, Water Resour Res., 13(2):1977 Cooley, R. L., 1979. A method of estimating parameters and assuming reliability for models of steady state groundwater flow 2, Application of statistical analysis, Water Resour. Res., 15(3):603-617 Cooley, R. L., 1982. Incorporation of prior information on parameters into nonlinear regression groundwater flow models, 1, Theory, WaterResour. Res., 18(4):965-976. Cooley R. L., 1983a Incorporation of prior information on paramters into nonlinear regression groundwater flow models, 2, Application, Water Resour. Res., 19(3):662676. Cooley, R. L. 1983b Some new procedures for numerical solution of variably saturated flow problems. Water Resour. Res., 19: 1271-1285. Cooper, H. H. And M. I. Rorabaugh. 1963. Groundwater movements and bank storage due to flood stages in surface streams. U.S. Geol. Surv. Water Supply Paper 1536-J. Cowardin, L. M ., Virginia, Golet, Francis C. and LaRoe, Edward T. 1979. Classification of wetlands and deepwater habitat of the United States. U.S. Fish and Wildlife Service Washington, DC. Cox, C. L., Jones, W. F., Quisenberry, V. L., Yo, F., 1994. One dimensional infiltration with moving finite elements and improved soil diffusivity, Water Resour. Res., 30: 14311438. Czarnecki, John B ., and Waddell, Richard K. 1984. Finite-Element Simulation of GroundWater Flow in the Vicinity of Yucca Mountain, Nevada-California, U.S. Department of Energy Washington, DC. Dagan G. 1985. Stochastic modeling of groundwater flow by unconditional and conditional probabilities: The inverse problem, Water Resour. Res., 21(1) : 65-72. Dean, Robert G. and Dalrymple Robert A. 1991. Water wave mechanics for engineers and scientists, World Scientific Singapore. Distefano N. and Rath A. 1975. An identification approach to subsurlace hydrological systems, Water Resour. Res., 11(6):1005 1012. PAGE 132 116 Doughty, Christine, 2000. Numerical model of water flow in a fractured basalt vadose zone: Box Canyon site, Idaho. Water Resour. Res. 36(12): 3521-3534. Erskine, A. D. 1991. The effect of tidal fluctuation on a coastal aquifer in the U.K. Ground Water, 29(4): 556-65. Finn, W. D. Liam, and Varoglu Erol. 1976. Variable domain finite element analysis of free surface flow problems. The first International Conference on Finite Elements in Water Resources, Princeton University NJ. Fleck, W. B., and McDonald, M. G. 1978 Three-dimensional finite-difference model of groundwter system underlying the Muskegon County wastewater disposal system, Michigan U.S. Geological Survey Journal of Research, 6(3) : 307-318 Fraisse Clyde W., and Campell, Kenneth L. 1997. FHANTM (field hydrologic and nutrient transport model), Version 2.0 User's Manual. University of Florida, Gainesville. France, P. W. 1976. Finite element analysis of two and three dimensional unconfined seepage problems, Proceedings of the First International Conference on Finite Elements in Water Resources, Princeton University, NJ. Freeze, R. A. 1971a Three dimensional transient, saturated-unsaturated flow m a groundwater basin. Water Resour. Res. 7: 347-366. Freeze, R. A 1971b. Influence of the unsaturated flow domain on seepage through earth dams. Water Resour. Res. 7: 929-941. Freeze, R. A., and Witherspoon, P A 1966, Theoretical analysis of regional groundwater flow; 1. Analytical and numerical solutions to the mathematical model. Water Resources Research 2: 641-656. Fretwell J. D., Williams J. S. and Redman P. J (Compliers). 1996. National water summary on wetland resources. U.S. Geological Survey Water Supply Paper 2425 431 p. Gambolati, G., Gatto, P. and Freeze, R. A. 1973 Mathematical simulation of the subsidence of Venice 2 results. Water Resources Research, 9: 721-733. Gleason P. J. 1989. Northern Martin County wellfield wetlands impact study aquifer performance test. James M. Montgomery Consulting Engineers Lake Worth, FL. Gonthier Gerard J ., Kleiss And Barbara A. 1996 Groundwater flow patterns and water budget of a bottomland forested wetland Black swamp, eastern Arkansas U.S. Geological Survey Water-Resources Investigations Report 95 -4192 PAGE 133 117 Gregg D 0 1966. An analysis of groundwater fluctuation caused by ocean tides in Glynn County Georgia, Ground Water, 4: 211-232. Groundwater Modeling System. 1996. Brigham Young University. Provo, UT. Gupta SumantK., Tanji, Kenneth K. andLuthin, James N. 1975. A three-Dimensional finite element ground water model California Water Resources Center, Davis Gureghian, A. B., Ward D.S. and Cleary R. W. 1979 Simultaneous transport of water and reacting solutes through multilayered soils under transient unsaturated flow conditions. Journal of Hydrology, 41: 253-278 Guymon G. L., Sco tt, V. H. and Hermann, L. R. 1970. A general numerical solution of the two-dimensional diffusion-convection equation by finite element method. Water Resources Research 6 ( 6): 1611-1617 Hoeksema, R. J. and Kitanidis, 1984 An application of the geostatistical approach to the inverse problem in two-dimensional groundwater modeling, Water Resour Res. 20(7): 1003-1020. Hughes, T. J. R. 1987. The Finite Element Method Prentice Hall Englewood Cliffs NJ. Hunt R. J 1996 Do created wetlands replace the wetlands that are destroyed? U.S. Geological Survey Fact Sheet FS 96 -2 46, 4 p Hunt, R. J., Bullen T. D. Krabbenhoft D. P. and Kendall C. 1998. Using isotropes of water and strontium to investigate the hydrology of a natural and a constructed wetland : Ground Water, 36(3): 434-443. Hunt, R. J., Krabbenhoft, D P ., and Anderson, M. P. 1996. Groundwate r inflow measurements in wetland systems: Water Resources Research 32(3): 495-507 Hunt R. J. Krabbenhoft D. P. and Anderson M.P., Assessing hydrogeochemical heterogeneity in natural and constructed wetlands Bio ge ochernist ry, 39: 2 712 93 Huyakorn, P. S., Guv a nasen V and Wadsworth, T. D. 1986. Three-Dimensional finite techniques for simulating unconfined flow with seepage faces. Proceeding of the 6th International Conference Lisbon Portugal. Huyakorn P. S. and Pinder, G F. 1983. Computational Methods in Subsurface Flow Academic Press New York Huyakorn, P. S. and Thomas, S. D. 1984. Techniqu es for making finite e l eme nt s competitive in modeling flow in v ariably saturated porous media Water Resou rces Research 20(8) PAGE 134 118 Istok, Jonathan 1989. Groundwater Modeling by the finite element method, American Geophysical Union, Washjngton, DC. Jain, S. K., Storm B. Bathurst, J.C., Refsgaard, J.C., and Singh R. D., 1992. Application of the SHE to catchments in India-Part 2: Field experiments and simulation studies on the Kolar Subcatchment of the Narmada River. Journal of Hydrology vol. 140 p. 25-47. Javandel, I., and Witherspoon, P. A. 1968 Applciation of the finite element method to transient flow in porous media. Society of Petroleum Engineers Journal 8(3) : 241252. Keidser, A., and Rosbjerg, D. 1991. A comparison of four inverse approaches to groundwater flow and transport parameter identification, Water Resour. Res. 27(9): 2219-2232. Kirkland M R. Hills R. G. and Wierenga P. J 1992. Algorithms for solving Richards equation for variably saturated soils. Water Resour. Res., 12(1): 57-64. Kirkner, D J., Theis, T L. and Jennings, A. A. 1984. Multicomponent solute transport with sorption and soluble complexation. Advances in Water Resources, 7:120-125. Knopman Debra S. and Voss Clifford I. 1989. Mul tiobjecti ve sampling design for paramter estimation and model discrimination in groundwater solute transport Water Resour. Res., 25(10):2245-2258. Konikow L. F. and Bredehoeft J. D. 1978. Computer model of two-dimensional solute transport and dispersion in groudnwater. U.S. Geological Investigation, Book 7. Krabbenhoft, D. P ., Anderson, Mary P. and Bowser Carl J. 1990a. Estimating groundwater exchange with lakes 2. Calibration of a threedimensional solute transport model to a stable isotope plume. Water Resources Research 26(10): 2455-2462. Krabbenhoft, D. P ., Bowser Carl J., Anderson, Mary P and Valley, John W 1990b Estimating groundwater exchange with lakes 1. The stable isotope mass balance method. Water Resources Research 2 6(10): 2445 2453. Krabbenhoft D P ., and Webster, A.W. 1995. Transient hydro ge olo gica l controls on the chemistry of a seepage lake: Water Resources Research, 31(9): 2295-2305 Lal A. M. Wasantha 2000. Numerical errors in groundwater and overland flow models Water Resour. Res. 36(5): 1237-1247 PAGE 135 119 Lamorey, Gregg and Jacobson Elizabeth 1998. Incorporation of constraints on hydraulic head gradients near no-flow boundary conditions in the determination of spatial drift and their use in an inverse groundwater flow model, Water Resour. Res. 34(11):2889-2910 Lal A. M. Wasantha 2000. Numerical errors in groundwater and overland flow models Water Resour. Res. 36(5): 1237-1247. Lapidus, L. and Pinder G. F. 1982. Numerical solution of partial differential equations in science and engineering John Wiley & Sons, New York. Lebedev N N 1965. Special functions and their applications Prentice-hall Inc. Englewood Cliffs, NJ. Lee, Terrie M., 2000. Effects of nearshore recharge on groundwater interactions with a lake in mantled karst terrain. Water Resour. Res. 36(8): 2167-2182. Loaiciga Hugo and Marino Miguel A 1987 The inverse problem for confined aquifer flow identification and estimation with extensions, 23(1):92-104. Mayer, A. S. and Miller, C T 1988 A Three-Dimensional Finite Element-Finite Difference Model for Simulating Confined and Unconfined Groundwater Flow, Proceedings of the VIl International Conference MIT, Cambridge MA. McDonald, M. G., and Harbaugh A. W. 1988 A modular three dimensional finite-difference groundwater flow model : US Geological Survey Techniques of Water-Resources Investigations of the U.S. Geological Survey, book 6 Chap Al, not sequentially paginated. Michalak Ann M., and Kitanidis, Peter K., 2004. Estimation of historical groundwater contaminant distribution using the adjoint state method applied to geostatistical inverse modeling. Water Resour. Res Vol.40, W08302. Middleton Beth. 1998. Wetland restoration, flood pulsing and disturbance dynamics. John Wiley & Sons Inc. New York. Mills J G. and Zwarich M A. 1986 Transient groundwater flow surrounding a recharge slough in a till plain : Canadian Journal of Soil Science 66(1) : 1 2 1 134 Mishra S. and Parker J C. 1989 Parameter estimation for coupled unsaturated flow and transport Water Resour. Res. 25(3) : 385 396. Mitchell A. R. and Griffiths D F. 1980. The finite difference method in partial d i fference equations. Wiley, New York PAGE 136 120 Mitsch, W. J. and Gosselink J. G. 1986 Wetlands. Van Nostrand Reinhold, New York. Mitsch, W. J and Gosselink, J G. 1993. Wetlands-2nd ed. Van Nostrand Reinhold, New York. Narasimhan T. N. and Witherspoon, P.A. 1976. An integrated finite difference method for analyzing fluid flow in porous media. Water Resour. Res., 12(1): 57-64. Narasimhan, T. N ., Witherspoon, P. A ., and Edwards, A L. 1978. Numerical model for saturated-unsaturated flow i n deformable porous media 2. The algorithm Water Resour. Res. 14:225-261. Neuman, S. P., 1973a Calibration of distributed parameter groundwater flow models viewd as a multiple-objective decision process under uncertainty, Water Resour. Res., 9( 4): 1006-1021. Neuman S P 1973b Saturated unsaturated seepage by finite elements Proceedings American Soci ety of Civil Engineers, 99(HY12): 2233-2250. Neuman, S. P., and Witherspoon P A. 1970 Finite-element method of analyzing steady seepage with free surface. Water Resources Research 6(3): 889-897 Nielsen, P. 1990. Tidal dynamics of the water table in Beaches Water Resour. Res. 26: 2127-2134. Norero, A. L. 1969. A formula to express evapotranspiration as a function o f soil moisture and evaporative demands of the atmosphere. PhD Diss. Utah State University, Logan. Novitzki, R. P. 1979 Hydrologic characteristics of Wisconsin s wetlands and their influence on floods, stream flows and sediment. In wetland functions and values: The state of our understanding P. E. Greeson, J R. Clark J. E. Clark eds. American Water Resources Association Minneapolis, pp 377-388 Orlob G. T., and Woods P C 1967 Water-quality management in irrigation systems Journal of the Irrigation and Drainage Di vision American society of Civil Engineers 93 : 49 66. Oster C. A., Sonnichsen J. C. and Jaske, P T. 1970. Numerical calculation o f the convective diffusion equation. Water Resources Research 6 : 1746 175 2. Peaceman, D. W and Rachford H. H. 196 2, Numerical calculation of multidimen sional miscible displacement. Society of Petroleum Engineers Journal 2: 327-339 PAGE 137 121 Philip, J. R. 1957 The theory of infiltration: 1. The infiltration equation and its solution Soil Science, 83: 345-357. Pickens J. F., Gillham, R. W. and Cameron, D R. 1979. Finite-element analysis of the transport of water and solutes in tile-drained soils Journal of Hydrology, 40: 243264 Pinder, G. F. 1973 A Galerkin finite element simulation of groundwater contamination in Long Island New York. Water Resources Research, 9(6): 1657-1669. Pinder, G. F., and Bredehoeft J. D 1968. Application of the digital computer for aquifer analysis. Water Resources Research, 8(1):108-120. Pinder G. F., andFrind, E. 0. 1972. Application of Galerkin's procedure to aquifer analysis. Water Resources Research 8(1):108-120 Pinder, George F. and Gray, William G. 1977. Finite element simulation in surface and subsurface hydrology, Academic Press, New York. Price, H. S. Cavendish J.C. and Varga, R. W 1968. Numerical methods of higher order accuracy for diffusion convection equations. Society of Petroleum Engineers Journal 1: 293-303 Pricket T. A. 1975. ModeLing techniques for groundwater evaluation. In: Advances in Hydro science, Vol. 10. Academic Press, New York pp. 1-143. Reeves M. andDuguid, J. 0 1975 Water movement through saturated-unsaturated porous media: A finite-element Galerkin model. ORNL-4927 Oak Ridge National Laboratory Oak Ridge TN. Remson I., Appel, C. A. and Webster, R. A. 1965 Groundwater models solved by digital computer. Journal of the Hydraulics Division American Society of Civil Engineers 91(HY3): 133-147 Rubin J. 1968 Theoretical analysis of two-dimensional transient flow of water in unsaturated and partly saturated soils Soil Sci. Soc. Am. Proc. 32: 607 615. Ruddy Barbara C. and Williams Robert S. Jr. 1991. Hydrologic relations between streamflow and subalpine wetlands in grand county, Colorado. Water-Resou rces Investigations Report 90-4129, U.S. Geological Survey, Washington DC. Rushton B. 1996 Hydrologic budget for a freshwater marsh in Florida Water Resources Bulletin 32(1), 13-21. PAGE 138 122 Russo, D., 2003. Simulation of transport in three-dimensional heterogeneous unsaturated soils with upscaled dispersion coefficients. Water Resour. Res. 39(12), 1358. Sacks, L.A., Herman, J. S., Konikow, L. F., and Vela, A. L. 1992 Seasonal dynamics of groundwater-lake interactions at Donana National Park, Spain J. Hydro!., 136: 123154 Sego!, G., Pinder G. F. and Gray, W. G. 1975. A Galerkin-finite element technique for calculating the transition position of the salt-water front. Water Resources Research, 11(2): 343-347 Senger, R. K., and Fogg, G E. 1987. Regional under pressuring in deep brine aquifers, Palo Dura Basin, Texas I, Effects of hydrostratigraphy and topography. Water Resources Research, 23(8): 1481-1493. Serfes, M E. 1992. Determining the mean hydraulic gradient of ground water affected by tidal fluctuations. Ground Water, 29:549-555. Serrano, Sergio E 2003. Modeling Groundwater flow under transient nonlinear free surface. Journal of Hydrologic Engineering, Vol. 8, No. 3. Shaffranek, Raymond W. 1996. Coupling models for cane! and wetland interactions in the south Florida ecosystem. U S. Dept. of the Interior, U.S. Geological Survey. Smith, A J., and Townley, L. R., 2002. The influence of regional setting on the interaction between shallow lakes and aquifers, 38(9), 1170 doi: 10.1029/2001 wr000781. Sonenshein, Roy S., and Hofstetter, Ro n ald H 1990. Hydrologic effects of well-field operations in a wetland Dade County, Florida. U.S Geological Survey, Water Resources Investigations Report 90 4143 Stone, H L., and Brain, P. T. L. 1963. Numerical solution of convective transport problems AlCHE J., 9: 5205-5213. Sun Hongbing. 1997 A two-dimensional analytical solution of groundwater response to tidal loading in an estuary, Water Resour. Res., 33(6): 1429-1435. Sun, Ne-Zheng, 1994 Inverse problems in groundwater modeling Kluwer Academic Publishers, Boston. Sun, Ne Zheng, Jeng Ming Chin, and Yeh William W. G. 1995. A proposed geological parameterization method for parameter identification in three-dimensional groundwater modeling, Water Resour. Res. 31(1):89-102. PAGE 139 123 Sun Ne-Zheng Yeh W W-G. 1985. Identification of parameter structure in groundwater inverse problem Water Resour. Res., 21(6):869-883. Swamee, Prabhata K. and Singh, Sushil K., 2003. Estimation of aquifer diffusivity from stream stage variation. Journal of Hydrologic Engineering, Vol. 8 No. 1. Switt, R. S., Annable, M. D. Wise W.R., and Walser J. A. E 1998. Hydrologic impacts of groundwater fluctuations on wetlands system Proceeding of the 1998 International Water Resources Engineering Conference and Mini-Symposia, Memphis, TN. Tanji, K. K., Dutt G R ., Paul J L. and Doneen L. D. 1867. A computer method for predicting salt concentrations in soils at variable moisture contents Hilgardia 38 : 307-318. Tannehill John C. Anderson Dale A. and Pletcher, Richard H. 1997 Computational fluid mechanics and heat transfer. Taylor & Francis, Washington DC. Therrien, R. and Sudicky E. A. 1996. Three-dimensional analysis of variably-saturated flow and solute transport in discretely-fractured porous media, Journal of Contaminant Hydrology 23: Throrenz Carsten Kosakowski Georg Kolditz Olaf, and Berkowitz Brian 2 002. An experimental and numerical investigation of saltwater movement in coupled saturated-partially saturated systems. Water Resour. Res 38(6) Townley, Lloyd R. and Trefry, Michael G ., 2000. Surface water-groundwater interaction near shallow circular lakes: Flow geometry in three dimensions. Water Resour Res 36(4): 935-949. Trescott, P C. and Larson, S. P. 1977. Comparison of iterative methods for solving two dimensional groundwater flow equations. Water Resources Research 13(1): 1 2 5-136. Trescott P. C. Pinder G. F. and Larson S. P 1976 Finite difference model for aquifer simulation in two dimensions with results of numerical experiments. U. S. Geological Survey Techniques of Water Resources Investigation Book 7 Automated Data Processing and Computations, 116pp. Tritscher, Peter and Read W. Wayne, and Broadbridge Philip, 2000. Specific yield for a two dimensional flow Water Resour. Res 36(6): 1393-1402. Truesdell B P 1995 The use of a wetland-groundwater interaction test to investi gate the hydrologic c onnection between isolated wetland and the surrounding groundwater. Master s Project University of Florida, Gainesvill e. PAGE 140 124 U.S. Geological Survey. 1996. Ground-water flow pattern and water buget of a bottomland forested wetland Black Swamp, eastern Arkansas. Van Genuchten, M. Th., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soil. Soil Sci Soc Am. J., 44: 892-898 Van Genuchten M.Th., Pinder, G. F. and Frind, E. 0. 1977. Simulation of two-dimesional contaminant transport with isoparametric-Herrnitain finite elements. Water Resources Research 13(2): 451 458. Vauclin, M., Khanhi D and Vachaud, G 1979. Experimental an d numerical study of a transient, two-dimensional unsaturated-saturated water table recharge problem. Water Resour Res ., 15(5): 1089-1101. Volker Ray E. and Zhang Qi 2001. Comment on "Analytical solut ion o f groundwater response to tidal fluctuation in a leaky confined aquifer by Jiu Jimmy Jiao and Zhonghua Tang. Water Resour. Res 37(1): 185-186. Voss Clifford I. 1984 A finite-element simulation model for saturated-unsaturated fluid density-dependent groundwater flow with energy transport or chemically reactive single-species solute transort U.S. Geological Survey and U S. Air Force Engineering and Services Center Tyndall A. F. B ., FL. Walser J.E., Annable, M. D., Wise W.R., and Switt R S 1998. Investigating Isolated Wetlands-Aquifer Interaction Through Modeling Proceeding of the 1998 International Water Resources Engineering Conference and Mini-Symposia, Memphis, TN. Wierenga, P. J. 1977. Solute distribution profiles computed with steady state and transient water movement models. Soil Science Society of America Journal 41: 1050-1055 Winter T C. 1978 Numerical simulation of steady -s tate three-dimensional groundwater flow near lakes. Water Resour. Res ., 14: 245-254. Winter T C. 1983 The interaction of l akes with variably saturated porous media. Water Resour. Res., 19 : 1203 1218. Winter, T C. 1986. Effect of groundwater recharge on configuration of the water table beneath sand dunes and on seepage in lakes in the Sandhills of Nebraska, U.S A. J. Hydrology 1986. Winter T C. and Woo, M. K. Eds. 1990. Hydrology of Lakes and Wetlands. In Surface water hydrology Geological Society of American, Bounder, CO. PAGE 141 125 Wise, W R., Annable, M. D. Walser, J. A. E. Switt, R. S. and Shaw, D. T. A Wetland aquifer interaction test, 2000 Journal of Hydrology, 227 (2000) 257-272. Yeh G. T., 1992. Users' Manual: A finite element model of water flow through saturated unsaturated porous media Department of Ci vii Engineering The Pennsylvania State University Ye, S. Xue, Y., and Xie, C., 2004. Application of the mutltiscaJe finite element method to flow in heterogeneous porous media. Water Resour Res. Vol. 40, W09202. Yeh, William W-G 1986 Review of parameter identification procedure in groundwater hydrology: The inverse problem, Water Resour. Res. 22(2):95 108. Zhang Jinqi and Yeh T.-C. Jim 1997. An iterative geostatisticaJ inverse method for steady flow in the vadose zone, Water Resour. Res., 33(1) : 63-71. Zhang Xiaoxian, and Ewen, John 2000. Efficient method for simulating gravity-dominated water flow in unsaturated soils. 36(9): 2777-2780. Zienkiewicz 0. C. 1971. The finite element method in engineering science McGraw-Hill, London Zienkiewicz, 0 C ., Meyer, P. and Cheung, Y. K. 1966. Solution of anisotropic seepage problems by finite elements. Proceedings American Society of Civil Engineers, 92 EMI: 111-120 Zienkiewicz, 0. C. and Parekh C J. 1970 Transient field problems : Two-dimensional and three dimensional analysis by isoparametric finite elements. International Journal of Numerical Methods in Engineering 2: 61 71. PAGE 142 BIOGRAPI-IlCAL SKETCH Jinhua (Jay ) Zhong was born on August 1 1966 in a small village of Jiangsu Province China He received his Bachelor of Science degree in civil eng i neering with special emphasis in hydrology and water resources engineering; in 1988 ; and the Master o f Science degree in 1991 from the Department of Hydrology and Water Resources Engineering at Hohai University Nanjing P R. China. Soon after earning his master s degree he began to work as an engineer in the Institute of Environmental Hydraulics at Hohai University In 1995, he came to the United States and entered the Graduate School of the Un iv ersity of Florida. In 2001 he obtained his MS degree from Electrical and Computer Engineerin g Department of the Univ ersity of Florida. He earned the Doctor of Philosophy degree from the Civil and Coastal Engineering Department of the University of Florida in 2004 126 PAGE 143 I certify that I have read this study and tha t in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate; in scope and quality as a dissertation for the degree of Doctor of Philosophy. Louis H. Motz, Chairman Associate Professor of Civil and Coastal Engineering I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate; in scope and quality as a dissertation for the degree of Doctor of Philosophy. Kirk Hatfield, Cochair Associate Professor of Civil and Coastal Engineering I certify that I have read this st udy and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate; in scope and quality as a dissertation for the degree of Doctor of Philosophy. I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate; in scope and quality as a dissertation for the degree of Doctor of Philosophy. PAGE 144 I certify that I have read this study and that in my opinion it conforms to acceptable standards of scholarly presentation and is fully adequate; in scope and quality as a dissertation for the degree of Doctor of Philosophy. 11P~ Michael D. Annable Professor of Environmental Engineering Sciences This dissertation was submitted to the Graduate Faculty of the College of Engineering and to the Graduate School and was accepted as partial fulfillment of the requirements for the degree of Doctor of Philosophy. December, 2004 Pramod Khargonekar Dean, College of Engineering Kenneth J. Gerhardt Interim Dean, Graduate School PAGE 145 , L 17 20 o'j J 3 ,z_lo UNIVERSITY OF F LOR I DA II Ill lllllllllll l lllll llllll lllllll l lll lll l l lll llllllll lllllll I I 3 1262 08554 5 712 xml version 1.0 encoding UTF-8 REPORT xmlns http:www.fcla.edudlsmddaitss xmlns:xsi http:www.w3.org2001XMLSchema-instance xsi:schemaLocation http:www.fcla.edudlsmddaitssdaitssReport.xsd INGEST IEID EL3NROQEB_SOTAHS INGEST_TIME 2012-03-13T14:42:47Z PACKAGE AA00008997_00001 AGREEMENT_INFO ACCOUNT UF PROJECT UFDC FILES |