Two-dimensional analytical and three-dimensional finite-element method modeling of the interactions between wetlands and...

MISSING IMAGE

Material Information

Title:
Two-dimensional analytical and three-dimensional finite-element method modeling of the interactions between wetlands and groundwater
Physical Description:
xvi, 126 leaves : ill. ; 29 cm.
Language:
English
Creator:
Zhong, Jinhua
Publication Date:

Subjects

Subjects / Keywords:
Civil and Coastal Engineering thesis, Ph. D   ( lcsh )
Dissertations, Academic -- Civil and Coastal Engineering -- UF   ( lcsh )
Genre:
bibliography   ( marcgt )
theses   ( marcgt )
non-fiction   ( marcgt )

Notes

Thesis:
Thesis (Ph. D.)--University of Florida, 2004.
Bibliography:
Includes bibliographical references.
Statement of Responsibility:
by Jinhua Zhong.
General Note:
Printout.
General Note:
Vita.

Record Information

Source Institution:
University of Florida
Rights Management:
All applicable rights reserved by the source institution and holding location.
Resource Identifier:
aleph - 003163679
System ID:
AA00008997:00001

Full Text










TWO-DIMENSIONAL ANALYTICAL AND THREE-DIMENSIONAL
FINITE-ELEMENT METHOD MODELING OF THE INTERACTIONS BETWEEN
WETLANDS AND GROUNDWATER


















By

JINHUA ZHONG


A DISSERTATION PRESENTED TO THE GRADUATE SCHOOL
OF THE UNIVERSITY OF FLORIDA IN PARTIAL FULFILLMENT
OF THE REQUIREMENTS FOR THE DEGREE OF
DOCTOR OF PHILOSOPHY

UNIVERSITY OF FLORIDA


2004















TABLE OF CONTENTS

ACKNOWLEDGMENT ................................................. ii

LIST OF TABLES ....................................................... v

LIST OF FIGURES ..................................................... vi

LIST OF SYM BOLS .................................................... ix

A B STR A CT .......................................................... xiii

CHAPTER ........................................................... .

1 INTRODUCTION ...................................................... 1

1.1 W etland ...................................................... 1
1.2 Wetland Water Budget ........................................... 2
1.3 Goals and Objectives ............................................ 4
2 ANALYTICAL GROUNDWATER MODEL RESPONSE TO THE PERIODIC
OR NONPERIODIC WATER STAGE FLUCTUATION IN A CIRCULAR
LAKE OR WETLAND ............................................. 7
2.1 Introduction ................................................... 7
2.2 Problem Formulation: Groundwater Flow Around a Circular Wetland ..... 8
2.3 Analytical Solution for the Periodic Boundary Condition ................ 9
2.4 Analytical Solution for the Non-periodic Boundary Condition ........... 12
2.5 Application .................................................. 14
2.6 Comparison of Analytical Solution and Numerical Solutions ............ 17
2.7 Sum m ary .................................................... 18
3 TRANSIENT THREE-DIMENSIONAL DEFORMABLE FINITE-ELEMENT
SATURATED GROUNDWATER MODEL .......................... 19
3.1 Introduction .................................................. 19









3.2 Groundwater Flow Models Based on Finite-Difference
and Finite-Element Approximation .................... .. ..... 21
3.3 Three-Dimensional Deformable Finite-Element
Saturated Groundwater Model ................................. 24
3.3.1 Governing Equation and Boundary Conditions ............... 27
3.3.2 Finite Element Model of Saturated Groundwater Flow ......... 28
3.3.3 Adaptation of the Saturated Finite-Element Model to
Describe the Free Surface .............................. 36
3.4 Validating the Code of the Three-Dimensional Deformable
Finite-Element Groundwater Flow Model ....................... 37
3.4.1 Theis Solution ........................................ 37
3.4.2 Groundwater Flow around a Circular Wetland ............... 41
3.5 A Numerical Model Application to Simulate the Interaction Between
a W etland and an Aquifer .................................... 42
3.6 Application of Three-Dimensional Deformable Saturated
Groundwater Flow Model to SV5 Wetland ...................... 50
3.6.1 Site Description-SV5 ................................... 50
3.6.2 Field Methods and Wetland Aquifer Interaction Test ........... 52
3.6.3 Stage-Volume Relationships .............................. 54
3.6.4 Creation of Finite-Element Model Mesh .................... 56
3.6.5 Solution Procedure ..................................... 57
3.6.6 Model Calibration and Results ..................... .... 58
3.7 Comparison with Results of Wise and others ........................ 69
3.8 Sum m ary .................................................... 71
4 INVERSE MODEL OF THREE-DIMENSIONAL FINITE-ELEMENT
METHOD SATURATED GROUNDWATER MODEL IN SEARCHING
FOR THE HYDRAULIC CONDUCTIVITY OF THE PEAT LAYER ....... 69
4.1 Introduction .................................................. 69
4.2 Adjoint Problem for Three-Dimensional Finite-Element
Saturated Groundwater Model ............................... 73
4.3 Searching Objective Function .................................... 80









4.4 Adjoint State Controlling Equations ............................... 80
4.5 Searching Method and Stop Criterion .............................. 82
4.6 Field Application on SV5 ....................................... 82
4.8 Summary ................................ .................... 85
5 THREE-DIMENSIONAL FINITE-ELEMENT VARIABLY-SATURATED
GROUNDWATER MODEL ....................................... 87
5.1 Introduction ................................ .................. 87
5.2 Three-Dimensional Finite-Element Variably-Saturated
Groundwater Model ......................................... 88
5.2.1 Governing Equation .................................... 88
5.2.2 Linearizing the Governing Equations ....................... 89
5.2.3 Boundary Conditions for the Variably-Saturated Model ........ 91
5.3 Finite-Element Method for Modeling Three-Dimensional
Variably-Saturated Flow ..................................... 94
5.4 Transient, Variably-Saturated Water-Table Recharge Example Test ...... 98
5.5 Application on SV5 W etland .................................... 101
5.5.1 Finite-Element Mesh for Three-Dimensional
Variably-Saturated Flow Model ......................... 101
5.5.2 Modeling Results and Summary .......................... 103
6 SUMMARY AND CONCLUSIONS ..................................... 111

REFERENCES ....................................................... 113

BIOGRAPHICAL SKETCH ............................................. 126














LIST OF TABLES


Figure page

2-1 Parameters and their values ............................................ 17

3-1 Selected early references for the finite-difference method .................... 22

3-2 Selected early references for the finite-element application to model groundwater 25

3-3 Parameters and their values .................................. .......... 42

3-4 Data for calculating the discharge to the wetland from the aquifer at 10 hours .... 63

3-5 Data for calculating the discharge to the wetland from the aquifer at 20 hours .... 64

3-6 Data for calculating the discharge to the wetland from the aquifer at 30 hours .... 64

3-7 Data for calculating the discharge to the wetland from the aquifer at 40 hours .... 65

3-8 Data for calculating the discharge to the wetland from the aquifer at 50 hours .... 65

3-9 Sensitivities calculated by saturated groundwater flow model ................. 66

4-1 Error and sensitivities based on the adjoint method model .................... 84

5-1 Data for calculating the discharge to the wetland from the aquifer at 10 hours ... 107

5-2 Data for calculating the discharge to the wetland from the aquifer at 20 hours ... 108

5-3 Data for calculating the discharge to the wetland from the aquifer at 30 hours ... 108

5-4 Data for calculating the discharge to the wetland from the aquifer at 40 hours ... 109

5-5 Data for calculating the discharge to the wetland from the aquifer at 50 hours ... 109

v









5-6 Sensitivities calculated by variably saturated groundwater flow model ......... 110














LIST OF FIGURES


Figure page


2-1 Setup of the groundwater problem around a circular wetland ................... 9

2-2 Transient groundwater level changes at radial distances 1,2,3 and 4 meters
away from the wetland shoreline ..................................... 15

2-3 Groundwater distribution over the distance to the wetland shoreline
at 6,7,8,9,10 h .................................................... 15

2-4 Three-dimensional view of the groundwater table around
a circular wetland at t=3hr .......................................... 16

2-5 Comparison of the analytical solution and numerical solution
at 2 and 5 meters from the wetland shore line ........................... 18

3-1 Conceptual profile model of a wetland/groundwater system .................. 26

3-2 Comparison of finite element mesh and finite difference mesh ................ 26

3-3 Setup of the saturated groundwater problem around a wetland ................. 27

3-4 Flowchart illustrating the iteration procedure of finding
the coupling boundary condition ..................................... 29

3-5 Unsteady flow to a well in phreatic aquifer at time t ......................... 38

3-6 Finite-Element Method mesh of the wedge used to simulate
of a well pumping in phreatic aquifer ................................. 40

3-7 Comparison of the numerical and analytical head for the example problem
of a transient pumping test in a phreatic aquifer ......................... 40

3-8 Comparison of the numerical and analytical simulated groundwater table ........ 41

vii









3-9 Three-dimensional finite-element mesh for the simulation of wetland and aquifer 43

3-10 Relationship between the area of wetland water surface ..................... 44

3-11 Relationship between the deepest water depth and the volume of wetland ....... 45

3-12 Illustration of the initial condition where it was assumed a wetland
surface-water stage and the phreatic surface elevation were the same ........ 45

3-13 Plan view of the numerical mesh at time zero. ............................ 46

3-14 Three-Dimensional Finite-Element mesh and simulated head contour at hour 5 .. 46

3-15 Simulated head contour at hour 5 ...................................... 47

3-16 Simulated iso-surface plot of 10.6 meter head at hour 5 ..................... 47

3-17 Three-dimensional Finite-Element mesh and head contour
on a cross section on a cross section through wetland ..................... 48

3-18 Simulated head contour on the vertical cross section through wetland .......... 48

3-19 Modeled iso-surface head plot of 11.1 m,l1 m and 10.5 m ................... 49

3-20 Mesh and modeled head contour with vertical cross section
through the pumping well .......................................... 49

3-21 Head contour (m) on the vertical cross section through the pumping well ....... 50

3-22 W ell location map (Switt et al. 1998) ................................... 51

3-23 Top surface of peat layer ............................................. 54

3-24 Bottom surface of peat layer .......................................... 55

3-25 Relationship between wetland volume and water stage ..................... 55

3-26 Finite-element mesh of wedge of wetland and aquifer
for the saturated flow model ........................................ 57

3-27 Vertical cross-section view of three-dimensional finite-element mesh
and head from the saturated flow model ............................... 61









3-28 Piezometric head at positions beneath wetland with time .................... 61

3-29 Variation in the inundated area of wetland during simulation ................ 62

3-30 Discharge between wetland and aquifer with respect to time ................. 62

3-31 Comparison between calculated head and observed head .................... 63

3-32 Comparison results using Walser's parameter values ....................... 67

3-33 3-D numerical model simulated results over specific storage ................. 67

4-1 Sketch showing setup of the saturated groundwater
inverse problem around a wetland .................................... 74

4-2 Relation between error and peat layer hydraulic conductivity .................. 84

4-3 Flowchart illustrating the iteration procedure of inverse problem .............. 85

4-4 Vertical cross-section of the adjoint state distribution at T = 10 hours ........... 86

5-1 Sketch showing setup of the variably saturated groundwater
flow problem around a wetland ...................................... 92

5-2 Schematic of relative evapotranspiration (ET/PET),
as affected by soil water potential,j .................................. 93

5-3 Schematic diagram of the flow domain ................................... 99

5-4 Water moisture distribution (x:y=2:1) at 4 hours .......................... 100

5-5 Water moisture distribution at 8 hours (x:y=2:1) .......................... 101

5-6 Finite-element mesh of the wedge part for
variably saturated flow model (Unite is m) ............................ 102

5-7 Calculated piezometric head variation beneath wetland using
the three-dimensional variably saturated groundwater flow model .......... 104

5-8 Model predicted discharge from the aquifer to the overlaying wetland ......... 104

5-9 Inundated area of wetland during simulation ............................. 105









5-10 Simulated and observed heads in wetland and in aquifer ................... 105

5-11 Simulated head from the saturated flow model and the variably
saturated flow model ............................................. 106

5-12 Horizontal view of the three-dimensional finite-element mesh and piezometric
head distribution from the variably saturated flow model at time of 9 hours .. 106

5-12 Enlarged top and middle part of the mesh and head contour at time of 9 hours .. 107














LIST OF SYMBOLS


C = the conductance of the wetland confining unit

Co = coefficient

di(t) = the unknown values of hydraulic head and at time t

ET = evapotranspiration [ULT]

F(t) = the Fourier coefficient

G = ground-water level [L]

Gi = groundwater inflow [L3/T]

Go = groundwater outflow [L3 /T]

GWNs = the net flow or seepage of groundwater to the wetland surface [L3fT]

h = piezometric head [L]

h(e) = the approximate solution for hydraulic head within element e [L]

Ho = initial saturated aquifer thickness [L]

H' = Sobolev space

Ho(t) = the Hankel function order 0

hoo = the regional average water table elevation, [L]

i = reference integer spatial interval

Jo = the Bessel function of the first kind order 0

k = a constant parameter

K = the saturated hydraulic conductivity tensor [ IT]

Kh/KI = anisotropic ratio

Ks = hydraulic conductivity [LIT]

Kxx = horizontal hydraulic conductivity in x direction [L/T]









Kyy = horizontal hydraulic conductivity in y direction [L/T]
K= vertical hydraulic conductivity [L/T]

Ki. = the saturated hydraulic conductivity tensor [L/T]

K = the relative permeability while the degree of water saturation is Sw
L = average peat thickness [L]

n = current integer time level
Nie) = the interpolation functions for each node within element e
n = a unit outward normal vector to the boundary

No = the number of observation wells

N, = the total number of nodes used in the numerical solution

