UFDC Home  Search all Groups  UF Institutional Repository  UF Institutional Repository  UF Theses & Dissertations  Vendor Digitized Files   Help 
Material Information
Subjects
Notes
Record Information

Full Text 
TWODIMENSIONAL ANALYTICAL AND THREEDIMENSIONAL FINITEELEMENT 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 Nonperiodic 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 THREEDIMENSIONAL DEFORMABLE FINITEELEMENT SATURATED GROUNDWATER MODEL .......................... 19 3.1 Introduction .................................................. 19 3.2 Groundwater Flow Models Based on FiniteDifference and FiniteElement Approximation .................... .. ..... 21 3.3 ThreeDimensional Deformable FiniteElement 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 FiniteElement Model to Describe the Free Surface .............................. 36 3.4 Validating the Code of the ThreeDimensional Deformable FiniteElement 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 ThreeDimensional Deformable Saturated Groundwater Flow Model to SV5 Wetland ...................... 50 3.6.1 Site DescriptionSV5 ................................... 50 3.6.2 Field Methods and Wetland Aquifer Interaction Test ........... 52 3.6.3 StageVolume Relationships .............................. 54 3.6.4 Creation of FiniteElement 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 THREEDIMENSIONAL FINITEELEMENT METHOD SATURATED GROUNDWATER MODEL IN SEARCHING FOR THE HYDRAULIC CONDUCTIVITY OF THE PEAT LAYER ....... 69 4.1 Introduction .................................................. 69 4.2 Adjoint Problem for ThreeDimensional FiniteElement 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 THREEDIMENSIONAL FINITEELEMENT VARIABLYSATURATED GROUNDWATER MODEL ....................................... 87 5.1 Introduction ................................ .................. 87 5.2 ThreeDimensional FiniteElement VariablySaturated Groundwater Model ......................................... 88 5.2.1 Governing Equation .................................... 88 5.2.2 Linearizing the Governing Equations ....................... 89 5.2.3 Boundary Conditions for the VariablySaturated Model ........ 91 5.3 FiniteElement Method for Modeling ThreeDimensional VariablySaturated Flow ..................................... 94 5.4 Transient, VariablySaturated WaterTable Recharge Example Test ...... 98 5.5 Application on SV5 W etland .................................... 101 5.5.1 FiniteElement Mesh for ThreeDimensional VariablySaturated 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 21 Parameters and their values ............................................ 17 31 Selected early references for the finitedifference method .................... 22 32 Selected early references for the finiteelement application to model groundwater 25 33 Parameters and their values .................................. .......... 42 34 Data for calculating the discharge to the wetland from the aquifer at 10 hours .... 63 35 Data for calculating the discharge to the wetland from the aquifer at 20 hours .... 64 36 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 38 Data for calculating the discharge to the wetland from the aquifer at 50 hours .... 65 39 Sensitivities calculated by saturated groundwater flow model ................. 66 41 Error and sensitivities based on the adjoint method model .................... 84 51 Data for calculating the discharge to the wetland from the aquifer at 10 hours ... 107 52 Data for calculating the discharge to the wetland from the aquifer at 20 hours ... 108 53 Data for calculating the discharge to the wetland from the aquifer at 30 hours ... 108 54 Data for calculating the discharge to the wetland from the aquifer at 40 hours ... 109 55 Data for calculating the discharge to the wetland from the aquifer at 50 hours ... 109 v 56 Sensitivities calculated by variably saturated groundwater flow model ......... 110 LIST OF FIGURES Figure page 21 Setup of the groundwater problem around a circular wetland ................... 9 22 Transient groundwater level changes at radial distances 1,2,3 and 4 meters away from the wetland shoreline ..................................... 15 23 Groundwater distribution over the distance to the wetland shoreline at 6,7,8,9,10 h .................................................... 15 24 Threedimensional view of the groundwater table around a circular wetland at t=3hr .......................................... 16 25 Comparison of the analytical solution and numerical solution at 2 and 5 meters from the wetland shore line ........................... 18 31 Conceptual profile model of a wetland/groundwater system .................. 26 32 Comparison of finite element mesh and finite difference mesh ................ 26 33 Setup of the saturated groundwater problem around a wetland ................. 27 34 Flowchart illustrating the iteration procedure of finding the coupling boundary condition ..................................... 29 35 Unsteady flow to a well in phreatic aquifer at time t ......................... 38 36 FiniteElement Method mesh of the wedge used to simulate of a well pumping in phreatic aquifer ................................. 40 37 Comparison of the numerical and analytical head for the example problem of a transient pumping test in a phreatic aquifer ......................... 40 38 Comparison of the numerical and analytical simulated groundwater table ........ 41 vii 39 Threedimensional finiteelement mesh for the simulation of wetland and aquifer 43 310 Relationship between the area of wetland water surface ..................... 44 311 Relationship between the deepest water depth and the volume of wetland ....... 45 312 Illustration of the initial condition where it was assumed a wetland surfacewater stage and the phreatic surface elevation were the same ........ 45 313 Plan view of the numerical mesh at time zero. ............................ 46 314 ThreeDimensional FiniteElement mesh and simulated head contour at hour 5 .. 46 315 Simulated head contour at hour 5 ...................................... 47 316 Simulated isosurface plot of 10.6 meter head at hour 5 ..................... 47 317 Threedimensional FiniteElement mesh and head contour on a cross section on a cross section through wetland ..................... 48 318 Simulated head contour on the vertical cross section through wetland .......... 48 319 Modeled isosurface head plot of 11.1 m,l1 m and 10.5 m ................... 49 320 Mesh and modeled head contour with vertical cross section through the pumping well .......................................... 49 321 Head contour (m) on the vertical cross section through the pumping well ....... 50 322 W ell location map (Switt et al. 1998) ................................... 51 323 Top surface of peat layer ............................................. 54 324 Bottom surface of peat layer .......................................... 55 325 Relationship between wetland volume and water stage ..................... 55 326 Finiteelement mesh of wedge of wetland and aquifer for the saturated flow model ........................................ 57 327 Vertical crosssection view of threedimensional finiteelement mesh and head from the saturated flow model ............................... 61 328 Piezometric head at positions beneath wetland with time .................... 61 329 Variation in the inundated area of wetland during simulation ................ 62 330 Discharge between wetland and aquifer with respect to time ................. 62 331 Comparison between calculated head and observed head .................... 63 332 Comparison results using Walser's parameter values ....................... 67 333 3D numerical model simulated results over specific storage ................. 67 41 Sketch showing setup of the saturated groundwater inverse problem around a wetland .................................... 74 42 Relation between error and peat layer hydraulic conductivity .................. 84 43 Flowchart illustrating the iteration procedure of inverse problem .............. 85 44 Vertical crosssection of the adjoint state distribution at T = 10 hours ........... 86 51 Sketch showing setup of the variably saturated groundwater flow problem around a wetland ...................................... 92 52 Schematic of relative evapotranspiration (ET/PET), as affected by soil water potential,j .................................. 93 53 Schematic diagram of the flow domain ................................... 99 54 Water moisture distribution (x:y=2:1) at 4 hours .......................... 100 55 Water moisture distribution at 8 hours (x:y=2:1) .......................... 101 56 Finiteelement mesh of the wedge part for variably saturated flow model (Unite is m) ............................ 102 57 Calculated piezometric head variation beneath wetland using the threedimensional variably saturated groundwater flow model .......... 104 58 Model predicted discharge from the aquifer to the overlaying wetland ......... 104 59 Inundated area of wetland during simulation ............................. 105 510 Simulated and observed heads in wetland and in aquifer ................... 105 511 Simulated head from the saturated flow model and the variably saturated flow model ............................................. 106 512 Horizontal view of the threedimensional finiteelement mesh and piezometric head distribution from the variably saturated flow model at time of 9 hours .. 106 512 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 = groundwater 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 xdirection [L] Az = the mesh space step in ydirection [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 TWODIMENSIONAL ANALYTICAL AND THREEDIMENSIONAL 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 surfacewater 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 finitedifference 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 wetlandaquifer interactions. The first model is a saturated groundwater flow model that incorporates adaptive simulation technologies to permit realtime 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 realtime. 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 lifesupport 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,ETSoG,T (12) 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 crosssectional area through which the flow occurs. Periodic or nonperiodic fluctuations of the wetland stage in a flowthrough 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 onedimensional 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 twodimensional analytical solution was obtained by considering a shift in phase and a change in amplitude along a straight coastline (Sun 1997). Threedimensional 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 flowthrough 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 onedimensional study by Cooper and Rorabaugh (1963), most investigations addressed surfacewater 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 surfacewater flood waves. This goal was achieved by meeting the single objective of developing an analytical model using Fourier transforms. A limited number of threedimensional 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 threedimensional finite difference groundwater model and studied the dynamic coupling between wetlands and the subsurface. Lee (2000) used the HYDRUS2D flow model and simulated vertical twodimensional groundwater flow through Lake Barco in Florida. Neither FEMWATER (1996) nor MODFLOW can simulate or incorporate a moving realtime wetland boundary caused by a change in surface water stage. Furthermore, as a saturated threedimensional finitedifference 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) threedimensional flow, 3) the movement of the phreatic surface, and 4) the coupling that occurs between the transient surfacewater 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 threedimensional saturated groundwater flow model that employs deformable hexahedral elements. The finiteelement 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 threedimensional 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 shortterm pumping from a phreatic aquifer affects wetlands and lake water stages. CHAPTER 2 ANALYTICAL GROUNDWATER MODEL RESPONSE TO THE PERIODIC OR NONPERIODIC 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 onedimensional 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 onedimensional relation between tidal fluctuations and groundwater table perturbations. Nielsen (1990), Serfes (1992) and Carr (1969) also published papers investigating the similar problems. A twodimensional analytical solution was obtained by considering a phase shift in watertable perturbations and the associated change of amplitude along a straight coastline (Sun 1997). An analytical solution representing vertical twodimensional groundwatersurface water interaction was obtained by Anderson (2003). Threedimensional 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 surfacewater 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 onedimensional 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 twodimensional groundwater flow in a homogeneous isotropic aquifer can be expressed in radial coordinates as: 2h 1 Dh S Oh 2 h (21) 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 21) in which K is the hydraulic conductivity and ho is 9 the average unconfined aquifer thickness equal to the regional average watertable elevation above an underlying impermeable confining unit. The boundary condition along the shore line of the wetland is Eq. 22, where r o is the radius of the circular wetland. It is also assumed in Eq. 23 that as ro, the groundwater table elevation approaches the regional average water table elevation, hoo. B.C: h(r,t) ,=r, = hoo + f(r ,t) @ r= r. (22) h(r,t) ,r= = hoo @ r= o (23) wetland or lake _ h_ //777/7 //7777//77///77/7 77777/ Figure 21. 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 24 is obtained: B.C: h(r,t) r=ro = hoo + h ei't @ r= r, (24) h(r,t) , = hoo @ r= co (25) where ho is the half height of the periodic constituent and to is the angular frequency. No initial condition is necessary, since here the steadystate 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 (26) (26) Substitution of Eq. 26 into Eq. 21, Eq. 22, and Eq. 23 yields d2g +l.O g S.ig=0 (27) Br2 r ar T where the boundary conditions are: g(r)=] @ r=ro (28) g(r)=O @ r= (29) Equation 27 can now be rewritten to obtain: v2 +vg +v2g= av2 Ov (210) v 2= T (211) The solution to Equation 210 subject to the boundary conditions in Equation 28 and 29 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() (212) H' (vC) and when wo is less than 0: g(r) = C iY) = CH (u)= (213) 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. 212 and Eq. 213 into Eq. 26 gives the following solution for the groundwater flow around a circular lake affected by periodic function: h,(rt)=hoo+h0o )e H.(.vo) (214) for 0 > 0 and: h.(r,t)=h0 +h0 e Nt( S (o2)(o) (215) for c < 0. The asymptotic approximations of the functions Ho')(z) and Ho(2)(z) are: H (z)= 2 i*) [ 132**(2k 1)2(2 iz)k+O(z 1) (216) S7 z k= 22kk! I i 2 /) 1/2 (i e( )~ 132..(2k 1)2 2)(z)= 2 / 4 ( iz)k+O(zi n1) VZ k=o 22kM! uo= '*(1O'i)ro V0 = *(i1)ro (217) u= .*(1i)r 2v= T F2T 2.4 Analytical Solution for the NonPeriodic 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" (218) 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 (219) H2)(u) h,(r,t)=hoo + F(o) "*eu ,, k) ( o < 0 ) (220) as given by Eq. 214 and Eq. 215. The inverse Fourier transformation of Equation 219 and Eq. 220 leads to: h(r,t)=h oo+ + 2n =h0 2x h' H'(v) e e Hf ~(v0) J H (u0) (221) 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 22 is then h(r,t)=h o+ In 2x (222) (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 V2n G (223) The analytical solution obtained when Eq. 223 is combined with Eq. 221 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) (224) The transient changes of groundwater table described by Eq. 224 at four different distances from the wetland are shown in Figure 22. The groundwater distributions along the radius at four different times are shown in Figure 23. From both figures, it is evident that groundwater fluctuations are greatest near the wetland. From Figure 23, 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 22. 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 23. Groundwater distribution over the distance to the wetland/lake shoreline at 6,7,8,9, and 10 h. 16 Figure 24 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 21. 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 24. Threedimensional 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. 21Eq. 23. An implicit finitedifference scheme was used to represent the radial form of the groundwater flow equation Eq. 21 as: h, 2h" i+h,' I h11 hiI _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 21. By considering the different boundary conditions (Eq. 23 and Eq. 24), the corresponding numerical results are found [Figure 25]. Analytical and numerical results illustrate time variations in groundwater levels. Figure 25 shows that the analytical solution and numerical solution compare very well. Table 21. 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 25. 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 systemwide water balance. CHAPTER 3 TRANSIENT THREEDIMENSIONAL DEFORMABLE FINITEELEMENT SATURATED GROUNDWATER MODEL 3.1 Introduction Recent research has focused on the details of various wetland functions such as waterquality treatment, floodwave 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 sideslope flow (including overland and smallchannel 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, groundwater flow paths and rates, and soilprofile development. The location of recharge and discharge areas can change in response to seasonal or weatherrelated factors and can affect the groundwater flow direction. As a result, groundwater and surface water quality and the presence of various wetlandplant 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 threedimensional deformable finiteelement 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 threedimensional 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 GroundWater Flow Models Based on FiniteDifference and FiniteElement Approximation Several types of numerical methods have been used to solve groundwater flow problems, the two principal ones being the finitedifference method and the finiteelement method. The finitedifference method was initially applied to the flow of fluids in petroleum reservoirs, and it wasn't until the mid1960'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 onedimensional, steadystate groundwater flow in an isotropic and homogeneous aquifer, mathematical formulation and computer implementation are easily understood; and (2) the accuracy of solutions to steadystate and transient groundwater flow problems is generally quite good. The finitedifference 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 finiteelement method. Table 3.1 lists selected early references where finitedifference methods were applied to model groundwater flow. More recent applications of hydrologic modeling using finitedifference 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 31. Selected early references for the finitedifference 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 finitedifference meshes. For lakes, a finitedifference mesh may be suitable; however, because the bed of most wetlands exhibits significant slope, a poorly designed fixed finitedifference 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 wetlandflow 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 studysite 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 finiteelement 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 finiteelement method include: (1) for simple problems, the finite element method requires a greater amount of mathematical and computer programming sophistication than does the finitedifference method; and (2) there are fewer well documented computer programs. 24 The finiteelement method provides approximations of much higher order than does finitedifference 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 higherorder 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 finiteelement applications to model groundwater flow are listed in Table 32. 3.3 ThreeDimensional Deformable FiniteElement Saturated Groundwater Model To simulate threedimensional 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 31). This transient boundary suggests that a deformable finiteelement mesh may be needed in the model to characterize adequately the moving groundwater/surfacewater interface. 25 Table 32. Selected references of early finiteelement 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 surfacewater and groundwater systems, where the groundwater/surfacewater interface moves, must use welldefined 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 31. Conceptual profile model of a wetland/groundwater system Finally, typical finitedifference modeling is not able to trace the curvature of the bottom of a wetland (see Figure 32). Finite elements, however, are easily configured to fit nonlinear 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/surfacewater interface, and tracing the curvature of a wetland bottom, a special threedimensional finiteelement model is needed. FiniteElement Method FiniteDifference Method Figure 32. Comparison of finiteelement mesh and finitedifference mesh 27 3.3.1 Governing Equation and Boundary Conditions Eq. 31 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 (31) 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 33, the boundary of volume 0 is characterized through a combination of specified flux boundaries r, and specified head boundaries P2. Pertinent specific noflux 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 noflow boundaries F EG, FH, ACDB, and CD are specified head boundaries rF Figure 33. 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 34 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. 31, two spaces, S and V, are defined (Hughes 1987): S= {uIuEH1, u(xyz)2=qr2 } (32) V= uuEHuH, u(x,y,z)r2=O} (33) where H' is Sobolev space comprised of functions which are squareintegrable 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 34. 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. 31 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. 34 expresses the governing groundwater flow equation integrated over the modeling domain. Manipulation of the first term on the left hand side of Eq. 34 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 (35) (35) 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. 34 leads to: f w[(KVh)n]dP+ f w[(KVh)n]ddf (Vw)(KVh)dQ j wQdQ = wS dQ a 11 adt (36) 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. 36 equals zero. If it is further assumed that Q equals zero, then Eq. 36 reduces to: w[(KVh)n]dr (Vw)(KVh)dQ= wS, dQ (37) 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 (38) 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. 38 and Eq. 39) are substituted into Eq. 37, 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 (310) [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 (311) aN.(e) ax aN (e) ay aN(e) az (312) and: N (e) C (e). e) ... e)] (313) V(e) NP where V(e) is the volume of element e. Eq. 310 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] (314) at By setting the residuals equal to zero, we have: ad,(t) [K] + + [C] I =[F] (315) at Eq. 315 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 finiteelement mesh. 34 Although several methods are available for solving this system of equations, it has become standard practice to use the finitedifference method: ah d(t+At)d(t) at At (316) If a variable, co is defined, such that: At At Then: (318) which can be extended to the vector of [d]: (317) [d(e)] = (1 )[d(t)]+ [d(+At)] (319) If Eq. 319 and Eq. 318 are substituted into Eq. 315, the finitedifference 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((1o)[F(t)]+([F(t+At)] (320) d(e) = (1o)d(t) +a d(t+AOt) 35 Eq. 320 represents a system of linear equations with unknown variables of d,(t), and where N,(x,y,z) are the basis functions in a trilinear 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. 320 begins by specifying the initial values of [d] (i.e., the values of head at time to= 0). Next, the system Eq. 320 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 finitedifference formulations are defined. When o=0, it is the forwarddifference method is obtained: ( [C] )[d(t+At)] = ( [C] At[K])[d(t)] +At[F(t)] (321) When w=1/2, the centraldifference, or CrankNicholson, method is obtained: ([C] + 1 At[K] [d(t+At)] = [C] 1 1 At [K] [d(t)] + (F(t)]+[F(t)]) (322) When (=1, it is the backwarddifference method is obatined: ( [C] + At[K] )[d(t+At)] = [C] [d(t)] +At[F(t)] (323) 36 3.3.3 Adaptation of the Saturated FiniteElement Model to Describe the Free Surface Several papers have been written on the groundwater freesurface problem (e.g.,Taylor 1967, Neuman 1970, Finn 1976, France 1976, and Huyakom 1986). Most of these papers discussed the steadystate freesurface 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 + (324) 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 (325) For the steadystate 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. 324 at the free surface, there is: ah h h  Ki + Kj + K k = q. (326) 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 (327) 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 NewtonRaphson iterative method. Thus, a new finiteelement mesh is generated by modifying the coordinates of the nodes in the new flow domain. 38 3.4 Validating the Code of ThreeDimensional Deformable FiniteElement Groundwater Flow Model 3.4.1 Theis Solution To validate the code of the threedimensional deformable finiteelement 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 35 shows a schematic of unsteady flow to a well in a phreatic aquifer. Ho 7/ ////777 //// Figure 35. 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 32, 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 (328) 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 (329) 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 threedimensional finiteelement model, the specific storativity, Ss, is 0.0001 m1; 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 [LT2M1]. Because of problem symetry, modeling was limited to simulating hydraulic head variations over a wedge section of the aquifer using a 21x11x3 mesh (Figure 36). 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 36. Finiteelement 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 overrelaxation (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 threedimensional finiteelement 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. 330: 1 2o2 (330) 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 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 38 are the threedimensional finiteelement model simulated results. The triangle points and the diamond points are the groundwater analytical solutions from Eq. 224 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 33. 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 finiteelement 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 39 illustrates the numerical mesh corresponding to the 100.9 100.8  3D num d2m 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 38. Comparison of the numerically and analytically simulated groundwater table. Table 33. 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 realtime 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 Upperextent 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+01 7.99 999E01 "E+0I aquiclude 9 9999E+01 9.9999E+01 Figure 39. Threedimensional finiteelement 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 39). 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 310 and Figure 311 were used to define the relationship between the wetland water stage, the surface water area of wetland, and the corresponding wetland volume. Figure 312 illustrates the initial condition where it was assumed that the wetland surfacewater stage and phreatic surface elevation were the same. Figure 313 illustrates a 44 plan view of the numerical mesh at time zero. Figure 314 through Figure 321 represent various depictions of head contours and isosurface 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 317 and Figure 318), 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 320 and Figure 321). 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 310. 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 (sm) Figure 311. Relationship between the deepest water depth and the volume of wetland. N Figure 312. 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.2005E04 Figure 313. Plan view of the numerical mesh at time zero. x , Figure 314. Threedimensional finiteelement mesh and simulated head contours at hour 5. 2 4999EOI I 2005E04 49E 9999E 01 7 4999E,01 5 9999E+01 7 9999E+01 9 9999E+01 9 9999E+01 Figure 315. simulated head contour at hour 5. zi Figure 316. Simulated isosurface plot of 10.6 meter head at hour 5. 48 52005E04 2005E04 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 317. Threedimensional finiteelement mesh and head contour on a cross section through wetland 0.046 04 .5 200SE.04 95+2005E.04 Figure 318. Simulated head contour on the vertical cross section through the wetland x Y 5 2005E04 5 2005E04 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 319. Modeled isosurface plot of head equal to 11.1 m, 11 m and 10.5 m Figure 320. 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 5hour 1 901< E.01 59999 oo 9 9S99 99 9EI 9 (Si + Figure 321. 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 5hour 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 finitedifference grids, the finiteelement 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 finiteelement method and finitedifference method could be an effective choice for a large scale practical problem. 3.6 Application of ThreeDimensional Deformable Saturated Groundwater Flow Model to SV5 Wetland 3.6.1 Site DescriptionSV5 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 watersupply wells located from a few hundred meters to four kilometers from SV5 (Figure 322). 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 322. 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 1218m below the surface, consists of white, gray and brown, largely fine to coarsegrained quartz sand intermixed with shell beds. It receives recharge directly from precipitation. The second layer, which is 36 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 3743 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 subtropical, 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 threemeter 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 SelfLeveling 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 threedimensional graphical interpolation program, to develop threedimensional 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 stagevolume 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 fivegallon 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 StageVolume Relationships Water and peat depth data were used to produce a threedimensional contour plot of SV5 (Figure 323 and Figure 324). The regression equation that describes the volume of water above the peat layer is V =3454S2148.45s (330) 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 VpV, volumee of saturated peat. The coefficient of determination equals unity for both of these regressions. Figure 323. Top surface of peat layer of wetland Figure 324. Bottom surface of peat layer Stagevolume 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 325. Relationship between wetland volume and water stage 56 Switt et al. (1998) used a similar formula (Eq. 332) to calculate the interaction between a wetland and the underlying groundwater: SPK GW) SA (332) 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 dewatered 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 38. Thus, seepage for the dewatered area of the wetland was SP=Ki( W SA, (333) 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 FiniteElement 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 326). 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 326. Finiteelement mesh of wedge of wetland and aquifer for the saturated flow model 3.6.5 Solution Procedure To simulate the wetlandaquifer 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 nonlinear 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 finiteelement 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 finiteelement 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 = 105 m1 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 threedimensional groundwater flow model near lakes. In this finiteelement 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 327 shows the head distribution beneath the wetland. Figure 328 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 330 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 331 shows excellent agreement between the simulated results given by the saturated flow model and observed groundwater levels. Table 34 to Table 38 show the modeled data used to calculate the interaction discharges at 5 different times. Table 38 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 327. Vertical crosssection view of threedimensional finiteelement 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 328. 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 329. 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 330. Discharge between wetland and aquifer with respect to time 0 20 40 60 Time (hour) WETL OBS WETL CAL GW OBS GW CAL 80 100 Figure 331. Comparison between calculated head and observed head Table 34. 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 35. 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 36. 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 37. 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 38. 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 39. Sensitivities calculated by saturated groundwater flow model Parameter (PR) Percentage varied APR I hih211 (m) II hih211/APR Ka(m/day) 0.10 0.85 0.9924E02 1.1675E02 Ka(m/day) +0.10 +0.85 0.9055E02 1.0653E02 Kp(m/day) 0.10 0.02 0.6410E02 0.3205 Kp(m/day) +0.10 +0.02 0.5564E02 0.2782 So 0.50 0.000005 0.2068E04 4.136 So +0.50 +0.000005 0.1265E04 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 threedimension finite element model in this dissertation, and the simulated results were plotted on the Figure 332. 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 threedimension 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 333. 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 332. Comparison results using Walser's parameter values. Figure 333. 3D numerical model simulated results over specific storage Table 310 shows the hydraulic conductivities in sand layer and peat layer. The variations are rather large. 68 Table 310. Hydraulic Conductivity values for sand and peat layers k in sand layer k in peat layer (m/day) (m/day) Switt et al.(1998) 712 0.03 Slug test 2 Rycroft et al. (1975) 107102 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 threedimensional finiteelement numerical saturated groundwater flow model using deformable mesh and realtime wetlandaquifer 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 THREEDIMENSIONAL FINITEELEMENT 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 inputoutput 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 timeconsuming. 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 firstorder 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 excitationresponse functions of the real system. Since the model structure for the wetland/groundwater system is assumed to follow Eq.5153, 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 threedimensional 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 finitedifference schemes. Sun and Yeh (1985) extended the method to the case of a finiteelement scheme. The sensitivity coefficients can be computed using the following equation: S q(xt) Vh(xy,) ddyd (41) 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 ThreeDimensional FiniteElement 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 (42) subject to the initial and boundary conditions: h(x,y,z,0) = fo (xy,z)EF (43) hir, =f, KVhnlr =f2 (44) 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 41 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 41. 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. 42 and subsidiary conditions, the following is obtained: v.(6 KVh) + V *( KV6h ) = S, (xvz)Q, Ot5tf (45) subject to the conditions: 6h(x,y,z)l o = 0 (xy,z)eQ (46) 6hIr2 = 0, [6KVh+KV8h]nr = 0 (47) Eq. 45Eq. 47 represent the variational problem of the primary problem in Eq. 42. The variational problem relates the variation 6k and 6h. Assuming h is an arbitrary function having continuous secondorder space derivatives on 0 and firstorder derivatives in [0, tf], then multiplying Eq. 45 by and integrating the result over timespace domains yields: (ST68h)I d f [S_6h + TV(8kVh) + TV'(kV6h) dOdt = 0 (48a) 0 00 If we define: V(x,y,z,tf) = 0 (x0y,z)eQ (48b) then the first term of Eq. 48a vanishes because of condition 47. Thus, Eq. 48b reduces to: h a + TV(6kVh) + V.(kV6h) ]ddt = 0 (49) 0 Q where: YV(8kVh)dQ = T(8kVh)nds  8kVhVWdV 'V(kV8h)dD = Y(kVbh)nd   kV8hVfd = I a a 8h[V(kVT)]ad  6hkVWTn ds Using identities of Eq. 410, Eq. 411 and Eq. 412, Eq. 49 can be written as: and: (410) (411) (412) ffS36h + h[V(kVT)]65kVhVIY d dt 0 0 + f If6kVh*nds+ fIkV6hnds f 8hkVWnds dt= 0 (413) 0 s s s Using conditions 47 and defining: Yr2 = 0, kV'nlr, = 0 then the second term of Eq. 413 is equal to zero. Eq. 413 produces: [s +V(kVT) dhd} dt (V .Vh)dkdOdt = 0 (414) 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 (415) 00 a Taking the variation of E(h, k) yields: I E(h, k)= 0 f 0 a The summation of Eq. 416 and Eq. 414 gives:  s5 V(kVW) dhdQdt + (417) S ++VTVh 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: (418) 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. 417 related to unknown 5h will vanish. Thus: f68k+ f 6h/d dt = 0 9k 9h ) (416) 8E = 0 0 t 8E= ff( +VYVh dkdQdt 00 From Eq. 419, the following partial derivative can be defined: 6E Ak I f f+VTVh 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. (419) (420) (421) (422) (422a) 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 inputoutput relation of the model fit any observed excitationresponse 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) (423) where: L E(k) = [hahobs (424) 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. 423). 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. 42) 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 (425) subject to: TPL=0=0, (x,y)eO 1r2=0=0 kV'nIr, =0 where h is the total head; T = T (x,y,tfT) 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(XXobs)8(yy obs) (ZZobsb)6(ttj) i=1 (428) where ( xobs, Yobs, Zobs) is the location of the observation well; and t, (j =1,2,,J) are the observation times. Integrating, Eq. 428 over time space domains produces: E(h,k)= f(h,k;x,y,z,t)dQdt (429) By solving the governing Equation 42 and the adjoint Equation 425 once, we may obtain the following functional partial derivative required to solve Equation 423 is obatined: and (426) (4.27) E f f 1Vhddt (430) 4.5 Searching Method and Stopping Criterion Since the gradient VE(k) is known, many searching methods including the gradient method, quasiNewton method, and the FletcherReeves method could be used to solve the inverse problem 423 (Sun 1985). Two stopping criteria used here are: IIK,K1il and: IIE(kn)E(k,1) 1 Using the adjoint method and the FletcherReeves 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(xxobs)(yyobs)6(zzbs)8(ttj) (433) Waterlevel data were measured every hour for a total duration, of 87 hours. An example Waterlevel 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.429 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, ) (434) 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 42 and Table 41 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 surfacewater 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 finiteelement 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 43 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 42. Relation between error and peat hydraulic conductivity Table 41. 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 