p = the number of nodes within element e
Pn = net precipitation ( total precipitation minus intercepted
precipitation[ L3/T]
Pc = the capillary pressure ( Pc = 1) [MILT2]
PET = the daily potential evapotranspiration [L/T]

Q = pumping discharge [L3T]
r = the radial coordinate [L] radial distance from the well [L]

ro = the radius of the circular lake [L]
s = drawdown [L]

S = the specific storage
Si = surface inflows, including flooding stream [ L/T]

So = surface outflows [L3T]
Ss = the specific storativity [L-']
SA, = wetland wetted surface area [L2]
SAv = wetland dry surface area [L2]
SP = exchanging flux [L3/T]
S = the degree of saturation, equal to 06/6









S = the residual saturation

Se = effective saturation

t = time [T]

tf = the time duration of the simulation [T]

t = the observation times [T]

T = the transmissivity in the aquifer [L2IT]

T, = tidal inflow (+) or outflow (-) [L3/T]

u = space function

Vw = the wetland volume [L3]

Vp = the volume above the bottom of the peat layer [L3]

W = wetland water level [L]

Wi( = the weighting function

x, y, z = the axes in the physical domain

X = the function of time

Yo = the Bessel function of the second kind order 0

z = the vertical coordinate or elevation head ( L)

V = the gradient operator

AH = the average difference in hydraulic head [L]

Ar = the radial step [L]

At = the time step [T]

AV = change in the volume of water storage in wetland [L3]

Ax = the mesh space step in x-direction [L]

Az = the mesh space step in y-direction [L]

60 = the angular frequency

S,r] = the axes in the computational domain
r = the boundary of volume 0









F, = flux boundary

F2 = head boundary
y = specific weight of a unit volume of water = pg [M/L2T2]
i = the pressure [M/LT2]
Q, = the exclusive subdomain of node i
8K = variation of K [L/T]
8h = variation of head h [L]

'(x,y,tf- = the adjoint state of h
T)
S= tf- t [T]
a and m = parameters obtained from a fit of above equations to experimental
results
6 = the water content, [L/L3]

6s = the saturated water content, [L3/L3]
Or = residual moisture content, [L3/L3]
= pressure head; [M/LT2]
l* = the value of i when ET = 0.5PET, [M/LT2]

wi = the i corresponds to PE/PET=0.05, [M/LT2]
,W, = the value of i when PE=0.95 PET [M/LT2]














Abstract of Dissertation Presented to the Graduate School
of the University of Florida in Partial Fulfillment of the Requirements
for the Degree of Doctor of Philosophy

TWO-DIMENSIONAL ANALYTICAL AND THREE-DIMENSIONAL
FINITE ELEMENT METHOD MODELING OF THE INTERACTIONS BETWEEN
WETLANDS AND GROUNDWATER

By
Jinhua Zhong
December, 2004
Chairman: Louis H. Motz
Cochair: Kirk Hatfield
Major Department: Civil and Coastal Engineering

Simulating saturated/unsaturated flow between wetlands and the contiguous subsurface

is complex, because it involves the modeling of wetland storage changes that in turn cause

changes in pertinent hydrologic boundary conditions. Existing groundwater models are

incapable of simulating the flow between a wetland and an aquifer, if small changes in

wetland storage cause significant vertical and lateral movements of the wetland boundary.

In addition to the simulation complexities introduced by a change in wetland geometry and

storage, periodic and nonperiodic fluctuations of the wetland stage induce corresponding

elevation changes in the phreatic surface of the contiguous aquifer. As these induced waves

propagate through an aquifer, friction causes a loss of energy, which is manifested as

spatially dampened potentiometric head perturbations along the inland direction. An

analytical model was developed that describes subsurface flows around a wetland induced

by any nonperiodic fluctuations in surface-water stage, such as a flood wave. Not considered,









however, are any changes in wetland storage, or wetland geometry as a function of surface

water stage. The analytical model was validated using a finite-difference numerical

groundwater model that is based on the governing equation expressed in radial coordinates.

Comparisons of analytical and numerical results show excellent agreement.

Two numerical models were specifically designed to simulate wetland-aquifer

interactions. The first model is a saturated groundwater flow model that incorporates

adaptive simulation technologies to permit real-time simulation of moving wetland

boundaries and their associated local influence on the phreatic surface. This model is

numerically efficient and uses deformable hexahedral finite elements to trace phreatic surface

changes in real-time. The second model also uses deformable hexahedral finite elements;

however, this model simulates variably saturated groundwater flow. For both models, the

numerical elements deform as required to characterize the horizontal and vertical extent of

the moving wetland boundary (which serves to define a location within the porous system,

where the pressure is known). Simulation results from both models were applied to a field

site where water was rapidly withdrawn from an isolated wetland to observe system response

and evaluate the wetland aquifer hydraulic connection.














CHAPTER 1
INTRODUCTION

1.1 Wetland

Wetlands are transitional lands between terrestrial and aquatic systems where the

water table is usually at or near the surface, or the land is covered by shallow water

(Cowardin et al. 1979). Wetlands are also called marsh, fen, peatland or water and can be

natural or artificial and permanent or temporary. The water can be static or flowing and fresh

or brackish, or it can include areas of marine water the depth of which at low tide does not

exceed 6 meters (Ramsar Convention 1971, Middleton 1998).

Wetlands perform various environmental functions that include increasing flood

mitigation, enhancing storm abatement (i.e., in the case of coastal wetlands), adding

groundwater recharge, providing water quality treatment, serving as wildlife habitats, and

performing aesthetic functions (Mitsch and Gosselink 1993). These functional attributes can

be classified into three basic types: hydrologic, water quality, and life support. Measures of

the life-support function must integrate several wetland attributes, including water quality,

hydroperiod and hydrodynamics, habitat structure, food chain support, biogeographic setting,

nutrient status, corridors for migration, and primary production (Preston and Bedford 1988).

In recent years, increasing interest has focused on the study and protection of

wetlands. Development, drainage, population pressure, and pollution have destroyed

extensive areas of these sensitive ecosystems. Another increasingly significance to wetlands

is municipal well fields supplying water to large urban areas. These well fields affect

1








2

wetlands by reducing surficial aquifer levels over large areas that encompass isolated and

interconnected wetlands (Sonenshein and Hofstetter 1990). As a result of depressing the

water table, historically wet ecosystems have been drained (Switt et al. 1998).

1.2 Wetland Water Budget

Hydrologic conditions are extremely important for the maintenance of a wetland's

structure and function. A hydrologic budget is often used to quantify the volume and rate of

water flowing into and out of a wetland. The primary hydrologic inflows include

precipitation, streams, groundwater, and, in coastal wetlands, tides (Mitsch and Gosselink

1986, Rosenberry 2000). The wetland water budget equation as defined by Mitsch and

Gosselink (1993) is:



AV/At = Pn +S, +G,-ET-So-G,T (1-2)



where

AV = change in the volume of water storage in wetland [L3]

At = change in time [T]

Pn = net precipitation (total precipitation minus intercepted precipitation (I)) [ L3/T]

Si = surface inflows, including flooding stream [ L3T]

Gi = groundwater inflow [L3/T]

ET = evapotranspiration [L3/T]

So = surface outflows [L3 /T]

Go = groundwater outflows [L3/T]

T = tidal inflow (+) or outflow (-) [L3T]








3

Terms in the wetland water budget vary in importance depending on the wetland type,

and not all terms are required to characterize a specific type of wetland. In most wetland

systems, precipitation and evapotranspiration are major components of the water budget

(Winter and Woo 1990). Groundwater can also have a significant influence on some wetland

systems, but not always.

Three of the most common relationships that characterize a wetland are groundwater

discharge, groundwater inflow and outflow, and groundwater recharge. Groundwater inflow

occurs when the contiguous wetland surface water (or groundwater) stage is lower than the

phreatic surface of the surrounding aquifer. When the surface water stage or groundwater

levels inside a wetland are higher than the phreatic surface in the surrounding contiguous

aquifer, groundwater and surface water will flow out of wetland. When a wetland is situated

well above the groundwater, the wetland is identified as being perched. This type of

hydrologic system, referred to as a surface water depression wetland by Novitzki (1979),

loses water only through infiltration into the ground and through evapotranspiration.

The rate of water exchange between a wetland and local aquifer is a function of the

hydraulic gradient between the two hydrologic subsystems. The discharge, in accordance

with Darcy's law (Truesdell 1995), is proportional to the hydraulic gradient, the hydraulic

conductivity, and the cross-sectional area through which the flow occurs. Periodic or

nonperiodic fluctuations of the wetland stage in a flow-through wetland will induce

fluctuations in the elevation of the phreatic surface in the contiguous aquifer. As these

induced waves propagate through an aquifer, friction causes a loss of energy, which is

manifested as spatially dampened head perturbations. This damping effect progressively

diminishes the amplitude of piezometric head along the inland direction (Carr and Van Der








4

Kamp 1969). Many studies have examined the problem through one-dimensional analytical

solutions (Cooper 1963, Gregg 1966, Nielsen 1990, Serfes 1992, Li, Barry et al. 2000, 2001,

Boufadel and Peridier 2002, Seerano 2003, Swamee and Singh 2003). A two-dimensional

analytical solution was obtained by considering a shift in phase and a change in amplitude

along a straight coastline (Sun 1997). Three-dimensional flow regimes near a shallow

circular lake was studied using numerical and analytical mathematical models of aquifer flow

(Townley and Trefry 2000).

To further the understanding of the interactions between lakes or wetlands and

groundwater, without the restrictive assumptions typically associated with analytical

solutions, numerical models are usually used. For example, it is possible to obtain numerical

solutions for the case of anisotropic and nonhomogeneous conditions. The accuracy of

solutions obtained by numerical methods can be excellent; however, it depends on several

factors, including the type of numerical method used, the complexity of boundary and initial

conditions, and the computational precision of the computer used to implement the method.

Regardless of these advantages, few investigators have chosen to focus on the numerical

simulation of interactions between wetlands and groundwater ( McDonald and Harbaugh

1988, Schaffranek 1996, Lee 2000).

1.3 Goals and Objectives

Analytical models to simulate the complex groundwater surface water interactions

that occur with isolated and flow-through wetlands would be valuable tools for developing

hydrologic budgets for wetlands. Given that most surface water fluctuations are nonperiodic,

there is a need to develop analytical groundwater models for simulating groundwater table

fluctuations induced by both periodic and nonperiodic surface water flood waves. However,








5

with the exception of the one-dimensional study by Cooper and Rorabaugh (1963), most

investigations addressed surface-water fluctuations that were periodic. Thus, the first goal

of this dissertation was to develop an analytical method for simulating groundwater table

fluctuations induced by both periodic and nonperiodic surface-water flood waves. This goal

was achieved by meeting the single objective of developing an analytical model using

Fourier transforms.

A limited number of three-dimensional numerical groundwater models have been

developed to simulate groundwater wetland interactions; however, these models are generally

inadequate tools for simulating the complex interactions between a wetland and a contiguous

unconfined aquifer. For example, Walser, Annable and Wise (1998) used FEMWATER

(Yeh, 1992), a variably saturated flow model, to simulate the exchange of water between a

wetland and an aquifer. McDonald and Harbaugh (1988) used MODFLOW and developed

a three-dimensional finite difference groundwater model and studied the dynamic coupling

between wetlands and the subsurface. Lee (2000) used the HYDRUS-2D flow model and

simulated vertical two-dimensional groundwater flow through Lake Barco in Florida. Neither

FEMWATER (1996) nor MODFLOW can simulate or incorporate a moving real-time

wetland boundary caused by a change in surface water stage. Furthermore, as a saturated

three-dimensional finite-difference model, MODFLOW can not simulate fluctuations of the

phreatic surface as a moving transient boundary as required to define fully the three-

dimensional flow. Finally, the water stage of a wetland itself is a function of the transient

exchange of water between a surface and subsurface systems; which is not fully modeled in

any existing groundwater flow codes. Thus, to simulate groundwater surface water

interaction without the above identified shortcomings of existing models, a new model is










needed.

The second goal of the dissertation was to develop two numerical models capable

of simulating 1) a transient and moving wetland boundary, 2) three-dimensional flow, 3) the

movement of the phreatic surface, and 4) the coupling that occurs between the transient

surface-water stage and the phreatic surface. To achieve this goal, two objectives were

pursued. The first objective was to develop two numerical groundwater models using finite

elements. The second objective was to apply both numerical models to simulate groundwater

and surface water fluctuations over an isolated wetland in southern Florida. Of the two

models developed, the first is a three-dimensional saturated groundwater flow model that

employs deformable hexahedral elements. The finite-element mesh deforms in space and

time to trace the location of the phreatic surface of an unconfined aquifer and the stage of the

contiguous surface water body. Based on the above model and the adjoint state method, an

inverse model was developed to estimate the hydraulic conductivity of peat layers which

typically exist as a wetland bottom that serves to isolate surface wetlands from subsurface

groundwater systems. The second numerical model, a three-dimensional variable saturated

flow model, was also developed using deformable mesh. This model was used to investigate

prediction differences between the saturated flow code and the variable saturated model.

From the groundwater flow pattern given by both models, an improved understanding was

obtained for the interactions occurring between surface waters and groundwaters. In addition,

an improved understanding was obtained of how short-term pumping from a phreatic aquifer

affects wetlands and lake water stages.














CHAPTER 2
ANALYTICAL GROUNDWATER MODEL RESPONSE TO THE PERIODIC
OR NON-PERIODIC WATER STAGE FLUCTUATIONS IN A CIRCULAR WETLAND


2.1 Introduction

Periodic or nonperiodic fluctuations of the surface water stage in wetlands or lakes

will induce fluctuations in the elevation of the phreatic surface of contiguous unconfined

aquifers. As these induced waves propagate through the aquifer, friction causes a loss of

energy, which is manifested as spatially dampened head perturbations. This damping effect

progressively diminishes the amplitude of fluctuations in the piezometric head along the

inland direction (Carr and Van Kamp 1969). Many studies have examined the problem

through one-dimensional analytical solutions (e.g.,Cooper and Rorabaugh 1963, Gregg 1966,

Nielsen 1990, and Serfes 1992, Li, Barry et al. 2000, 2001, Boufadel and Peridier 2002,

Seerano 2003, Swamee and Singh 2003). Cooper and Rorabaugh (1963) obtained the

analytical solution for groundwater movement and bank storage due to flood stages in surface

streams. Gregg (1966) found a one-dimensional relation between tidal fluctuations and

groundwater table perturbations. Nielsen (1990), Serfes (1992) and Carr (1969) also

published papers investigating the similar problems. A two-dimensional analytical solution

was obtained by considering a phase shift in water-table perturbations and the associated

change of amplitude along a straight coastline (Sun 1997). An analytical solution

representing vertical two-dimensional groundwater-surface water interaction was obtained

by Anderson (2003). Three-dimensional flow regimes near a shallow circular lake was

7








8

studied using numerical and analytical mathematical models of aquifer flow (Townley and

Trefry 2000, Smith and Townley 2002).

With the exception of Cooper and Rorabaugh (1963), all of the above studies

addressed periodic surface-water fluctuations, in spite of the fact that most stage variations

are not periodic i.e., consider, for example, the pervasive nonperiodic flooding of rivers,

lakes, and wetlands.

Cooper's (1963) nonperiodic analytical model for a river/groundwater system was

a one-dimensional solution along an axis perpendicular to the river. However, this analytical

model is not suitable for simulating the nonperiodic interactions between a wetland and an

aquifer. In this chapter, two analytical groundwater models are developed to describe

transient groundwater levels affected by periodic and nonperiodic fluctuations of the surface-

water stage in a wetland or a lake.

2.2 Problem Formulation: Groundwater Flow Around a Circular Wetland

Considering a circular wetland that is hydraulically contiguous to a surfacial aquifer

and ignoring recharge from rainfall, the governing equation for two-dimensional groundwater

flow in a homogeneous isotropic aquifer can be expressed in radial coordinates as:


2h 1 Dh S Oh
2 h (2-1)
ar 2 r r Tat


where h is the piezometric head, r is the radial coordinate, t is time, S is the specific storage

and Tis the transmissivity in the aquifer. When a thick unconfined aquifer is considered, the

gradient and small fluctuations relative to saturated thickness are assumed to be small, T is

approximately equal to Khoo (Figure 2-1) in which K is the hydraulic conductivity and ho is








9

the average unconfined aquifer thickness equal to the regional average water-table elevation

above an underlying impermeable confining unit.

The boundary condition along the shore line of the wetland is Eq. 2-2, where r o is

the radius of the circular wetland. It is also assumed in Eq. 2-3 that as r--o, the groundwater

table elevation approaches the regional average water table elevation, hoo.



B.C: h(r,t) ,=r, = hoo + f(r ,t) @ r= r. (2-2)



h(r,t) ,r= = hoo @ r= o (2-3)


wetland or lake _



h_

//777/7- //7777//77///77/7 77777/


Figure 2-1. Setup of the groundwater problem around a circular wetland

2.3 Analytical Solution for the Periodic Boundary Condition

Since the boundary condition is a periodic function of time, substituting forfl(rt),

Equation 2-4 is obtained:

B.C: h(r,t) r=ro = hoo + h ei't @ r= r, (2-4)

h(r,t)| ,- = hoo @ r= co (2-5)

where ho is the half height of the periodic constituent and to is the angular frequency. No

initial condition is necessary, since here the steady-state oscillatory solution is sought as t -,.

Because of the oscillatory boundary condition, a solution is sought in the form











h(r,t) =ho +hog(r) e (2-6)
(2-6)

Substitution of Eq. 2-6 into Eq. 2-1, Eq. 2-2, and Eq. 2-3 yields



d2g +l.O g S.i-g=0 (2-7)
Br2 r ar T


where the boundary conditions are:

g(r)=] @ r=ro (2-8)

g(r)=O @ r=- (2-9)

Equation 2-7 can now be rewritten to obtain:


v2 +vg +v2g=
av2 Ov (2-10)





v 2= T (2-11)



The solution to Equation 2-10 subject to the boundary conditions in Equation 2-8 and 2-9

depends on the value of w. Considering the boundary condition along the shoreline,

groundwater around a circular water body will have the following forms.

When (o is greater than 0:



g(r)=Co( +iY)= CoH v)- H() (2-12)
H' (vC)









and when wo is less than 0:



g(r) = C -iY) = CH (u)= (2-13)



where Jo is the Bessel function of the first kind order 0; Yo is the Bessel function of the

second kind order 0; Ho(t) is the Hankel function order 0; and CO is a coefficient (Lebedev

1965). Substitution of Eq. 2-12 and Eq. 2-13 into Eq. 2-6 gives the following solution for the

groundwater flow around a circular lake affected by periodic function:




h,(rt)=hoo+h0o )e
H.(.vo) (2-14)



for 0 > 0 and:


h.(r,t)=h0 +h0 -e Nt(
S (o2)(o) (2-15)



for c < 0.

The asymptotic approximations of the functions Ho')(z) and Ho(2)(z) are:




H (z)= 2 i*) [ 1-32**(2k- 1)2(2 iz)-k+O(z --1) (2-16)
S7 z k= 22k-k! I












i 2 /) 1/2 (-i e( )~ 132..-(2k 1)2
2)(z)= 2 / 4 ( iz)-k+O(|zi -n-1)
VZ k=o 22k-M!


uo= '*(-1-O'i)-ro





V0 = -*(i-1)ro


(2-17)


u= .*(-1-i)-r
2v= T





F2T


2.4 Analytical Solution for the Non-Periodic Boundary Condition

A flood wave will typically induce a nonperiodic variation in the wetland stage. Here,

f(ro,t), which is a nonperiodic function of time, is assumed to describe transient variations

in the wetland stage. Using the Fourier transformation, f(ro ,t) can be represented as:


J(ro,t)= fF()e s"'d"


(2-18)


where F(wo) is the Fourier coefficient, and F(o)dw is the amplitude of oscillation for the

component with a frequency in the range of (o,(+do). For this F(x), the corresponding

fluctuations in the groundwater table are described by:


(W >0)


HO 2)(v))
h.(r,t)=hoo + F(w) *e w


(2-19)














H2)(u)
h,(r,t)=hoo + F(o) "*eu
,, k)


( o < 0 ) (2-20)


as given by Eq. 2-14 and Eq. 2-15. The inverse Fourier transformation of Equation 2-19 and

Eq. 2-20 leads to:


h(r,t)=h oo+ +
2n


=h0 2x
-h'


H'(v)
e e
Hf ~(v0)


J H (u0)


(2-21)


wl.dj


From the above analytical model, it is clear that the unsteady groundwater table around a

circular wetland depends not only on the instantaneous values of the wetland water stage but

also on the accumulative effect represented by the integration. However, physically it is not

possible for the wetland water stage at a later time to affect the groundwater table at an

earlier time. Hence, the phreatic surface depends on its current stage and historical stage. The


_na


=ooh +-








14

final analytical solution to the nonperiodic boundary condition 2-2 is then


h(r,t)=h o+-
In 2x


(2-22)


(1),.a f (v)
H()(o)
0


2.5 Application

If it is assumed that wetland stage fluctuations can be described by a nonperiodic

Gaussian function of time, then:


1 202
(ro, t) =- e 2a
V2-n G


(2-23)


The analytical solution obtained when Eq. 2-23 is combined with Eq. 2-21 is:


1 r 1 ( t H(2)
h(r,t)=h o+ e 202 2 e) e(t-)i*d
In 2 f H2)(Uo)


(2-24)

The transient changes of groundwater table described by Eq. 2-24 at four different

distances from the wetland are shown in Figure 2-2. The groundwater distributions along the

radius at four different times are shown in Figure 2-3. From both figures, it is evident that

groundwater fluctuations are greatest near the wetland. From Figure 2-3, it can be seen that

the variation of the groundwater table over time is no longer symmetrical although the


" '(-)'-d*o dr












variation of the wetland water stage is a Gaussian distribution which is symmetrical over

time.


101 I

100.8 | I

100.6

100.4

100.2

100

99.8
0 5 10 15 20 25
Time (hr)


lake

d=1m

d=2m

d=3m

d=4m


Figure 2-2. Transient groundwater level changes at radial distances 1,2,3, and 4 meters
away from the wetland shoreline.


100.3 -
H00.25
:100.2
3-
00.15
100.1
900.05
(o 100 t
99.95
50


t=6hr

t=7hr

t=8hr

t=9hr

t=10hr


55 60 65
Distance to lake shoreline (m)


Figure 2-3. Groundwater distribution over the distance to
the wetland/lake shoreline at 6,7,8,9, and 10 h.








16

Figure 2-4 shows the groundwater table elevations around a circular wetland after 3

hours. The radius of the wetland is 50 meters. All other parameters are listed in Table 2-1.

Theoretically, when o is very large, it means that the frequency of the corresponding

boundary condition is very large and the period is very short. When the absolute value of a

is very large, its effect on the groundwater around the wetland is very small. This is because

the period of the wave component is too small to have any impact on the groundwater around

the wetland. That is to say, when the integration of a is taken, a certain negative number

could be used to replace the negative infinite number, and a certain positive number could

also be used to replace the positive infinite number. Accurate model results are obtained

when w takes the range of [-5, 5]. Thus, it may be assumed that the boundary condition has

no effect on the ground water table if the frequency of the boundary water level is very large.






H (M)













Figure 2-4. Three-dimensional view of the groundwater table
around a circular wetland at t=3 h.








17

2.6 Comparison of Analytical Solution and Numerical Solutions

A comparison was made between the analytical and the numerical solutions to the

radial flow problem specified through Eq. 2-1-Eq. 2-3. An implicit finite-difference scheme

was used to represent the radial form of the groundwater flow equation Eq. 2-1 as:


h,- -2h" i+h,' I h11 -hi-I _S h, -h(-5)
(Ar)2 r, 2Ar T At


Here i is the reference integer spatial interval; n is the current integer time level; At is the

time step; and Ar is the radial step.

A circular wetland which has a radius of 50 meters was assumed. All other system

parameters are given in the Table 2-1. By considering the different boundary conditions (Eq.

2-3 and Eq. 2-4), the corresponding numerical results are found [Figure 2-5]. Analytical and

numerical results illustrate time variations in groundwater levels. Figure 2-5 shows that the

analytical solution and numerical solution compare very well.

Table 2-1. Parameters and their values

Parameter Value Parameter Value

Hoo (m) 100.0 T (m2/h) 0.635

Ar (m) 0.25 S 0.2

At (hr) 0.02 to (h) 5.0
k (m/hr) 0.00635 o 0.4












0.4

S0.3 analy(d=2m)

0.2
0.2 ............. analy(d=5m)
0.1
-o 0num (d=2m)

2. num (d=5m)
-0.1 7 __
0 5 10 15 20 25
Time (hr)



Figure 2-5. Comparison of analytical solution and numerical solution
at 2 and 5 meters from the wetland shoreline.


2.7 Summary

The analytical model presented in this chapter can be used to simulate aquifer surface

fluctuations caused by the transient wetland stages regardless of the temporal function

assumed to describe variations of water stage. In other words, the model can deal with any

general transient of the boundary function. The analytical model was validated using a finite-

difference numerical groundwater model for a homogenous aquifer. This first numerical

method is not recommended for simulations involving heterogenous aquifers. Rather, a new

method is presented for simulating wetland/aquifer interactions in such aquifers, which

requires a simple transformation of the equations. The excellent agreement between two

models shows that the analytical solution presented is valid. This model is also useful in a

situation where groundwater and the surface water components are being considered in a

system-wide water balance.














CHAPTER 3 TRANSIENT THREE-DIMENSIONAL DEFORMABLE
FINITE-ELEMENT SATURATED GROUNDWATER MODEL


3.1 Introduction

Recent research has focused on the details of various wetland functions such as

water-quality treatment, flood-wave suppression, and groundwater/surface water interactions

(Wise et al. 2000, Krabbenhoft and Webster 1995, Anderson and Cheng 1993, Ruddy and

Williams 1991, and Mills and Zwarich 1986). For example, Ruddy and Williams (1991)

investigated the hydrologic relation between the South Fork Williams Fork stream in the

Colorado River basin and four adjacent subalpine wetlands, and they also evaluated the

potential effects on these wetlands when stream flow is diverted. They concluded that stream

flow was the primary source of water for two of the four wetlands, whereas precipitation,

snowmelt, valley side-slope flow (including overland and small-channel flow), and

groundwater inflow served to define the transient hydrological behavior of all wetlands.

The above studies showed that wetland functions (i.e., water quality treatment,

floodwave suppression, etc.) are correlated with recharge and discharge areas, ground-water-

flow paths and rates, and soil-profile development. The location of recharge and discharge

areas can change in response to seasonal or weather-related factors and can affect the

groundwater flow direction. As a result, groundwater and surface water quality and the

presence of various wetland-plant communities may all vary in response to these changing

conditions.








20

Analytical models such as the one presented in Chapter 2 can be used to simulate

deep wetlands or lakes. However, the primary limitation of analytical methods is that

solutions can only be obtained by imposing severely restrictive assumptions on aquifer

properties and boundary and initial conditions. For example, an assumption commonly made

to obtain analytical solutions is that the aquifer is isotropic and homogeneous for hydraulic

conductivity. These assumptions, however, are not valid and often not adequate to describe

groundwater flow or solute transport in most field situations.

Numerical methods have the primary advantage over analytical approaches in that

they do not require such restrictive assumptions. For example, it is possible to obtain

numerical solutions for the case of anisotropic and nonhomogeneous conditions. The

accuracy of solutions obtained by numerical methods can be excellent. However, this

accuracy depends on several factors, including the type of numerical method used, the

complexity of boundary and initial conditions, and the computational precision of the

computer used to implement the method. Regardless of these advantage, only a few

investigators have chosen to focus on the numerical simulation of interactions between

wetlands and groundwater (e.g., Schaffranek 1996 and McDonald and Harbaugh 1988, Lee

2000, Rosenberry 2000). In this chapter, a general discussion of the

advantages/disadvantages of two numerical methods commonly used to model hydrologic

systems is briefly discussed. Next, the mathematical development of a three-dimensional

deformable finite-element saturated groundwater model is presented. This model and the

associated code are then validated by comparing model predictions with analytical modeling

results. A more complicated application of the three-dimensional model is presented to

illustrate the type of information provided through the numerical approach taken. Finally, the








21

model is used to simulate an actual wetland using field data collected when the hydrologic

system was responding to induced stress.

3.2 Ground-Water Flow Models Based on Finite-Difference
and Finite-Element Approximation

Several types of numerical methods have been used to solve groundwater flow

problems, the two principal ones being the finite-difference method and the finite-element

method. The finite-difference method was initially applied to the flow of fluids in petroleum

reservoirs, and it wasn't until the mid-1960's that problems of groundwater flow and solute

transport were examined. The method has a number of advantages that contribute to its

continued widespread use and popularity: (1) for simple problems such as one-dimensional,

steady-state groundwater flow in an isotropic and homogeneous aquifer, mathematical

formulation and computer implementation are easily understood; and (2) the accuracy of

solutions to steady-state and transient groundwater flow problems is generally quite good.

The finite-difference method also has disadvantages: (1) the method works best for

rectangular or prismatic aquifers of uniform composition (for example, it is difficult to

incorporate irregular or curved aquifer boundaries, anisotropic and heterogeneous aquifer

properties, or sloping soil and rock layers into the numerical model without introducing

numerous mathematical and computer programming complexities); and (2) the accuracy of

solutions to solute transport problems is less than can be obtained through application of the

finite-element method. Table 3.1 lists selected early references where finite-difference

methods were applied to model groundwater flow. More recent applications of hydrologic

modeling using finite-difference for simulating of surface water /groundwater interactions

include studies by Mc Donald and Harbaugh (1998), Winter (1983, 1986, and 1989), Sacks








22

(1992) and Krabbenhoft and Anderson (1990).

Table 3-1. Selected early references for the finite-difference method.

Topic Reference
Early Developments in Petroleum Bruce et al. (1953), Peace and
Reservoir Modeling Rachford (1962)

Saturated Groundwater Flow Remson et al. (1965), Freeze and
Whitherspoon (1966), Pinder and
Bredehoeft (1968).

Unsaturated Groundwater Flow Philip (1957), Ashroft et al. (1962),
Freeze (1971), Brutsaert (1973),
Healy (1990), Kirkland et al. (1992).

Solute Transport Stone and Brian (1963), Oster et al.
(1970), Tanji et al. (1967), Wierenga
(1977)

Application to Field Problems Orlob and Woods (1967), Gambolati
et al. (1973), Fleck and McDonald
(1978), Jain 1992

Comprehensive Reference Trescott and Larson (1977), Ames
(1977), Mitchell and Griffiths
(1980), Lapidus and Pinder (1982)

Computer Program Trescott et al. (1976), Konikow and
Bredehoft (1978), Healy (1990)


Winter (1983, 1986, 1989), Sacks (1992) and Krabbenhoft and Anderson (1990)

conducted numerical investigations on the interactions between groundwaters and lakes.

Both Winter and Krabbenhoft and Anderson used fixed finite-difference meshes. For lakes,

a finite-difference mesh may be suitable; however, because the bed of most wetlands exhibits

significant slope, a poorly designed fixed finite-difference mesh may not adequately

characterize water fluxes across a bed situated between an overlying body of water and an

underlying aquifer. The issue of properly modeling fluxes across a sloping bed may not be








23

as important as the fact that accurate simulations of wetland-flow systems require

measurements of hydrologic properties at finer scales and frequencies than are commonly

made (Hunt 1996, Hunt and others 1996 1997).

Assuming the hydrological model is correct, the accuracy of any computer simulated

hydrologic system, such as a wetland, depends on the quality of the data used to calibrate

model. Through a well instrumented study-site in the Kankakee River Valley in northwest

Indiana, the USGS and the USEPA conducted detailed measurements of many various

variables required to characterize a wetland hydrologic system. McDonald and Harbaugh

(1998) developed a groundwater model using the available data and conducted the

simulations of observed hydrologic process using the computer code commonly known as

MODFLOW. The focus of the research effort was to study the effects of various wetland

restoration activities on local groundwater levels. Computer simulations were used to predict

the effects of wetland restoration on the existing hydrology of the wetlands and adjacent

farmlands and to suggest residence times and path ways for surface water and groundwaters.

The finite-element method was first used to solve groundwater flow and solute

transport problems in the early 1970s. The method has several advantages: (1) irregular or

curved aquifer boundaries, anisotropic and heterogeneous aquifer properties, and sloping soil

and rock layers can be incorporated easily into the numerical model; and (2) the accuracy of

solution to groundwater flow and solute transport problems is quite good. The main

disadvantages of the finite-element method include: (1) for simple problems, the finite

element method requires a greater amount of mathematical and computer programming

sophistication than does the finite-difference method; and (2) there are fewer well-

documented computer programs.








24

The finite-element method provides approximations of much higher order than does

finite-difference methodology. The whole modeling domain can be subdivided into irregular

elements as dictated by the physical geometry of the problem. The size of each element can

vary. Small elements can be used in areas where rapid change occurs with model parameters

and predicted state variables, whereas the large elements can be used where such variations

are less severe. Inhomogeneities and anisotropy are easily accommodated with the finite-

element method. The advantage of irregular elements and higher-order approximations is

certainly important when the "space" domain is complex. However, this advantage is less

obvious in problems involving time (evolution) (Gupta 1975). Selected references of early

finite-element applications to model groundwater flow are listed in Table 3-2.

3.3 Three-Dimensional Deformable Finite-Element
Saturated Groundwater Model

To simulate three-dimensional groundwater flow in phreatic aquifer systems and the

exchange of water between a wetland and the contiguous aquifer, a numerical model is

needed that is capable of modeling the free groundwater surface. This requires first that the

model assume a tentative location of the free surface and then solve for the distribution of

the dependent variable, *ii, over a fixed domain. All this requires that the flux boundary

condition across the surface be satisfied and that the pressure be atmospheric. To model the

exchange of water between a wetland and an aquifer, it is necessary to recognize that

significant vertical and lateral movements of the wetland boundary typically occur with

changes in wetland storage (see Figure 3-1). This transient boundary suggests that a

deformable finite-element mesh may be needed in the model to characterize adequately the

moving groundwater/surface-water interface.








25

Table 3-2. Selected references of early finite-element applications to model groundwater.

Topic Reference
Early Developments in Petroleum Price et al. (1968)
Reservoir Modeling
Saturated Groundwater Flow Zienkiewicz et al. (1966), Javandel
and Witherspoon (1968),
Zienkiewicz and Parekh (1970),
Pinder and Frind (1972).
Unsaturated Groundwater Flow Neuman (1973), Gureghian et al.
(1979), Pickens and Gillham (1980)
Yeh and Ward (1981),Huyakom
(1984), Yeh (1992).
Solute Transport Price et al. (1968), Guymon et al.
(1970), Neuman (1973), Van
Genuchten et al. (1977), Kirkner et
al. (1984).

Application to Field Problems Pinder (1973), Gupta and Tinji
(1976), Senger and Fogg (1987).
Comprehensive References Ziekiewicz (1971), Pinder and Gray
(1977), Lapidus and Pinder (1982),
Huyakorn and Pinder (1983).
Computer Programs Neuman and Witherspoon (1970),
Reeves and Duguid (1975), Segol et
al. (1975), Pickens et al. (1979)
Huyakom (1984), Yeh (1992).


Numerical groundwater modeling that is focused on simulating the interaction

between surface-water and groundwater systems, where the groundwater/surface-water

interface moves, must use well-defined linking conditions to describe the exchange of water

across the interface. These linking conditions are unknown in the beginning of the modeling

time step and found through a flux type interaction that then forecasts the linking boundary

conditions at the end of each time step. Thus, the model has transient boundary conditions.












Wetland
SWater Surface at Time t,
SWater Surface at Tim --


Phreatic Surface at Time t, Phreatic Surface at Time t2
Groundwater Aquifer
Impermeable Aquacludc




Figure 3-1. Conceptual profile model of a wetland/groundwater system

Finally, typical finite-difference modeling is not able to trace the curvature of the bottom of

a wetland (see Figure 3-2). Finite elements, however, are easily configured to fit non-linear

spatial geometries. In conclusion, to cope with the above described three technical problems

of modeling, i.e., the free surface location, incorporating dynamic changes in the vertical and

horizontal extent of the ground/surface-water interface, and tracing the curvature of a

wetland bottom, a special three-dimensional finite-element model is needed.










Finite-Element Method Finite-Difference Method


Figure 3-2. Comparison of finite-element mesh and finite-difference mesh








27

3.3.1 Governing Equation and Boundary Conditions

Eq. 3-1 below can be used to describe incompressible, saturated groundwater flow

beneath a wetland and in a phreatic aquifer:



V.(KVh) Q = S, ah ij= 1,2,3 (3-1)
at


where K is the saturated hydraulic conductivity tensor ( LT -'); h is equal to z+* (L); i is

pressure head (L); z is elevation head (L); and S, is the specific storativity (L'1). From Figure

3-3, the boundary of volume 0 is characterized through a combination of specified flux

boundaries r, and specified head boundaries P2. Pertinent specific no-flux boundaries include

the groundwater phreatic surface (EC and DF) and the bottom impermeable boundary (GH).

Relevant specified head boundaries r2 include the wetland head boundary (AB and CD) and

lateral boundaries (EG and FH).



E F

Peat Layer
Groundwater Aquifer
G HI



GH, EC, and DF are no-flow boundaries F

EG, FH, ACDB, and CD are specified head boundaries rF


Figure 3-3. Setup of the saturated groundwater problem around a wetland








28

As to the boundary condition controlled by the wetland water table, the initial water

table of the wetland and the pumping discharge were taken as the only inputs of the coupling

model. During pumping, the pumping discharge lowered the wetland water table, and it also

triggered an inflow from the aquifer to the wetland through the peat layer. In each simulation

time step, the pumping discharge, estimated inflow from the aquifer to the wetland, initial

water table of the wetland, and relation between the volume of the wetland and wetland

water table were used as inputs in a mass balance equation that was solved to determine the

water table at the end of the time step. The estimated water table at the end of the time step

was put back into the saturated groundwater model to adjust the inflow from the aquifer to

the wetland estimated initially. Therefore, several iterations in each time step were necessary

to get an accurate estimated inflow from the aquifer to the wetland and a convergent solution

of the water table in the wetland at the end of the time step. When pumping was stopped, the

pumping discharge was taken as zero. Figure 3-4 shows the procedure of finding the inflow

from the aquifer to the wetland and the coupling boundary.

3.3.2 Finite Element Model of Saturated Groundwater Flow

To obtain a finite element approximation of Eq. 3-1, two spaces, S and V, are defined

(Hughes 1987):



S= {uIuEH1, u(xyz)2=qr2 } (3-2)



V= uuEHuH, u(x,y,z)r2=O} (3-3)


where H' is Sobolev space comprised of functions which are square-integrable generalized












Predict a wetland water table hl(t+At)



Simulate groundwater flow based
on the predicted wetland water table


Find the groundwater inflow to the wetland

Wetland
pumping Based on the water mass balance equation of the
discharge Q wetland to find the wetland water table h2(t+At)







Yes

Next time step



Figure 3-4. Flowchart illustrating the iteration procedure of finding
the coupling boundary condition

derivatives through order 1; u is a space function that is member of Sobolev space; qr` is the

vertical flux along P2; and x, y, z are three spatial coordinates. Any function from S space is

equal to the flux boundary value along specified head boundary (F2). Also, any function u

from V space is equal to zero along head boundary (F2). Assuming a basis function w is

chosen from V space, then the following integration of Eq. 3-1 over space domain Q is

always true:










w[V-(KVh)]dQ wQdO= wS, -'dQ, ij= 1,2,3 (3.4)
J fa f at


where Eq. 3-4 expresses the governing groundwater flow equation integrated over the
modeling domain.
Manipulation of the first term on the left hand side of Eq. 3-4 begins with the
application of Green's theorem, which is an expanded form of the divergence theorem. This
converts the integral into two terms, one of which is evaluated only on the surface of the
region to be simulated. Green's theorem is:



(v.-P)AdQ= j ( -.)AdT- (-VA)d (3-5)
(3-5)
o r a

where A is a scalar, and P is a vector quantity. The boundary of volume 0 is defined by F,
and n is a unit outward normal vector to the boundary. Application of the Galerkin criterion
(Hughes 1986) and Green's theorem to Eq. 3-4 leads to:



f w[(KVh)-n]dP+ f w[(KVh)-n]dd-f (Vw)-(KVh)dQ




j w-QdQ = wS- dQ
a 11 adt (3-6)








31

Because w is a member of V space, it is equal to zero at all specified head boundaries ("2);

thus, the second term in Eq. 3-6 equals zero. If it is further assumed that Q equals zero, then

Eq. 3-6 reduces to:




w[(KVh)-n]dr- (Vw)-(KVh)dQ= wS, -dQ (3-7)





For an approximate solution of h(x, y, z, t), it is assumed that:



h(e)(x,y,z,t)= Nie)(x,y,z) d,(t) N,e V heS (3-8)




i=1



where h(e) is the approximate solution for hydraulic head within element e, Ni(e) are

interpolation functions for each node within element e, p is the number of nodes within

element e, and d,(t) are the unknown values of hydraulic head for each node within element

e and at time t. At the specified head boundaries, F2, di(t) = h(x,y,z)I2. wi(e) is the weighting

function for node i, and the limits of the integration are chosen to represent the volume of the

elements. In Galerkin's method, the weighting function for each node in the element is set

equal to the interpolation function for the node, wi(e) = Ne).

When the approximate solutions (Eq. 3-8 and Eq. 3-9) are substituted into Eq. 3-7,

the resulting equation is not satisfied exactly; consequently, an error or residual occurs at

every point in the problem domain. However, if it is also assumed that the aquifer is







32
isotropic, that is to say, values of saturated hydraulic conductivity in the three coordinate

directions are constant within an element (but can vary from one element to the next), the

contribution of any element e to the residual at a node i associated with the element can be

described through the following system of equations:


R (e)(t) dite)It W
=[K(e) [~)
Re)(t) d(e)(t)


adi (e)(t)
at


ad(e)(t)
at


(3-10)


-[Fe]


where p = the number of nodes in each element (i.e., triangle element: p = 3), and Ri(e) = the

residual error at node i in element e:


(e ) =fff
V(9


aN(ie)
ax


aN(e)
ax


aN(e)
ay


aN(e)
ay


aN (e)
az

a N(e)
az


K= 0 0
0 K 0
0 0 K


and:


[F(e) ffrace

1 p J


aN (e)
ax
aN (e)
ay
aNle)
az


(3-11)


aN.(e)
ax
aN (e)
ay
aN(e)
az


(3-12)









and:



N (e)
C (e). e) ... e)] (3-13)

V(e) NP



where V(e) is the volume of element e. Eq. 3-10 represents a system of equations pertinent

to a single element, and there is an equation for every node in the element. When equations

for all elements in the mesh are combined, the global equations are obtained for the system:



ad,(t)
RI(t) dd(t) t
=[K] + [C] -[F] (3-14)


at



By setting the residuals equal to zero, we have:


ad,(t)

[K] + + [C] I =[F] (3-15)


at


Eq. 3-15 is a system of ordinary differential equations, whose solution provides values of

head d and the derivative with respect to time, ad/at, at each node in the finite-element mesh.








34

Although several methods are available for solving this system of equations, it has become

standard practice to use the finite-difference method:


ah d(t+At)-d(t)
at At


(3-16)


If a variable, co is defined, such that:


A-t
At


Then:


(3-18)


which can be extended to the vector of [d]:


(3-17)


[d(e)] = (1 -)[d(t)]+ [d(+At)]


(3-19)


If Eq. 3-19 and Eq. 3-18 are substituted into Eq. 3-15, the finite-difference formulation for

the transient saturated flow equation is as follows:



[C] + w At[K] )[d(t+At)] =( [C] (1 -w) At [K]) [d(t)]+At((1-o)[F(t)]+([F(t+At)]


(3-20)


d(e) = (1-o)d(t) +a d(t+AOt)







35
Eq. 3-20 represents a system of linear equations with unknown variables of d,(t), and

where N,(x,y,z) are the basis functions in a tri-linear hexahedral element domain. Ni(x,y,z)

does not change with time. The most common basis functions used are hat functions, which

are linear functions of space variables.

The solution of Eq. 3-20 begins by specifying the initial values of [d] (i.e., the values

of head at time to= 0). Next, the system Eq. 3-20 is solved to obtain values of {d(t+At)} at

the end of the first time step 1. This solution process is repeated for the next time step, and

so on until the specified time duration of simulation is satisfied. Depending on the choice of

w, several different subsets of the finite-difference formulations are defined.

When o=0, it is the forward-difference method is obtained:



( [C] )[d(t+At)] = ( [C] At[K])[d(t)] +At[F(t)] (3-21)



When w=1/2, the central-difference, or Crank-Nicholson, method is obtained:



([C] + 1 At[K] [d(t+At)] = [C] 1 -1 At [K] [d(t)] + (F(t)]+[F(t)])


(3-22)

When (=1, it is the backward-difference method is obatined:


( [C] + At[K] )[d(t+At)] = [C] [d(t)] +At[F(t)] (3-23)








36

3.3.3 Adaptation of the Saturated Finite-Element Model to Describe the Free Surface

Several papers have been written on the groundwater free-surface problem

(e.g.,Taylor 1967, Neuman 1970, Finn 1976, France 1976, and Huyakom 1986). Most of

these papers discussed the steady-state free-surface problem except for France (1976). This

thesis examines the transient problem starting from the basic ideas that Neuman (1970)

presented for steady flow. Basically, both France (1976) and Neuman (1970) tried to satisfy

the two boundary conditions and are consistent. This chapter is also based on the two

conditions and adapts the method they used to locate the free surface.

The primary dependent variable to be used in modeling the free surface presented

here is the hydraulic head h defined as:



h = z + (3-24)
Y


where z is the vertical coordinate; y is specific weight of a unit volume of water pg and ip

is the pressure. Along the free surface, there are two conditions which must be satisfied.

First, qr=0 (pressure atmospheric), and hence for a given position, h is prescribed on it. That

is:



h =z (3-25)

For the steady-state scenario, the second condition is simply that there is a specified

flux across the free surface boundary (i.e., a zero flux or a flow imposed equal to net

recharge). Thus, in addition to Eq. 3-24 at the free surface, there is:












ah h h -
K-i + K-j + K- k = q. (3-26)
ax 'y zz


If steady state has not been reached and the boundary is moving with a velocity

normal to its instantaneous configuration of V., then the quantity of liquid entering its unit

area is modified, and:


Ah- h- 9h-
K.,hax i+ K qj + K- k = q ,+V =q (3-27)



where i, j, k are unit orthogonal vectors; S is the specific yield coefficient (or the effective

porosity) relating the total volume of material to the quantity of liquid which can be drained

from it, and q, equals the net flux across the free surface.

The transient problem is solved by considering the solution to be a number of steady-

state problems at small intervals of time At. For a given time level, the governing equation

is solved to obtain values of nodal heads and velocities of free surface. The components of

actual velocities are determined by dividing by the specific yield coefficient. By repeating

this process for all modal points on the free surface, a new configuration and a new set of

boundary conditions are defined. The actual nodal coordinates of the free surface are found

using the Newton-Raphson iterative method. Thus, a new finite-element mesh is generated

by modifying the coordinates of the nodes in the new flow domain.








38

3.4 Validating the Code of Three-Dimensional Deformable
Finite-Element Groundwater Flow Model

3.4.1 Theis Solution

To validate the code of the three-dimensional deformable finite-element groundwater

flow model, a comparison was made of results generated by the numerical model and those

created with an analytical model describing unsteady flow to a well in a phreatic aquifer.

Figure 3-5 shows a schematic of unsteady flow to a well in a phreatic aquifer.







Ho

7/ ////777 ////



Figure 3-5. Unsteady flow to a well in a phreatic aquifer at time t.

where s = drawdown(L), Ho= initial saturated aquifer thickness (L), Q = pumping discharge

(L3/T), and r = radial distance from the well (L).

As shown in Figure 3-2, the boundary condition along the phreatic surface is

nonlinear and the shape and position of the phreatic surface are unknown. As a result, an

exact analytical solution that describes groundwater flow is not known. However, by the

employing the Dupuit assumptions and assuming recharge is zero, the unsteady flow to a

fully penetrating well in a phreatic aquifer can be described in radial coordinates by:


2h 2 1 9h2 2S 9h
ar 2 r ar K at (3-28)










Initial condition: h(to) = Ho

Boundary condition: h(t) (p) = Ho

Here, S is the phreatic aquifer specific yield (or effective porosity, ne). If drawdowns are

sufficiently small compared to the average depth of flow, (say 0.02Ho), then the approximate

solution to this problem is the Theis solution (Bear 1976):


s W(u)
41fT (3-29)

where u = r2S/(4T t); W(u) is the well function, and T=KHo= initial transmissivity.

The conditions that define the test problem used to validate the numerical

groundwater flow model include: a pumping discharge, Q, of 1.5 m3/hr; an initial saturated

phreatic aquifer thickness, Ho, of 20 m; a specific yield S (or ne), of 0.1; and a hydraulic

conductivity, K, of 15 m/day. For the three-dimensional finite-element model, the specific

storativity, Ss, is 0.0001 m-1; where S,=pg(a+n3); n equals porosity; a equals solid matrix

compressibility [LT2M-']; g is acceleration of gravity [LT2]; and 3 is fluid compressibility

[LT2M-1]. Because of problem symetry, modeling was limited to simulating hydraulic head

variations over a wedge section of the aquifer using a 21x11x3 mesh (Figure 3-6).

The boundary condition controlled by well water table, the initial water table of well

and the pumping discharge were taken as the only input. The pumping pulled water out of

the well, on the other hand, the groundwater supplied well. In each simulation time step,

based on the known water table in the beginning of time step, and an estimated water table

of well, a new discharge supplied from groundwater to well was calculated from the

groundwater modeling. Based this new calculated supplying discharge to the well, the

pumping discharge, and the circular area of the well, a more accurate well water table could









40

be estimated again. The groundwater flow model was solved again. After many iterations,

a convergent and accurate well water table could be obtained before the modeling of next

time step.
z




20




1 35
5






Figure 3-6. Finite-element mesh of the wedge used to simulate of
a well pumping in phreatic aquifer


Because the mesh is deformable, the known pressure at the free surface and flux were

used to locate the free surface. During calculation, a slice successive over-relaxation (SSOR)

method was used as an iterative approach to locate the free surface. The results illustrated

in Figure 3-7 validate the coding of the numerical model.

3.4.2 Groundwater Flow around a Circular Wetland

To validate the three-dimensional finite-element saturated groundwater flow model

again, the numerical model was used to simulate the groundwater flow around a circular

wetland and induced by a flood wave in the wetland. It is assumed that the wetland stage

fluctuations can be described by a nonperiodic Gaussian fucntion of time as Eq. 3-30:



1 2o2 (3-30)
V~lO













0.1 i T.......

0.08 .. ..
.0 .- FEM (t=5.0hr)
0.06 -
3 .ana (t=5.0hr)
0.04 .
0.04 FEM (t=2.5hr)
0.02 ana (t=2.5hr)

0
0 5 10 15 20 25 30 35
The radius to well




Figure 3-7. Comparison of the numerical and analytically simulated head
for the example problem of a transient pump test in a phreatic aquifer


Unlike the coupling boundary condition that was used in the last case, a flood

wave boundary condition used in this numerical model here. The solid line and the dash

line in Figure 3-8 are the three-dimensional finite-element model simulated results. The

triangle points and the diamond points are the groundwater analytical solutions from Eq.

2-24 obtained in Chapter 2. The circular wetland which has a radius of 50 meter was

assumed. All other system parameters are given in the Table 3-3. The numerical model

and analytical solution matched very well.

3.5 A Numerical Model Application to Simulate the Interaction
Between a Wetland and an Aquifer

The finite-element model was applied to simulate the interaction between a wetland

and a contiguous phreatic aquifer stressed by an active productive well. A horizontal,

irregular quadrilateral mesh was used to define the irregular boundaries of the wetland and

the area beyond the wetland. Figure 3-9 illustrates the numerical mesh corresponding to the













100.9
100.8 -- 3D num d-2m
100.7 ....... 3D num d=5m
3 100.6 A analytical d=2m
100.5 analytical d=5m
100.4
100.3 --
J 100.2
100.1 -s
100 77 7
99.9
2 7 12 17
Time (hr)


Figure 3-8. Comparison of the numerically and analytically simulated groundwater table.

Table 3-3. Parameters and their values


Parameter Value Parameter Value

Hoo (m) 100.0 T (m2/h) 0.635

Ar (m) 0.25 S 0.01

At (hr) 0.02 to (h) 5.0

k (m/hr) 0.00635 a 0.4



horizontal and vertical extent of the phreatic aquifer. Beyond the wetland, the vertical extent

of the mesh defines the thickness of the phreatic aquifer. Over the inundated portions of

wetland, the thickness of the mesh defines the vertical distance between the wetland bottom

and the underlying aquiclude. The wetland is shown in the middle of this area. The model

domain is 104x104x10 meters. The four lateral boundaries and bottom boundary are

assumed to be impermeable. An active production well is located within the model domain.

It was assumed that the initial groundwater table and the water stage of wetland were the

same. The deepest water depth in the wetland is 1.0 meter. Model runs were conducted to









43

simulate the exchange of water between wetland and aquifer. Such interactions between the

wetland and the aquifer were examined under transient pumping conditions and precipitation.

Since the real-time stage of the wetland is used as a boundary condition over the wetted part

of wetland, the horizontal and vertical extent of the mesh must be adjusted to trace the

vertical and horizontal movement of the transient wetland water surface.


z
Upper-extent of the phreatic aquifer defined
x vY by the bottom elevation inundated wetland

Top entent of the phreatic aquifer P Mg well









3.9999E+01 .9999E+01
S 9999E0i1 5 9999E-+0-1
7.99 999E01 "E+0I
aquiclude 9 9999E+01 9.9999E+01


Figure 3-9. Three-dimensional finite-element mesh for the simulation of wetland and aquifer

The modeled area was divided into a 21x21x10 hexahedron mesh. In the middle, a

9x9 quadrilateral mesh was used to represent the wetland area (Figure 3-9). Hydrologic

conditions that define the test problem include: a net precipitation rate of 0.02 m/hour; a

specific yield, S, of 0.075. a hydraulic conductivity, K, of 10 m/day; and a pumping

discharge, Q, of 1.554 m3/hr. Finally, data from Figures 3-10 and Figure 3-11 were used to

define the relationship between the wetland water stage, the surface water area of wetland,

and the corresponding wetland volume.

Figure 3-12 illustrates the initial condition where it was assumed that the wetland

surface-water stage and phreatic surface elevation were the same. Figure 3-13 illustrates a









44

plan view of the numerical mesh at time zero. Figure 3-14 through Figure 3-21 represent

various depictions of head contours and iso-surface plots of head predicted at hour 5 of the

simulation. These figures show that water moves laterally into the wetland and pumping well

because the same net rainfall forces the water table to increase faster than in the water stage

in the wetland. The directions of flow also depend on the active production well. If

production is large, more water flows to the well instead of flowing to the wetland.

Otherwise, some water flows towards and into the wetland as groundwater inflow. Beneath

the wetland, the hydraulic head increases with time and with increasing depth (Figure 3-17

and Figure 3-18), which shows that water flows from the aquifer into the wetland. As to the

pumping well, the contours are almost vertical and parallel to each other near the well

(Figure 3-20 and Figure 3-21). Groundwater moves towards the well in a horizontal direction

because the well fully penetrates the aquifer.




,700
n 600
... ........ .
CU 500 ....
400
.. ..........
CU 300





0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6
The water depth of wetland D (m)


Figure 3-10. Relationship between the area of wetland water surface.












500 --- ,----------j

450 ........ .. .. .

?400 -F -

E 3 5 0 -- !- --"
> 300 -

1- 250 __......_ _-

200
400 450 500 550 600 650
The surface area of wetland S (s-m)






Figure 3-11. Relationship between the deepest water depth and the volume of wetland.


N


Figure 3-12. Illustration of the initial condition where it was assumed the wetland surface
water stage and the phreatic surface elevation were the same


























5 9999E+-01





-9.9999E+01



9.9999E+01 4.9999E+01 -5.2005E-04







Figure 3-13. Plan view of the numerical mesh at time zero.


x ,


Figure 3-14. Three-dimensional finite-element mesh
and simulated head contours at hour 5.


































2 4999EOI I 2005E-04
49E 9999E 01
7 4999E,01 5 9999E+01
7 9999E+01
9 9999E+0-1 9 9999E+01





Figure 3-15. simulated head contour at hour 5.













zi


Figure 3-16. Simulated iso-surface plot of 10.6 meter head at hour 5.











48



















-52005E-04 -2005E-04
19999E4 0 1 9999E+01
3.9999E+01 39999ME 1
5.9999E+01 5.9999E+01
7.9999E+01 7.9999E+01
9.9999E+01 9.9999E+01






Figure 3-17. Three-dimensional finite-element mesh and head contour

on a cross section through wetland













0.-046 04



.5 200SE.04 95+2005E.04


Figure 3-18. Simulated head contour on the vertical cross section through the wetland





















x Y


-5 2005E-04 --5 2005E-04
9999E+01 L9999E.o
3.9999E+01 3,9999E+01
5.9999E+01 '.9999EE1I
7.9999E+01 7 9999E+01
9.9999E+01 p9f9999E401





Figure 3-19. Modeled iso-surface plot of head equal to 11.1 m, 11 m and 10.5 m


Figure 3-20. Mesh and modeled head contour with vertical cross section
through the pumping well




















A computer with a Pentium 1.4GHz processor took 4 hours to complete the 5-hour
1 901< E.01 59999 oo 9
9S99 99 9EI 9 (Si +


Figure 3-21. Head contour (m) on the vertical cross section through the pumping well

A computer with a Pentium 1.4GHz processor took 4 hours to complete the 5-hour

imulation usimulng a time step of 0.1 hours and 3 iterations per time step. To use this model

over a real system, a finer mesh would be desirable. Unlike most rectangular finite-difference

grids, the finite-element grid used here precludes the use of more efficient solution

algorithms that divide simulation domain into three individual directions. In order to use this

model to investigate wetland and aquifer interaction on a large scale, it would be desirable

to introduce a simplified mesh that would in turn be more amenable to efficient solvers. With

computer hardware developing, the simulation speed can be accelerated consequently. To

improve the simulation speed of the model from fundamental, a kind of hybrid method which

combines finite-element method and finite-difference method could be an effective choice

for a large scale practical problem.

3.6 Application of Three-Dimensional Deformable Saturated
Groundwater Flow Model to SV5 Wetland

3.6.1 Site Description-SV5

SV5 is an isolated wetland found in pine flatwoods located in southeast Florida.

It is part of the Savannas wetland system located south of the city of Jensen Beach in Martin








51

County. SV5 is located approximately two km west of the Atlantic Ocean and is adjacent to

a large municipal well field. The South Florida Water Management District (SFWMD)

monitors this wetland system. There are 13 surficial aquifer municipal water-supply wells

located from a few hundred meters to four kilometers from SV5 (Figure 3-22). The wetland

is roughly circular in shape, with a diameter of approximately 60 meters (Wise et al. 2000,

Walser 1998).




N A


SV5

A A A A AA 0* 0 0 AA

SFWMD Well \ /
A Monitor Well

o Interior Well




Figure 3-22. Well Locations (Switt et al. 1998)

The surficial aquifer that underlies the wetland is composed of three layers. The first

layer, which is located from land surface to 12-18m below the surface, consists of white, gray

and brown, largely fine- to coarse-grained quartz sand intermixed with shell beds. It receives

recharge directly from precipitation. The second layer, which is 3-6 m thick, is comprised of

fine to very fine sand with small amounts of shell and clay. The third layer, which extends

to a depth of between 37-43 m, consists of limestone intermixed with sand and shell. The

upper confining layer lies beneath the surficial aquifer. This layer consists of low-








52

permeability rocks that are primarily clastic. The Floridian aquifer system lies under the

upper confining layer. The system is quite thick and consists of limestone of variable

permeability (Fetter 1994). The climate in the area is considered humid sub-tropical, and the

area receives an average of 1.5 m of rainfall per year (Walser 1998).

3.6.2 Field Methods and Wetland Aquifer Interaction Test

The SV5 wetland aquifer interaction test was conducted (Wise et al. 2000). Eighteen

monitoring wells were installed. Wells were placed along transects running from east to west

and north to south. Six wells were placed outside the wetland on the western transect and two

on each of the other three transects. The wells were identified by the transect designation

followed by the distance location from the center of the wetland, i.e., well 37 is the well

located 37 m from the center of the wetland along the west transect. SFWMD installed two

wells and a staff gage at the site. The wells measure groundwater. The groundwater well is

screened at a depth of six meters and grouted from the surface to a depth of five meters.

The water and peat depths in the wetland were measured along four transects through

the center of the wetland at 450 angles from one another. Peat and water depth measurements

were taken at three-meter increments. The peat depth was determined by pushing a piece of

1 cm rebar into the peat at a steady rate until significant resistance was felt. This method is

most easily applied to sites such as SV5 that are free of significant tree roots. The depth of

the water table was determined by measuring the water level in the exterior monitoring wells.

Elevation data for the sites were measured utilizing a Topcon Self-Leveling Rotating laser

level. The South Florida Water Management District staff gage on the site was used as a

relative reference for the elevation data, and all elevation data were expressed in terms of

NGVD of 1929.








53

The water and peat depth information was entered into Surfer, a three-dimensional

graphical interpolation program, to develop three-dimensional contours of the wetland. A

stage -volume curve to describe the wetland was developed by utilizing Surfer to determine

the wetland volume at varying water levels. A regression curve was fitted to this data to

determine the equation that describes the stage-volume curve (Switt et al. 1998, Wise et al.

2000). This curve describes the top of the peat layer. The process was repeated to develop

a curve that describes the bottom of the peat layer.

The goal of the test was to stress hydraulically a wetland by pumping the standing

water to lower the water level a predetermined amount, which varies from wetland to

wetland, and then monitoring the response in wells and subsequent recovery of the water in

the wetland. For the wetland test to be effective, the wetland must be a discharge wetland at

the end of the pumping phase. The wetland test was conducted utilizing a 10.2 cm centrifugal

trash pump to pump the water from the wetland until the water level in the wetland dropped

30.5 cm. The water was pumped into a drainage ditch located approximately 100 m from the

wetland. This distance was chosen to ensure that the water pumped from the wetland did not

reenter the system. The discharge rate of the pump was approximately 70 m3/hr. The

discharge rate was measured by determining the amount of time it took to fill a five-gallon

plastic bucket at the end of the discharge hose. It took nine hours to lower the wetland water

level 30.5 cm.

During the testing, water levels in all of wells were manually measured every hour,

with the exception of the surface water well, which was measured every 30 minutes, utilizing

a Slope Indicator water tape. Once the pump was turned off, water levels in all of the wells

were manually measured every hour for 22 hours. Exterior and interior well water levels








54

were recorded for the wells in the wetland area.

3.6.3 Stage-Volume Relationships

Water and peat depth data were used to produce a three-dimensional contour plot of

SV5 (Figure 3-23 and Figure 3-24). The regression equation that describes the volume of

water above the peat layer is


V =3454S2-148.45s
(3-30)

where Vw is the wetland volume (m3) at s, the stage level (m) (Switt et al. 1998, Wise et al.

2000). The bottom of the peat layer is described by the follow regression equation (Figure

3.25):



V =-116s4 +607s3 +938s2+108s (3,31)



where Vp is the volume above the bottom of the peat layer (m3); and Vp-V, volumee of

saturated peat. The coefficient of determination equals unity for both of these regressions.


Figure 3-23. Top surface of peat layer of wetland































Figure 3-24. Bottom surface of peat layer







Stage-volume curve


1600
1400 -- -
1200 -
1000
800
600
400
200
0
0 0.1 0.2 0.3 0.4
Stage


Vw~S

Vp-~S


0.5 0.6 0.7


Figure 3-25. Relationship between wetland volume and water stage








56

Switt et al. (1998) used a similar formula (Eq. 3-32) to calculate the interaction

between a wetland and the underlying groundwater:



SP-K G-W) SA (3-32)




where SP= exchanging flux (m3/day), K,= hydraulic conductivity (m/day), G = groundwater

level (m), W = wetland water level (m), L = average peat thickness (m), and SA, = wetland

wetted surface area (m2). Over the de-watered area of the wetland, seepage decreases from

the above value at the edge of the water to zero at the edge of the wetland. Switt et al. (1998)

assumed a standing water head difference equal to 1/3 of the value use in Equation 3-8. Thus,

seepage for the de-watered area of the wetland was



SP=Ki( W SA, (3-33)



where SArw= wetland dry surface area (m2).

The hydraulic conductivity of the organic soil in the peat layer is critical in wetland

systems because it controls the amount of groundwater that flows from the aquifer into the

wetland or from the wetland into the aquifer (Truesdell 1995). Determining a value of

hydraulic conductivity for the soils in wetland system is particularly challenging.

3.6.4 Creation of Finite-Element Mesh

The model domain was approximately three times the size of the wetland. Because

the wetland is approximately radically symmetrical it was possible to develop a highly

refined 4.00 numerical wedge to represent the wetland aquifer system. The aquifer was








57

divided into 10 layers to obtain a 11x26x3 mesh (Figure 3-26). The peat layer was divided

into three layers producing a 3x11x3 mesh. The mesh was divided and fitted along the shape

of peat layer.
z


-Wtland Pel Layer
15



5


40 20 0 0




Figure 3-26. Finite-element mesh of wedge of wetland
and aquifer for the saturated flow model


3.6.5 Solution Procedure

To simulate the wetland-aquifer interaction observed at SV5, the model must

correctly simulate the movement of the free surface as changes occur with the wetland water

stage. Two boundary conditions must be satisfied on free surfaces; first, there can be the flux

across the surface must be represented, and, second, the pressure must be atmospheric. The

initial fixed domain solution usually will not satisfy both boundary conditions

simultaneously. Therefore, the trial free surface is adjusted and a solution obtained for the

revised fixed domain. This process is repeated until the required adjustment in the free

surface is considered negligible. Thus, the solution of the non-linear free surface problem is

achieved by a series of solutions to linear problems with prescribed boundaries. Finally,

because the wetted part of wetland may change significantly with variations of wetland water

stage, the finite-element mesh must be adjusted horizontally and vertically between time










steps.

Before the pumping test began, the SV5 wetland water surface was about 60 meters

in diameter. However, after pumping stopped it was less than 30 meters. Thus, the horizontal

inundated area changed from 2750 m2 to 750 m2. To model the influence of a changing

horizontal free surface required that the finite-element mesh be adjusted horizontally. The

zone of mesh adjustment would be within the peat layer and include an upper boundary

representing the free surface or the surface water stages and a lower boundary corresponding

to the vertical limit of the peat layer. Adjusting the mesh horizontally keeps the same number

of nodes on the boundary. Thus, the physical domain of the mesh may change horizontally

and vertically. However, the calculation domain does not. After each simulation, information

defining the location of every node is updated.

Since the wetland boundary was assumed to vary over time, the total head was taken

as a boundary condition. Thus, at any time, the wetland boundary can be viewed as an

equipotential line. Water moves vertically to and from the equipotential line represented by

the free surface boundary. The vertical hydraulic gradient at the wetland boundary times the

product of the wetland area and hydraulic conductivity equals the volume of water exchanged

between the wetland and the aquifer.

3.6.6 Model Calibration and Results

Because the aquifer underlying SV5 is composed of sand, a specific storativity of

S = 10-5 m-1 was used (Bear 1979). The saturated water content and residual water content

in the peat layer were 0s = 0.7 and Or = 0.6. In the aquifer, 6, = 0.34 and Or = 0.18, based on

Walser (1998). Water stage level observation data and groundwater monitoring data collect

during the pump test and many hours after pumping stopped were used to calibrate the










numerical model.

In Krabbenhoft's (1996) interaction model between groundwater and lakes, an

anisotropy ratio of 100 was obtained after several trial and error simulations. Also an

anisotropy ratio of horizontal to vertical hydraulic conductivity of 10 was obtained by Winter

(1978) in his numerical simulation of steady state three-dimensional groundwater flow model

near lakes.

In this finite-element model, it was assumed that the hydraulic conductivity of the

peat layer was different from the hydraulic conductivity of the underlying aquifer, and the

aquifer was also assumed anisotropic. Furthermore, it was assumed that the hydraulic

conductivities and anisotropies of both layers were spatially uniform within each layer

throughout the aquifer zone. By adjusting the anisotropy ratio (Kh/Kv) and checking for the

best agreement between modeled and observed data, it was determined that the best results

were obtained when Kh/Kv= 1.0, Kxx= Kyy = 8.7 m/day, and Kzz = 8.7 m/day in the aquifer,

while Kp = 0.2 m/day in the peat layer. This number is different from the value of 0.03m/day

obtained by Switt (1998) and 0.087m/day obtained by Wise (2000). The possible difference

may be caused by several reasons. First, the peat layer depth used in the numerical model is

not a constant. When the water stage in the wetland is lowered, the peat layer depth in the

inundated area is much greater than the average value of 0.25 m used in the analytical model.

The variation of groundwater head immediately beneath the peat layer may not be the same

value estimated using a groundwater head variation at a depth of 6 meters beneath of the

wetland. Thus, the second reason may be that the head distribution beneath the wetland

obtained from the numerical model is different from the values used by the analytical model.

Since the aquifer beneath the SV5 wetland is composed of sand, the specific








60

storativity is extremely small. The heads along the depth are very different. Figure 3-27

shows the head distribution beneath the wetland. Figure 3-28 shows the head variation with

time at different depths beneath the wetland. The piezometric head changes relatively quickly

in the deep zone, but the head does not change much near the wetland. Figure 3.29 shows the

variation of the inundated area of the wetland. The area changes from 2,600 m2 to 700 m2just

after pumping. Figure 3-30 shows the variation of the interaction discharge. The discharge

supplied from the aquifer to the wetland increases quickly and then decreases slowly after

pumping was stopped. Figure 3-31 shows excellent agreement between the simulated results

given by the saturated flow model and observed groundwater levels. Table 3-4 to Table 3-8

show the modeled data used to calculate the interaction discharges at 5 different times.

Table 3-8 shows the sensitivities of the head relative to changes in the hydraulic

conductivities in the pear layer and the aquifer and to changes in the specific storativity. The

hydraulic conductivities in both the peat layer and the aquifer play the most important part

in the change of head distribution and interaction discharge. The sensitivity of the specific

storativity is relatively small.









61



Lx


-50 0 50


Figure 3-27. Vertical cross-section view of three-dimensional finite-element
mesh and head from the saturated flow model at a time of 87 hours.


._ 15.9

E
. 15.8



15.7
0 20 40 60
Time (hour)


80 100


Figure 3-28. Piezometric head at positions beneath the wetland with time















2.5
2.5 ----------------------



S1.5.





0
S0.5 .......-- .. ........ .. 4.-^-.......--



0 20 40 60 80
Time (hour)


Figure 3-29. Variation in the inundated area of wetland during simulation





3000
< .... .. ...... ............... i .. ..
2500 ---

S2000 -

S1500

1000 -

500 I-


0 20 40
Time (hour)


60 80


Figure 3-30. Discharge between wetland and aquifer with respect to time





























0 20 40 60
Time (hour)


WETL OBS

WETL CAL

G-W OBS

G-W CAL


80 100


Figure 3-31. Comparison between calculated head and observed head





Table 3-4. Data for calculating the discharge to the wetland from the aquifer at 10 hours


Elements Area Head in Peat Head on Peat Distance of Two
of Peat Layer (m2) Layer (m) Surface (m) Nodes (m)

1 6.95828 15.7810 15.7223 0.301929

2 20.8748 15.7813 15.7223 0.305786

3 34.7914 15.7819 15.7223 0.306822

4 48.7075 15.7827 15.7223 0.303675

5 62.6240 15.7843 15.7223 0.297326

6 76.5414 15.7868 15.7223 0.289907

7 90.4590 15.7898 15.7223 0.281774

8 104.377 15.7936 15.7223 0.271725

9 118.296 15.7979 15.7223 0.260863

10 132.215 15.8028 15.7223 0.248279

Q to wetland at 10 hr 1.48466 mn/hr











Table 3-5. Data for calculating the discharge to the wetland from the aquifer at 20 hours


Elements Area Head in Peat Head on Peat Distance of Two
of Peat Layer (m2) Layer (m) Surface (m) Nodes (m)

1 8.36449 15.7953 15.7427 0.302115

2 25.0935 15.7956 15.7427 0.306344

3 41.8223 15.7960 15.7427 0.306027

4 58.5508 15.7964 15.7427 0.302183

5 75.2802 15.7982 15.7427 0.294043

6 92.0102 15.8003 15.7427 0.287151

7 108.741 15.8036 15.7427 0.275314

8 125.472 15.8070 15.7427 0.264839

9 142.204 15.8114 15.7427 0.250437

10 158.937 15.8162 15.7427 0.236900

Q to wetland at 20hr 1.65871 n3/hr


Table 3-6. Data for calculating the discharge to the wetland from the aquifer at 40 hours


Elements Area Head in Peat Head on Peat Distance of Two
of Peat Layer (m2) Layer (m) Surface (m) Nodes (m)

1 10.9288 15.8222 15.7797 0.302417

2 32.7864 15.8225 15.7797 0.307252

3 54.6436 15.8225 15.7797 0.304707

4 76.5010 15.8229 15.7797 0.298378

5 98.3597 15.8243 15.7797 0.289097

6 120.219 15.8264 15.7797 0.277635

7 142.080 15.8288 15.7797 0.266032

8 163.941 15.8321 15.7797 0.249581

9 185.804 15.8359 15.7797 0.233460

10 207.669 15.8406 15.7797 0.216548

Q to wetland at 40hr 1.86649 m'/hr












Table 3-7. Data for calculating the discharge to the wetland from the aquifer at 60 hours


Elements Area Head in Peat Head on Peat Distance of Two
of Peat Layer (m2) Layer (m) Surface (m) Nodes (m)

1 13.1219 15.8463 15.8115 0.302649

2 39.3657 15.8466 15.8115 0.307800

3 65.6089 15.8464 15.8115 0.303651

4 91.8530 15.8469 15.8115 0.295196

5 118.098 15.8480 15.8115 0.285844

6 144.345 15.8498 15.8115 0.271268

7 170.593 15.8519 15.8115 0.255422

8 196.843 15.8544 15.8115 0.239015

9 223.095 15.8579 15.8115 0.219572

10 249.349 15.8626 15.8115 0.199788

Q to wetland at 60hr 1.94458 m'/hr


Table 3-8. Data for calculating the discharge to the wetland from the aquifer at 80 hours


Elements Area Head in Peat Head on Peat Distance of Two
of Peat Layer (m2) Layer (m) Surface (m) Nodes (m)

1 15.0334 15.8681 15.8391 0.302835

2 45.1003 15.8683 15.8391 0.307522

3 75.1664 15.8680 15.8391 0.302700

4 105.234 15.8685 15.8391 0.292596

5 135.303 15.8695 15.8391 0.280887

6 165.374 15.8709 15.8391 0.266576

7 195.446 15.8727 15.8391 0.247805

8 225.521 15.8749 15.8391 0.228349

9 255.599 15.8777 15.8391 0.208494

10 285.679 15.8819 15.8391 0.186883

Q to wetland at 80hr 1.93622 m3/hr








66

Table 3-9. Sensitivities calculated by saturated groundwater flow model

Parameter (PR) Percentage varied APR I hi-h211 (m) II hi-h211/APR
Ka(m/day) 0.10 -0.85 0.9924E-02 -1.1675E-02
Ka(m/day) +0.10 +0.85 0.9055E-02 1.0653E-02
Kp(m/day) -0.10 -0.02 0.6410E-02 -0.3205
Kp(m/day) +0.10 +0.02 0.5564E-02 0.2782
So -0.50 0.000005 0.2068E-04 -4.136
So +0.50 +0.000005 0.1265E-04 2.53

3.7 Comparison with Results of Wise's (2000) and others

Parameter values obtained by Wise (2000), kpeat= 0.087 (m/day), by Walser (1998),

kand = 0.75 (m/day), kpeat = 0.03 (m/day) for SV5 were put in the three-dimension finite-

element model in this dissertation, and the simulated results were plotted on the Figure 3-32.

When the specific storage was set on 0.01 (1/m), both the simulated wetland water table and

groundwater are pretty close to the observed conditions. The simulated wetland water table

and the groundwater table corresponding to the specific storage of 0.001 (1/m) are much

different from the observed results.

There are two main differences between Walser's numerical simulation and the

models in this dissertation. First, unlike a coupling boundary condition used in this

dissertation, the fixed observed water table process over time was used in Walser's numerical

simulation. The other difference is that a fixed three-dimension mesh was used in Walser's

numerical simulation.

When kaqu at 8m/day and kpeacat 0.2m/day were used, and specific storage was varied,

the simulated results are plotted on the Figure 3-33. The simulated wetland water tables for

both specific storage values are very close to the observed wetland water table results. The









67

groundwater table for the specific storage of 0.01(1/m), however, is above the observed

results.


4.3 -
4.25
S 4.2 --. -- gwobs
S4.15 .- ----gw sim(Ss=0.01)
.) 4.1 ..- ''' gw sim(Ss=0.001)
S4.05 -- wetlobs
...-----wetilsim(Ss=0.01)
4
S- -... wet sim(Ss=0.001)
0 3.95 ,,-_ _
3.9 __............. ... -
3.85 -
0 10 20 30 40 50 60 70 80 90
Time (hr)



Figure 3-32. Comparison results using Walser's parameter values.


Figure 3-33. 3-D numerical model simulated results over specific storage

Table 3-10 shows the hydraulic conductivities in sand layer and peat layer. The

variations are rather large.








68

Table 3-10. Hydraulic Conductivity values for sand and peat layers

k in sand layer k in peat layer
(m/day) (m/day)
Switt et al.(1998) 7-12 0.03

Slug test 2

Rycroft et al. (1975) 10-7-102

Wise et al. (000) 0.083

Walser Numerical Model (1998) 0.75 0.03/0.3

Models in this dissertation 8 0.2

3.8 Summary

A three-dimensional finite-element numerical saturated groundwater flow model

using deformable mesh and real-time wetland-aquifer coupling boundary is developed in this

chapter. Considering the area of the water surface changes quickly because most wetlands

have a shallow bank slope, the model has the capability of tracking the interaction the

interaction boundary between groundwater and a wetland. The model is verified using two

analytical solutions, one is Theis solution, the other is the analytical solution found in

Chapter 2 of this dissertation, which describes groundwater flow around a circular wetland.

The model is further used to simulate the interaction problem of groundwater and the SV5

wetland. The pumping from the SV5 wetland and the recovery of the SV5 wetland surface

afterwards are closely simulated by this model, the interaction boundaries are effectively

traced by the deformable mesh.














CHAPTER 4
INVERSE MODEL OF THREE-DIMENSIONAL FINITE-ELEMENT METHOD
SATURATED GROUNDWATER MODEL IN SEARCHING
FOR THE HYDRAULIC CONDUCTIVITY OF THE PEAT LAYER


4.1 Introduction

4.1.1 Inverse Problem

During the past twenty years, mathematical modeling has been extensively used in

the study of groundwater problems. However, mathematical models require a number of

parameters to characterize the aquifer. Usually, prior knowledge of the values of these

parameters is limited, and they must be quantified by calibration which requires observations

of hydraulic heads or solute concentrations.

Given a parametric model of the physical system and values of the model parameters,

the prediction of system output (response) to any input is referred to as the direct problem,

or the forward problem (simulation). The forward problem predicts unknown system states

by solving appropriate governing equations.

Conversely, the estimation of model parameters given the parameteric system model

and an input-output relationship is known as the inverse problem (calibration). The typical

inverse problem involves estimating the values of model parameters, solving the direct

problem with this set of parameters to calculate the predicted system response, and then

adjusting model parameters until deviations between observed and model predicted

responses are in some specified sense minimized (Sun 1994).








70

Trial and error and graphical matching techniques are the most basic and traditional

methods to solve inverse problem. They are widely used because they are simple to

implement. Because the number and combinations of parameter adjustments are not

bounded, "trial and error" is rather flexible but also time-consuming. Additionally, the

solution is strongly dependent on the skill of the practitioner (Keidser 1991).

4.1.2 Parameterization

Most groundwater inverse algorithms adopt either a blocked or a geostatistical

(random field) description of spatial variability. The first approach divides the region of

interest into a number of discrete blocks or zones that are believed to correspond to distinct

geological units (Yeh 1986, Carrera and Neuman 1986a,b,c, and Cooley

1977,1979,1982,1983). Each block is characterized by a set of spatially uniform

hydrogeologic properties which are treated as parameters in an appropriate inverse problem.

Thus, the unknown distributed paramters can be represented by a number of constants. The

weakness of this method is obvious, i.e., although we have some of the geological

information is available, it is still difficult to divide the flow region into a limited number of

zones. If the zonation pattern is incorrect, the identified paramter values also will be incorrect

(Sun and Yeh 1985).

The geostatistical alternative views the properties of interest as stationary random

fields that vary relatively smoothly over space (Hoeksma and Kitanidis 1984, Dagan

1985,Yeh 1996, and Zhang 1997, Michalak 2004). In the geostatistical interpolation method

of kriging, the unknown parameter is regarded as a stochastic field described by some

statistical parameters, such as the mean, variance, and correlation distance. Then the

identified parameter vector consists of only a few statistical parameters if the aquifer is








71

basically homogeneous. The problem of overparameterization may be avoided, and the

inverse solution obtained by the maximum likelihood estimate and cokriging is generally

stable. The correctness of the inverse solution, however, depends on the correctness of some

statistical assumptions and the structure of covariance functions.

Although the two approaches are based on different parameterizations, they both treat

hydrogeologic properties as spatial functions. This suggests that it should be possible to

formulate a general inverse theory which encompasses both the blocked and geostatistical

alternatives, as well as hybird methods which lie between these extremes.

4.1.3 Parameter Identification Methods

Many methods are available for the solution of the common groundwater modeling

inverse problem whereby the conductivity or transmissivity of a porous medium are

estimated using available observations of head or pressure and other information (Distefano

1975, Chang 1976, Loaiciga 1987, Mishra 1989, Knopman 1989, Keidser 1991, Sun 1995,

Yeh 1996, Zhang 1997, Lamorey 1998, Anderman 1999, Michalak 2004). Neuman (1973)

classified inverse modeling solution techniques as either "direct" or "indirect".

The direct method or equation error method is used when head variations and

derivatives (usually estimated) are known over the entire flow region. If the measurement

and model errors are negligible, the original governing equation becomes a linear first-order

partial differential equation of the hyperbolic type in terms of the unknown parameters. With

the aid of boundary conditions and flow data, a direct solution for the unknown parameters

may be possible. In practice, observation wells are sparsely distributed in the flow region and

only a limited number of observation wells are available. To formulate the inverse problem

by the equation error criterion, missing data have to be estimated by interpolation. The








72

interpolation data contain errors resulting from the interpolation, and, thus, cause additional

errors in the results of parameter identification.

The solution approach of the indirect method to solve inverse problems is to

formulate the inverse problem in term of an optimization framework, based on the most

common output error criterion, i.e., output error least squares criterion. Using this approach,

a set of parameters is identified which best describes observed system states.

With inverse modeling, the structure and parameters of a model are adjusted

simultaneously or sequentially so as to obtain a parameter input system and output relations

that characterize observed excitation-response functions of the real system. Since the model

structure for the wetland/groundwater system is assumed to follow Eq.5-1-5-3, the problem

of determining the hydraulic conductivity of the peat layer from the observed water stage and

total head in the aquifer may be taken as a typical inverse problem. The investigation of the

inverse problem based on the three-dimensional saturated groundwater flow model and the

adjoint state method is presented in this chapter.

4.1.4 Computation of Sensitivity Coefficients

Sensitivity coefficients, i.e., the partial derivatives of head with respect to each of the

parameters, play an important role in the solution of the inverse problem. In the Gauss-

Newton algorithm, elements of the Jacobian matrix are represented by the sensitivity

coefficients, ahi/aTI, i=1,....,M, 1=1,...,L. If h is the head vector, the sensitivity coefficients

are 8h,1/8T, 1=1 ,...,L. There are three methods used in calculation of sensitivity coefficients.

There are the influence coefficient method, sensitivity equation method, and variational

method. The influence coefficient method determines the sensitivity coefficients by

perturbing each parameter once at a time.








73

In the sensitivity equation method, a set of sensitivity equations is obtained by taking

the partial derivatives with respect to each parameter in the governing equation and in the

initial and boundary conditions. The sensitivity coefficients are obtained by solving the new

governing equations.

The variational method was first used for solving the inverse problem of parameter

identification by Jacquard and Jain (1965) and then by Carter et al. (1974, 1982) associated

with finite-difference schemes. Sun and Yeh (1985) extended the method to the case of a

finite-element scheme. The sensitivity coefficients can be computed using the following

equation:



S q(xt-) Vh(xy,) ddyd (4-1)

0 0,



j=1,2,..., No i=1,2,...,N,

where Q is the exclusive subdomain of node i; V is the gradient operator; h(x,y,t) is the

solution of the governing equation; No is the number of observation wells; Nn is the total

number of nodes used in the numerical solution; and q(x,y,t) is the time derivative of q(x,y,t),

which is the solution of the set of adjoint equations.

4.2 Adjoint Problem for Three-Dimensional Finite-Element
Saturated Groundwater Model

The focus of this section is the formulation of the adjoint problem which defines

sensitivities to be used in determining the hydraulic characteristics of the peat layer and the

underlying aquifer. It is assumed that the thickness of the peat and aquifer are known. As








74

discussed in Chapter 3, three dimensional groundwater flow in an unconfined aquifer is

governed by:



V (KVh)+ Q = S,- (xyz)eO, Ot: tf (4-2)



subject to the initial and boundary conditions:


h(x,y,z,0) = fo (xy,z)EF (4-3)


hir, =f, KVh-nlr =f2 (4-4)



fo is the initial head; f, is the specific head condition; f2 is the specified flux condition; and

t f is the time duration. The boundary of volume Q shown in Figure 4-1 is denoted by r

including wetland head boundary (r2), no flow boundary (rF), and lateral head boundary

(r2).



E A B F
-- I
Peat Layer
Groundwater Aquifer
G I HI



AB, EG, FH or CD are head boundaries r2

GH or EC, DF are no flow boundaries r


Figure 4-1. Sketch showing setup of the saturated groundwater
inverse problem around a wetland








75

The hydraulic conductivities of the peat layer and the underlying aquifer are the

unknown parameters. Any variation 6K of K where K=hydraulic conductivity must cause a

variation 6h of head h, because h is dependent on K through the governing equation and

subsidiary conditions. Taking the variation of Eq. 4-2 and subsidiary conditions, the

following is obtained:



v.(6 KVh) + V *( KV6h ) = S, (xvz)Q, Ot5tf (4-5)



subject to the conditions:



6h(x,y,z)l o = 0 (xy,z)eQ (4-6)



6hIr2 = 0, [6KVh+KV8h]-nr = 0 (4-7)



Eq. 4-5-Eq. 4-7 represent the variational problem of the primary problem in Eq. 4-2. The

variational problem relates the variation 6k and 6h. Assuming h is an arbitrary function

having continuous second-order space derivatives on 0 and first-order derivatives in [0, tf],

then multiplying Eq. 4-5 by and integrating the result over time-space domains yields:





|(ST68h)I d f [S_6h + TV-(8kVh) + TV'(kV6h) dOdt = 0 (4-8a)
0 00










If we define:


V(x,y,z,tf) = 0


(x0y,z)eQ


(4-8b)


then the first term of Eq. 4-8a vanishes because of condition 4-7. Thus, Eq. 4-8b reduces to:





h a + TV-(6kVh) + V.-(kV6h) ]ddt = 0 (4-9)
0 Q


where:


YV-(8kVh)dQ =


T(8kVh)-nds -


8kVh-VWdV


'V-(kV8h)dD = Y(kVbh)nd -


- kV8h-Vfd = I
a a


8h[V-(kVT)]ad -


6hkVWTn ds


Using identities of Eq. 4-10, Eq. 4-11 and Eq. 4-12, Eq. 4-9 can be written as:


and:


(4-10)


(4-11)


(4-12)












ffS36h + h[V-(kVT)]-65kVh-VIY d dt
0 0





+ f If6kVh*nds+ fIkV6h-nds f 8hkVW-nds dt= 0 (4-13)
0 s s s


Using conditions 4-7 and defining:



Yr2 = 0, kV'-nlr, = 0


then the second term of Eq. 4-13 is equal to zero. Eq. 4-13 produces:



[s +V-(kVT) dhd} dt (V .-Vh)dkdOdt = 0 (4-14)
00 0 0

In order to identify the hydraulic conductivity k(x,y,z) that would lead to the best prediction

of observed heads, a performance criterion must be defined. A general form of the

performance function may be written as:



E(h,k)= fh,k;x,y,z,t)dDdt = 0 (4-15)
00 a










Taking the variation of E(h, k) yields:


I

E(h, k)= 0 f
0 a


The summation of Eq. 4-16 and Eq. 4-14 gives:


- s5 -V(kVW) dhdQdt +


(4-17)


S ++VT-Vh dkdQdt

0 a


If W(x,y,z,t) is selected to satisfy the following PDE:



S -t +V-(kVT)- -a=0
a T h



subject to conditions:


(4-18)


t ((x,y,z,t ) = 0,


r = 0,
T Ir, "


kVT'nlr =0
2*


then, the first integral on the right hand side of Eq. 4-17 related to unknown 5h will vanish.


Thus:


f68k+ f 6h/d dt = 0
9k 9h )


(4-16)


8E =
0 0










t

8E= ff( +VY-Vh dkdQdt

00

From Eq. 4-19, the following partial derivative can be defined:


6E
Ak


I

f f-+VT-Vh d(dt
(8k )
00


By introducing the following time transformation:


T = tf- t



the adjoint problem can be expressed as


i -_ V-(kVY) +
at


subject to:


I lr=o=0


kV'nr, =0


where T = J(x,y,z,tf -T). Solving each equation once, all components of the gradient VE(k)


can be obtained.


(4-19)


(4-20)


(4-21)


(4-22)


(4-22a)


S1|=o=0, (xy)EQ








80

4.3 Searching Objective Function

To find the best parameter value, the model structure and model parameters are

solved simultaneously or sequentially so as to make the input-output relation of the model

fit any observed excitation-response relationship of the real system. Since the model structure

of the problem in this study is determined, the problem of determining the hydraulic

conductivity of the peat layer and the underlying aquifer from the observed water stage and

total head in the aquifer may be taken as a typical inverse problem:



E(Kest) =min E(K) (4-23)




where:



L
E(k) = [ha-hobs (4-24)
1=1


The traditional method to solve this kind of problems is using the trial and error method

based on the output least squares criterion (Eq. 4-23).

4.4 Adjoint State Controlling Equations

In this study the adjoint state method was used to solve the inverse problem. From

the governing equation (Eq. 4-2) and the associated boundary conditions, the adjoint

formulation of the primary problem may be found using Green's first and second theorems

and other transforms:


S -V(k VT)+ -f =0
s rT h (4-25)










subject to:


TPL=0=0, (x,y)eO


1r2=0=0


kV'nIr, =0


where h is the total head; T = T (x,y,tf-T) is the adjoint state of h; tf is the time duration of

the simulation criterion and f(h,k;x,y,t) expresses the 'goodness' of fit between the predicted

head, h, and the observed head, hobs expressed through the ordinary least squares (OLS):


f(h,k; x,y,z,t) = 11 h(k) -hob o6(X-Xobs)8(y-y obs) (Z-Zobsb)6(t-tj)
i=1


(4-28)


where ( xobs, Yobs, Zobs) is the location of the observation well; and t, (j =1,2,---,J) are the

observation times. Integrating, Eq. 4-28 over time space domains produces:


E(h,k)= f(h,k;x,y,z,t)dQdt


(4-29)


By solving the governing Equation 4-2 and the adjoint Equation 4-25 once, we may obtain

the following functional partial derivative required to solve Equation 4-23 is obatined:


and


(4-26)


(4.27)












E f f 1Vhddt (4-30)



4.5 Searching Method and Stopping Criterion

Since the gradient VE(k) is known, many searching methods including the gradient

method, quasi-Newton method, and the Fletcher-Reeves method could be used to solve the

inverse problem 4-23 (Sun 1985). Two stopping criteria used here are:



IIK,-K-1il
and:



IIE(kn)-E(k,1) 1


Using the adjoint method and the Fletcher-Reeves algorithm to update the search sequence,

a highly efficient numerical scheme for solving the inverse problem can be formulated.

4.6 Field Application on SV5

The inverse model was used to estimate the hydraulic conductivity of the peat layer

underlying the SV5 wetland (Wise et al. 2000). Since the piezometric head at a depth of 5

meters beneath the wetland was observed, f is defined as:


J
f(h,k;x,y,z,t) = I h (k) hobI 8(x-xobs)(y-yobs)6(z-zbs)8(t-tj) (4-33)
Water-level data were measured every hour for a total duration, of 87 hours. An example


Water-level data were measured every hour for a total duration, tf, of 87 hours. An example








83

of how the hydraulic conductivity, K, is updated between iterations is as follows. Assuming

an initial estimate of the vertical hydraulic conductivity K of the peat layer was

K,=1.00m/day, then aE/9K = -0.0125 by Eq.4-29 and E(K,)= 0.48998. Because it is not

possible to match the modeling result with the observed result exactly, the final E(k) is not

zero. The sensitivity was used to decide the direction of searching. Assuming E(K) is

supposed to be the smallest value, the updated value of K=K2 is assumed to be 0.6 m/day as

generated from the formula:



E(K) = E(K) + K, (K2 K, ) (4-34)



Based on the new K = K2 = 0.6 m/day, the updated gradient 9E/dK = 0.0188, and E(K2)=

0.41467. Then a new K is assumed until E(K) is minimized. The final K is 0.2 m/day.

Figure 4-2 and Table 4-1 are the calculated values from the adjoint method model. In this

model, the observed groundwater data were used in the output least squares criterion. The

observed water stage in the wetland was not used in the model. This is the major reason that

the final result from the inverse model is a little bit different from the result obtained in

Chapter 3. The saturated flow model in Chapter 3 uses both the observation of groundwater

head beneath the wetland and the surface water stage of wetland to calibrate the model. In

the inverse model, however, it is hard to use the surface-water observation result, since it is

a boundary of the model.

Because the wetted area of surface water associated with the wetland changes with

time, the space location of each finite-element point tied to the location of the water surface

must be saved before the adjoint state equations are solved. Thus, for each assumed hydraulic









84

conductivity in the peat layer, the forward problem, and adjoint state problem must be

solved first, and then aE/ak for this assumed value of K could be found. Figure 4-3 shows

the procedure for solving an inverse problem using adjoint method.


0 0.2 0.4 0.6
K in peat layer (m/day)


0.8 1


Figure 4-2. Relation between error and peat hydraulic conductivity


Table 4-1.


Error and sensitivity based on the adjoint method model


K (m/day) E(k) aE/aK

1 1.0 10.5698 0.26965

2 0.6 8.94509 0.40555

3 0.4 7.17941 1.03330

4 0.2 3.47742 0.31495

5 0.1 4.82392 3.72980

6 0.08 5.59932 3.14520