»:«

NJ

'Mflm :Mfl

 

A .. “tr

 

 

,

,y.. w

 

 

4

 

 

 

95.1

4I.-t5]-H7hpuhu

’2

 

"J. 7 .'-

 

 

     
   
  

:41..."
.1 23:3:
3:23; 3...4

.

  

 

1'.

 

”nu“

‘wm-Im

'H
lath

 
  

   
 
  

   

‘42.!
1A A»:

     
 
  

 
 

 
 

.
O ' ’. -

9“ ,‘ n}—f.9h" >- .: - ‘

32*. Marni. #39? #35,, r‘ :14,” :_

Q93?" " 1:?" ‘

._ . - . 0

- ‘ V4

    
 
 
   
 

4
I

 

fi

 

33’”
\3I In
21:34

   
 
  

     
   
       

_ 7 - ‘ 5.44, a3 a. ”45:44.:
cm." .3 .~..‘L,.:.““"L..-.. 4.4..»
' 3 .4 4

  
 
 

l ‘
._, {Egvd‘

E—\‘

 

lb

MICHIGAN STATE UN RSI AR

n I it Im In:“ismwiliiirinmiil

3 1293 00896 7352

 

 

 

 

 

 

 

 

 

 

 

 

I!

 

This is to certify that the

dissertation entitled

Interactions Between Stream-Wetland Boundaries
and an Aquifer Stressed by a Pumping Well

presented by

Yakup Darama

has been accepted towards fulfillment
of the requirements for

PhD Civil and Environmental

degree in

 

 

Engineering

Date 7/é /¢ /
/ I

MSU is an Affirmative Action/Equal Opportunity Institution 0- 12771

 

f

.. LEERARV
Michigan Stet.)
University

m m-“

F."

 

*1

PLACE IN RETURN BOX to remove this checkout from your record.
TO AVOID FINES return on or before date due.

DATE DUE DATE DUE DATE DUE

 

P 25130?
i

i-

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

L______.

"Tfl  

MSU Is An Affirmative Action/Equal Opportunity Inaitution
cMRchna-o.‘

 

 

 

 

 

 

 

 

 

 

 

 

_____.

INTERACTIONS BETWEEN STREAM-WETLAND BOUNDARIES AND AN
BQUIEER STRESS!” BY 3 FUNDING WELL

BY

Yakup Darama

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Civil and Environmental Engineering

1991

ABSTRACT

INTERACTIONS BETWEEN STREAN-WETLAND BOUNDARIES AND AN
AQUIEER STRESSED BY A PUNPING WELL

BY
Yakup Darama

The rates of stream depletion (q,) and captured wetland
evapotranspiration (q,) caused by a pumping well located in a
connecting water table aquifer are studied.

An analytical model is developed to predict q, caused by
cyclic pumping from a deep phreatic aquifer. Equations and
dimensionless plots are developed for the volume of stream
depletion over a pumping cycle (vac) and for the maximum rate
of stream depletion at dynamic equilibrium. These plots
provide a way to determine the time required to reach
equilibrium and the error in q, produced by neglecting the
nonuniform pumping rate within a pumping cycle. Analysis
showed that vac ending at time t, is equal to the volume of
stream depletion between the start of pumping and time t by
continuous pumping and under some circumstances, approximating
the effects of cyclic pumping by a continuous pumping at the
equivalent cycle-average rate is not adequate.

A second analytical model is developed to predict q, and
q, caused by a continuously pumping well located in a shallow
phreatic aquifer. Linear variation of q, as a function of the
drawdown of the water table is incorporated to the model .

Scaled plots of q, are developed for transient and for steady

ii
state conditions. Analysis showed that as the pumping
distance increases, q, increases, while q, decreases, and the
system reaches equilibrium faster. Analysis also showed that
linear ET assumption can produce errors in q,, and q,.

A third model is developed to predict the steady state g,
by a pumping well including the effects of the nonlinear
variation of q; as a function of water table depth.
Dimensional analysis is used to determine the relationship
between the scaled q,, the scaled pumping distance (a/hl) , the
scaled hydraulic conductivity (Ks/PET) , and the scaled initial
depth to the water table ((d,/h,,) ") . Families of scaled graphs
are developed for a wide range of these parameters. Analysis
showed that (1) (CL/1),)" has an increased effect on g, as KS/PET
increases, and (2) q, decreases as Ks/PET decreases while

(d,/h,,)", and a/h, increase.

This dissertation is dedicated to my wife Hanks-e for her
continuous support and deep love.

iii

ACKNOWLEDGEMENTS

The author wishes to express his deep appreciation and
gratitude to his advisor Dr. Roger B. Wallace for his
continuous encouragement, assistance and useful suggestions
throughout this study.

He also wishes to thank for interest and cooperation of
the other committee members, Dr. David C. Wiggert of the
Department of Civil and Environmental Engineering, Dr. George
Merva of the Department of Agricultural Engineering and Dr.
Raymond Kunze of the Department of Crop and Soil Sciences.

He is particularly grateful to the Ministry of Education
of the Republic of Turkey for awarding him a National Science
Scholarship. Also, Division of Engineering Research and
Institute of Water Research are acknowledged for the partial
financial support of this research.

The author also wishes to express his appreciation to his
colleague and friend Michael D. Annable for his valuable
suggestions during this study.

Valuable suggestions of Dr. Arthur T. Corey in
development. of the last stage. of 'this dissertation. are
appreciated.

The author also wishes to express his appreciation to the

reviewers of Water Resources Research, AGU who reviewed

iv

v
earlier versions of the third and the fourth chapters and
provided helpful comments.

Most of all, the author would like to express his sincere
gratitude to his beloved parents for whom he has deep respect

and love.

TABLE OF CONTENTS

CHAPTER ............... ............. ....... ........ .... PAGE
LIST OF FIGURES ....................................... ix
LIST OF TABLES ... ......... ...... ........ ... ........ ... xii
NOMENCLATURE ............ ..... . ................. . ...... xiii
CHAPTER I. INTRODUCTION ..... . .................. . ..... l
1.1. Objective and Scope of the Study ... ...... .. 3
1.1.1. Objective ....... ............ ......... 3
1.1.2. Scope ........ .............. .......... 4
CHAPTER II. LITERATURE REVIEW .. .............. ........ 6

2.1. Introduction ............................... 6
2.2. Stream Depletion by Pumping Wells .......... 6

2.3. Steady State Evaporation from the
Water Table ................................ 11

2.4. summary ......OIOOOOIOOIOOOIO ..... 0.0.0.0... 15

CHAPTER III. STREAM DEPLETION BY CYCLIC
PUMPING OF WELLS ...... ........... . ...... 17

3.1. Introduction ............................... 17

3.2. Stream Depletion During a Single
Pumping CYCle .OOOOOOOOOOOOOOOO00.0.00...... 19

3.2.1. Stream Depletion During Pumping ...... 19
3.2.2. Stream Depletion After Pumping
Stops ......OOOOOOOOOOOOOOOOO ......... 21

3.3. Rate of Stream Depletion During
Cyclic Pumping ................ ..... .. ..... . 23

3.3.1. Conceptual Description ............... 23

3.3.2. Analytical Solution for Stream
Depletion ......OOOOO......OOOOOOOOOOO 23

vi

vii

3.4. Volumetric Stream Depletion During
One Cycle of Pumping .............. .........

3.4.1. Graphical Interpretation .............
3.4.2. Analytical Solution ........ ..........

3.5. The Approach to and the Conditions at
Dynamic Equilibrium 000......00.0.000.00....

3.5.1. Conceptual Description ...............
3.5.2. Dependence of r on the Pumping
Cycle Characteristics, Aquifer
Response Time, and Time ..............
3.5.3. Accuracy of Stream Depletion
Rates Calculated with the
Cycle-Average Pumping Rate ...... .....

3.6. Summary and Conclusions ....................
CHAPTER IV. STREAM DEPLETION BY A PUMPING WELL
INCLUDING THE EFFECTS OF LINEAR
VARIATION OF CAPTURED WETLAND
EVAPOTRMSPIRATION 00000....00. ..........
4.10 IntrOduCtion ...0.0.0000000000.0 ...... 0.00..

4.2. Analytical Model Development ...............

4.2.1. Model Assumptions ......... ...........
4.2.2. Flow Equation ........................
4.2.3. The Relationship Between do,

PET, and Soil Properties .............
4.2.4. Solutions for Drawdown,

Stream Depletion Mined Water,

and Captured ET ................ ......
4.2.5. Discussion of the Solution ...........

4.3. Solution of Partial Differential
Equation Using A 2-D Finite Element

"Odel 0.000000000000000...0000.00.000..000O.

4.3.1. Formulation and Description of
the Hypothetical Aquifer for 2-D
Finite Element Model Application .....
4.3.2. Comparison of the Analytical
Model and 2-D Finite Element
Model Results ........................

4.4. Summary and Conclusions ......... ...........

27

27

29

31

31

31

40

42

47

47

49

49
51
57

62
72

75

75

80

81

CHAPTER V.

5.1.

5020

viii
STREAM DEPLETION BY A PUMPING WELL
INCLUDING THE EFFECTS OF NONLINEAR
VARIATION OF CAPTURED WETLAND
EVAPOTRANSPIRATION ............. ..... ......
Introduction ........ ...... . ....... . ........

Model Development ..........................

5 0 2 0 1 0 nadel Assumptions 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
5 0 2 0 2 0 Flow Equation 0 0 0 0 0 0 0 0 0 0 . 0 0 0 0 0 0 0 0 0 . 0 0 0
5.2.3. The Nonlinear Relationship for q;

as a function of the depth to
water table 00000000000000.0000. ......

5.2.3. Scaling the Flow Equation ....... .....

5.3.

Numerical Solution of Flow Equation for
Hypothetical Situation ................ .....

5.3.1. The Hypothetical Aquifer .............
5.3.2. Accuracy of the Solution .............
5.3.3. Discussion of Scaled Families

5.4.
CHAPTER VI.
APPENDICES

APPENDIX A:

APPENDIX B:

APPENDIX C:

LIST OF REFERENCES 00.000.00.00 .................. 00.00

of Curves for Stream Depletion .......
Summary and Conclusions ....................

SUMMARY AND CONCLUDING REMARKS ...........

Solution of Equation 4.8 by Laplace
TranSformations 00000000000000.0000000.00.0

Listing of Computer Program for Stream
Depletion Caused by Cyclic Pumping of

wells 00.0.0.0....000000000000.000000.00000

Listing of Computer Program for Stream
Depletion by a Pumping Well Including

the Effects of Linear Variation of

Wetland ET ................................

87
87
90
90
91

96

98

105

106
108

118

124

128

136

137

141

151

160

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

LIST OF FIGURES

Cross section of the aquifer described

by analytical methOd 00.00.000.0000000.0.0..

Stream depletion during constant-steady

pumping ....................................

Stream depletion during and after

pumping ....................................

Construction of total stream depletion
with the method of superposition for

series of cyclic pumping ....... ..... .......

Example calculation of the rate of
stream depletion Equation 3.6 at

t = 450 days .............. ..... . ...........

(a) Construction of total stream

depletion for cyclic pumping when t,

and t,I are constant and (b) detail of

Figure 3.6a between t-t, and t

(a) Total stream depletion during the
cycle ending at t and (b) alternate

representation of Figure 3.7a ..............

The dependence of r, the dimensionless
volume of stream depletion during one

cycle on the scaled times a and 7 ..... .....

Comparison of stream depletion rate
caused by cyclic pumping and the equivalent
cycle-average rate of steady-continuous

pumping when 7 = 900 .......................

3.10. The accuracy of stream depletion rates
calculated with the cycle-average

pumping rate ....... ..... ..... .... ....... ..

3.11. Comparison of stream depletion rate

caused by cyclic pumping and the equivalent

cycle-average rate of steady-continuous
pumping when 7 = 0.5 and r =

ix

0095 00000000.

18

21

.22

24

26

28

30

33

39

43

44

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

4.1.

408.

4.10.

4.11.

4.12.

X

Cross section of the semi-infinite
stream-wetland aquifer system during
pumping 000.000.000.000.000.00.000.000000000

Cross section of the leaky water
table aquifer during pumping
(After Hantush, 1964a)

Cross section of the semi-infinite
stream-wetland aquifer system
before pumping 00.00000000000000000000000000

Dependence of evapotranspiration on
depth to water table when
Ks = 3000 cm/day ....... . ..............

'Linear and nonlinear dependence of

captured ET, qg, on the depth to
water table 00.000000000000000000 00000000000

Transient behaviour of stream depletion
and captured ET rates when pumping is
continuous ............................ .....

The dependence of r“ and an during
continuous pumping on e, the scaled
pumping distance, and a, the scaled

time 0.00.00000.00.00.000.0000000000000000.0

Dependence of r“, the scaled steady state
stream depletion on e, the scaled pumping
distance when ET varies as a linear function
of water table depth .......................

Comparison of steady state stream depletion
rates obtained from the nonlinear ET model
(AQUIFEM) when act-0.151 m, with the
depletion rates obtained from the linear ET
model for different values of do (Analytical
model) .....................................
(a) Plan view of the hypothetical aquifer,
(b) cross section II - II', (c) cross
section I - I' ............................

Comparison of the boundary fluxes produced
by the analytical model and by the 2-D
finite element model ......................

Water table profile computed by the
2-D finite element model along the cross
section II - II' at t = 2000 days .........

48

54

55

59

61

65

68

71

74

77

79

82

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

Figure

xi

4.13. Water table profile computed by the

2-D finite element model along the
cross section I - I’at t - 2000 days ......

Cross section of the semi-infinite
stream-wetland aquifer system described
by Equation 501 0000000000000000000000000000

Cross section of the semi-infinite
stream-wetland aquifer system
prior to pumping ...........................

(a) Dependence of steady state ET flux

on the depth to water table,

(b) Dependence of steady state captured

ET flux on the depth to water table ........

Comparison of the stream depletion rates
obtained from the linear ET, the nonlinear
ET, the step ET, and no ET models when

a = 500 m, K5431 m/d, and d,= 0.01 m

Comparison of the stream depletion rates
obtained from the linear ET, the nonlinear
ET, the step ET, and no ET models when

a = 500 m, Ks=8.64 m/d, and d,= 0.01 m

Comparison of the stream depletion rates
obtained from the linear ET, the nonlinear
ET, the step ET, and no ET models when

a = 500 m, K,=.864 m/d, and d,= 0.01 m

Variation of captured boundary fluxes
obtained from the nonlinear ET model
when a=1000 m, K,=8.64 m/d, d,=0.1 m

Dependence of the scaled steady state
stream depletion on the scaled pumping
distance and the scaled initial depth

to water table when Its/PET = 10‘

Dependence of the scaled steady state
stream depletion on the scaled pumping
distance and the scaled initial depth

to water table when Ks/PET = 103

5.10. Dependence of the scaled steady state

stream depletion on the scaled pumping
distance and the scaled initial depth
to water table when Its/PET = 102 . . . . . . . . . . .

83

92

94

97

110

111

112

116

121

122

123

Table

Table

Table

Table

Table

Table

Table

Table

LIST or IABLBB

Data for dimensionless volume of stream
depletion, 1, corresponding to selected
values of a and 7 ...... ..... ................

Range of values of hydraulic conductivity,
specific yield, aquifer thickness,
transmissivity, and 7 .......................

Values of time to practical state of
dynamic equilibrium, t" at 1=0.95 for
corresponding values of 7 ...................

Values of B corresponding to selected
values of A and 7 ............. . .............

Characteristics of the hypothetical
aquifer 000000000000000000000000000 0000000000

Characteristics of the grid configurations
used in the 2-D finite element model ........

Range of values of aquifer parameters

used in the nonlinear ET model . ...... .......
Comparison of the scaled steady state q,
obtained from configurations 2 and 3 ........

xii

34

36

38

41

78

78

107

118

(1’!be

U‘

NOMENCLATURE

perpendicular distance from the pumped well to the stream,
[L].
the first term in the summation of infinite series in
Equation 4.5, [1].
soil parameter in Gardner's unsaturated hydraulic
conductivity function.
soil parameter function in Warrick's equation, [1].
aquifer thickness, [L].
average saturated thickness, [L].
thickness of the semi permeable layer, [L].
soil parameter in Gardner's unsaturated hydraulic
conductivity function.
evapotranspiration capture coefficient, [L].
coefficient to determine the value of do in Equation 4.14,
[1] .
depth to the water table at time t, [L].
depth where evapotranspiration is less than potential
evapotranspiration, [L].
average critical depth where evapotranspiration flux is
zero, [L].

initial depth to water table, [L].

ET

N

PC

xiv
water table depth where linear ET is approximately zero,
[L].
evapotranspiration rate, [L/T].
flux per unit area that leaves or enters the aquifer,
[L/T].
piezometric head in the aquifer at t > 0, [L].
depth of water in a piezometer that is open at the base
of the aquifer, [L].
displacement head at imbibition cycle, [L].
required elevation of the piezometric surface for
evapotranspiration to be equal to the potential rate, [L].
initial head in the aquifer at t = 0, [L].
average piezometric head in the aquifer at t > 0, [L].
index that counts the number (less one) of pumping periods
each of duration t,.
index that counts the number of terms in the summation of
the terms inflWarrick's equation for steady state ET flux.
hydraulic conductivity of the semi-permeable layer, [L/T] .
unsaturated hydraulic conductivity, [L/T].
saturated hydraulic conductivity, [L/T].
pore-size distribution parameter in Anat et al.'s and
Warrick's equation.
total number of pumping periods in time t.

capillary pressure, [M/DTU.

PET potential evapotranspiration rate, [L/T].

xv

q steady state evaporation flux in Gardner's equation,
[L/T].

go stream depletion rate produced by the first pumping cycle

at t > o, [L3/T].

q, stream depletion rate produced by the second pumping
cycle at t > 0, [L3/T].

g, groundwater flow to the stream before pumping, [L3/T] .

qb, groundwater flow to the stream during pumping, [L3/T] .

qc rate of stream depletion caused by continuous steady
pumping at dynamic equilibrium where the effects of
captured evapotranspiration are not accounted, [L3/ T] .

q, evapotranspiration flux at a point calculated by Warricks
and Anat et a1 relationships, [L/T].

q; captured evapotranspiration flux at a point, [L/T],
or captured evapotranspiration rate from the entire
surface area at time t > 0, [L3/T].

q, rate of leakage from the lower confined aquifer to the
phreatic aquifer, [L3/T].

q; rate of mined water from the aquifer storage at time
t > o, [L3/T] .

q” maximum rate of stream depletion caused by cyclic pumping
at dynamic equilibrium, [L3/T].

q, rate of stream depletion at time t > 0, [L3/T].

qs horizontal boundary flux, [L/T].

Q pumping rate, [L3/T].

Q“ cycle-average pumping rate, [L3/T] .

rl

‘Qw

Vet

xvi
radial distance between the pumping well and an arbitrary
point P, [L].
radial distance between the image well and an arbitrary
point P, [L].
drawdown in the aquifer produced by pumping well, [L].
Jenkins' stream depletion function, [T].
suction head, [L]
specific yield of the aquifer.
time since pumping started, [T].
aquifer response time, [T].
time to practical condition of dynamic equilibrium, [T].
time for one complete cycle of pumping (on and off), [T].
time for period of pumping during one cycle, [T].
t-t;i,‘time parameter to simplify Equation 3.6, [T].
t-t,-t,i, time parameter to simplify Equation 3.6, [T].

aquifer transmissivity, [IR/T]. ‘

 

Jt,/[4 (t-tdi) f, scaled time to simplify Equation 3.6.

 

.[ta/[4 (t-tP-tfifi, scaled time to simplify Equation 3.6,
volume of stream depletion by continuous pumping from a
deep water table aquifer starting time at time zero, [1?].
volume of stream depletion by continuous pumping from a
deep water table aquifer starting time at time t,, Lnfi.
volume of stream depletion by cyclic pumping during the
period t-t, to t, [L3].

volume of evapotranspiration captured by continuous

pumping during time t, [IF].

Vi

WU

xvii
volume of water mined from aquifer storage by continuous
pumping during time t, [1?].
volume of stream depletion produced by continuous pumping
during time t, [1?].
well function for leaky aquifers.
rectangular coordinate system perpendicular to the river,
[L]
rectangular coordinate system along the river boundary,
[L] -
depth to water table in Gardner's equation for steady
state evaporation flux, [L]

parameter to simplfy Equations 4.5 and 5.5, [1?].

3:00k L0tt028

6()

t/t scaled time

error produced by approximating the effects of cyclic
pumping by continuous-steady pumping at the equivalent
cycle-average rate.

(0.5/I—-6I;), scaled aquifer response time with inclusion
of linear variation of wetland ET.

tP/tu scaled pumping period

unit step function.

scaled distance between the well and the stream.

1/4a, scaled time.

a/7, scaled time without inclusion of linear variation of

wetland ET.

fl.

xviii
(0.5/I346IE), scaled aquifer response time with inclusion
of linear variation of wetland ET.
td/tp, scaled cycle length.
1/[4(a-7)], scaled time without inclusion of linear
variation of wetland ET.
scaled stream depletion volume during cyclic pumping.
criteria for the practical state of dynamic equilibrium
during cyclic pumping.
scaled stream depletion with inclusion of wetland ET.
scaled mined water from the aquifer storage.
criteria for steady state condition for continuous pumping

from a shallow water table aquifer.

CHAPTER I.
INTRODUCTION

Groundwater is one of the major water resources in the
United States as well as around the world. In the past and
present, groundwater has been a major source of fresh water
for many municipalities, industries and irrigation. In the
United States alone, groundwater is estimated to supply water
for about half of the population and about one-third of all
irrigation water (Bouwer, 1978). Reports (William et al.,
1973; Henningsen, 1977; Bedell and Van Til, 1979; Fulcher et
a1. , 1986) indicate that groundwater usage for irrigation
purposes has been increasing tremendously over the past two
decades and will be increasing in the future. Since groundwa-
ter supplies are not unlimited, adequate management of
groundwater is essential to keep the groundwater resources
usable for this and future generations.

It is well known that surface water features are
hydraulically connected to the groundwater in the adjacent
aquifers. Because of the hydraulic connection between the
surface water features and the groundwater in the aquifer,
water removed from the aquifer to irrigate crops reduces the

flow available to support habitats in streams, lakes, and

1

2

wetlands. Reduction of streamflow and/or evapotranspiration
(ET) for phreatophyte use occurs because the irrigation water
is returned to the atmosphere rather than to the aquifer.
This depletion rate is often referred to as capture
(Bredehoeft et al., 1982) that can be considered as reduced
groundwater flow to surface water features or increased flow
from the surface water features toward the pumping well. This
capture often adversely effects the quantity and the quality
of water in ‘wetlands, lakes and streams during' drought
periods. The magnitude and timing of the capture caused by
pumpage are transmitted through the aquifer to the surface
water features and depend upon the aquifer properties
(storativity or specific yield and transmissivity), the
hydraulic connection between the aquifer and stream and/or
wetland, and the distance between the pumping well and the
stream. Numerous problems, past and present, have
illustrated the excessive reduction in the quantity of
streamflow and wetland evapotranspiration caused by excessive
pumping for irrigation purpose. In order to reduce the
adverse impacts of pumping on streamflow, legislation was
created in Massachusetts (Mueller and Male, 1990) to regulate
pumping rates at certain times.

To determine the reduction in wetland ET due to drawdown
of the water table by a pumping well, it is necessary to
recognize the importance of the unsaturated zone. The
unsaturated zone is near the ground surface and plays a

critical role in the determination of evaporation or the ET

3
rate from the ground surface to the atmosphere. Several
investigators (Gardner, 1958; Anat et al., 1965, Ripple et
al., 1972; Skaggs, 1978; Markar and Mein, 1987; Warrick, 1988)
indicated that the rate of actual ET is limited by the soil
water distribution in the unsaturated zone when the ET rate is
less than the potential rate. They concluded that the
distribution of soil water in the unsaturated zone strongly
depended on the unsaturated soil properties and the depth to

water table.

1.1. Objective and Scope of the Study

1.1.1. Objective

The overall objective of this research is to gain an
improved understanding of the transient and steady state
behavior of stream depletion caused by wells that pump water
(1) from a deep water table aquifer where there is no ET loss
and (2) from a shallow water table aquifer where there is ET
loss from 'the ground surface» The specific objectives
involved in this study are:

1. To develop an analytical model to predict stream depletion
caused by nonuniform cyclic pumping during the period when
flow in the aquifer establishes a condition of dynamic
equilibrium.

2. To develop an analytical solution that predicts stream

depletion and captured wetland ET that varies as a linear

4

function of water table depth in a semi-infinite wetland
area where water is pumped from a shallow phreatic aquifer.

3. To investigate the adequacy of the analytical solution for
stream depletion where the effects of linear variation of
captured wetland evapotranspiration are included.

4. To develop a method that can approximate the nonlinear
variation of wetland ET accurately.

5. To develop a method that can be used to predict steady
state stream depletion which accounts for the effects of

nonlinear variation of captured wetland ET.
1.1.2. Scope

This study is carried out in three phases. In the first
phase of the study, two analytical models are developed for
stream depletion rate and volume produced by a pumping well.
The first model describes the development of the analytical
solution for stream depletion produced by nonuniform cyclic
pumping from a deep semi-infinite phreatic aquifer. The
solutions are obtained by the application of the superposition
principle to Jenkins' solution for stream depletion produced
by steady continuous pumping of a well (Jenkins, 1968). The
second model describes the development of the analytical
solution for stream depletion and captured wetland ET caused
by continuous pumping of a well located in a shallow semi-
infinite phreatic aquifer. The effects of ET are incorporated

into the solutions as a linear function of water tablezdepth.

5

The second phase of the research consists of two stages.
The first stage is to develop’atmethod.that could describe the
nonlinear behavior of wetland ET accurately. Anat et al.’s
(1965) method to predict the steady state ET from the water
table is used to develop this method. The dimensional
analysis technique is used. to combine this method. with
Hantush’s equation for stream depletion so that the scaled
parameters that influence the stream depletion and captured
wetland ET could be determined. Dimensional analysis showed
that the scaled steady state stream depletion depends on three
independently scaled parameters: the pumping distance, the
hydraulic conductivity, and the initial depth to water table.

The second stage of the phase is a numerical study of
steady state stream depletion in a hypothetical stream-aquifer
system by using the two dimensional finite element model. The
numerical model is modified to incorporate the nonlinear
variation of wetland ET as a function of water table depth.
Anat et al.'s method is used to calculate ET losses in the
model. The results of the numerical simulations are evaluated
by the scaled parameters developed in the first stage. Three
families of dimensionless graphs are developed that can be
used to predict steady state stream depletion rates for a wide
range of aquifer properties. These graphs are valuable tools
that eliminate the necessity of the numerical model.

In the third phase, the study is. summarized and the

conclusions are presented.

CHAPTER II

LITERATURE REVIEI

2.1. Introduction

The objective of this study is to gain an improved
understanding of the transient and steady state behavior of
stream depletion and/ or captured wetland ET caused by a
pumping well located in a shallow water table aquifer. A
review of the literature for stream.depletion and evaporation
(or ET) from the water table was necessary to meet this
objective. This section is devoted to a review of the
previous research in these areas. This review will be divided
into two sections. The first section summarizes the studies
related to stream depletion produced by continuous pumping
wells“ The second section. reviews :models that. predict
evaporation or ET as a function of the depth to the water

table and unsaturated soil properties.

2.2. Stream Depletion by Pumping Ielle

The study of the interactions between surface water

features and aquifers started in the early 1940's. Theis

6

7

(1941) was the first investigator to derive an analytical
expression that determines the percentage of the pumped water
being diverted from the stream by a pumping well, provided
there is a hydraulic connection between the aquifer and the
stream. Later, Theis and Conover (1963) developed the first
dimensionless graph for calculating the percentage of stream
depletion during steady continuous pumping of a well.

Glover (1960) developed charts based on the theoretical
equations of Glover and Balmer (1954) to determine the capture
from the river which was computed in terms of the distance of
the well from the river, the properties of the aquifer
(storativity and transmissivity), and time.

Hantush (1955, 1964a, 1964b) derived the rate and volume
of stream depletion for a leaky water table aquifer
hydraulically connected to a straight stream. He developed
the stream depletion rate and volume in terms of the
complementary error function of the dimensionless pumping time
and dimensionless leakage coefficient. He developed these
equations by combining the drawdown equation for a leaky water
table aquifer (Hantush and Jacob, 1955) with the generalized
form of Darcy's equation.

Hantush (1965) also developed an analytical equation for
depletion rate and volume from streams with semipervious beds.
He replaced the resistance to flow due to the semiperviousness
of the stream bed which is determined by use of a pumping test
technique (Kazmann, 1948; Rorabaugh, 1956; and Hantush, 1959)

by an equivalent resistance. This resistance is caused by a

8

horizontal flow through a semipervious layer of insignificant
storage capacity which is lying between the aquifer and the
streambed. He introduced this resistance in his previous
analytical solutions as an additional length between the
pumping well and the stream. Hantush (1967) made further
refinements to his analysis and eliminated the assumption of
a straight river of infinite length. He developed an
analytical solution for stream depletion by pumping wells that
are located near streams that border a quadrant aquifer.

Jenkins (1968) defined a "lumped” parameter known as the
"sdf" which is a time scale that uniquely characterizes the
aquifer response curve. The "sdf" is based on the aquifer
parameters (transmissivity and specific yield) and the
perpendicular distance from the stream to the well. He
developed dimensionless plots of rate and volume of stream
depletion based on the equations developed by Glover and
Balmer (1954), and Hantush (1964a, 1965) so that calculation
of the stream depletion rate and volume as a function of sdf
and time would be easier for practical purposes. Furthermore,
Jenkins (1968) described the residual impacts after pumping
stopped and briefly discussed cyclic pumping. Jenkins (1968b)
and Moulder and Jenkins (1969) also modeled a non-straight
river of infinite length to determine stream depletion by
using electric-analog and digital computer.

Recently, Burns (1983) used Jenkins' method to describe
how to develop discrete unit response functions for unit

periods of stress and showed how the shape of the response

9
function depended on the value of sdf. He used the unit
response function and time series of aquifer stresses with a
discretized version of the convolution equation to compute the
time series of stream depletion.

None of the theoretical studies reviewed above indicate
that the effect of a discharging well on a nearby stream is
independent of the length of the reach of a stream-aquifer
system. Taylor (1971) confirmed this by using a digital
computer model of short and long reaches of a stream-aquifer
system. From the results of digital computations, he
concluded that short reaches of stream-aquifer systems gave
comparable results to long reaches with considerably less
effort and expense.

All of the analytical models briefly reviewed above rely
on the linearity of the governing flow equation of the
unconfined aquifer. The basic assumption required to have a
linear aquifer system is that the drawdowns are small compared
to the thickness of the aquifers. These models and the
analytical solutions developed in the present study considered
an average hydraulic conductivity and specific yield for the
entire aquifer. They did not consider how sensitive the
stream depletion was to stochastic variations of hydraulic
conductivity and specific yield of the aquifer.

In a recent study, Hantush and Harino (1989) developed a
stochastic water management model for the stream-aquifer
system. The model can be used to optimize pumping patterns

from the aquifer while minimizing the stream depletion rate

10

under specified system performance probability requirements.
The model considers the distribution of the hydraulic
conductivity and the specific yield of the aquifer. These
parameters were assumed to be log-normally distributed. The
stream depletion rate was introduced into the model as a
function of random aquifer parameters with a known probability
distribution. Hantush and Marino (1989) applied their model
to a hypothetical situation so that they could examine the
sensitivity of the pumping patterns and stream depletion by
varying the statistical properties of the aquifer parameters.
They concluded that the results obtained from the model
indicated that the stream depletion rate was insensitive to
the coefficient of variation of the log-specific yield and
slightly sensitive to the probability levels, whereas, the
depletion rate was highly sensitive to the hydraulic
conductivity statistics.

Analytical solutions that are briefly reviewed above for
computing streamflow depletion rates neglect conditions that
exist in typical stream aquifer systems. Some of these
conditions are (1) partial penetration of the stream to the
aquifer, (2) the presence of the semipermeable layer which
covers the streambed, (3) aquifer storage available to the
pumping well from areas beyond the stream and (4) hydraulic
disconnection between the stream and the aquifer. Recently,
Spalding and Khaleel (1991) analyzed these conditions by using
the 2-D finite element.model, AQUIFEH-l, (Townley and Wilson,

1980). From ‘the ‘results of numerical experiments they

11
concluded that the analytical solutions can be in error by 20
percent for neglecting the partial penetration of the stream,
by 45 percent for neglecting the existence of a semipermeable
layer between the stream and the aquifer, and by 21 percent
for neglecting the storage available in areas beyond the
stream. They also found that neglecting hydraulic
disconnection had only a minor effect. They indicated that the

combined effects of these errors can be as high as 85 percent.
2.3. Steady State Evaporation From the later Table

The steady state upward flow of water from the water
table through the soil profile for evaporation from the soil
surface was first studied by Moore (1939) . He introduced
water tables at the bottom of soil columns and allowed the
soil to absorb water. The soil surface was subjected to
evaporation and tensiometers were placed along the length of
the column. By using these tensiometers, he observed the
relationship between water pressure, moisture content, and
rate of water loss from the water table. As a result of his
experiments, he concluded that the finer soils supplied higher
evaporation rates even at the greater depths to the water
table.

Theoretical solutions of the flow equations for
evaporation processes from the water table were given by
several investigators, including Philip (1957), Gardner

(1958), Staley (1957), Anat et al. (1965), Ripple et a1.

12
(1972) and Warrick (1988).

Gardner (1958) analyzed upward flow rates with the water
table at any particular depth.and any suction head at the soil
surface. For homogeneous soils, he started with the
"diffusivity equation" that describes unsteady flow in
partially saturated porous media. He solved this equation for
the special case of steady one dimensional flow in a vertical

(upward) direction.

 

3
d3
2 =
{1+}? (2.1)

C

Here 2 is the depth to water table, S is the suction head, q
is the steady state evaporation flux, and K5 is the
unsaturated hydraulic conductivity (it is sometimes called
"capillary conductivity"). He defined an empirical equation
for the unsaturated hydraulic conductivity as:

K = a. (2 2)
‘ b,+S” .

 

Where a,, 1:3, and n are constants for particular soils. He
substituted Equation 2.2 into Equation 2.1 and obtained
analytical expressions for the steady state evaporation flux
when n = 1, 3/2, 2, 3, and 4. Gardner also showed maximum
values of evaporation as an inverse function of the depth to
the water table when the suction head approaches infinity.
Gardner and Fireman (1958) conducted laboratory studies

of evaporation from soil columns. The columns were saturated

13

and drained a few times in order to develop or simulate a
structure of the soil in the columns as close as possible to
soils in the field. The water was supplied from a reservoir
through porous cups at the lower end of the column. They
simulated different ranges of water table depths by varying
the vertical position of the inflow reservoir. In addition
they measured evaporation fluxes for two soils: Chino clay and
Pachappa sandy loam. They found agreement between the theory
and the experiment.

Staley (1957) conducted wind tunnel experiments to
determine how the evaporation rate from a fine sand column was
affected by the wind velocity and the depth of the water
table. He also derived functional relationships between the
evaporation rates under specific ambient conditions and the
depth of the water table. Staley’s approach was different
from that employed by Gardner. He employed Darcy's law for
the case of one dimensional flow in the vertical direction and
rearranged this equation to give an explicit expression
describing the rate of change of capillary pressure with the
depth to the water table as a function of evaporation flux and
capillary conductivity. He solved this expression for 2 which
is identical in form to Gardner's equation (Equation 2.1).

Schleusener and Corey (1959) studied the role of
hysteresis in reducing evaporation from the soils in the
presence of a water table. They concluded that the analysis
of upward movement of water from a water table in the absence

of hysteresis effects does not provide a satisfactory

14
explanation for the inverse relation that Gardner (1958)
found.

Anat et al., (1965) modified Duke's (1965) equation for
maximum rate of upward flow for evaporation from the water
table. They expanded the term in the infinite integral of
Duke's equation into a convergent series and integrated it
term by term. They obtained the first dimensionless explicit
equation for maximum upward flux from the water table in terms
of depth to water table, and the hydraulic properties of soil.

Ripple et al. (1972) also developed a set of
dimensionless equations and graphs that can be used to
estimate maximum steady state evaporation from the water
table. Their treatment of evaporation was more flexible than
previous investigators, because they considered the role of
meteorological conditions. Their procedure makes it possible
to estimate the steady state evaporation from soils including
layered soils, with a high water table. The field data
required include soil-moisture characteristic curves, water
table depth, and a record of air temperature, air humidity,
and wind velocity at one elevation (Hillel, 1980; Ripple et
al., 1972).

Recently, Warrick (1988) developed additional analytical
solutions for steady state evaporation from a shallow water
table. He incorporated Gardner’s unsaturated conductivity
equation (Equation 2.2) into Richard's equation for steady
state upward flow. His solution is valid for all n > 1,

including fractional values and bs= 0, therefore is more

 

ger-

is

 

 

the

th
th
sh
St

we

as

b:

We

th

Do

an

15
general than Gardner's solution. .However, Warrick's equation
is an implicit equation and requires iteration to determine

the value of steady state evaporation.

2.4. Summary

The literature review in the preceding sections shows
that analytical solutions used to predict stream depletion in
the case of nonuniform cyclic pumping are limited. It also
showed that there is no analytical solution available for
stream depletion which includes the effects of captured
wetland ET.

The present study will go beyond these studies in several
aspects. For example, in the first phase of the research, the
first analytical model will focus on stream depletion caused
by nonuniform cyclic pumping during the period when flow in
the aquifer establishes a condition of dynamic equilibrium.
It will show the magnitude of the error in depletion rates
estimated by assuming steady-continuous pumping at the cycle-
average rate. The second analytical model will focus on
stream depletion and/or linear variation of captured wetland
ET caused by the continuous pumping of a well from a shallow
water table aquifer.

The literature review in this chapter showed that ET from
the ground surface varies as a nonlinear function of the
position of the water table. This indicates that the second

analytical solution for stream depletion and captured wetland

16
ET has limitations. Therefore, the accuracy of the linear
variation of captured wetland ET will be investigated.

In order to obtain a more accurate solution for stream
depletion and captured wetland ET, numerical simulations will
be carried out by using a 2-D finite element model. The model
will be modified to incorporate the nonlinear variation of
wetland ET as a function of the position of the water table.
The numerical model will be applied to a hypothetical
situation to simulate steady state stream depletion and
captured wetland ET for a large range of aquifer parameters.
The present study will also make use of the dimensional
analysis technique to determine the scaled aquifer parameters.
The result of the numerical simulations will be plotted in
dimensionless graphs for a wide range of the scaled aquifer
parameters. These dimensionless graphs will be useful tools
for determining steady state stream depletion rates and
nonlinear captured wetland ET rates without employing the

numerical model.

CHAPTER III

BTREAH DEPLETION BY CYCLIC POHPIHG OP WELLS

3.1. Introduction

Groundwater withdrawals for irrigation have been
increasing over the past two decades in humid regions of the
0.8. These withdrawals occur primarily during summer months
and impact wetlands, streams, and lakes. For example, a
portion of the discharge from a pumping well located in an
aquifer that is hydraulically connected to a stream consists
of water captured from that stream (Figure 3.1) . The
groundwater withdrawal reduces the streamf low. This can occur
as a reduction of groundwater flow to the stream or an
increased flow from the stream to the aquifer. Besides
reducing the streamflow, secondary effects include the
possible reduction in quality of instream flow or the need to
improve the quality of point discharges in an attempt to
maintain instream water quality. This could be especially
pronounced during drought periods. Assuming that the demand
for consumptive water use will increase,it is important to
continue improving the methods available for estimating the

impacts of future groundwater use.

17

18

.Uocuws Honeymamsm an omnwuommp Houfisua

 

0:»

no coauomm mmouu .H.n ounces

 

19

Previous works lead to analytical solutions describing
the depletion of stream flow during a period of constant
pumping. Superposition was used to obtain similar solutions
for the period following a single episode of pumping and it
was recognized that other solutions could be obtained by this
method. In this chapter, an analytical method for stream
depletion as it approaches dynamic equilibrium in response to
cyclic pumping was derived.

This chapter focuses on stream depletion caused by
nonuniform cyclic pumping during the period when flow in the

aquifer establishes a condition of dynamic equilibrium.
3.2. Stream Depletion During a Single Pumping Cycle
3.2.1. Stream Depletion During Pumping

Previous investigators (Jenkins, 1968; Theis, 1941; Theis
and Conover, 1963; Glover and Balmer,1954; Hantush,1964a)
derived equations for rate and volume of stream depletion
caused by a constant, steady pumping of water from.an aquifer.
This stream depletion is also referred to as capture
(Bredehoeft et al.,1982). According to Jenkins (1968) the

dimensionless stream depletion rate at time t is

qr , _t_e on (3.1)
-c—) erfc («I 4t) Osts

and the dimensionless volume removed from the stream in

time t is

t
Vt t. q C. 2 '14!)
—= —+1—- ——e C 0st: (3.2)
or: (21: )0 “:7; °°
In Equations 3.1 and 3.2
aZS
t.=—1I'_z (3e3)

The term glis the response time of the aquifer which is used
in this study in preference to Jenkins' (1968) stream
depletion factor. In this stream depletion problem, t; is a
measure of the time required for the aquifer to reach a new
equilibrium after a pumping stress Q is first introduced at
the well located a distance "a" away from the stream boundary
where water is captured.

To develop Equations 3.1 and 3.2, the following
assumptions are required for the aquifer properties: (1) T is
constant with time and space, therefore, drawdown is
negligible compared to the thickness of the phreatic aquifer;
(2) there is no recharge to the aquifer under natural
conditions, so that the water-table is initially horizontal;
(3) the aquifer is isotropic, homogeneous and.semi-infinite in
areal extent; (4) the stream that forms the boundary is
straight and fully penetrates the aquifer; (5) water is
instantaneously released from aquifer storage; the well fully
penetrates the aquifer; (6) the pumping rate is steady during
any period.of pumping; (7) the temperature of the water in the
stream is constant and is the same as the temperature of the

water in the aquifer.

21

Equation 3.1 is plotted in Figure 3.2 where the solid
line shows the stream depletion rate as it approaches the
pumping rate which is shown by the broken line. As time goes
to infinity the depletion rate approaches the pumping rate.
The area under the q curve between zero and t gives the
accumulated volume of stream depletion (V) at time t and the
area between the solid line and the broken line in the same
period is the accumulated volume of water mined from aquifer

storage.

A
qt

   

0 ——__MWLML_
Wanna .. . 1
Water

   

 

o . t

Figure 3.2 Stream depletion during constant-steady pumping.
3.2.2. Stream Depletion After Pumping Stops

Stream depletion continues after pumping stops (Jenkins,

1968). If tP is the time the pump was on and t is the time

22
when conditions are assessed, the volume of water pumped from
the well during this time is Qt,. .Assuming t z t), then as t
approaches infinity, the volume of stream depletion v,
approaches QtP, if the stream is the only source of captured
water. The rate and volume of stream depletion at any time
after pumping ends can.be computed using Equations 3.1 and 3.2
and the method of superposition. According to Jenkins (1968) ,

the equation for rate of stream depletion after pumping ends

= - ‘3; _ _ t.. .. (3.4)
Q: 0(1 erf 4t) 0 [1 erf\J ——4(t-tp)] cpsts

Equation 3.4 is plotted in Figure 3.3.

is

 

 

   
   

Eq. 3.1 for time t- t,
[I]. 3.1 for time I.

~ -

‘I l. ' Q \v

., . .1"? 735.???”1} - - --

--. .1 'wés‘ "4.1., .4fi«.5~-2$; , -. , \ — .7 .
‘ ,. ...‘.. _. , . .'-" 4"“ » 31. .-. ..

0 t t

 

Figure 3.3. Stream depletion during and after pumping

23

3.3. Rate of Stream Depletion During Cyclic Pumping

3.3.1. Conceptual Description

The rate of stream depletion during cyclic pumping can be
calculated using the principle of superposition and Equation
3.4. A detailed example is given below. Suppose t, is the
length of the pumping period as shown in Figure 3.4.
Furthermore, suppose, this pattern repeats itself at intervals
in as shown in Figure 3.4. The final stream depletion rate is
obtained by adding the stream.depletion rate produced by each
individual pumping period. In the following section, an
equation is derived to calculate the rate of stream depletion

at any time after pumping starts.
3.3.2. Analytical Solution for Stream Depletion

An analytical expression for the rate of stream depletion
during cyclic pumping can be written for the conceptual
problem stated above. For constant t9 and t,, if ta, t“, U“

and U; are defined as follows:

 

 

t01=t-tdi (3.5a)

U , t. (3.5b)
01 4t°1

tuft-tp-tdi (3.50)

‘5 3.5d

U111: . ( )

41:m

24

.ocaficn 030%. no moaned .uOu cofiuamomuonsm
no vogue..— osu cue; :Owuoaep S355 150» no cofiuosuumcou .e.n ensur—

  

 

 

 

cote—ac:
:52; we 23. .23 ...

 

25

and the rate of stream depletion is

N-l
(I, = opus”) [l-erflUoiH-Mtni) [1-erf(Uni)] 0stsco (3.6)
-0

Here 6 is the unit step function which has a value of 1 when
its argument is greater than zero, and a value of zero when
its argument is less than or equal to zero. The arguments are
the terms (t-tdi) and (t-tP-tdi) . The summation over 1 starts
at zero in order to include the impact of the first pumping
cycle, and ends at (N-l) where N is the number of times the
pump is turned on in the time t.

The detailed example in Figure 3.5 is useful to show the
nature of Equation 3.6. If t, is 90 days and t, is 360 days,
then the stream depletion rate g, at t = 450 days can be
calculated as follows: At t -= 450 days, N-l = 1. When i = 0,
the term in the first square bracket of Equation .3.6
calculates the depletion rate q,' at t’= 450 days and the term
in the second square bracket of Equation 3.6 calculates the
depletion rate g" (point e in Figure 3.5) at time t"= 360 days
(t"= 450-90) which is associated with the recharge well; each
term is based on an application of Equation 3.1. The depletion
rate contribution g, caused by the first pumping cycle (point
f in Figure 3.5) can then be calculated by taking the
difference of these terms according to Equation 3.6 (q0=g0'-
go") . For i = 1 the term in the first bracket of Equation 3.6
calculates depletion rate q,’ at time t'= 90 days (t'=t-t,.i =

450-360) using Equation 3.1. Since the argument (t-tp-tdi)=0,

.msao one u u an o.n cofiuasom Sh
cofiuoaneo soouum no much on» no sawucnso~c0 cacaoxm .m.n ousofim

. a: 8m 2. e

 

 

27

the unit step function of the second term has a value of zero.
Therefore, the second term in Equation 3.6 is q,"-0. Thus, the
depletion rate contribution q, of the second pumping cycle
(point d in Figure 3.5) is according to Equation 3.6 g,=q,'-g,".
The total impact rate¢;.(point c in Figure 3.5) is calculated

by summing go and q, according to Equation 3.6.

3.4. Volumetric Stream Depletion During One Cycle of Pumping

3.4.1. Graphical Interpretation

Figure 3.6, a simplified version of Figure 3.4, shows the
total stream depletion during one cycle produced by a series
of pumping cycles. The total volume of stream depletion
between the time t, which need not be a time of maximum g"
and t-tg is the area under the curve ab in Figure 3.6b. This
volume is the same as the summation of the areas below the
individual curves cd, st, and gh in Figure 3.6b. For the
example shown in Figure 3.6 where N=3, at is the recession
limb from the first pumping cycle, cd is the recession limb
from the second pumping cycle, and gh is the rising limb of
the third pumping cycle. Since t,, Q, and t4 are constant, the
area under curve cd is the same as the area under curve c'e;
also, the area under curve gh is the same as the area under
curve 0c'. Therefore, the total volume of stream depletion

between the time t-t, and t is the same as the volume of

28

 

 

 

(b)

Figure 3.6. (a) Construction of total stream depletion for
cyclic pumping when t, and t, are constant and
(b) detail of Figure 3.6a between t - to and t.

29
stream depletion produced between zero and t by a single
period of pumping. Recognizing this fact reduces a
potentially complex mathematical problem to a very simple one,
and an analytical expression for volumetric stream depletion

over one cycle can easily be derived.
3.4.2. Analytical Solution

0n the basis of the above explanation, Figures 3.7a and
3.7b can be plotted. The shaded area in Figure 3.7a is the
same as the shaded area in Figure 3.7b. Each gives the total
volume of stream depletion between time t-tg,.and t discussed
in relation to Figure 3.6. This volume can be calculated
using Equation 3.2.

If V) is the volume represented by the area under curve
I and v2 is the volume represented by the area under curve II
(Figure 3.7b), then the shaded area represents the volume of
stream depletion during one cycle, vac, for the case of cyclic
pumping. Here

V

M: ‘,'1-v2 (3.7)

where application of Equation 3.2 gives

.. t. l t. _ '3; 75 (3.8a)
" Qt [(32+1) erfC{ TE) 4: fie t]

t. '-‘—‘ (3 8b)
v2: Qt’( (——‘,+1) )erfc{ 4:]- 4217-9 “ '

where, t'-t-tp

 

 

30

 

  

curvoll

 

Eq .3.4
I'is EqJLl

 

 

 

(b)

Figure 3.7. (a) Total stream«depletion.during the cycle ending
at t and (b) alternate representation of
Figure 3.7a

31

3.5. The Approach to and the Conditions at Dynamic Equilibrium
3.5.1. Conceptual Description

Dynamic equilibrium occurs when, during one cycle, the
volume of stream depletion vac (Equation 3.7) is equal to the
volume of water pumped from the well. Let r be the ratio of
var to the volume of water pumped during one cycle,

V

at,

Dynamic equilibrium occurs when 1 = 1. According to Equation
3.7, chc -. Qt, only as t -. on. Thus, dynamic equilibrium
occurs at t 81w. For practical purposes, dynamic equilibrium
occurs at a finite time if it is defined as occurring when

1 = r’. In this study the value of 1' was taken as 0.95. This

practical state of dynamic equilibrium occurs at t = tr

3.5.2. Dependence of 1 on the Pumping Cycle Characteristics,

Aquifer Response Time, and Time

Substituting Equations 3.8a and 3.8b into Equation 3.9

and simplifying produces the following equation for 1:

r=n((2c+1)erfcfl-ifle‘<)-(n-1)((2E+1>erchE-is/Ee") (340)
J5? J1?

32
where,

= an
1/4a

new:
II

= 1/[4(a-7)] (3-10a)
a-t/t

'Y = tP/td

Equation 3.10 is plotted in Figure 3.8 for a range of
dimensionless times a and 7. Table 3.1 contains data for the
plot. Note that 'r is independent of td/t,I even though the
stream depletion rate q, (Equation 3.6) shows such a
dependence. The range of values of dimensionless time used in
the figure were selected on the basis of the range of values
of hydraulic conductivity, specific yield, and aquifer
thicknesses given by Freeze and Cherry (1979) and summarized
in Table 3.2. A typical value for t5 was taken as 90 days.
The range of values for the distance from the pumping well to
the river was 30 to 6,000 meters.

Figure 3.8 shows that lines of 7 versus a collapse onto
a single line as r d 1’. This feature reveals a significant
aspect of the time required for the aquifer-stream system to
reach dynamic equilibrium. This may be seen most easily by
using Figure 3.8 to determine the time required to reach
equilibrium as the length of the pumping period decreases.
That is, fix'q" let r= 0.95, and use Figure 3.8 to determine
how long it takes, in each situation, for the aquifer

stream system to come into equilibrium. Table 3.3 shows the

33

.> one o moawu beacon can so oaoao 0:0 unease so«uoacoc
scouum no oas~o> mmoHsofimcoswo on» .e no oocooconoo one .m.m ousmfim

 

 

 

6

oo+ue ooo. ooe or F
btbL b F L pbbbbb hiEil'P hLPbLLb b: in —bbLLPb b L L omeo
T ...
r loto
I. Tomo
1 P II." r
I now u .6 No_. ".6 For n... .6 T.ood
I \ {oho
I To. u .o Tooo
. Amuoeunloc u .o .23
I lllllll llll < llllllllllmmumuwdl

» ILL P bbbbb p b bbTbPh b L P bbbbbb P P P

 

 

 

 

00.?

=.L

dm/OKQA

34

Table 3.1 Data for dimensionless volume of stream depletion

1, corresponding to selected values of a and 7.

 

 

 

7 10'7 10'1 1 101 102 103
a Dimensionless volume of stream depletion (r)

1.000 0.4795 0.4680 0.3068

1.500 0.5644 0.5570 0.4671

2.000 0.6171 0.6125 0.5586

2.500 0.6545 0.6514 0.6144

3.000 0.6828 0.6805 0.6531

3.500 0.7052 0.7034 0.6820

4.000 0.7235 0.7219 0.7047

4.500 0.7387 0.7374 0.7231

5.000 0.7518 0.7506 0.7385

5.100 0.7523 0.7510 0.7406 0.3018

6.000 0.7727 0.7719 0.7627 0.3685
7.000 0.7891 0.7885 0.7813 0.4436
8.000 0.8025 0.8019 0.7961 0.5250
9.000 0.8136 0.8131 0.8082 0.6034
10.000 0.8230 0.8226 0.8184 0.6902
15.000 0.8551 0.8548 0.8526 0.8170
20.000 0.8743 0.8742 0.8727 V0.8530
25.000 0.8875 0.8874 0.8863 0.8733
30.000 0.8972 0.8971 0.8964 0.8869

35.540
44.430
50.000
66.660
77.770
88.880
100.00
105.50
111.10
122.20
133.30
166.60
200.00
300.00
311.10
388.80
444.40
500.00
555.50
611.10
666.60
777.70
888.80
1000.00
1111.10

1222.20

Table 3.1 (cont'd.)

0.9055
0.9155
0.9203
0.9307
0.9358
0.9400
0.9434
0.9449
0.9463
0.9488
0.9510
0.9562
0.9600
0.9674
0.9685
0.9713
0.9732
0.9747
0.9760
0.9771
0.9781
0.9797
0.9810
0.9821
0.9830

0.9838

0.9055
0.9155
0.9203
0.9307
0.9358
0.9400
0.9434
0.9449
0.9463
0.9488
0.9510
0.9562
0.9600
0.9674
0.9685
0.9713
0.9732
0.9747
0.9760
0.9771
0.9781
0.9797
0.9810
0.9821
0.9830

0.9838

35

0.9047
0.9150
0.9199
0.9307
0.9358
0.9400
0.9434
0.9449
0.9463
0.9488
0.9510
0.9562
0.9600
0.9674
0.9685
0.9713
0.9732
0.9747
0-9760
0.9771
0.9781
0.9797
0.9810
0.9821
0.9830

0.9838

0.8985
0.9101
0.9159
0.9281
0.9339
0.9384
0.9421
0.9437
0.9452
0.9479
0.9497
0.9556
0.9596
0.9671
0.9784
0.9712
0.9731
0.9746
0.9759
0.9771
0.9781
0.9797
0.9810
0.9821
0.9830

0.9838

0.2931
0.3740
0.4298
0.5794
0.6831
0.7874
0.8920
0.9109
0.9188
0.9285
0.9349
0.9465
0.9532
0.9639
0.9649
0.9692
0.9715
0.9733
0.9748
0.9761
0.9772
0.9790
0.9805
0.9816
0.9826

0.9835

0.2917
0.3671
0.4211
0.4752
0.5294
0.5837
0.6380
0.7468
0.8557
0.9648
0.9742

0.9773

Table 3.1 (cont'd.)

36

1666.60 0.9861 0.9861 0.9861 0.9861 0.9859 0.9830
2222.20 0.9880 0.9880 0.9880 0.9880 0.9878 0.9862
3333.30 0.9902 0.9902 0.9902 0.9902 0.9901 0.9893
4000.00 0.9910 0.9910 0.9910 0.9910 0.9910 0.9904
5000.00 0.9920 0.9920 0.9920 0.9920 0.9920 0.9919

 

Table 3.2 Range of Values of Hydraulic Conductivity, Specific

Yield, Aquifer Thickness, Transmissivity and 7

 

Parameter Range

 

Hydraulic Conductivity, Ks (m/day) 0.009 - 8550
Specific Yield, S, 0.01 - 0.3
Aquifer Thickness, b (m) 6 - 100
Transmissivity, T=K,*b (mzlday) 0.06 - 8.5(10)5

7 10'7 - 106

 

 

 

 

37
results of these calculations. It can be seen from Table 3.3
that for small values of 7 in the range 10’757510‘, the time to
equilibrium is constant multiple of t,. This is revealed at
point A in Figure 3.8 which shows that 1 = 0.95 for values of
7 in the range 10'7575101 when the dimensionless time
a=1.27(10)2. Within this large range of a, t, depends on t,
only! The details of the dependence are given by the value of
a at point A which produces a simple expression for

determining time to equilibrium
t9= 127 t. for 104575101 (3.11)

On the other hand, large values of 7 imply situations
where the aquifer response time is small compared to the
length of one pumping period i.e. t, >> t“. In this situation
wells are close to the stream and stream depletion rates are
approximately the same as the pumping rate. Figure 3.9 shows
how the stream depletion rate (Equation 3.6) varies with time
when 7 is large. Here 7 = 900, t,= 0.1 days and the well is
pumped at a rate of 300 m3/day for 90 days of each year.

For large 7 the rate of stream depletion should not be
approximated by assuming pumping occurs at the equivalent
cycle-average rate. Figure 3.9 shows how poorly the capture
rate produced by pumping the same volume at the equivalent
cycle average rate approximates the capture produced by cyclic

pumping. Here t,= 360 days so, Of 300*90/360 = 75 m3/day.

38

Table 3.3. Values of time to practical state of dynamic
equilibrium t, at r = 0.95 for corresponding

values of 7

 

 

1 t.
103 9.80(10)2 t,
102 1.85(10)2 t,
101 1.27(10)2 1:.
10° 1.27(10)2 t,
10°l 1.27(10)2 t,
10'2 1.27410)2 t.
10'3 1.27(10)2 t.
104 1.27(10)2 t.
10" 1.27(10)2 t.
1045 1.27(10)2 t,
10'7 1.2“)(10)2 t,

 

 

 

5‘“.-.

.000 u > can: osqnfisa
wsossfiucOOI>Uooum mo ouch ousuo>cn0~0>o ucOHc>wsvo ecu use

osfimscm chaoao >0 pmmcco ouch scauoflacp somuuu uo somwumcaoo .m.n ousofim

 

 

 

 

 

 

 

 

 

 

 

Amxoov oEP
0.000 0.m0¢ 0.00m 0.0
fl D o
I o4
9 I. 100
3 T r
1 row.
T
,. om:
h 300 0:6an 000.86 00m
20.02160 .3 8:33 I
I 39. 0:363 son 00 one ova
. 30.. 9.652.
1 >00 00 ‘3 0.3300 xlx 00w
_ . a. . r . »»»»»»»»

 

1D
1D
0
D
D
D

 

‘1‘“““‘1‘

(flop/Cw) 3103

40

3.5.3. Accuracy of Stream Depletion Rates calculated with the

Cycle-Average Pumping Rate

The previous example illustrates two rather obvious
points: First, the stream depletion rate is essentially the
pumping rate for wells located close to the stream. Second
and more important, some measure is needed to characterize the
error in a stream depletion rate calculated by assuming
pumping at the cycle-average rate when infact the pumping is
nonuniform and cyclic. Let the error criterion be defined by

the equation

 

(3.12)

Here 3 is the error caused by approximating the effects of
cyclic pumping by assuming continuous, steady pumping at the
equivalent cycle-average rate; qc is the rate of stream
depletion caused by continuous, steady pumping; and gm is the
maximum rate of stream depletion caused by cyclic pumping.
Values of B have been determined at time t, when dynamic
equilibrium is first established. Table 3.4 shows the data
developed to define the variation in B; time to dynamic
equilibrium was calculated implicitly using Equation 3.10 for
'r = 0.95 for a range of values of 7. This calculated value of
t; was used in Equation 3.1 to calculate qt, assuming pumping
at the equivalent cycle-average rate, Q“. The maximum rate of

stream depletion at dynamic equilibrium, gm, was determined

41

Table 3.4. Values of 6 corresponding to selected values

 

 

 

of A and 7.
A 2 4 8 12
7 B
1000. 0.5000 0.7500 0.8750 0.9167
100. 0.4921 0.7460 0.8730 0.9154
33.3 0.4792 0.7360 0.8680 0.9120
10.0 0.4528 0.7164 0.8562 0.9039
5.00 0.4237 0.6938 0.8434 0.8950
3.33 - 0.3998 0.6744 0.8319 0.8871
2.00 0.3644 0.6433 0.8130 0.8739
1.33 0.3266 0.6073 0.7903 0.8578
1.00 0.2943 0.5740 0.7684 0.8420
0.700 0.2465 0.5213 0.7318 0.8154
0.500 0.1989 0.4596 0.6852 0.7805
0.333 0.1420 0.3699 0.6085 0.7206
0.200 0.0777 0.2459 0.4778 0.6095
0.133 0.0419 0.1569 0.3574 0.4942
0.100 0.0247 0.1062 0.2732 0.4037
0.075 0.0135 0.0678 0.1978 0.3140
0.050 0.0049 0.0324 0.1143 0.2012
0.033 0.0014 0.0134 0.0589 0.1156
0.020 0.0002 0.0034 0.0216 0.0493
0.014 0.0000 0.0011 0.0093 0.0242
0.010 '0.0000 0.0003 0.0037 0.0112

 

 

 

 

42
by a search.procedure using Equation 3.6 between the time t;4g
to t;. Value of 3 were then calculated using Equation 3.12
and are shown in Figure 3.10 as a function of the independent
variables 7 and A = td/tp.

Figure 3.10 can be used to determine 6. For example,
assume a well located 1272.8 meters away from the stream,
pumping 90 days in one year at the rate 300 m3/day. If the
aquifer transmissivity is 900 mz/day, and the specific yield
is 0.1, then the aquifer response time t, = 180 days, the
dimensionless pumping period 7 = 0.5 and the dimensionless
cycle length.A = 4. ‘With these parameters, the error 6 z 0.46
can be found from Figure 3.10. Figure 3.11 shows stream
depletion rates during one complete cycle at the equilibrium
time of 22933 days for the above conditions. This figure
illustrates how inadequate the cycle-average approximation.q¢
is, especially if the maximum rate of stream depletion gin
occurs during a critical time of the year.

The listing of the computer program that calculates 1,

gm, and B given in the Appendix B.

3.6. Summary and Conclusions

Stream depletion produced by nonuniform, cyclic pumping
from a well located in a hydraulically connected aquifer was
studied. The mathematical treatment required standard use of

superposition theory and the existing analytical solutions for

43

on» saws voucHsOHco nouns scauoamoo Soouum uo xoousoom one .oH.n ousofim

.ousu ocuussn oucuo>ouo~o>o

x.

 

 

 

 

 

0.00 P 0.0? O._. ...o Nolmoé
bbbbhirb [— —bbhbt b b —bbPLbbL b hb-hbbb h L rhb e
o o
w v
I 1N6
1 4o
N H K 1
I 10.0
r .V "n K .-
.. Imd
m ..u A
. moo u .... .
NP ...... K
thh P vP h _PPbPPh h b —bh;b P b LP —P:PIPL|h L —bbbPPh P P

 

 

 

0.?

Xow/3b_L = g

44

mucosaucoonaccoum no one» eocuo>ou
ocficssq ofiaomo an pomsoo ouch sofiuoH

.me.o u 5 use m.o u r con: ocfldscn
oH0>o uso~o>asvo on» use
doc soouum uo somfiuoasou .Ha.n gunman

0.0mnNN

 

 

 

Amxoov oEc.
0.0m0NN 0.0mnum 0.0mnum
uo
tn .3.
o.n om
. .... ...
ono A000 .. v
‘1

 

 

 

goo nnouu n o.
p _

 

 

(App/Cw) 9103

45
steady, continuous pumping. The solutions obtained are
appropriate only when the actual field conditions approach the
assumed conditions. Study focused on the approach to and
conditions at dynamic equilibrium.

The 1 criterion was developed for determining the time
required to achieve a practical state of dynamic equilibrium.
An analytical expression for 1 was obtained (Equation 3.10)
and plotted (Figure 3.8) for a practical range of the
independent variables. Equation 3.10 was obtained by
recognizing that the volume of stream depletion over one
cycle, from t-t, to t, is the same as the volume of stream
depletion between the start of pumping and the time t by a
single period of pumping (Figure 3.6). Equation 3.8 shows
that 1 depends on the values of a and 7 but that it is
independent of the value of t,/t, even though the stream
depletion rate shows such a dependence (Equation 3 . 6).
Analysis of the 1 relationship produced Equation 3.11 which
shows that time to equilibrium depends on t, alone for small
values of 7. Inspection of conditions when 7 is large lead to
the recognition that such wells are close enough to the stream
to reach equilibrium within the first pumping period so that
the nonuniform pumping pattern is an accurate representation
of the stream depletion pattern.

Having obtained a consistent basis for determining the
occurrence of dynamic equilibrium permitted study of stream

depletion rates at dynamic equilibrium. The B criterion

46
(Figure 3.10) was developed to characterize the error in
depletion rates estimated by assuming steady, continuous
pumping at the cycle-average rate. Figure 3.10 showed that
under some circumstances, the cycle-average approximation is

not an adequate representation of the pumping pattern.

CHAPTER IV

BTREAN DEPLETION BY A PUNPINO HELL INCLUDING THE EFFECTS OP

LINEAR VARIATION OP CAPTURE!) NETLAND EVAPOTRANSPIRATION

4.1. Introduction

A pumping well located in a wetland-stream aquifer
system, where the water table is close to the ground surface,
captures some of the well discharge from wetland evapotranspi-
ration (ET) and the rest from the bounding stream (Figure
4.1). Pumping lowers the water table in the wetland so that
water is captured by reducing ET losses. The water captured
from ET will not be captured from streamf low. Therefore,
estimates of streamflow reduction based on methods that do not
account for captured evapotranspiration are conservative and
over estimate the reduction in streamf low. On the other hand,
neglecting the fact that a lowered water table in the wetland
may reduce the flux of groundwater up to the wetland surface
may result in a failure to identify circumstances that can
threaten survival of wetland biology.

Analytical solutions describing the stream depletion
rates and volume during periods of steady continuous pumping

were summarized in the previous chapter. In addition to these

47

48

Ht “IIHH

—_‘\ 1; /v - ———-

>0

 

 

 

 

qs

 

'5 ho

 

HHH

 

 

 

 

 

 

 

777777] [’1 lili/I/I/I/I if I [III/I’llllillllilln

q,(x +Ar)

 

 

Figure 4.1. Cross section of the semi-infinite stream-wetland
aquifer system during pumping.

49

studies, Hantush (1955, 1964a) developed analytical
expressions for stream depletion, for induced leakage from the
lower confined aquifer and for mined water from the aquifer
storage by continuous pumping wells from a leaky water table
aquifer. He also developed solutions to estimate stream
depletion, induced leakage from the lower aquifer, and mined
water from the aquifer storage by gravity wells in sloping
leaky water table aquifers (Hantush, 1964b).

This chapter focuses on the development of analytical
solutions that predict captured wetland evapotranspiration and
stream depletion in a semi-infinite wetland area where water
is continuously pumped from a shallow phreatic aquifer that is
bounded by a stream. The model developed in this chapter is
based on the Hantush's solutions for pumping from a leaky
water table aquifer which is hydraulically connected to a

bounding stream.

4.2. Analytical Hodel Development

4 . 2 . 1 . Hodel Assumptions

In addition to the assumptions stated in the third
chapter, the following two groups of assumptions are required
in order to obtain the analytical model developed in this
chapter. The first group of assumptions for the aquifer
properties is: (1) there is no recharge to the aquifer from

the aquifer surface; however, there is a constant horizontal

50

boundary flux to the aquifer and this flux supplies the
groundwater for evapotranspiration.and.groundwater discharged
to the stream (see Figure 4.1); (2) the boundary flux is not
influenced by pumping; (3) water in the unsaturated zone,
except for the region close to the well, is always drained to
an equilibrium condition; that is, whenever drawdown occurs,
soil water in the unsaturated zone is redistributed
instantaneously. This assumption was justified by Skaggs and
Tang (1976, 1978) and Skaggs (1978) for shallow water table
aquifers with high hydraulic conductivity soils.

The second group of assumptions for wetland ET is as
follows: (1) Initially the water table is at the ground
surface, so that evapotranspiration from the wetland is at the
potential rate before pumping; (2) the ground surface is
shifted to the bottom of the root zone so that steady state
evaporation flux from the water table can be replaced by
steady state ET flux from the water table. This assumption
was also used by Skaggs (1978) to account for ET from a
shallow water table aquifer; (3) evapotranspiration captured
at the wetland surface varies linearly with water table depth
for drawdown sSd . Here do is the depth at which the
evapotranspiration flux becomes zero.

One of the assumptions in the previous chapter is that
the specific yield does not vary with the position of the
water table and time. This is a valid assumption for deep

water table aquifers. Duke (1972) and Gillham (1984) showed

51
that for shallow water table aquifers with fine textured
soils, the specific yield depend on the position of the water
table and time. Duke (1972) showed that for Touchet silt loam
the relative specific yield changes greatly for depths between
0 meter to 8 meter below the ground surface. This indicates
that the unsteady stream depletion and/or captured wetland ET
obtained from the analytical solution developed later in this
chapter (Equations 4.15 to 4.21) might be in error especially
when the value of hydraulic conductivity is low. Duke also
showed that for fine sand the relative specific yield values
are the same for two different water table depths when the
bubbling pressure head is 105 cm. This indicates that, as the
values. of hydraulic. conductivity' of soil increases, the
bubbling pressure head decreases; in addition a small drawdown
in the water table will result in a maximum yield. Therefore,
specific yield can be assumed constant for high conductivity
soils. This assumption was also used by Skaggs (1975, 1978),

and Bredeoheft et a1. (1982).

4.2.2. Plow Equation

The governing equation for flow in the control volume,
shown as a shaded area in Figure 4.1, can be described
approximately by the following differential equation (Hantush,

1964a):

52

3333+32F+2qu a: (EZ)§E. (4.1)
311’ By? K, T at

Here, T is the transmissivity, Sy the specific yield, K5 the
saturated hydraulic conductivity, go the captured
evapotranspiration flux, and fi(xhy,t) the depth of water in
the observation well screened throughout the aquifer during
pumping. This depth is the average hydraulic head that is
intercepted by the observation well and approximately equal to
both, h(x, y, t) , the depth of water above the base of the
aquifer during pumping and 11,01, y, t) , the depth of water in a
piezometer that is open at the base of the aquifer (Hantush,
1963, 1964a) .

As stated above, the captured evapotranspiration flux is
assumed to be a linear function of the water table depth in
the aquifer. Although Gardner (1958), Anat et al. (1965),
Ripple et al. (1972) and Warrick (1988) showed that
evapotranspiration is a highly nonlinear function of water
table depth, it was necessary to make this assumption for the
following reasons: (1) to obtain a partial differential
equation that would yield an analytical solution, (2) to gain
a better understanding of the general behavior of the system
for a wide range of aquifer parameters during pumping. It is
recognized that the linear variation of ET assumption in the
analytical solution might produce results that are in error.
However, the magnitude of the error is not known a priori.

Therefore, it is valuable to develop an analytical solution so

53
that the magnitude of error in the results caused by this
assumption could be investigated. The analytical solution can
also be used.as.a qualitative tool to show the behavior of the
system. Prickett and Lonnquist (1971) , Trescott et al.
(1976), Bredeoheft et al. (1982), and McDonald and Harbaugh
(1988) also used this assumption in a hypothetical situation
to show how phreatophyte evapotranspiration was captured by
steady continuous pumping. With this assumption, captured ET

can be written as

13310:, —Ii)
get: d
0

 

for (121—Ii) sdo (4-2)

Here, PET’is the potential evapotranspiration rate, h,is
the elevation of the piezometric surface required in.order for
the groundwater flux that supplies evapotranspiration to be
equal to the potential ET rate and do is the depth where q, is
equal to the potential ET rate. Substitution of Equation 4.2

into Equation 4.1 gives the following:

 

 

 

3252+azfi2+2pzflhrfi> : (51);}? (4.3)
T

3x2 6y2 doK, 8E

Equation 4.3 is the approximate partial differential
equation for groundwater flow in the system shown in Figure
4.1. This equation is identical in form to Hantush's (1964a)
equation for flow to a well in a horizontal, leaky water table

aquifer (Figure 4.2). Thus, his procedure can be followed to

54

K3

 

--_---- {_—-———

 

 

=‘I

 

 

 

 

 

 

 

 

 

 

— - -
IIIIIIIIIIIIIIIIIIIIIIIIIII III III II 'III III IIIIII.' ' ' ' 'I.‘ - - 'IIIIIIIIIIIII I.
IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII III III IIIIIIII‘ I I- 'IIIIII. K"I. IIIIIIIIIIIII I.
3$$$$$$$$$$IIIIIIIIIIIIIIIIIIIIIIIIII III III III I IIIIIIII. III :IIIII'I. ’05. ’nggsgggggssg 4"-
. . . . . .45555544555554545‘5555155 . 554 544 .545 ”’ 45’54454‘ 5445 .555554555'. . .5' '. . . . . . . . . . . . . . . .'
. I

 

InmerbnflmuLAqfifln

hnpumwflflebbmmn

\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
'\ ‘I :‘I I \I \I :I\ I\ I I I I I I I I I I I I I I I I I I I I I I I I I I
\ \ \ \ \ \ \ \’\ \ \’\’\’\’\’\’\’\’\’\’\’\’\’\’\’\’\’\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ S \

 

 

 

Figure 4.2. Cross section of the leaky water table aquifer
during pumping (After Hantush, 1964a).

55

solve Equation 4.3. This leads to analytical solutions that
predict captured wetland ET and captured streamflow.

Following Hantush (1964a) , suppose ho(x, y, to) is the
distribution of water table depth that would have prevailed in
the flow system if the well had not been pumped. That is, ho
is the solution to the boundary value problem shown in Figure
4.3 when the discharge from the pumping well is equal to zero
at t = to. The governing equation for flow in this system

then can be written by replacing h with ho in Equation 4.3:

 

 

 

 

 

 

 

 

6‘12: 631:: 223102 -ho) . (s )ah:
ax3+ay2+ doK: .7! .5? (4.4)
Q-0 (t-O)

PET . f f t 1 t J

‘ = f ::
Stream Aquifer K98, «— q:

«m——-Qb h° ‘*—-
l +-

Willi/I’llliillllld . (II/11111111111111)

 

 

 

Figure 4.3. Cross section of the semi-infinite stream-wetland
aquifer system before pumping.

The superposition principle may be used to combine

Equations 4.3, and 4.4. This gives the following differential

equation for determining drawdown from the initial condition

56

in the aquifer

 

2 2 2
a2(h.-fi2>+3’<ho-fiz>-2PET(ho‘m .(Ez)___a‘h°'fiz’ (4.5)
6x2 3Y2 CTR} T at

Equation 4.5 is a nonlinear partial differential equation

which Hantush (1964) linearized by assuming (ho-fi)<ho; that
is, drawdown is negligible compared to the saturated thick-
ness. Then, if the variables 2 anij are defined as follows:
z = 1203-757 (4-6)

13 = 0.5mm) (4.7)

Equations 4.6 and.4.7 may be substituted into Equation 4.5 and

the resultant differential equation manipulated to give

 

 

222.11-; . (5)12 (.8,
5x2 ay2 3 T at
Here,
Td
, . (4.9)
B VPET

 

and the initial and boundary conditions are:

Z(X.y.o) = o
Z(°o,y, t) = 0
Z(0,y,t) = 0 (4-10)
Z(x.1~,t) = 0

The point sink that represents the pumping well requires that

 

57
Lim (r 3;)“
r-0 K,
where, r = [(x-a)2+y2]"’. Equation 4.8 and its solution are

valid as long as the drawdown (hi-7i) in the aquifer is less
than or equal to do. Here B is an "ET capture coefficient"
that serves the role of Hantush's "leakage coefficient". The
linear dependence of g, on drawdown (hi-75) results in a

constant B as in the leakage problem.
4.2.3. The relationship between do, PET, and soil properties

To obtain an analytical solution to Equation 4.8, do in
Equation 4.9 must be determined in terms of PET and the
unsaturated soil hydraulic properties. Warrick's equation was
linearized to obtain a relationship between do, PET, and the
unsaturated soil hydraulic properties. According to Warrick
(1988) the steady state ET flux, q,, is related to the depth
to the water. table, d, and the unsaturated soil properties by

the following equation:
l/n .
(1:12,, (“A”) - ii. g: aj (4.11)
qe (K,+Q'.)2 -0

Here ho is the displacement head during imbibition cycle, 11 is

 

the pore-size distribution parameter (n > 1) and

58

 

 

 

. l l ” .
An «12) cosec(n)) (4 12a)
_ n
ao- n+1 (4.12b)
. 61.
(1+1) (K.+q. 6!, . (4.12s)
aj,1- 1+j+(1/n) for 120

Equation 4.11 is an implicit nonlinear relationship
between depth to water table and the ET flux. The general
nature of this equation is shown in Figure 4.4 for the
situation where Ks= 3000 cm/day and PET = 0.5 cm/day. This
figure shows that the ET flux is constant and equal to PET
when d is less than or equal to do. Here do is the critical
depth where the ET flux is equal to PET. As the water table
drops below do, the ET rate decreases rapidly at first but
goes to zero only as the depth to water table approaches m.
If the position of the water table is above do, the soil is
capable of transmitting higher flux than PET. Since the upper
limit of q, is limited by atmospheric conditions, then g,= PET
when d s d . If the position of the water table is below do,
then the ET flux is limited by the flux that can be
transmitted by the soil profile. The value of do can be
determined from the intersection of these two curves as shown
in Figure 4.4. Here do: 21.2 cm was obtained for Ko=3000
cm/day and PET-0.5 cm/day.

To find a relationship between do, PET and the soil

properties, q, in Equation 4.11 is assumed to vary linearly

59

We. I1llilzlti -31.. - in I:

4335 ooon n .0
sons mama» nouns ou cameo so scauouaamscuuoao>o mo occoocomoo .o.o oucaem

080 .o 638 12o; 8 £95
.00, .00 .00 do .00 0.0 .0... .mn .0m .0. .0

 

       

 

 

. . L . L . DOO.O
1000.0 S
1|
. 0
..oo. o o.
u (A
low—.0 S
I 1.
nu
loomo m.
soomo .5
Ioono %
noono $1
- T w
I fioooo /
I >oo\Eo 0.0 H bun. . . m.
u as}; ooon "no. 13 «.8 ".....xoml. o M
- — n — p P b — P — P — n b p F— r — - 000.0

 

 

60

between PET and zero as the position of the water table varies

between the ground surface and do.

g, = PET (l-Bq) for d < do (4.13a)
0

go = CtPET for d = do (4-13b)

The depth do is calculated by substituting Equation 4.13b into

Equation 4.11.

 

 

 

K. An m C PET K. ' (4 14)
do: hd - aj e
C*PET (K.+C*PET) 3 -o
. C PET
a _ ‘3 *1) {—K.+c pm") 31 for .20 (4.14:1)
1’1 1+j+(1/n) 3

Here C is a fractional coefficient. For example, C = 0.0001
implies that the ET flux is assumed to be zero when g, is
C*PET. Conceptually Figure 4.5 shows how steady state
captured ET flux varies linearly between zero and

(1-C)*PET as the position of the water table varies between
the ground surface and do.

In Equations 4.11 and 4.14, the displacement head and
pore-size distribution parameter can be determined by
empirical relationships given by Campbell (1985) . However
these relationships are not reliable for high conductivity
soils. Therefore, the values of the displacement head and the
pore-size distribution must be determined experimentally

(Corey, 1990).

61

 

 

 

 

 

 

 

 

 

A
/ Nonlinear variation of q,
°‘ 4
'o
a
8
o.
8 .
Liam vaxiation of q.
l -
0 dc Depth to water table, d, [L] do

Figure 4.5. Linear and nonlinear dependence of captured ET,
g“, on the depth to water table.

62

4.2.4. Solutions for Drawdown, Stream.Depletion, Hined‘later,

and Captured ET

According to Hantush (1964a) , Equation 4.8 can be solved
to determine the following relations (see Appendix A):

For drawdown,

= ___Q .5 — ' .5; (4‘15)
3 (Wm. B) W(u . B ))

where,
u=rzs/4Tt
u'2=r'2S/4Tt

r'2=(.x+a)’+y2

For the rate of stream depletion,
(I, = O-SQ(e" erfcwo) + e‘ erfc(ne)) (4.16)

where

B.=(o.s/\/E)- «Ex/E
n.=(0.5/JE)+ 673
For the volume of stream depletion,

q’ 1 ‘ -¢ .17
v: a Qt (34.23“ erfcm‘) + e err-C(33)) (4 )

63

For the rate that water is mined from aquifer storage,

(1., = 0 e"“/‘" erf(0.S/J'E) (4.18)

and for the volume of mined water,

1- qt+qfl
T (4.19)

V’ = Qt
«£2

Then, following Hantush, the captured ET flux can be

calculated as

q.c = Q - (q,+q,,) (4.20)

and the volume of captured ET,

V“ -= Qt: - (vr+vm) (4.21)

During the steady state (as t awn) equations for the rate and

volume of stream depletion can be simplified as:
q,=oe" (4.22)
= -‘ -——aBS 423
v2 06 (t 21 (. )

At steady state, Equations 4.18 and 4.19 reduce to;

q... = 0 (4.24)
v, = Elf—2 (l-e“) (4.25)

and equations for the rate and volume of captured ET become:

64

qgg 8 Q (l-e-e) (4.26)

 

v... = o (t-e“ (t-§1-£(§+B))+ 53,5) (4.27)

In Equation 4.15, W(u,r/B) is a leaky aquifer well
function which is available in tabular form in many references
(Hantush, 1956; Walton, 1962; and Hantush, 1964a). The
variable B in Equations 4.15 to 4.27 must be defined by
Equation 4.9 in order to account for captured ET as linear
function of water table depth.

Equations 4.16 to 4.27 are valid for s<do. .According to
Equation 4.2, the maximum possible ET capture at a point is
PET.) Thus, the solution for the case where q, = 0 can be
obtained by taking PET = 0 in Equation 4.9. Then 6 = 0 and
Equations 4.16 and 4.17 reduce to Glover and Balmer's (1954)
and Jenkins's (1968) equations for a nonleaky water table
aquifer where evapotranspiration processes are neglected.
Note that go in Equations 4.20 and 4.26 is the total captured
ET flux i.e. the captured flux integrated over the entire
surface area effected.

Figure 4.6 shows conceptually how stream depletion and
captured ET rates vary with time if the pumping rate Q is
continuous. The dashed line in this figure shows the stream
depletion when the effects of ET are neglected. The lower
solid curve shows the stream depletion rate when the effects
of ET capture are included, and the upper solid curve which

starts from the pumping rate at t = 0 and approaches the

'21

65

 

Cl" q,
Pumping Rate

 

’.
."

/‘/ Captured ET Volume
.//§ . <—— Captured ET Rate

 

 

 

 

 

 

 

Mined Water /' Y
/
./ i Stream depletion without ET
//_ Stream depletion rate with ET
Stream depletion volume
v v 4.
0 t

Figure 4.6. Transient behaviour of stream depletion and
captured ET rates when pumping is continuous.

66

stream depletion rate as time increases, shows the captured
ET. The area under the lower solid curve (q, curve) between
zero and time t gives the accumulated volume of stream
depletion (v,) at time t. The area between the upper solid
curve (captured ET curve) and the dashed horizontal line
(pumping rate) in.the sameoperiod.is the accumulated volume of
captured ET. The area between the two solid curves in the
same period is the accumulated volume of water mined from
aquifer storage.

The temporal response of the system exhibits the
following characteristics that can be seen in Figure 4.6.
Pumping causes a drawdown cone to develop in the aquifer.
Initially, most of the water removed is mined from aquifer
storage. As pumping continues, the drawdown cone enlarges.
This generates a horizontal flow to the well by reducing the
groundwater flow to the stream as well as reducing the
vertical flux from the water table that supplies ET at the
ground surface. The system reaches steady state when mining
ends at which time the discharge from.the well must equal the
total rate of capture from the stream and ET boundaries.

Equations 4.16 and 4.18 can be written in dimensionless
form as:

- f(e,a) (4.28a)

d
I

(I

f(e,a) (4.28b)

0'!
ll

67
_and Equation 4.22 can be written in dimensionless form as:
1“ = f (e) (4.29)
where,
74‘ qr/QI
Tm= gnu/Q.

Calculation of stream depletion rate and water removed
from the aquifer storage by Equations 4.16 and 4.18 require
evaluation of the complementary error function and the error
function. These equations were evaluated numerically (see the
listing of the computer program in Appendix B) and plotted in
Figure 4.7 in terms of the dimensionless parameters suggested
by Equations 4.28a and 4.28b.

In Figure 4.7, for the solid lines that increase as the
dimensionless time, 0; increases represent dimensionless
stream depletion, 1“; and the dashed lines that decrease with

increasing 0 represent dimensionless water mined, 1 from

a"
aquifer storage. This figure can be used to determine the
rates of stream depletion, mined water, and ET capture at any
desired time for the range of values of the dimensionless
pumping distance, 6, shown in the figure. In order to make
the value of e = 0 in Figure 4.7, either 'a' or PET has to be
zero. The 'a'- 0 case indicates that the well is in the
stream. The PET = 0 case indicates that the ET prior to

pumping is zero so there is no ET to capture. Figure 4.7

shows that as the dimensionless pumping distance increases,

68

.osemssa

msosswusoo useusu .u .oossumeu assassn uOHoOm on» use

.8 .oaeu uodmom ms» so

 

 

 

 

...s use ...C. no oosmusemou use .53. 0.30:

JSZL

 

 

 

5
00.000 3 00.00 P 00.03 00. — 0 3.0 30.0 00.0
Ill ' ll 0 e I, I 10.0
r 0|otlll~l00 IIII II II ..
I 6 z .
I. l/I/ // I II I. P O
I N fl 0 I rd I A 1 I I
l/ b I I
l Iflny 3 .Q NAU
”0 3 3 _ ou Eoobm 7
«3.333 3 Lose; uoEE III .
l I I ’\. ’3 1..” O
3 u 0 . 3 3 I
I. 3 3 3 .|¢20
- 3 3 3 .
3 3 3 .
I. 3am 3 03 1.0A0
1 30 at 6 1
l . 3% 9 8 I .
In 0 " Q ’8 a 3 m o
x 33 3 3 u
I. 33 3 3 .umuo
I ,3 3 33 I
.Ilnwavn.m .it\\ 3 3 .
l 3 3 3 fiumau
I 3 3 x/ r
mull—.0 u u 333/ II 10.0
1 o n a .I I, II T
g - —---—-b P b —:—-h b _ t—hnhpn- — P -'/. I I -I I IO.—.

 

 

IJJSJL z?

69
the portion of the well discharge captured from the stream
decreases; whereas, the portion of the well discharge captured
from the wetland ET increases. Furthermore, when the distance
between the pumping well and the stream increases, the aquifer

response time , t

4!

increases (Jenkins, 1968; Burns, 1983).
Burns (1983) showed that as t, increases, the response of the
recharge boundary (stream) to the pumping is slowed which
results in less stream depletion at any instant of time.
Moreover, as the distance between the pumping well and the
stream increases, the area from which ET may be captured
increases. This results in a larger rate of ET captured from
the wetland and a smaller rate of stream depletion.

Figure 4.7 can be used to determine the time when the
practical steady state condition is reached. Theoretically,
steady state is reached when mining stops, that is, qm= 0.
According to Equation 4.18, _-» 0 in two cases. The first
case is e-m and the second case is a-nD. In the first case, 6
can approach infinity either as adm or as Boo. For a given
stream-wetland-aquifer system the value of B is constant, thus
6 can approach infinity only as a-no. This situation can be
seen in Figure 4.7; as 6 increases, the time to steady state
decreases. As explained above, as 'a' increases, the value of
to increases. This indicates that the major portion of the
pumped water is being captured from the wetland ET long before

the effects of the drawdown reached the stream. In the second

case, a can approach infinity either as t-«n or as taoo. If

70
'a' is fixed, then a can approach infinity only as t -v 00. For
practical purposes, steady state is reached at a finite time
if it is defined as occurring when 1m=1‘m. In this study the
value of 1;, was taken as 0.01. This practical steady state

is reached at t = t . Substituting the value of 1“,“ and t, into

0

Equation 4.18 and solving for t, results in the following:

 

t = -£§[-4.605-ln.[erf 421]] (4'30)
Equation 4.30 is an implicit equation for t, which can be
solved by trial and error methods. The dimensionless form of
Equation 4.30 could have been plotted for determining the
value of tr. This would have eliminates the need to solve
Equation 4.30.

The dimensionless rate of stream depletion at steady
state is plotted in Figure 4.8 in terms of the dimensionless
pumping distance suggested by Equation 4.29. This figure
shows that for a fixed value of 'B’ the rate of stream
depletion decreases as 'a’ increases. This figure also shows
that at steady state, as 'a’ increases, captured ET increases
too, since gg= 0 at steady state. From this analysis, it can
be deduced that the coefficient 13 is the radius of the
characteristic area required to capture Q from the wetland

area where ET is at the potential rate.

.3- _300._.-¢.x-‘ .e.~.1‘.mr- '. N
I'

71

.sueeu mane» 00003 no sowuossu usesua s we meeuo> 90
secs .0 .ooseumeu assassn uedoom es» so sceueanou
seesaw ououm aumeuw ueaocm 0s» .8. no eosousoaon .0.e unseen

m\o Hm
000.03 000.3 02.0 030.0 300.0

nP-brklp P F-nl-PFF - —_-hph_ - - —--P-L - — Coo

 

l 3.0
1N0
10.0
[4.0
1.0.0
10.0

O/‘Jb = J31

I50
1.0.0
I00

I

I

thh-_ b h —-b_:-- n - —hb-- ... n _ O.F

 

 

 

 

72

4.2.5. Discussion of the Solution

The results of the analytical model developed in this
chapter are valid when the actual field conditions approach
the assumed conditions. The actual field conditions can be
described as a stream-wetland-aquifer system where the water
table is shallow. In this system, Duke (1972) and Gillham,
(1984) have shown that the specific yield depends on the
position of the water table and time. The variable specific
yield assumption is discussed in section 4.2.1 in detail,
thus, it will not be discussed in this section again.
Furthermore, Gardner (1958), Anat et al., (1965), Ripple et
al., (1972), and Warrick (1988) have shown that ET from the
wetland surface varies as a nonlinear function of the position
of the water table depth and the unsaturated soil hydraulic
properties. The effects of a linear variation of ET
assumption on the accuracy of results obtained from the
analytical model is investigated in this section.

The analytical model results at steady state might
therefore deviate from the results obtained from a model that
describes the actual system in more detail. As explained in
section 4.2.4 that as long as 'd' is less than dc, there is no
reduction in wetland ET and as 'd’ drops below dc, the ET rate
decreases as a nonlinear function of 'd'. Contrary to the
assumption of a nonlinear variation of ET, Figure 4.5 shows

that even small change in 'd' produces a reduction in the

73

wetland ET rate in the case of a linear ET assumption.
Therefore, a linear dependence of captured ET on the depth to
water table can produce a substantial amount of error in
predicting the rate of stream depletion and capture wetland ET
at steady state. Figure 4.9 shows the magnitude of error in
steady state stream depletion rates obtained from the
analytical model as compared to results obtained from a model
that employee the nonlinear dependence of ET on the depth to
water table. The latter results are obtained with a numerical
solution to the nonlinear form of Equation 4.1 (Equation 5. 1).
Figure 4.9 shows the results of the numerical model and the
analytical model at steady state for pumping distance varying
from O m to 1600 m. This figure shows the comparison for
Ito-43.2 m/day. The comparison is also made for the lower
range of values of K5 and presented in the next chapter.

Figure 4.9 also shows the steady state stream depletion
rates obtained from the analytical model for three different
values of do. These values of do are computed with Equation
4.14 by adjusting the value of C. It is apparent from the
figure that when the value of do is small, the results
obtained from the analytical model are in substantial error.
The amount of error could be decreased if the value of C in
Equation 4.14 could be altered in order to increase the value
of do. It was deduced from this analysis that a theoretical
basis is needed to justify the above method that determines

the value of do that reduces the error. However, it is not

74

.3000! acofitnacss on. no accuser
ucmuouuau now H0603 an ucoswa on» loan nonwcuno nouns sawuoHaop
on» 5? .a finance :25 Axum—585 deco... as 53:35: on»

souu cocficuno nouns scauoucop scouum ououu upcoum no somwucnsoo .m.v ousmwm

AEV Eombm E0: 8:329 9:955

.000F .00¢P .00.NF .omop .oww 0.00 .om: .0_0N .0

 

    

 

 

 

o.o
. $.08
-98..
.. ..odoo
. “E wand no.3 E .5251 .
_ . no
I at mum . .. 3 .5. BEST... [Odom
as 32.. "03 E .9251
E 53 nos 5 BuszozI
L _ b _ . _ r p F _ p _ . b P O.OOOP

(flop/ow) Jb aims Kpoeis

75
possible to define such a method at this stage of the study.
Therefore, the analytical model developed in this chapter is
restricted significantly as a quantitative tool, but it is

still valuable as a qualitative tool.

4.3. Solution of Partial Differential Equation Using A 2-D

Finite Eleaent.lodel

A 2-D horizontal finite element model, AQUIFEM-l (Townley
and Wilson, 1980), was modified to incorporate linear and
nonlinear variation of ET as a function of water table depth.
The objectives of this modification were (1) to compare the
analytical model results with a numerical solution of the
equation and (2) to prepare the numerical model for further
investigations“ The governing flow equation of the numerical
model is a second order nonlinear parabolic partial
differential equation and it is solved iteratively by
employing the finite element method, based on the Galerkin
technique with linear finite elements (Townley and Wilson,

1980).

4.3.1. Formulation and Description of the Hypothetical Aquifer

for 2-D Finite Element Model Application

In order to compare the finite element results with the

analytical solution, the aquifer is located adjacent to a

76

flowing surface water feature (Figure 4.10a). The aquifer is
4400 meters wide and 6300 meters long. Initially, the water
table in the aquifer is horizontal and very close to the soil
surface so that ET is at the potential rate. The aquifer is
bounded by a fully-penetrating stream on the southern side and
bounded by no-flow boundaries on the other sides“ The well is
located 800 meters away from the stream as shown in Figure
4.10a and pumps water from the aquifer at a constant rate of
3000 m3/day. Table 4.1 shows the characteristics of the
hypothetical aquifer. The aquifer domain shown in Figure
4.10a is divided into triangular finite elements for the
application of the numerical model. Three grid configurations
were designed in order to determine the effects of spatial
discretization of the aquifer domain on the accuracy of the
numerical model results. Table 4.2 shows the total number of
elements, the total number of nodes and the dimensions of the
grid in the vicinity of the well and the stream in each
configuration.

The finite element and analytical models were applied to
the aquifer described above to predict the stream depletion
and captured wetland ET. .Figure 4.11 shows the results of the

simulation from both models for comparison.

77

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

.1
“4400111. . NoFlow
‘ I Boundary
6300m E
= Q
11.... umyfigT ----------------- 11’
Si
'1.
(a
Q
11 1 II’
1' \1 _____
12. _ K‘ "o
Sy
-‘-\\\\ \\\\, \ \\
(b)
ET=PET ET<PET
e ¢ 1‘ Q ‘
1 IT TTTMHTT I.
‘ \L/
K
h s: q1' Stream
a
(C)

Figure 4.10. (a) Plan view of the hypothetical aquifer,
(b) cross section II-II', (c) cross section I-I' .

78

Table 4.2. Characteristics of the hypothetical aquifer

 

(1)

 

Aquifer thickness, Do (meter)

Initial saturated thickness, ho (meter)
Saturated hydraulic conductivity, Ks (meter/day)
Specific yield, Sy

Maximum ET depth, do (meter)

Potential ET rate, (meter/day)

 

 

 

Table 4.2. Characteristics of the grid configurations used in

the finite element model

 

of nodes of elements the smallest elem.

(Base width X

Total number Total number Dimensions of thel

 

 

 

 

 

height)
Configuration 1 416 744 200 X 200
Configuration 2 1081 2024 100 X 100
Configuration 3 2808 5396 50 X 50

 

 

79

.Hooo- acoaoso unease au~ on» an can Hooc-
Hcoqu>~cso ecu >0 teaspoon noxcnu ahevsson ecu uo confines-00 .-.v ousonm

 

 

 

 

9.30:... 000 u 0

C83 mEF
.oow .00. .03 .0: .03 .00 F .00 .00 .0... .om .0
F . _ . _ t _ . _ . b . _ . _ p _ . _ . .O
B 5;, 8:038 88.5 W
1-1.1! n .8003 ...
. + .
. x .
H + n
I x .. coop
.. ... t
.. .93
x . X t
n. m 95sz potsbaco H.000—
. . xr
m. Wooom
H n .5228 p n .
.1 N .5328 x ...oomm
. F ...:mzcoo + 4
H .89: .65. ....
p — . _ . — . _ p p . — p _ . — . — - .Ooon

(flop/Cw)
saxn” Mopunoq paJnldDQ

80
4.3.2. Comparison of the Analytical Model and 2-D Finite

Element Eodel Results

Figure 4.11 shows stream depletion rates and captured ET
rates obtained from both models using 1 day time increments
for the computations during continuous pumping. This figure
shows that stream depletion rates obtained from both models
agreed well throughout the simulation for each configuration.
table 4.1 and table 4.2. Stream depletion rates obtained from
the finite element.model are observed to converge to the rates
obtained from the analytical model as Ateo.

Stream depletion rates agreed well for all the configura-
tions described in Table 4.2. However, the captured ET rates
obtained from the finite element model are considerably
different from those obtained with the analytical model during
early pumping for configurations 1 and 2. These results
agreed well throughout the simulation for configuration 3.
Figure 4.11 shows the simulation. results for the 'three
configurations. This clearly indicates that as spatial
discretization of the aquifer domain gets smaller in the.
vicinity of the well and the stream, captured ET rates
obtained from the numerical model converge to rates obtained
from the analytical model. At late time, when steady flow is
established, the two models converge for all configurations.

As stated earlier, the analytical model describes the

behavior of a semi-infinite aquifer while the numerical model

81
was applied to the finite aquifer domain shown in Figure 4.10.
It is apparent that there is no influence of the no-flow
boundaries in the results obtained with the analytical model,
whereas, the results produced by the numerical model could be
influenced if the drawdown cone reached the no-f low boundaries
of the aquifer. Figure 4.12 and Figure 4.13 show that the
water table is unaffected at the no-flow boundaries. Figure
4.12 shows the water table elevations at cross section II-II'
when a steady state condition is obtained. The well is
located 2200 meters from both the eastern and western no-flow
boundaries. The piezometric head at the boundary computed by
the finite element model is 30.0 meters which is the same as
the initial piezometric head and the drawdown at the well is
approximately equal to 1.6 meters. This amount of drawdown at
the well indicates that the analytical solution is in error in
the vicinity of the well since drawdown at the well is greater
than do = 0.7 meter. Figure 4.13 shows the water table
elevations at cross section I - I' at steady state. The water
table profile between the well and the northern no-flow
boundary of the aquifer indicates that this boundary of the

aquifer did not influence the results either.

4.4. Summary and Conclusions

An analytical model was developed that predicts the

stream depletion by a steady continuous pumping well where the

82

.naeu coca I u as .HH I an scuuoee esouo one occuo

sono- ucososo 32.3 sun 05 an sauna-co £385 .33 noun: .2; 933..

0325 82.55

.00NN..0£0N 0&0. damp .00.¢_ .OAWN— .onwo— .0w0 .0w0 .0“: .0wN .0

 

 

 

 

N 0:00:00 o 000.0N
— ...:00000 x
a
I 1000.0N
.0
I 1000.0N
D
I x 1000.0N
D
I p p x c o x 9 p x #00000
. L b p b p F L — n _ p b p _ oom.on

96.. coo
.__I__ .3003 30..

N
o

“O

a
o

ESQ 5500265

(renew) OVEH

83

.e>sc 000~ I u us.n I H scuuooe uuouo ecu anode
doves unusuao causuu nun 0:0 an 00000.00 suuuoum annoy have: .n~.v ousoum

0825 wozfima

PhD—blbbb Dbbh Pble Fur-_npbb— hhbb pnnbbbbhp DDDb—bL:

.0000 .0000 00.00 00.0.0 .00.? .0000 .omnn 00.0N .000N .00.0_ .000? .000

 

N 0:00:00 0
P0000000 x

E

D
I

h.
T u xox m xox u u u m u m m mouomomoupuomomomom

 

PbbbLhnnb—bhbb—nPr-bhh-PhLL-bbhbbkbb-lrbb.hb[—b-bLb-bub

00.0N

100.0N

100.0N

100.0N

100.00

 

 

98.. ooomu .
..I. 000000 00000 00 2100.5 0300305

.80».

(Jexew) Cl‘v’ElH

84

effects of captured ET are taken into account. The analytical
model is based on Hantush's solution for pumping from a leaky
phreatic aquifer connected to a bounding fully penetrating
stream. Hantush's leakage from the lower confined aquifer was
interpreted to describe a reduction in wetland ET resulting
from the drawdown of the water table in a shallow phreatic
aquifer. It was assumed that ET varies as a linear function
of the drawdown in the phreatic aquifer so that ET is at the
potential rate with zero drawdown and zero at a depth do. An
explicit relationship was derived that shows how potential ET
and do are related to Hantush's leakage coefficient. Proper-
ties of unsaturated soils were incorporated into the solution
by using Warrick's relationship between ET flux and depth to
water table to determine do.

Application of the analytical model showed that as the
dimensionless distance between the well and the stream
increases (increasing 6) , the portion of the pumping well
discharge captured from wetland ET increases; whereas, the
portion captured from the stream decreases. This resulted
from increased aquifer response time and increased area from
which ET could be captured. A dimensionless graph (Figure
4.7) was developed that shows the relationship between stream
depletion, water mined from aquifer storage, captured wetland
ET, and time. The analysis indicated that Figure 4.7 can be
used to determine the time, t,, when the system reaches

practical steady state. Analysis of Figure 4.7 also showed

85
that as scaled pumping distance (6) increases, the system
reaches steady state faster than the case where the effects of
captured ET was neglected.

A.dimensionless graph (Figure 4.8) for the rate of stream
depletion at steady state was also developed as a function of
the scaled pumping distance. This graph could easily be used
to determine the stream depletion rate and captured ET rate at
steady state for known value of the scaled pumping distance.

The solutions obtained in this study are appropriate only
when the actual field conditions approach the assumed
conditions. The accuracy of the steady state stream
depletion, based on the assumption that captured ET rate
linearly depends on the position of the water table, was
investigated. The analysis indicated that this assumption
could result in substantial error in stream depletion rate
when do was small. The analysis also indicated that the
magnitude of this error could be decreased as the value of do
was increased. However, there was no theoretical procedure to
determine the value of do which could minimize the error.

Additional error could be expected from the analytical
model results during the transient period. This error would
be produced as a result of the constant specific yield
assumption for soils that have low hydraulic conductivity
values. It was concluded from the above analysis that the
analytical model is restricted significantly as a quantitative

tool, but it is still valuable as a qualitative tool.

86
It was determined that the analytical model is also a
valuable tool for determining the optimum finite element grid
density which produces accurate stream depletion rates and

captured wetland ET rates during the transient period.

STREAM DEPLETION BY A PUNPING HELL INCLUDING THE EFFECTS OF

NONLINEAR VARIATION Ol' CAPTURE!) NBTLAND EVAPOTRANSPIRATION

5.1. Introduction

In the previous chapter, an analytical solution for
stream. depletion and captured ‘wetland ET jproduced by a
continuous pumping well from a shallow'water table aquifer is
derived. A linear variation of captured wetland ET is
incorporated into the flow equation in order to obtain the
analytical solution. The assumption of linear variation of
wetland ET as a function of water table depth is found to
overestimate the captured wetland ET and underestimate the
stream depletion for the values of do defined by Equation
4.14. When pumping starts and the drawdown cone develops and
propagates in the aquifer, reduction in wetland ET starts
immediately with the assumption of the previous chapter.
Contrary to the assumption of a linear variation of ET, it is
shown in previous chapter (section 4.2.3) and by Ripple et
al., (1972) that there is no change in ET until the depth to
the water table is equal to the critical depth, do. When the

depth to the water table exceeds do, ET is reduced according

87

88

to Warrick’s relationship (Equation 4.11) given in chapter 4.

The objective of this chapter is to develop a method that
can be used to predict the steady state stream depletion rate
produced by a pumping well which includes the effects of
captured wetland ET. .For'this method.a nonlinear variation of
wetland ET’as a function of the water table depth and the soil
properties is incorporated into the flow'equation. IBecause of
the nonlinear behavior of ET, the governing flow equation
becomes nonlinear, as shown in section 5.2.2. Thus, there is
no known analytical solution. To solve the flow equation for
stream depletion a 2-D finite element model (Townley and
Wilson, 1980) is employed for a hypothetical situation. To
generalize the numerical results for stream depletion,
dimensional analysis is employed to find the functional
relationship between the stream depletion rate, the aquifer
properties, and the pumping characteristics. As a result of
dimensional analysis, the scaled steady state stream depletion
will be shown to depend on five independent scaled variables
in the case where ET varied as a nonlinear function of the
water ‘table depth- Therefore, it is not practical to
generalize the simulation results in this case. In order to
reduce ‘the number of scaled independent (parameters, the
nonlinear variation of ET is approximated by a step variation
of ET as a function of water table depth. This assumption
reduced the number of scaled independent parameters to three;
thus, it. is practical to jpresent. the numerical results

graphically. The scaling of the flow equation in the case of

89
the nonlinear variation of ET as well as the step variation of
ET is explained in detail in section 5.2.4.

The accuracy of stream depletion rates obtained from the
numerical model by incorporating the step ET assumption is
determined by comparing stream depletion rates resulting from
the nonlinear ET assumption. The comparison is made during
the transient period because the numerical model results
oscillate in the step ET case. The oscillation is caused by
the spatial discretization of the aquifer domain and the
insufficiently'small time intervals that.are practical for the
transient computations. . The accuracy analysis of the
depletion rates is presented in section 5.3.2. Excellent
agreement is found between stream depletion rates that were
obtained from the numerical models. This indicated that
steady state stream depletion rates can be obtained
numerically by employing the nonlinear variation of ET and
that such results can be presented in the dimensionless form
developed by employing the step assumption.

From the results of the computer simulations, families of
dimensionless graphs were developed for a wide range of soil
hydraulic properties and pumping distances. Provided the
basic assumptions of the model are satisfied, these graphs
provide a means for easily determining the steady state stream
depletion and captured wetland ET rates if the following
parameters are known: the transmissivity of the aquifer, the
displacement head, the pore size distribution parameter (Anat

et al., 1965), the initial position of the water table, the

90
perpendicular distance between the river and the well, the

pumping rate, and the annual average PET rate.

5.2. model Development

5.2.1. lodel Assumptions

In addition to the assumptions for the aquifer properties
stated in chapters 3 and 4, the following assumptions are
required for wetland ET. These assumptions are necessary to
make use of the advantage of dimensional analysis technique
and to develop the families of dimensionless graphs from the
results of numerical simulations: (1) the ground surface is
located from the bottom of the root zone. With this
assumption, the steady state evaporation rate from the water
table can be replaced by the steady state ET rate from the
water table. This assumption was also used by Skaggs (1978)
to determine the reduction in steady state ET rate from a
shallow water table aquifer. Skaggs also incorporated the
effects of the storage in the root zone into his water
management model. However, the effects of the storage in the
root zone will not be included in this study in order to limit
the number of dimensionless parameters; (2) the rate of
evapotranspiration captured at the wetland surface is zero
when the depth to the water table is less than or equal to the
average critical depth, do, and occurs at the potential rate

when the position of the water table exceeds do (see Figure

-4 .. “a .eI' . _r 1-?‘n‘-n.~ulc

V53? "

. 91
5.1) . The concept of do as it relates to captured wetland ET
will be explained in more detail in section 5.2.4; (3) the
initial depth to the water table, d" is constant and parallel
to the ground surface; (4) the average critical depth is
always greater than the initial depth to the water table so
that the ET rate from the wetland surface is at the potential

rate before pumping.
5.2.2. Flow Equation

The governing differential equation for two-dimensional,
essentially horizontal, groundwater flow in a homogeneous

isotropic aquifer is (see Figure 5.1):

a

.2. 0 i .62 11.5242.
8x(hax)+dy(hay+Ko K,at (5'1)

In equation 5.1, K, is the saturated hydraulic conductivity,
5, specific yield, F flux per unit area that leaves or enters
the aquifer and, h(x,y,t) is the piezometric head during
pumping. If the bottom of the aquifer is coincident with the
horizontal datum, then, h(x, y, t) is the saturated aquifer
thickness. Because of the product (bah/61m , Equation 5.1 is
a nonlinear second order parabolic partial differential
equation. This equation is solved for h(x,y,t) by the 2-D
finite element groundwater flow model and the values of

h(x, y, t) are used to calculate the stream depletion rates and

92

 

 

 

 

 

 

 

 

 

Q
L
I’ET
1 351.1 1
‘0 _ ———- -—~- ""7... . u...‘
Stream K951 «4.! 4— qt
i?) h
<— 4.. «2.5.2.; h ° ‘—
Wi/I/l/III/l/i/l/ I I/IL/IILUW

 

P—a—‘i

PET

  
 

qr(l' + Ar) _1...

9.0)

 

Figure 5.1. Cross section of the semi-infinite stream-wetland
aquifer system described by Equation 5.1.

93
captured wetland ET rates as discussed in section 5.2.4.
To determine the most influential parameters on stream
depletion and captured wetland ET rates, Equation 5.1 is

manipulated as follows. The product (bah/6x9 can be written

 

as
an _ lam
32 2 6X1 (5.25)
and
K, at 2T t:

Substituting Equations 5.2a and 5.2b into Equation 5.1 and
manipulating gives
flj.fli.§. 32.113: (5.3)
8x2 ayz K; T at
Equation 5.3 is the partial differential equation for
groundwater flow in the system shown in Figure 5.1. This
equation looks similar to Equation 4.3 except for the third
term on the left hand side. This term represents the scaled
flux per unit surface area that leaves or enters the aquifer
and varies from PET/K5 to zero depending on the position of
the water table.

Suppose ho(x, y, to) describes the position of the water
table that would have prevailed in the flow system if the well
had not been pumped. That is, ho is the solution to the
boundary value problem shown in Figure 5.2, when the discharge

from the well is equal to zero at t = to. The flow equation

94
for this undisturbed system can be written by replacing h with

ho and ZIP/Ks with ZPET/Ks in Equation 5.3.

 

 

 

2 2 2
32ho+32ho+2pET= fgaho (5.4)

Equation 5.4 is linear in hf, since PET/Ks is constant.

0'0 (t-O)

rPEtrtH H ”J

 

 

 

 

 

 

g! -.-.-. ... .... ' ...... 0175....-. ......‘__.
‘ _ rm... 1 T" :
Stream Aquifer Kus, ‘— q.

<—%

 

 

 

 

 

'77,779,7777,177,111,1171177714I}iillllll/IIIIIIIIII

Figure 5.2. Cross section of the semi-infinite stream-wetland
aquifer system prior to pumping.

Combining Equations 5.3 and 5.4 by the superposition principle
yields

a2 (223-122) of a“ (123-122) + zq,= 5: 3013—122)

(5 . 5)
3x2 ay2 K, T at

 

In Equation 5.5 g,/Ko is the scaled ET flux leaving the aquifer

95
surface. This term is variable according to Anat et al.’s
(1965) equation 'that. is. presented in. the next section.
Therefore, Equation 5.5 is a nonlinear partial differential
equation and no analytical solution is known. Although,
Equation 5.5 is nonlinear, it can be manipulated so that the
most influential parameters that affect the stream depletion
and captured wetland ET could be obtained.‘ If the parameter

2 is defined as follows:

z = 1102-112 (5.6)

Substituting Equation 5.6 into Equation 5.5 yields

322.323.2qc. 1913—3 (5.7)
8x2 By2 R; 2'6:

The initial and boundary conditions are:

Z(x.y.0) = 0
Z(°°,y, t) = 0 (5.73)
Z(0.y. t) = o

Z(x, :I:°°, t) = O

The point sink that represents the discharging well requires

that

 

- .32 .. - 0
1:34 61’ nK, (5°7b)

where r’=(x-a) 2+y’. Equation 5.7 is still nonlinear because

the third term on the left hand side did not change.

96
5.2.3. The Nonlinear Relationship for q, as a Function of

Water Table Depth

In the previous chapter, Warrick's (1988) equation was
linearized to determine do in term of PET and the unsaturated
soil hydraulic properties. During the latter part of this
study, it was found that an explicit relationship had been
developed by Anat et a1. (1965) for estimating the steady
state ET flux in terms of the depth to the water table and the
soil hydraulic properties. An investigation was conducted to
determine if this relationship is as accurate as Warrick's
relationship. According to Anat et al. (1965), the steady
state ET flux is related to the depth to water table and the
unsaturated soil hydraulic properties by the following

equation:

 

. hold“) 1.886 n (5.8)
q’ a 0' (1+ n2+1 J)

The ET flux in Equation 5.8 may be limited by either of
two boundaries. The first is the atmospheric boundary at the
soil surface where the upper limit of the ET flux is PET, and
the second is the water table boundary where the ET flux is
negligible when the position of the water table is below do.
Figure 5.3a shows how the steady state ET flux varies between
the two boundaries and Figure 5.3b shows how the steady state

captured ET flux varies. The dashed lines in Figure 5.3a and

 

97

 

 

PET

 

Figure 5.3.

 

W

 

 

f d
d: (b)

(a) Dependence of steady state ET flux on the
depth to water table, (b) Dependence of steady
state captured ET flux on the depth to water
table.

98

Figure 5.3b depicts the nonlinear relation calculated by
Equation 5.8. The dashed line in Figure 5.3a shows that the
steady state ET flux is constant and equal to PET until the
depth to water table reaches do. .As the water table drops to
a depth that is greater than do, the steady state ET flux
starts to decline very rapidly, and theoretically, it becomes
zero only at depth d=w.

The critical depth where the ET rate is equal to the PET
rate can be found by replacing q, with PET and d with do in

Equation 5.8 and solving for dc;. this yields,

d = h ( K: )1,n 14.103132 (5.9)
c d PET 112-+1

 

5.2.4. Scaling the Flow Equation

An analytical solution of Equation 5.7 with respect to
given boundary conditions is not known. Therefore, a
numerical method (finite element method) is used to obtain the
solution. .

The numerical model used in this study requires
considerable time to operate and substantial expense for
computer time. Consequently, it is an objective of this
part of the study to obtain results in a form that may be
applied to as wide a range of field conditions as practical.

This is accomplished by expressing the governing equation and

99
its boundary conditions in terms of such dimensionless
variables that the model produces similar results within the
largest possible range of physical situations. This involves
choosing the most appropriate set of variables to describe the
physical system shown in Figure 5.1 and non-dimensionalizing
the variables with appropriate scaling variables.

In order to include the effects of a nonlinear variation
of q,/Ks on stream depletion, it is necessary to express this
term in terms of the most appropriate soil parameters
suggested by Equation 5.8, the initial conditions, and the
atmospheric bouhdary conditions. Equation 5.8 and the
discussion that followed it in the previous section shows

that,

 

22. .2 _d_1 PETES 51o
f(h.""h.' K.'h.) ( ’

Here d/ho is the scaled depth to water table, d,/ho is the
scaled initial depth to water table, do/ho is the scaled
critical depth. In Equation 5.10, the first two terms on the
right.hand side are introduced by Equation 5.8, the third term
is introduced by the initial condition, and.the last two terms
are introduced by the atmospheric boundary condition. In this
equation d,is the initial depth to water table, an important
physical parameter that influences the rate of stream
depletion and captured ET. For example, suppose dfi=05 that

is, the water table is at the ground surface before pumping.

100

In this case the water table has to be lowered more than do

before any wetland ET is captured. If initially d,=do, the.

reduction in wetland ET starts as soon as the water table
drops and this could produce a higher ET capture rate and
lower stream depletion rate.

As indicated above the objective of this section is to
find.the most.appropriate scaled parameters that influence the
magnitude of stream depletion. Therefore, dimensional
analysis was used to find these scaled physical parameters.
The mathematical relationship between the stream depletion

rate and Equation 5.7 is given by Hantush (1964b).

qr = 31911;?! dy (5.11)

According to Equation 5.11, the stream depletion rate can be
determined if the gradient of z in the vicinity of the stream
can be integrated along the stream boundary. The functional
relationship between the stream depletion rate and other
parameters suggested by Equations 5.7 and 5.11 can be written

as,

— . 1“ (Z) (5.12)

where,

..Qfa 32.52 5.13
Z K}. [ ,x,y,1gr Tt ( )

101
Equations 5.10 and 5.13 can be substituted into Equation 5.12.
The resulting equation can be manipulated by employing the
dimensional analysis technique, incorporating boundary
conditions of x-O and integrating from y=-co to y-oo as

suggested by Equation 5.11. This yields

 

3; = f £,n,$,£, K8 , Tt (5.14)
0 11,, do. do PET azsy

Equation 5.14 shows that the scaled stream depletion during a
transient period depends on six scaled independent parameters
when q, varies as a nonlinear function of the depth to the
water tableu The reason that the parameter d does not show up
in Equation 5.14 is because the parameter 2 is related to the
drawdown produced by the pumping well. In other words, the
solution of Equation 5.7 for z is equivalent to solving for
drawdown in the aquifer. Since d is the combination of the
drawdown and d,, it is redundant to include d in Equations 5.13
and 5.14. As tea, the system approaches steady state. Thus,
the scaled time can be eliminated from Equation 5.14.

g: _ a di a K.

_ .. _ ._ _ 5.15
o f(ho’n' dc' do.’ PET) ( )

 

Equation 5. 15 shows that even at steady state, the scaled
stream depletion depends on five independent scaled parameters
in the case of a nonlinear variation.of wetland ET. .Although,

the dimensional analysis technique revealed the appropriate

102
dimensionless parameters that influence the steady state
depletion rate, it is impractical to develop useful results in
this case. The impracticality stems from the extremely large
number of computer solutions that would be required and the
difficulty of using the information if obtained.

To reduce the number of scaled parameters introduced by
the assumption of a nonlinear variation of wetland ET and to
take advantage of dimensional analysis, the step variation of
ET assumption is adopted. With this assumption, qo/Ko is
either equal to PET/K5 or zero depending on whether the depth
to the water table is greater than do (See Figures 5.3a and
5.3b) . Moreover, Corey (1990) indicated that incorporation of
the step variation of wetland ET into the numerical model
would produce more accurate results than the linear variation
of wetland ET.

As shown in Figures 5.3a and 5.3b, it is necessary to
determine do so that the step variation of wetland ET can
accurately approximate the nonlinear variation of wetland ET.
The value of do shown in Figure 5.3 is determined by the
following procedure: (1) Equation 5.8 is evaluated numerically
for a wide range of soil types. (2) As a first guess, the
vertical line AB in Figures 5.3a and 5.3b is located on do.
(3) The shaded areas on the right and on the left side of AB
are computed. (4) The position of AB is adjusted by locating
it into a new position and then the shaded areas at the new
position are computed. (5) The depth do defines the position

of the vertical line AB so that the two areas are equal to

F“

.0 .______.. ”-..; .....- .. —_ we

.103
each other.

The procedure described above to find the value of do is
applied to soil types varying from fine clean sand to silty
sand. The soil types that are used to evaluate Equations 5.8
and 5.9 numerically are obtained from Anat et al.’s (1965)
report (page 23-26). Thus, the soil parameters, K5,)y, and n
are the same as theirs. The annual average PET rate used in
Equation 5.9 is taken as 0.005 m/day for humid.regions (Merva,
1990). The results of the numerical evaluations of Equations
5.8 and 5.9 show that the value of do is approximately 1.08
times the value of do for fine clean sand and 1.13 times for
silty sand. In this study, do= 1.1do is taken for practical
purposes. Thus, the value of do can be written in terms of
the displacement head, the saturated conductivity, the PET
rate and the pore-size distribution parameter using Equation

5.9 as

 

3,: 1 1h (( K, )1” (1+1.886)) (5.16)
c . d PET n2+1

For n26, the second term on the right hand side of Equation
5.16 is approximately equal to 1, thus Equation 5.16 can be

simplified further,

 

1
a; = 14(1)” “8 )3 (5.17)

The functional relationship between qe/Ks, the soil

. "...-av *‘u‘ul ‘ Wu.--M-—-ty

a. ‘

104
parameters, the atmospheric conditions and the initial
conditions for the block ET case can be obtained with the aid
of Equation 5.10 and the discussion that immediately follows.
As indicated before, the first two terms on the right hand
side of Equation 5.10 are introduced by Equation 5.8 which
determines the dashed curves in Figures 5.3a and 5.3b. Since
Equation 5.8 is not used to determine the ET flux as a
function of the depth to water table in the step ET case,
these two terms on the right hand side of Equation 5.10 can be
eliminated and also do can be replaced by do. Thus, Equation

5.10 in the step ET case simplifies to

 

& = _d_ _‘2 PET (5.18)
K, ac'ac’ K,

Equations 5.17 and 5.18 can be substituted into Equations 5.12
and 5.13. The resulting equation can be manipulated by
dimensional analysis and by the boundary conditions of x=0 and

by integrating from ys-co to y-uo suggested by Equation 5.11.

 

22.,53 3.!" K8 Tt (5.19)
Q ho' 1:, 'PET’ .323,

Equation 5.19 indicates that, in the case of step ET
variation, the solution of Equation 5.7 for the scaled stream
depletion rate, 9J0. depends on four independent scaled
variables. These variables are the scaled initial depth to

water table, (ch/ho)"; the scaled pumping distance, a/ho; the

 

105
scaled saturated hydraulic conductivity, Ks/PET; and the
scaled time, Tt/afig.
At steady state the scaled time in Equation 5.19 can be

eliminated. Thus, Equation 5.19 becomes,

 

3.! = f(£,($)n, Kn] (5.20)
ha PET

Since the number of independent scaled variables has been
reduced to three, it is practical to represent the results of

numerical simulations graphically.

5.3. numerical Solution of Flow Equation for Hypothetical

situation

“This section of the study is related to the development
of dimensionless graphs that can be used to predict the scaled
steady state stream depletion as expressed by Equation 5.20.
In order to develop families of dimensionless curves for
stream depletion rate, the 2-D finite element model Aquifem-l
(Townley and Wilson, 1980) is modified first to incorporate
the nonlinear behavior of wetland ET as a function of water
table depth and soil hydraulic properties, and second to
incorporate the step variation of wetland ET as a function of
the water table depth. For the remainder of this chapter, the
first model will be referred to as the nonlinear ET model

(NL—ET model) and the second model will be referred to as the

106
step ET model (ST—ET model). The governing flow equation of
the numerical model is given by Equation 5.1 which is
nonlinear and it is solved iteratively'by'employing'the finite
element method, based on the Galerkin technique with linear

finite elements.

5.3.1. The Hypothetical Aquifer

The hypothetical aquifer used for the numerical
simulations in this chapter is the same one (Figure 4.10) that
is described in chapter 4. In order to account for the
effects of scaled.pumping distance, the‘well is located at the
center line (I - I' in Figure 4.10) and the position of the
well varies from 100 meters to 2500 meters from the stream,
and pumps water from the aquifer at a constant rate. Table
5.1 summarizes the range of values of the hypothetical aquifer
properties. In this table, the range of values of the aquifer
thickness and the range of values of the initial depth to
water table are assumed. The range of the values of the
annual average PET rate was suggested by Merva (1990) for
humid regions. Also, the range of values of the saturated
hydraulic conductivity, the displacement head and the pore-
size distribution parameter are given by Brooks and Corey
(1964), Anat et al (1965), and Freeze and Cherry (1979).
These soil hydraulic properties are given for the upper range
of silt loam soils and for the upper to lower range of clean

sand. The range of the values of do is calculated according

.L“ g .'-'l.:.\. \_' ._

107
to Equation 5.10.

The aquifer domain shown in Figure 4.10a is divided into
triangular elements for the application of the numerical
model. The two grid configurations used in this part of the
study (configurations 2 and 3) are the same as those in
chapter 4. Table 4.2 summarizes the total number of elements,
the total number of nodes and the dimensions of the smallest
elements in the vicinity of the pumping well in each

configuration.

Table 5.1 Range of Values of the Aquifer Parameters

 

 

Parameter Range
Hydraulic conductivity“l” (m/day) 0.85 - 50
Displacement headflm (m) .10 - .70
Pore-size distribution parameter an, n 6 - 12
Average critical depth“h dg=1.1 dk (m) .20 - 1.50
initial depth to water table“), di (m) 0 - 1.20
Total aquifer thicknessm'(m) 30.00
Potential ET rate“ (m/day) 0.005 - 0.009

 

 

(1) Brooks and Corey (1964)
(2) Anat et al. (1965)

(3) Freeze and Cherry (1979)
(4) Computed by Eq. 5.10

(5) Assumed

(6) Merva (1990)

 

 

 

108

5.3.2. Accuracy of the Solution

The finite element model is applied to the aquifer
described above to predict the stream depletion produced by a
pumping well. The results of the numerical simulations are
used first to check the influence of the grid discretization
on the accuracy of stream depletion rates and captured ET
rates, and second to compare the accuracy of the depletion
rates that are obtained from the ST-ET model with the true
values of depletion rates. It is assumed that the true values
of depletion rates are obtained from NL-ET model. This
comparison is made within the range of aquifer parameters that
are given in Table 5.1. Finally, having established a method
for obtaining accurate solutions with reasonable computational
effort, steady state simulations with the NL-ET model are used
to develop dimensionless families of graphs for q, for the
range of aquifer parameters given in Table 5.1.

It is necessary to determine how configuration 3
influences the accuracy of stream depletion rates as well as
captured ET rates that are obtained from the numerical model.
This was done by comparing the numerical results with those
for two analytical models presented in chapter 3 and chapter
4, for steady continuous pumping. The numerical model is
applied to configuration 3 for two cases. In the first case,
the transient stream depletion is determined by neglecting the
effects of wetland ET. In the second case, transient stream

depletion rates and captured wetland ET rates are determined

109

by assuming a linear variation of ET as a function of water
table depth. The numerical model for each case is executed
six times for the three different values of Ks given in Table
5.1, and the two different values of 'a'. x; is chosen such
that the lowest, the mid, and the highest values of the range
are used. It is found that the depletion rates that are
obtained from the numerical model and from the analytical
models are in excellent agreement (the relative error < 1%)
for' configuration. 3. This. indicates 'that. spatial
discretization in configuration 3 is adequate.

It is concluded in section 5.2.4 that the depletion rates
can be presented graphically by using the ST-ET parameters if
the ST-ET model can accurately approximate the true values of
stream depletion. Therefore, it is necessary to determine the
accuracy of the depletion rates that are obtained from the
ST-ET model. This is determined by comparing transient q,
results for ST-ET and NL—ET generated using configuration 3.
The comparison of the results from both models is made for
1g=43.2 m/day, 8.64 m/day and 0.864 m/day when a=500 m and
cy=0.01 m and plotted in Figures 5.4, 5.5, and 5.6. Here the
values of ho and n are varied according to the values of Ks,
since these parameters depend on each other (Brooks and Corey,
1964; Anat et al., 1965; Corey 1986). Figures 5.4, 5.5, and
5.6 indicate that as the values of hydraulic conductivity
decrease the accuracy of the depletion rates obtained from the

ST-ET model increases.

".37).! attu‘fl'fi 73‘wofie. all? t v-

‘AIUM. ’. ‘5.

110

.a 3.6 «6 can .3... «alias .... can u a
c053 mamooa 9m 0: can .Bm noun 0:» .BE ucoswaco: 0:» .am Housed
on» scum Omcwmuno mmucu sawumammc Bumuum 0:» Ho comwuumaoo .v.m musofim

 

 

 

 

 

A963 Etc.
000— com _ com com com com 00¢ com com 00. o
_ _ . P . _ _ _ . _ . _ . _ . _ . _ . _ .
o o
I 0.2: %
u ........................... Edged ....................... e
1 Odom o
I . w
l i Iodon an”.
s \ . m.
I \ 10.8... m,
I T I.
\\ . o
I. \\\ \ ,IO DOW U
1 E 58:82 \\\\\u\\\ x , J
nllllIWnuudnlnuuunlhnunfl lllllllll D
... ...... E 8% m.
w m
C.
I /
I D.
o
.. B oz ..odom Mx
_ _ . _ _ _ _ _ . _ . _ . _ . _ _ _ _ 0.009

 

 

>on\E Nanny "mx .6me com. H o

111

.a 3.6 "6 can 6? 36".; .5 com
cons mamooa 8m 0: can .Bm noun 0:» .Bm newsflaso: 0:» .9m ummcwa
may deny omcflcuno mmumu :0wymanmo acmnum on» no comfiumnaoo .m.m musmflm

fl

 

 

A983 mEfi
000— com com com com com 00¢ com com 00, o
b — _ — _ _ _ b p — b — p — p — . — p .
o o

r T
.1. ........................................................ x 10.09 5

J
,. - m
I Odom
I \ a w
I \ [Odom P
I .... il'hflrllhllfhlil‘\ a. (ml
inhhffinhhffnhhlr 10.00... we:
1 l m.
I [Odom u
... I J
H. wqoom m.
I lodom \w/
1 pm Comet. Ix: . n;
.| PM QOHW ..l 1.0.Oom W
1 pm tomczcoz (I. - D
.| PM 02 l T0.00m r/AK

- — p _ p F b — p L! P P P — p — P h P 0.000—.

 

 

%.m umx BEE com u o

112

000.

owm 00m

 

 

.a 3.6 "6 can .3... Stuns .a com u
can: maoooa 9m 0: can .Bm noun ecu .Bm unusuacoc mnu .Bm nausea
mzu scum omcumuno mmumu sowucaomo acouum ecu no comuumoeou .o.m museum

A963 mEP

con 0mm

_ P . _ . _

owm

_

oov
_

_

p

_

5 52m

B 895.52

p

hm oz
_ .

 

em .605.- I ...

M

 

 

63E $2 nmx

..39: com H o

(flop/Cw) 910.1 uoneldep moans

113

Figure 5.4 shows that the relative error is approximately
3 % between the depletion rates obtained from the ST-ET model
and the NL-ET model when K,=43.2 m/day and d,=0.01 m. Since
the maximum relative error occurs when K5843.2 m/day, it is
necessary to investigate the accuracy of the ST-ET model
results for K‘s-43.2 m/day and for different range of values of
d, and 'a'. The ST—ET model is run when K,=43.2 m/day, d,=do,

=200 m and a=1000 m. The maximum relative error in these
cases is less than 6 %. This indicated that the ST-ET model
results are accurate for practical purposes, thus, the ST-ET
model parameterization (Equation 5.20) can be used to
represent the NL-ET model results (Equation 5.15).

Stream depletion rates obtained by incorporating the
linear variation of wetland ET and by neglecting the effects
of captured wetland ET (No ET model) are also plotted in
Figures 5.4, 5.5, and 5.6 for comparison. These figures
clearly document that the depletion rates obtained by
incorporating a step variation of wetland ET are very accurate
compared to depletion rates obtained by incorporating a linear
variation of wetland ET.

Although it is not readily apparent, Figure 5.4 shows
that stream depletion rates are fluctuating in the case of a
step variation of ET as the system approaches steady state.
These fluctuations occur due to the spatial discretization of
the aquifer domain and the insufficiently small time intervals
that are practical for the transient computations. It is

observed that these fluctuations could be eliminated if the

114
grid density were finer and time interval smaller. However,
this would make the computational cost as well as execution
time unfeasible.

In order to obtain ‘the scaled steady state stream
depletion suggested by Equation 5.20, the numerical model is
forced to predict the final steady state solution by using the
information at the final time step of the transient
calculations. This indicates that the accuracy of the final
steady state results depend. on 'the accuracy of results
obtained from the final time step of the transient
computation. Because of the oscillations of q, obtained by
ST—ET model during the transient period, predictions of steady
state q, with this model are inaccurate. It is shown that
during the transient period, stream depletion rates that are
obtained from the ST-ET model do accurately approximate the
depletion rates obtained from the NL-ET model. Therefore the
NL-ET model is used to predict the final steady state stream
depletion rate, and the results are presented by using the
ST-ET model parameters.

As indicated above, the final steady state stream
depletion is computed by using the results of the final time
step of the transient computations. Numerical simulations for
stream depletion in the case of nonlinear variation of ET
showed that the system approaches the final steady state
asymptotically. However, when the values of Roeare small and
'a' is large, the system approaches the steady state much

slower than in the case where R; is large and 'a' is small.

115

In these kinds of cases, the numerical model requires an
unreasonably long time for transient computations so that it
could predict the final steady state stream depletion rate.
Thus, the CPU time required to run the model would not be
practical. To reduce the CPU time and maintain the accuracy
in these cases, the steady state depletion rate is predicted
by extrapolating the lines for stream depletion and captured
ET until they intersect (dashed lines in Figure 5.7). The
error introduced by extrapolation is controlled by carrying
out the transient computations for a period long enough that
the rate that water is mined from aquifer storage (the rate
shown as BC in Figure 5.7) is less than 10% of the stream
depletion rate (rate AB shown in Figure 5.7). If the
extrapolation had not been used, the error in depletion rate
at this time would be 10 %. Thus, the error in steady state q,
is less than 10% after the extrapolation.

Figure 5.7 shows that for K5= 8.64 m/day and a=1000
meters, both stream depletion rate and captured ET rate
approach. their steady state ‘values very slowly as time
increases. It is observed from the numerical simulations that
this asymptotical approach to a steady state value is much
slower for lower range of values of R5 and higher values of
'a' than shown in Figure 5.7. Therefore, a longer period of
time is required to execute the numerical model for smaller
values of K51and larger values of 'a' for configuration 3 so
that the above criteria for extrapolation is achieved. For

example, it took 26 hours of CPU time to obtain the results at

116

.5 Toni .33 «Young .s 8on :95 Hmooa am

(89:13: wnu Eouu penumuno mmxsau unconson poucummo mo cauucuum> .56 museum

A283 mEz.

oooop oo_om oo_om coach oouom oopom oo_o¢ 00.3.. oo_o~.. coop

._.

o

_ ed

 

[T r T

I

 

4

3c. concave Ecobm

 

m
.llnunnHHHHHHHHHHHIw .32, 85:
u

Bot ...m ooaaucco

.—+~._»_._._.

16.2:
r.
rodow

‘—

Todon

10.00».
fiodom
r

Iodom

19.005
lodom

r-

fiodom

 

 

>8}: :3 aux 52: 8o— .. o

0.000—

(flop/ow)
sexnl; Mopunoq pelnidog

117
time t-5000 days for Ks=8.64 m/day and a=1000 meter shown in
Figure 5.7. In order to reduce this time and to control the
relative error within 10 %, the numerical model is applied to
configuration 2.

Since configuration 2 has coarser grid density than
configuration 3, it is important to check the influence of
configuration 2 on the accuracy of the depletion rate before
the development of the families of dimensionless graphs. The
scaled steady state stream depletion rates obtained by
applying the numerical model to configuration 2 are compared
with those obtained by applying the model to configuration 3
for three different values of Its/PET and a/ho. Table 5.2
shows these scaled depletion rates, the values of Ks/PET, a/hd,
and the relative error of configuration 2 depletion rates.
The relative error is determined by comparing the results of
configuration 2 with the results of configuration 3. Table
5.2 shows that the amount of relative error in q, obtained by
using configuration 2 increases as the value of a/ho increased
when Ks/PET - 10‘. In order to find the maximum amount of
relative error in steady state q,, the numerical model is

10‘. The

executed for a/h‘, = 20000 (a=2000 m.) and Ks/PET
relative error in this case is also less than 10% (9.6%).
This indicates that as the distance between the pumping well
and the stream increases, the growth rate in relative error in
steady state q, decreases. Therefore, at steady state,
configuration 2 produces approximately 90% accurate results.

This amount of accuracy in this study is assumed ‘to be

‘np.~._

118

acceptable for practical purposes, since using this

configuration reduced the CPU time about one order of

magnitude.

Table 5.2 Comparison of the scaled steady state q, obtained

from configurations 2 and 3

 

 

 

 

 

 

 

 

 

 

 

Config. 2 Config. 3
KS/PET a/h, q,/Q q,/Q Relative err. (1:)

2000 1.000 0.988 1

10‘ 6000 1.000 0.964 5
10000 0.972 0.902 8

800 0.649 0.670 3

103 2000 0.420 0.425 1
4000 0.287 0.291 1

286 0.325 0.326 .3

102 715 0.160 0.159 .6
1430 0.107 0.106 .9

5.3.3. Discussion of Scaled Families of Curves for Stream

Depletion

Dimensionless families of graphs that may be used to
predict the steady state stream depletion rates and captured

wetland ET rates are developed by using the results of

119

numerical simulations. Graphs (Figures 5.8, 5.9, and 5.10)
are plotted for values of Ks/PET equal to 10‘, 103, and 102 and
(ct/ho)" between 0 and 100. The symbols in Figures 5.8, 5.9,
and 5.10 show the values of steady state stream depletion
rates obtained from the numerical simulations. These data
points are connected to each other by fitting the best curve.
These figures provide a means for easily determining the
steady state stream depletion rates for values of Ks/PET,
(CL/ho)", and a/ho within the range of values presented.

Figures 5.8, 5.9, and 5.10 show that the scaled steady
state stream depletion, q,/Q, approaches zero as a/hd
approaches infinity; whereas, qr/Q approaches 1.0 as a/ho
approaches zero. This can be interpreted as follows: as the
pumping distance increases, the aquifer response time
increases and this produces a lower rate of stream depletion.
An increase in the value of a/ho also indicates that the area
from which ET may be captured increases, thus, producing
higher ET capture. This trend can be observed from these
figures for all the range of values of Ko/PET and (cf/ho)".

Figure 5.8 shows that when the value of the scaled
conductivity is high, the steady state stream depletion rate
can not be neglected even for large values of a/ho. This
figure also confirms the importance of the initial depth to
water table as explained in section 5.2.4. For instance, as
the value of ((1)/ho)” approaches zero, steady state stream
depletion rates increase, and the major portion of the well

discharge is captured from the stream, although the value of

120
a/hd is large. As the value of d, approaches zero, captured ET
from the wetland decreases. This occurs because the sum (di-i-s)
approaches the drawdown, s, and ET capture begins only when
s>do instead of (d,+s) >do, and this causes a lower ET capture
rate and a higher stream depletion rate.

In Figures 5.8, 5.9 and 5.10, the lowest curve for (CL/ho)’I
is the case where the product d,/do =5 1. In this case, ET
capture from the wetland starts immediately after the pumping
starts, and as the value of a/ho increases, the major portion
of the well discharge is being captured from wetland ET. The
upper horizontal line in Figures 5.8, 5.9, and 5.10 (q,/Q =
1.0) is the case where d, >do. In this case, the ET rate from
the wetland surface is zero, thus, all the pumped water is
captured from the stream for all the values of a/ho.

Figures 5.9 and 5.10 show the variation of the steady
state stream depletion rates for lower values of Ks/PET. The
curves in these figures are closer to each other than the
curves in Figure 5.8 for the given range of (CL/ho)". This can
be interpreted.as follows: the effects of the initial position
of the water table is a more influential parameter for
determining steady state stream depletion rate when K5 is
higher. Although d, is more effective when K5 is higher,
Figures 5.9 and 5.10 also show the effects of (CL/ho)" on the
magnitude of q,/Q distinctively. This indicates that the
initial position of the water table can influence the stream
depletion and/or captured ET rates even when the value of

saturated hydraulic conductivity is low.

121

..S a mafia. can: 033
ueums ou cudec Hmuuucu peacee ecu one eeccueuo acucaco peacee
ecu co scuueameo Sceuue eucuu acceue ceaeom ecu no eeceocemeo .m.m eusmum

Pic .eocBmfi mEaEna 028m

 

.ooooe .ooom .ooo_
e _ e e s p — b p b — p _ 0.0
1.2 u Ea\mx 8. u :03}; I Ito

1. E u :31}; III

. ... swig elm [8
T? u swig ole IS
0 if}; I

 

l¢.0
[0.0
[0.0
1.5.0
[0.0

 

[0.0

o/Jb ‘uoneldap moans paloog

 

 

 

 

122

.noH u 33% secs e33
ueums ou cuaec Hmuuwcu peaeom ecu 0cm eoceumwu acumasa ceaeee
ecu co scuueumec Emeuue eueum acceum peacem ecu mo eocencemeo .m.m eucmum

nz\c .eucccmfi mEaEza 028m

 

 

 

.ooooe .ooom .ooo, .oom
\— - LP h h P L it P — b b b n — 00.0

. x

l: o m.

. a0

ION O D.

(one B

1

m
13.0

w

I load W

_ r000 .mul

I .61

Elm robe m.

U

I 00.0 ...D

I o n causes 010 n? u fining. om o W
— P . b p — F . » PL . p L _ - OO.P

 

 

123

..3 u .533. :95 383
ueums ou cudeu Heuuucu ceamoe ecu 0cm eocmumuc ocusaso peacom
ecu co cauueameo Sceuum eucum aoceue ceaeom ecu no eeceocemeo .oa.m eucmum

Uc\c .eoccbmfi ocfiEsa ne_oom

 

 

 

 

 

 

.ooo_ .oon .002
r p _ b p L L — L - p 00.0
. x
12 o m.
. ad
lONO D.
Tome E
J
m
rote
w
[one may
_ _ 88 01
I W
.. «L: n 923 I 12.6 m.
I To. ... c333 elm. [8.6 .0
7.2 u :03an I m0
I . o .... c333 Ole N? n Emmy. road W
p — r b L P P p n p OO.F

 

124

As indicated above, Figures 5.8, 5.9, and 5.10 can also
be used to predict the captured wetland ET since the system is
at steady state. An estimate could be made based on the
assumption that all pumped water is obtained from captured ET.
That would result in an estimate that is at.most 20 % in error
when the value of a/h, is greater than 5000 (see Figure 5.9) .
The scaled pumping distance decreases for the same amount of
relative error caused by the assumption.that all of the pumped
water that is captured from the wetland ET decreases as Ks
decreases, as seen ianigure 5.10. 'This figure shows that the
maximuerelative<error at steady state is less than.10% if one
assumes that all of the pumped water is captured from the
wetland ET when a/h, 2 2000. This clearly indicates that as
the value of R; decreases while 'a’ increases, the aquifer
response time increases. Thus, steady state stream depletion
will be low, and the major portion of the well discharge will

be captured from the wetland ET.

5.4. Summary and Conclusions

In this part of the study, an investigation was conducted
to develop a model that can be used to predict the steady
state stream depletion rates caused by a continuous pumping
wellm The effects of nonlinear variation.of wetland ET on the
steady state stream depletion rate was investigated in the

model. The model was based on combining the governing flow

125
equation for a water table aquifer and Anat et al's equation
for nonlinear ET from the water table.

Dimensional analysis was performed to find the functional
relationship between the scaled steady state stream depletion
and the physical aquifer parameters. It was found that the
scaled steady state stream depletion depended on five
independent scaled aquifer parameters in the case where ET
varies as a nonlinear function of the water table depth. It
was also found that these independent parameters were reduced
to three when the step ET assumption was employed to
approximate the nonlinear ET assumption. The scaled
parameters in the case of step ET assumption were: the scaled
hydraulic conductivity, the scaled pumping distance from the
stream, and the scaled initial depth to water table.

To incorporate step variation of wetland ET into the
model, an explicit relationship was found for the average
critical ET depth, do, by combining Anat et al’s equation with
the atmospheric boundary condition. In the model it was
assumed that captured ET from wetland is zero when the sum of
the initial depth to water table and.the drawdown is less than
or equal to do; whereas, captured ET is equal to PET when the
sum of the initial depth to water table and drawdown is
greater than do. It was found that stream depletion rates
obtained by employing this assumption were very accurate as
compared to stream depletion rates obtained by incorporating
the linear variation of the ET assumption.

A 2-D finite element model for saturated groundwater,

126

AQUIFEM-l, was modified to incorporate the nonlinear behavior
of wetland ET as a function of the position of the water table
and the soil hydraulic properties. The numerical model was
applied to a hypothetical situation for a wide range of scaled
independent aquifer parameters. The results of the numerical
simulations were used to develop the families of dimensionless
graphs that could be used to estimate the rate of steady state
stream depletion for a practical range of scaled independent
parameters. Three families of dimensionless graphs were
developed between the range of values of the scaled hydraulic
conductivity of 10‘ - 102, the scaled initial depth to water
table of 0-100, and the scaled pumping distance of 100-25000.
The following conclusions were made from the analysis of
the dimensionless graphs: the scaled steady state stream
depletion rate approached zero as the scaled pumping distance
approached infinity; whereas, the scaled steady state stream
depletion approached 1 as the scaled pumping distance
approached zero. When the value of the saturated hydraulic
conductivity was high, the error caused by the assumption that
all of the pumped water is captured from wetland ET is
substantial even for large values of a/ho. The dimensionless
graphs showed that steady state stream depletion rates
strongly depended on the initial position of the water table.
Analysis showed that as the saturated conductivity increased,
the effect of the initial position of the water table on the

magnitude of stream depletion rate was more influential.

127
Analysis also showed, that as the value of saturated
conductivity decreased, the error caused by the assumption

that all of the pumped water is captured from wetland ET,

decreased.

CHAPTER VI
BOMB! m CONCLUDING me

In the present study, three different mathematical models
were developed. These models can be used to predict the
stream depletion and/or captured wetland ET rates produced by
a pumping well from the phreatic aquifer which is
hydraulically connected to an adjacent stream. The
mathematical treatment required standard use of superposition
theory and the existing analytical solutions for steady,
continuous pumping. The solutions obtained are appropriate
only when the actual field conditions approach the assumed
conditions. The study focused on the approach to dynamic
equilibrium as well as conditions at dynamic equilibrium.

The first model is an analytical solution for stream
depletion produced by nonuniform, cyclic pumping from a well
located in a deep water table aquifer. The effects of
captured wetland ET were not included in this model. An
analytical expression for dimensionless volume of stream
depletion (r) was developed. This expression can be used to
determine the time required to achieve a practical state of
dynamic equilibrium. An analytical expression for 'r was

obtained (Equation 3.10) and plotted (Figure 3.8) for a

128

129

practical range of the independent variables. Equation 3.10
was obtained by recognizing that the volume of stream
depletion over one cycle, from t-tg‘to t, is the same as the
volume of stream depletion between the start of pumping and
the time t by a single period of pumping (Figure 3.6).
Equation 3.8 shows that 1 depends on the values of a and 7 but
that it is independent of the value of tad; even though the
stream depletion rate shows such a dependence (Equation 3.6).
Analysis of the 1 relationship produced Equation 3.11 which
shows that time to equilibrium depends on t; alone for small
values of 7. An inspection of conditions when 7 is large
leads to the recognition that such wells are close enough to
the stream to reach equilibrium within the first pumping
period so that the nonuniform pumping pattern is an accurate
representation of the stream depletion pattern.

Having obtained a consistent basis for determining the
occurrence of dynamic equilibrium permitted the study of
stream depletion rates at dynamic equilibrium. The 6
criterion (Figure 3.10) was developed to characterize the
error in depletion rates estimated by assuming steady,
continuous pumping at the cycle-average rate. Figure 3.10
showed that under some circumstances, the cycle-average
approximation is not an adequate representation of the pumping
pattern.

The second model that was developed is also an analytical
model to predict the stream depletion by a steady continuous

pumping well where the effects of linear variation of captured

<= ~':,. ”1]
02.14%

-r'm‘s‘ a.

r99.

fiz- '

130

ET were taken into account. The analytical model is based on
Hantush's solution for pumping from a leaky phreatic aquifer
connected to a bounding fully penetrating stream. Hantush's
leakage from the lower confined aquifer was interpreted to
describe a reduction in wetland ET resulting from the drawdown
of the water table in a shallow phreatic aquifer. It was
assumed that ET varies as a linear function of the drawdown in
the phreatic aquifer so that ET is at the potential rate with
zero drawdown and zero at a depth do. An explicit
relationship was derived that showed how potential ET and d0
were related to Hantush's leakage coefficient. Properties of
unsaturated soils were incorporated into the solution by using
Warrick's relationship between ET flux and depth to water
table to determine do.

Application of the analytical model to a hypothetical
situation showed that as the dimensionless distance between
the well and the stream increases (increasing 6) , the portion
of the pumping ' well discharge captured from wetland ET
increases; whereas, the portion captured from the stream
decreases. This resulted from increased aquifer response time
and increased area from which ET could be captured. A
dimensionless graph (Figure 4.7) was developed that shows the
relationship between stream depletion, water mined from
aquifer storage, captured wetland ET, and time. The analysis
showed that Figure 4.7 could be used to determine the time,

t when the system reaches a practical steady state. The

(I

analysis of Figure 4.7 also showed that as scaled pumping

131
distance (6) increased, the system reached steady state faster
than the case where the effects of captured ET was neglected.

A dimensionless graph (Figure 4.8) for the rate of stream
depletion at steady state was also developed as a function of
the scaled pumping distance. This graph was used to determine
the stream depletion rate and captured ET rate at steady state
for known value of the scaled pumping distance. T h e
solutions obtained in this study are appropriate only when the
actual field conditions approach the assumed conditions. The
accuracy of the steady state stream depletion based on the
assumption that the captured ET rate linearly depended on the
position of the water table was investigated. The analysis
indicated that this assumption could result in substantial
error in stream depletion rate when (10 was small. The
analysis also indicated that the magnitude of this error could
be decreased as the value of (10 was increased. However, there
was no theoretical procedure to determine the value of d0
which could minimize the error.

Additional error could be expected from the analytical
model results during the transient period as a result of a
constant specific yield assumption for soils that have low
hydraulic conductivity values. It was concluded from the
above analysis that the analytical model is restricted
significantly as a quantitative tool, but it is still valuable
as a qualitative tool.

An existing finite element model was modified extensively

in order to incorporate captured wetland ET as a linear

132

function of the depth to water table in the aquifer.
Application of the finite element model and the analytical
model to a hypothetical situation showed that the results
obtained from the finite element model could be made to
converge to the results obtained from the analytical model if
the spatial discritization of the aquifer domain and time
interval were made sufficiently small. Therefore, the
analytical model is a valuable tool for determining the
optimum finite element grid density which could produce
accurate results during the transient period.

In the third model, the effects of a nonlinear variation
of wetland ET was investigated on steady state stream
depletion rates caused by a continuous pumping well. The
model was based on combining the governing flow equation for
water table aquifer and.Anat et.al's equation for nonlinear ET
from the water table.

Dimensional analysis was performed to find the functional
relationship between the scaled steady state stream depletion
..and the physical aquifer parameters. It was found that the
scaled steady state stream depletion depended on five
independent scaled parameters in the case where ET varied as
a nonlinear function of the water table depth. It was also
found that the number of these independent parameters was
reduced to three when the step ET assumption was employed to
approximate the nonlinear ET assumption. These scaled
parameters were the saturated hydraulic conductivity, the

pumping distance from the stream, and the initial depth to

e. Ime~.-.L- - -uc.‘-1 ...L‘1

133
water table.

To incorporate the step variation of wetland ET into the
model, an explicit relationship was found for the average
critical depth, do, by combining Anat et al’s equation with
the atmospheric boundary condition. In the model it was
assumed that captured ET from the wetland is zero when the sum
of the initial depth to water table and the drawdown is less
than or equal to do; whereas, the captured ET is equal to PET
when the sum of the initial depth to water table and drawdown
is greater than do. It was found that stream depletion rates
obtained by employing this assumption were very accurate
compared to stream depletion rates obtained by incorporating
the linear variation of ET assumption.

A 2-D finite element model for saturated groundwater,
AQUIFEH-1, was modified to incorporate the nonlinear behavior
of wetland ET as a function of water table depth and soil
hydraulic properties. The numerical model was applied to a
hypothetical situation for a wide range of scaled independent
aquifer parameters. The results of the numerical simulations
were used to develop families of dimensionless graphs that
could be used to estimate the rate of steady state stream
depletion for a practical range of scaled independent
parameters. Three families of dimensionless graphs were
developed between the range of values of the scaled hydraulic
conductivity of 10‘ - 10’, the scaled initial depth to water
table of 0 - 100, and the scaled pumping distance of 100 -

25000.

134

The following conclusions were made from an analysis of
the dimensionless graphs: the scaled steady state stream
depletion rate approached zero as the scaled pumping distance
approached infinity whereas, the scaled steady state stream
depletion approached 1 as the scaled pumping distance
approached zero. When the value of the saturated hydraulic
conductivity was high, the error caused by the assumption that
all the water is captured from wetland ET is substantial even
for large values of a/ho. The dimensionless graphs showed
that the steady state stream depletion rates strongly depended
on the initial position of the water table. Analysis showed
that as the saturated conductivity increased, the effect of
the initial position of the water table on the magnitude of
stream depletion rate was more influential. Analysis also
showed, that as the valuezof the scaled saturated conductivity
decreased, the error caused by the assumption that all the
pumped water was captured from wetland ET, decreased.

At this point some suggestions are made for future
analytical and numerical work on stream depletion and captured
wetland ET rates caused by either continuous or cyclic pumping
of the well.

For a more complete understanding of the physics of the
problem, it is important to modify the numerical model further
by incorporating the 1-D unsteady Richards' Equation between
the water table and the ground surface. This would allow the
specific yield to vary with the position of the water table

and time.

135

Partial penetration of the aquifer by a stream including
the presence of the clogging layer between the streambed and
the aquifer is more realistic than a fully penetrating stream
with no clogging layer. Although this option was already
available in the numerical model, it was not used in the
numerical part of the present study. The effects of these
parameters on stream depletion and captured wetland ET rates
can be investigated in the future study.

It is also interesting to study the effects of the
partial penetration of stream with the presence of a streambed
clogging layer on stream depletion caused by cyclic pumping of
a well. This can be accomplished by combining Hantush's
theory (1965) with the superposition principle that was used
in the present study.

Also as a natural extension to the second analytical
model, a theoretical development is suggested for determining
the value of do that improves the accuracy of the analytical
model as a quantitative tool.

Finally, it is of equal importance to investigate further
the possibility of obtaining an analytical solution for stream
depletion by a pumping well where the effect of step variation

of captured ET is accounted.

APPENDICES

APPENDIX A

SOLUTION OP EQUATION 4.8 E! LAPLACE TRANSPORNLIIONB

Analytical solution of Equation 4.8 can be obtained by

applying Laplace transformations (Hantush, 1964a).

%+%§-_BZE = %% (A.1a)
Z(x,y,0) = 0 (A.1b)
z(oo,y, t) = 0 (A.1c)
Z(x, 10°, t) = 0 (A.1d)

Equations A.1a, A.1b, A.1c, A.1d can be written in radial

 

coordinates.
.52_Z.+l.a£-£ a £22 (A.2a)
r3 r r 32 T t
z<:,0) - 0 (A.2b)
Z(°°,t) - 0 (A-ZC)
£10m (1.3% = 33(- (A.2d)

137

138
Here r is the radial distance from the center of a pumping
well located at (xo,yo) to any point in the surrounding area
and is equal to [(x-xo)2+(y-yo)2]m. Applying ‘ Laplace
transformations to Equations A.2a, A.2b, A.2c, and A.2d one

can obtain,

 

fl+l§Z-(£§+i) 2 = o (A.3a)
3:2 rar T 32

7(co'p) = O (A.3b)
lim (Ia—2) = - Q for t > 0 (A-3C)
r-0 8: «pk"

In Equations A.3a, A.3b, and A.3c, p is the transformation
variable. Equation A.3a is the zero order modified Bessel

Equation and the general solution is (Hantush, 1964a)

2 = C1K0(Nr) +C2Io (Nr) (A- 4)

where N’- pS/T + 1/B’. In Equation A.4, Io and Ko are the zero
order modified Bessel functions of the first and the second
kinds, respectively. Since Ko(oo)=0, Io(co)=oo, and

6Ko(Nr)/6r=-N K,(Nr) and as xeo, lim x K,(x)-1, (Carslaw, and
Jaeger, 1959) the constants Cl and C2 can be obtained by
evaluating Equations A.3c and A.4. K, in the above expression
is the first order modified Bessel function of the second

kind.

 

:1 Q
1 npK, (A.5a)
CZ = 0 (AeSb)

Substituting Equations A.5a and A.5b into Equation A.4 yields,

NI
u

 

Q 1
— K NI 0
"P199 o( ) (A 6)

Substituting the value of N into Equation A.6 and taking the
inverse Laplace transformation of the resultant equation gives
the analytical solution of boundary value problem defined by

Equation A.1 (Hantush, 1964a)

 

2
z= 0 W(L§,£) (A.7)
21:16, 42': B

where W is the well function for leaky aquifers which is
available in tabular form in the literature (Hantush, 1956;
and Walton, 1962). Equation A.7 is identical to Equation 4.15
that is given in chapter 4.

The rate of stream depletion per unit with of an aquifer
can be obtained by combining Equation A.7 with Darcy’s Law
(Hantush, 1964a)

qzw=0.5K, (g3) - (A,3)

x-O

The total rate of stream depletion then can be obtained by

140
integrating Equation A.8 along the stream boundary (Hantush,

1964b)

= ” a_z A.9
q, O'SK'£(aXx-ody ( )

and the volume of stream depletion at time t can be obtained

by integrating Equation A.9 from t=0 to t (Hantush, 1964b)

c-co

V .-.- [qr dt (A010)

2‘
8.0

The mathematical reduction of Equations A.9 and A.10 gives
Equations 4.17 and 4.18. 'To evaluate the infinite integral in
Equation A.9 and the definite integral in Equation A.10,

Laplace transformations were used (Hantush, 1964b).

APPENDIX B

LISTING OF COMPUTER PROGRAM

FOR

STREAM DEPLETION CAUSED BY CYCLIC PUMPING OF WELLS

141

()0OOOOOOOOOOOOOOOOOOOOOOOO

142

********************************************************

THIS PROGRAM CALCULATES THE IMPACT OF IRRIGATION
WATER ON INSTERAM FLOW QUANTITY BY USING UNSTEADY
GROUNDWATER FLOW EQUATIONS. THE PROGRAM COMPUTES

THE IMPACTS FOR CYCLIC PUMPING. IT IS ASSUMED THAT
THE PUMPING RATIO Q AND AQUIFER CHARACTERISTICS

ARE CONSTANT WITH TIME AND SPACE.

WRITTEN AND DEVELOPED
BY
Yakup Darama Ph.d Candidate

Michigan State University
College of Engineering

Department of Civil and Env. Engineering

East Lansing, Michigan 48824

 

 

INTEGER ANSWER,MUNIT,EUNIT,OPTION,FLAG

IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION RATEQ(4000)

CHARACTER FLNAME*10

WRITE(*,*) ' OUTPUT FILE NAME FOR BETAl AND BETAz ?'
READ(*,'(A)')FLNAME
OPEN(UNIT=5,FILE=FLNAME,STATUS='NEW',FORM='FORMATTED')
WRITE(*,*) ' OUTPUT FILE NAME FOR SDF ?'
READ(*,'(A)')FLNAME
OPEN(UNIT=6,FILE=FLNAME,STATUS='NEW',FORM='FORMATTED')
WRITE(*,*) ' OUTPUT FILE NAME FOR TAU ?'
READ(*,'(A)')FLNAME
OPEN(UNIT=7,FILE=FLNAME,STATUS='NEW',FORM='FORMATTED')

************* INPUT VARIABLES *********

143

DX
PI

1.0
4.0*ATAN(DX)

WRITE(*,7)
FORMAT(/,1OX,'UNIT SYSTEM FOR INPUT DATA',/
* ,10x,'0 = METRIC UNIT ',/
* ,on,'1 = ENGLISH UNIT ',//
* ,lOX,'PLEASE SPECIFY --> '/)
READ(*,*)IO
IF(IO.EQ.O) THEN

MUNIT = o
ELSE

EUNIT = 1
ENDIF
WRITE(*,*) ' STORATIVITY Sy ?'
READ(*,*) STOR

IF(EUNIT.EQ.1) THEN

WRITE(*,*) ' TRANSMISSIVITY IN ft‘2/day ?'

READ(*,*) TRANS

WRITE(*,*) ' DISTANCE BETWEEN WELL AND RIVER IN ft ?'
READ(*,*) DIST

WRITE(*,*) ' PUMPING RATE Q IN ft‘3/day ?'

READ(*.*) Q

ELSEIF(MUNIT.EQ.0) THEN
WRITE(*,*) ' TRANSMISSIVITY IN m‘2/day ?'
READ(*,*) TRANS
WRITE(*,*) ' DISTANCE BETWEEN WELL AND RIVER IN m. ?'
READ(*,*) DIST
WRITE(*,*) ' PUMPING RATE Q IN m‘3/day ?'
READ(*.*) Q
ENDIF

WRITE(*,*)' PUMPING PERIOD DURING ONE CYCLE TP ?'
READ(*,*) TP

WRITE(*,*)' CYCLE LENGTH TD in days ?'
READ(*,*) TD

WRITE(*,*)' TIME INCREMENT DT in days ?'
READ(*,*) DT

WRITE(*,*)' DLT in days FOR CONV TO TAU=0.95 ?'
READ(*,*) DLTl

WRITE(*,*)' NUMBER OF MAX ITERATION MAXIT ?'
READ(*,*) MAXIT

WRITE(*,*)' TOLERANCE FOR TAUSTR=0.95+TOL ?'
READ(*,*) TOL

WRITE(*,*)' THE VALUE OF T1 in days ?'

 

144

READ(*,*) TI

DDT DLTl

0
TA = ( DIST*DIST*STOR )/TRANS
c
WRITE(5,8)TA
WRITE(*,8)TA
8 FORMAT(10X,'AQUIFER RESPONSE TIME = ',DZO.12)
c
I = o
FLAG = 1
c
C ----------------------------------------------------

c
50 x1 = DSQRT((TA)/(4.0*TI))
C ----------------------
CALL ERRFNC (X1,ERFCX)
C ----------------------
xx1 = x1*x1
EXPTRl = l.0/(EXP(XX1))
A = ((((TA)/(2.0*TI))+l.0)*ERFCX)
B = (2.0*X1*EXPTR1)/(DSQRT(PI))
VC1 = (TI/TP)*(A-B)
C
IF(TI.GT.TP)THEN
c _
TMDIF =TI-TP
x2 = DSQRT((TA)/(4.0*TMDIF))
C ------------,. .........
CALL ERRFNC (X2,ERFCX)
C ----------------------
xxz = x2*x2
EXPTRZ = 1.0/(EXP(xx2))
AA = ((((TA)/(2.0*TMDIF))+l.0)*ERFCX)
BB = (2.0*X2*EXPTR2)/(DSQRT(PI))
v02 = (TMDIF/TP)*(AA-BB)
TAU = VC1-VC2
c
ELSEIF(TI.LE.TP)THEN
TAU = VC1
ENDIF
c
c -------------------- - ----------------------------

C ***** FIND THE VALUE OF te WHEN TAU = 0.95 ******

145

TAUSTR = 0.95+TOL
IF(I.LE.MAXIT) THEN
I a 1+1
IF((TAU.GE.0.95).AND.(TAU.LE.TAUSTR)) THEN
TE = TI
TTA = TE/TA
WRITE(*,*) ' CONVERGENCE REACHED'
GO TO 45
ELSEIF(TAU.LT.0.95) THEN
IF(FLAG.EQ.2) THEN
DLT = DLT/2.0
ENDIF
FLAG = 1
T1 = TI+DLT
WRITE(*,51)TAU,TAUSTR
51 FORMAT(5x,'TAU =',F20.10,5X,'TAUSTR =',FlO.7)
GO TO 50
ELSEIF(TAU.GT.TAUSTR) THEN
IF(FLAG.EQ.1) THEN
DLT = DLT/2.0
ENDIF
FLAG = 2
TI = TI-DLT
GO TO 50
WRITE(*,52)TAU,TAUSTR
52 FORMAT(5X,'TAU =',F20.7,5X,'TAUSTR =',F10.6)
ENDIF

ELSE
WRITE(*,56)I,MAXIT
56 FORMAT(5X,I4,1X,'CONVERGENCE NOT SUCCEEDED',/,
1 5X,'INCREASE THE VALUE OF MAXIT >',15)
GO TO 999
ENDIF

C **** COMPUTATION OF IMPACT AT ANY FUTURE TIME ******

45 WRITE(*,46)
46 FORMAT(/,10X,'CALCULATE qmax ?',/
*,1ox,'(enter 0PTION)',/

*,10X,'TYPE 0 < ----------- YES',/
*,lOX,'TYPE 1 < ----------- NO',//
*,10X,'PLEASE SPECIFY ------- >',/)
READ(*,*)OPTION
IF(OPTION.EQ.1) THEN
GO TO 190

ELSEIF(OPTION.EQ.0) THEN
WRITE(*,91)

146

91 FORMAT(5X,' calculating the max. q')

N = INT(TE/TD)+1
NN = INT(TD/DT)+5
DO 4000 I=1,NN
RATEQ(I) = 0.00
4000 CONTINUE
C
K = NN
5000 IF(K.GE.1) THEN
C
TI = TE-((NN-K)*DT)+DT
C
write(*,*)K
SUM1 = 0.00
SUM2 = 0.00
I = 0
100 IF(I.LE.N) THEN

TI-I*TD
TI-TP-I*TD

ARGl
ARG2

IF(ARG1.GT.0) THEN
x = DSQRT((TA)/(4*ARGl))
CALL ERRFNC (X,ERFCX)
SUM1 = Q*ERFCX

ENDIF

IF(ARG2.GT.0) THEN
x = DSQRT((TA)/(4*ARGZ))
CALL ERRFNC (X,ERFCX)
SUM2 = Q*ERFCX

ENDIF

IF((ARG1.GT.0).AND.(ARG2.LE.0)) THEN
RATEQ(K) = RATEQ(K)+SUM1
ELSEIF((ARGl.GT.0).AND.(ARG2.GT.0)) THEN
RATEQ(K) = RATEQ(K)+SUMl-SUM2
ELSEIF((ARG1.LE.O).AND.(ARGZ.LE.O))THEN
SUM1 = 0.D0
SUM2 = 0.00
RATEQ(K) = RATEQ(K)+SUM1-SUM2
ENDIF

I = I+1
GO TO 100

147

ENDIF
C
WRITE(*,21)K,RATEQ(K),TI
21 FORMAT(5X,'RATEQ(',I3,')=',F12.5,5X,'t =',F10.2)
C
C —————————————————————————————————————————————————————————
C **** FIND MAX IMPACT DURING THE CYCLE WHEN TAU=0.95 *****
C ——————————————————————————————————————————————————————————
c
IF(K.LE.(NN-2)) THEN
IF((RATEQ(K+1).GT.RATEQ(K)).AND.
1 (RATEQ(K+1).GT.RATEQ(K+2)))THEN
QMAx = RATEQ(K+1)
TMAx = TI+DT
GO TO 450
ENDIF
ENDIF
c
C ------------------------------------------------

C ****** CHECK THE DERIVATIVE OF EQ. 3.1 *********

IF(K.EQ.(NN-1)) THEN
SLOPE = (RATEQ(K)-RATEQ(K+1))/DT
WRITE(*,22)SLOPE
22 FORMAT(5X,'SLOPE =',F15.7)
IF(SLOPE.LT.0.DO) THEN
K = ((NN-5)/2)+5
GO TO 5000
ELSE
K = K-l
GO TO 5000
ENDIF
ENDIF

K = K-l
GO TO 5000

ENDIF

C ***** CALCULATING THE STREAM DEPLETION ********
C ***** FOR CYCLE-AVEREGE CONTINUOUS PUMPING *******

C
450 X = DSQRT((TA)/(4.0*TMAX))
C
CALL ERRFNC (X,ERFCX)
C

QTD = Q*TP/TD
QCl = ERFCX*QTD

00000 0 0 0

184

185

186

187

255

256

****** THE ERROR BETA IN EQUATION 11. ******

148

x = DSQRT((TA)/(4.0*TE))
CALL ERRFNC (X,ERFCX)

QCZ = ERFCX*QTD

BETA1 = 1.0-QC1/QMAx
BETA2 = 1.0-QC2/QMAX
TPTA = TP/TA
TDTP = TD/TP

ENDIF

CONTINUE

IF(EUNIT.EQ.1) THEN
WRITE(S,184)TMAX,RATEQ,QC1
WRITE(*,184)TMAX,RATEQ,QC1
FORMAT(10X,'At t=',F18.0,' days',2X,
1'qmax =',D15.10,2X,'cu.ft/day'
2,/,40X,'qc =',D15.10,2X,'cu.ft/day')
WRITE(5,185)TE,QC2
WRITE(*,185)TE,QC2

FORMAT(10X,'At te=',F18.0,' days',2X,
1'qc =',D15.10,2X,'cu.ft/day')
ELSEIF(MUNIT.EQ.0) THEN
WRITE(S,186)TMAX,QMAX,QC1
WRITE(*,186)TMAX,QMAX,QC1
FORMAT(10X,'At t=',F18.0,' days',2X,
l'qmax =',D15.10,2X,'cu.m/day',
2/,40X,'qc =',D15.10,2X,'cu.m/day')
WRITE(5,187)TE,QC2
WRITE(*,187)TE,QC2

FORMAT(10X,'At te=',F18.0,' days',2X,
1'qc =',D15.10,2X,'cu.m/day')

ENDIF

WRITE(5,255)BETA1,TPTA,TDTP
WRITE(*,255)BETA1,TPTA,TDTP
FORMAT(10X,'BETA1 =',F10.6,5x,'tp/ta =',F10.3,
15X,'TD/TP =',F10.3)
WRITE(5,256)BETA2,TPTA,TDTP
WRITE(*,256)BETA2,TPTA,TDTP
FORMAT(10X,'BETA2 =',FlO.6,5X,'tp/ta =',F10.3,

149

15X,'TD/TP =',F10.3)

WRITE(6,257)TPTA,BETA2
257 FORMAT(1OX,F10.4,5X,F10.5)

WRITE(7,191)TP,TTA,TAU
WRITE(*,191)TP,TTA,TAU

191 FORMAT(10X,'Tp =',F5.1,1X,'days',5x,
l't/ta =',D20.12,5X,'TAU =',020.12)

WRITE(*,888)
888 FORMAT(/,10X,'CALCULATE THE qmax',/
*,10X,'FOR ANOTHER TIME (enter ANSWER) ?',/

*,10X,'TYPE O < ----------- YES',/
*,10X,'TYPE 1 < ----------- NO',//
*.10X.'PLEASE SPECIFY .............. >',/)
READ(*,*)ANSWER

IF(ANSWER.EQ.0) THEN
WRITE(*,*)' NEW DISTANCE BETWEEN WELL AND RIVER ?'
READ(*,*) DIST
GO To 99
ELSEIF(ANSWER.EQ.1) THEN
WRITE(*,*)' PROGRAM TERMINATED'

GO TO 999

ENDIF

C

999 CONTINUE
STOP
END

C

C

C **********************************************

SUBROUTINE ERRFNC (X,ERFCX)
***************************

***********************************************

THIS SUBROUTINE COMPUTES THE VALUE OF ERROR
FUNCTION FOR A GIVEN ARGUMENT X BY USING
THE RATIONAL POLYNOMIAL APPROXIMATIONS

 

 

00000000

IMPLICIT REAL*8 (A-H,O-Z)
DIMENSION CONST(5)

CONSTANT TERMS IN THE APPROXIMATING FORMULA

--’---------------------------------. .........

NCONST = 5

PT = 0.3275911

CONST(1) = 0.254829592

CONST(2) = -0.284496736
CONST(3) = 1.421413741

CONST(4) = -1.453152027
CONST(5) = 1.061405429

APPROXIMATION OF ERROR FUNCTION X (erf X) BY
RATIONAL APPROXIMATIONS

SUM = 0.D0
PARM = l.0/(l.0+PT*X)
COEFF = EXP(-(X*X))

Do 500 M=1,NCONST
AI = FLOAT(M)
SUM = SUM+CONST(M)*(PARM**(AI))

500 CONTINUE

ERFX = 1.0-(SUM*COEFF)
ERFCX = 1.0-ERFX
RETURN

END

APPENDIX C

LISTING OF COMPUTER PROGRAM
FOR
STREAM DEPLETION BY A PUMPING WELL INCLUDING THE EFFECTS

OF LINEAR VARIATION OF WETLAND ET

151

OOOOOOOOOOOOOOOOOOOO

152

PROGRAM TO CALCULATE STREAM DEPLETION AND ET REDUCTION PRODUCED
BY CONTINUOUS PUMPING WELLS LOCATED IN HOMOCENEOUS, AND
ISOTROPIC AQUIFER.

WRITTEN AND DEVELOPED

l l
| I
| l
| l
l l
| l
I By I
| Yakup Darama Ph.d Candidate |
I l
| Michigan State University |
| College Of Engineering |
l l
l |
l l
I l
| I

Department Of Civil and Environmental Engineering

East Lansing, Michigan 48824

......................................................................

IMPLICIT REAL*8(A-H,O-Z)

CHARACTER FLNAME*10,FLIN*IO

WRITE(*,*) ' INPUT FILE NAME - ?'

READ(*,'(A)')FLIN

OPEN(UNIT-5,FILE-FLIN,STATUS-'OLD')

WRITE(*,*) ' OUTPUT FILE NAME FOR LEAKY DEPL. - ?'
READ(*,'(A)')FLNAME
OPEN(UNIT-6,FILE-FLNAME,STATUS-'NEW',FORM-'FORMATTED')
WRITE(*,*) ' ENTER THE OUTPUT FILE NAME FOR CAPTURED ET - ?'
READ(*,'(A)')FLNAME
OPEN(UNIT-7,FILE-FLNAME,STATUS-'NEW',FORM-‘FORMATTED')
WRITE(*,*) ' ENTER THE OUTPUT FILE NAME FOR MINED WATER - ?'
READ(*,'(A)')FLNAME
OPEN(UNIT-8,FILE-FLNAME,STATUS-'NEW',FORM-'FORMATTED')

. INPUT VARIABLES

000

DX - 1.0
PI - 4.0*DATAN(DX)

CALL INPUT(STOR,TRANS,DIST,Q,PET,DT,SOILN,RHOB,SANDM,
* SATK,DTPARM,EPS)
CALL DEPTH(PI,SOILN,RHOB,SANDM,SATK,PET,EPS,DO)

15

16

153

WRITE(*,*) 'D0 -',D0

TI - 0.00

RATEQ - 0.00

RATEQL - 0.00

RATEQS - 0.00

ETRATE - 0.00

QRQET - RATEQL + ETRATE

TOL - O.9999*Q

ENDTIM - 17S

IF(QRQET.LE.TOL) THEN
0T - 0T * DTPARM
TI - TI + 0T
00 - SDF(DIST,STOR,TRANS,TI)
TA - ATA(DIST,STOR,TRANS)
ALPHA - DIST/ALEAK(TRANS,DO,PET)

CALL BFLUX1(Q,TI,TA,U0,RATEQ)

CALL BFLUX2(Q,TI,TA,UO,ALPHA,RATEQL)
CALL REDST0(Q,TI,U0,ALFHA,RATEQS)
CALL ETCAP(TI,Q,RATEQS,RATEQL,ETRATE)

QRQET - RATEQL + ETRATE
WRITE(*,16)TI,QRQET
F0RMAT(5X,F10.3,5x,F12.6)
Go TO 15

ENDIF

STOP
END

C WWWMW”***************

DOUBLE PRECISION FUNCTION ALEAK(TRANS,DO,PET)
*********************************************

IMPLICIT REAL*8(A-H,O-Z)

. CALCULATION OF PARAMETER ALEAK FOR LINEAR ET CAPTURE

ALEAK - DSQRT((TRANS*DO)/PET)

RETURN
END

154

C *******************************************************************

SUBROUTINE BFLUX1(Q,TI,TA,UO,RATEQ)
c ***********************************

C SUBROUTINE TO CALCULATE BOUNDARY FLUXES WHEN ET CAPTURE IS NOT
C CONCERNED.

c ....................................................................
C
IMPLICIT REAL*8(A-H,O-Z)
C
XIN - U0
CALL ERRFNC (XIN,ERFCX)
RATEQ - Q * ERFCX
TITAl - TI/TA
C WRITE(*,22)TI,RATEQ
C WRITE(7,22)TI,RATEQ
C22 FORMAT(5X,F10.3,5X,F10.4)
C
RETURN
END
C
C
C .
C **********************************************************************
SUBROUTINE BFLUX2(Q,TI,TA,U0,ALPHA,RATEQL)
C ******************************************
c ......................................................................
C SUBROUTINE TO CALCULATE BOUNDARY FLUXES WHEN ET CAPTURE IS CONCERNED
c ......................................................................
C

IMPLICIT REAL*8(A-H,O-Z)

BETA - no - ((ALPHA/2.DO)/UO)
CALL ERRFNC (BETA,ERFCX)
TERMl - (1.00/EXP(ALPHA)) * ERFCX
ETA - U0 + ((ALPHA/2.DO)/U0)
CALL ERRFNC (ETA,ERFCX)
TERM2 - EXF(ALFHA) * ERch
RATEQL - (Q/2.DO) * (TERMl + TERMZ)
TITA2 - TI/TA
DIMQRT - RATEQL/Q

c WRITE(*,21)TI,RATEQL
WRITE(6,21)TI,RATEQL .

21 F0RMAT(sx,F10.3,5x,F12.6)

RETURN
END

155

c t**************************************************************
SUBROUTINE REDSTO(Q,TI,UO,ALPHA,RATEQS)

C SUBROUTINE TO CALCULATE THE RATE OF MINED WATER FROM AQUIFER
C STORAGE

IMPLICIT REAL*8(A-H,O-Z)

XVAL - U0
CALL ERRFNC (XVAL,ERFCX)
ERFUO - 1.00-ERFCX
HALPUO - - (0.5 * ALPHA/UO)**2
RATEQS - Q * EXP(HALPUO) * ERFUO
c WRITE(*,24)TI,RATEQS
' WRITE(8,24)TI,RATEQS
24 F0RMAT(5x,F10.3,5x,F12.6)

C
RETURN
END
C
C
C
c ***************************************************************

SUBROUTINE ETCAP(TI,Q,RATEQS,RATEQL,ETRATE)
*******************************************

IMPLICIT REAL$8(A-H,O-Z)

C
ETRATE - Q - (RATEQL + RATEQS)
C AVET - Q - ETRATE
C
C WRITE(*,26)TI,ETRATE

WRITE(7,26)TI,ETRATE
26 FORMAT(5X,F10.3,5X,F12.6)

RETURN
END

 

156

C HMWMMWWWWWW**

SUBROUTINE DEPTH(PI,SOILN,RHOB,SANDM,SATK,PET,EPS,DO)
c *****************************************************

C SUBROUTINE TO CALCULATE DEPTH do WHERE EVAPOTRANSPIRATION FLUX IS
C NEGLIGIBLE

C ......................................................................
C
C
IMPLICIT REAL*8(A-H,O-Z)
C
C .... Soil parameters
0
DC - SANDM
SIGMAG - 1
C
C .... Air Entry Suction SE in cm of water
C

PSIES - -0.S/DSQRT(DG)
SLOPE - -2.0 * PSIES + 0.2 * SIGMAG
POWl - 0.67 * SLOPE
PSI - PSIES * (RHOB/1.3)**POW1
SE - -10.2 * PSI
WRITE(*,121)SE
121 FORMAT(SX,’ SE - ',F10.5)
APP - PET * EPS

R - APP/(SATK + APP)

...Computing the depth do where ET - 0.

00C)

PIN - PI/SOILN

COSEK - 1.00/SIN(PIN)

AN - (PIN * COSEK)**SOILN

ANOM - sATx * AN

00 - SE * ((ANOM/APP)**(1.0/SOILN) - R)
160 WRITE(*.113)DO
113 F0RMAT(SX,'dO -',FlO.3)

RETURN

END

157

C WMW****H***H******************

200

SUBROUTINE ERRFNC (X,ERFCX)
***************************

IMPLICIT REAL*8(A-H,O-Z)
DIMENSION CONST(5)

REAL AI

. CONSTANT TERMS IN THE APPROXIMATING FORMULA

NCONST - 5

PT - 0.3275911

CONST(1) - 0.254829592
CONST(2) - -O.284496736
CONST(3) - 1.421413741
CONST(h) - -1.453152027
CONST(5) - 1.061405429

.. APPROXIMATION OF ERROR FUNCTION X BY RATIONAL APPROXIMATIONS

IF(X.LI.0.DO)THEN
Y - DABS(X)
ELSE
Y - x
ENDIF
SUM - 0.00
EARN - 1.0/(1.0+PT*Y)
COEFF - EXP(-(Y*Y))
DO 200 Mp1,Nc0NST

AI - FLOAT(M)

SUM - SUM+CONST(M)*(PARM**(AI))
CONTINUE
IF(X.LI.O.DO)THEN

EREX - 1-(SUM*COEFF)
ERFCX - 1.0+EREx
ELSE
EREx - 1.0-(SUM*COEFF)
ERFCX - 1.0-EREX
ENDIF

RETURN
END

158

c *****************************************************************
SUBROUTINE INPUT(STOR,TRANS,DIST,Q,PET,DT,SOILN,RHOB,SANDM,

* SATK,DTPARM,EPS)
c *********************************
c ......................................................................
C SUBROUTINE TO READ THE INPUT DATA FROM THE FILE
0 ......................................................................
C

IMPLICIT REAL*8(A-H,O-Z)
C
C .... AQUIFER PROPERTIES
C

READ(5,*) STOR

READ(S,*) TRANS

READ(5,*) DIST

READ(5.*) Q

READ(S,*) DT

READ(5,*) PET
C
C . .. SOIL (MATERIAL) PROPERTIES
C

READ(5,*) SOILN

READ(5.*) RHOB

READ(5,*) SANDM

READ(5,*) SATK
C
C .... CONTROL PARAMETERS
C

READ(5,*) EPS

READ(5,*) DTPARM
C

RETURN
END

159

c *************************************************************
DOUBLE PRECISION FUNCTION SDF(DIST,STOR,TRANS,TI)

c *************************************************
C .............................................................
C FUNCTION TO CALCULATE AQUIFER RESPONSE TIME
C .............................................................
C

IMPLICIT REAL*8(A-H,O-Z)
C
C

SDF - DSQRT((DIST * DIST * STOR)/(4.0 * TI * TRANS))
C .

RETURN

END
C

c *************************************************************
DOUBLE PRECISION FUNCTION ATA(DIST,STOR,TRANS)

c *********************************************
C IMPLICIT REAL$8(A-H,O-Z)

g .... CALCULATION OF AQUIFER RESPONSE TIME TA

: ATA - DIST*DIST*STOR/TRANS'

RETURN
END

LIST 01' REFERENCES

Anat, A., H. R. Duke, and A. T. Corey, 1965, "Steady Upward
Flow from Water Table", Hydrology Paper No. 7, Colorado
State University, Fort Collins, Colorado.

Bedell, D. J., and R. L. van Til, 1979, "Irrigation in
Michigan, 1977", Water Management Division, Department of
Natural Resources, Lansing, Michigan, 41 pp.

Bower. H., 1978. MW. Chapter 1. p- 1-11.

McGraw-Hill Inc. New York.

Bredehoeft, J. D., S. S. Papadopulus, and H. H. Cooper, Jr.,
1982, "Groundwater: The Water Budget Myth", In Scientific
Basis of Water Resources Management, p. 51-57, National
Research Council, Geophysics Study Committee, National
Academy Press, Washington D. C.

Brooks, R. H., and A. T. Corey, 1964, "Hydraulic Properties of
Porous Media", Hydrology Papers, No. 3, pp. 27, Colorado
State University, Fort Collins, Colorado.

Burns, A. W., 1983, "Simulated.Hydrologic Effects of Possible
Ground-Water and Surface-Water Management Alternatives in
and Near the Platte River, South-Central Nebraska", U.S.
Geologic Survey Professional Paper, 1277-6.

Campbell. 6. 8.. 1985. 521W. Transvort
Models for Soil-Plant Systems, Development in Soil

Science 14, Elsevier Science Pub. Co. Inc. New York.

Carslaw, H. S., and J. C. Jaeger, 1959, ggngyggign_gf_ng§;_in
Solids, Second Edition, Oxford University Press, London.

Chung. S 0.. 198$.W1Wmmemm
Wm Ph D Dissertation. Iowa

State University, Ames, Iowa.

Corey. A- T-. 1986. WWW
Media, Chapter 4, p. 131-162, Water Resources

Publications, Littleton, Colorado.

Corey, A. T., 1990, ”Personal communication".

160

161

Duke. H- R., 1965. uaxiEum_Bate_9f_nnuard_£leg_fren_fleter
Tables, Masters Thesis, Department of Agricultural
Engineering, Colorado State University, Fort Collins,
Colorado.

Duke, H. R., 1972, "Capillary Properties of Soils-Influence
Upon Specific Yield", Transactions of the ASAE, 15, p.
688-691

Freeze, R..A., and.J. A. Cherry, 1979, figggggggggr, Chapter 2,
p. 29-61, Prentice-Hall, Englewood Cliffs, N.J.

Fulcher, G. W., S. A. Miller, and R. V. Til, 1986, "Effects of
Consumptive Water Uses on Drough Flows on the River
Raisin", Engineering-WaterfiManagement.Division, Michigan
Department of Natural Resources, Lansing, MI.

Gardner, W. R., 1958, "Some Steady-State Solutions of the
Unsaturated Moisture Flow Equation with Application to
Evaporation from a Water Table", Soil Science Vol. 85
p. 223-232.

Gardner, W. R. , and Fireman, Milton, 1958, "Laboratory Studies
of Evaporation from Soil Columns in the Presence of a
Water Table", Soil Science, Vol. 85, No. 5.

Gillham, R. W., 1984, "The Capillary Fringe and Its Effect on
Water-Table Response", Journal of Hydrology, 67: 307-324

Glover, R. E., and Balmer, C. G., 1954, "River Depletion
Resulting from Pumping a Well Near a River”, American
Geophysical Union Trans., Vol. 35, pt. 3, p. 468-470.

Glover, R. E. , 1960, "Groundwater Surface Water Relationship",
A Paper at Groundwater Section of Western Resources
Conference, Boulder, Colorado, Colorado State University,
Paper CER60REG, 8 pp.

Hantush, M. S. , 1955, ”Discussion of River Depletion Resulting
from Pumping a Well Near a River", Trans. American
Geophysical Union, Vol. 36, p. 345-347

Hantush, M. S., and Jacob, C. E, 1955, "Non-Steady Radial
Flow in an Infinite Leaky Aquifer", Trans. American
Geophysical Union, Vol. 36, No. 1, p. 95-100.

Hantush, M. S., 1956, "Analysis of Data from Pumping Test in
Leaky.Aguifers", Trans. American Geophysical Union, Vol.
37, No. 6, p. 702-714.

Hantush M. S. , 1959, "Analysis of Data from Pumping Wells Near
a River", Journal of Geophysical Research, Vol. 64, No.

162
11, p. 1921-1932.

Hantush, M.S., 1963, ”Hydraulics of Gravity Wells in Sloping
Sands", Transactions of ASCE, v. 128, part 1, pp. 1423-
1442.

Hantush, M.S., 1964a, u s ' o
95L, Advances in Hydroscience, Vol. 1, New York, Academic
Press, p. 2829-2838

Hantush, M. S., 1964b, "Depletion of Storage, Leakage, and
River Flow by Gravity Wells in Sloping Sands", Journal of
Geophysical Research, Vol. 69, No. 12., p. 2551-2560.

Hantush, M. S., 1965, "Wells Near Streams with Semipervious
Beds”, Journal of Geophysical Research, Vol. 70, No. 12,
p. 2829-2838.

Hantush, M. S., 1967, "Depletion of Flow in Rigth-Angle Stream
Bends” , Water Resources Research, Vol . 3 , No . 1 p . 23 5-240 .

Hantush, Mohamed M. S. and Marino, M. A., 1989, "Chance-
Constrained Model for Management of Stream-Aquifer
Systems", ASCE Journal of Water Resources Planning and
Management, Vol. 115, No. 3, p. 259-277.

Henningsen, F. J., 1977, "Irrigation Expansion in Michigan",
Michigan State University Cooperative Extension Service,
East Lansing, Michigan.

Hillel. D- 1980. Wise. Academic Press.

New York.

Jenkins, C. T., 1968, ”Computation of Rate and Volume of
Stream Depletion by Wells", Techniques of Water-Resources
Investigations of the U.S. Geological Survey, Book 4,
Chapter D1, 17 p.

Jenkins, C. T., 1968b, "Electric-Analog and Digital-Computer
Analysis of Stream Depletion by Wells", Ground Water
v. 6, n. 6, pp. 27-34.

Kazmann, R. G. , 1948, "The Induced Infiltration of River Water
to Wells", Trans. American Geophysical Union, Vol. 29, p.
85-92.

Markar, M. S., and R. G. Mein, 1987, "Modeling of
Evapotranspiration from Homogeneous Soils", Water
Resources Research, Vol. 23, No. 10, p. 2001-2007.

McDonald, M. G., and A. W. Harbaugh, 1988, "A Modular Three-
Dimensional Finite-Difference Groundwater Flow Model",

163

Techniques of Water-Resources Investigations of the USGS,
Chapter Al, Book 6, Washington.

Merva, G., 1990, "Personal communication"

Moore, R. E. , 1939, "Water Conduction from Shallow Water
Tables", Hilgardia, Vol. 12, No. 6. p. 383-426

Moulder, E. A., and C. T. Jenkins, 1969, "Analog-Digital
Models of Stream-Aquifer Systems", Ground Water, v. 7,
n. 5, pp. 19-24.

Mueller, F. A., and J. W. Male, 1990, "A.Management Model for
Specification of Groundwater Withdrawal Permits. ", Paper
submitted to Water Resources Research, Department of
Civil Engineering, University of Massachusettes, Amherst.

Philip, J. R., 1957, "Evaporation, Moisture and Heat Fields in
the Soils", Journal of Meteorology, Vol. 14, p. 354-366.

Prickett, T. A. and C. G. Lonnquist., 1971, ”Selected Digital
Computer Techniques for Groundwater Resource
Evaluation," Illinois State Water Survey, Urbana,
Bulletin 55.

Ripple, C. D., J. Rubin, and T. E. A. van Hylckama, 1972,
"Estimating Steady State Evaporation Rates from Bare
Soils Under Conditions of High Water Table", U.S.
Geological Survey Water Supply Paper, 2019-A.

Rorabaugh, M. I. , 1956, "Groundwater in Northeastern
Louisville and Kentucky with Reference to Induced
Infiltration", U.S. Geological Survey Water Supply Paper,
1360-B, 169 pp.

Schleusener, R. A., and Corey, A. T., 1959, "The Role of
Hysteresis in Reducing Evaporation from Soils in Contact
with a Water Table", Journal of Geophysical Research,
Vol. 64, No. 4, p. 469-475

Skaggs, R. W., 1975, "Drawdown Solutions for Simultaneous
Drainage and ET”, Journal of Irrigation and Drainage
Division, ASCE, Vol. 101(IR4), p. 279-291.

Skaggs, R. W. and Y. K. Tang, 1976, "Saturated and Unsaturated
Flow to Parallel Drains", Journal of Irrigation and
Drainage Division, ASCE, Vol. 102(IR2), p. 221-238.

Skaggs, R. W. and Y. K. Tang, 1979, "Effects of Drain
Diameter, Openings and Envelopes on Water Table
Drawdown", Transactions of the ASAE, Vol. 22, n. 2 pp.
326-333.

164

Skaggs. R- We: 1978. W
W, Water Resources Research Institute of

the University of North Carolina, Report No. 134. 178 pp.

Spalding, C. P., and R. Khaleel, 1991, "An Evaluation of
Analytical Solutions to Estimate Drawdowns and Stream
Depletion by Wells", Water Resources Research, Vol. 27,
No. 4, p. 597-609.

Staley, R. W., 1957, t W e on
W. Masters ThesiS. Colorado

State University, Fort Collins, Colorado.

Taylor, 0. J., 1971, "A Shortcut for Computing Stream
Depletion by Wells Using Analog or Digital Models" ,
Groundwater, Vol. 9, No. 2, p. 9-11.

Theis, C. V., 1941, "The Effect of a Well on the Flow of a
Nearby Stream", Trans. American Geophysical Union, Vol.
22, pt. 3, p. 734-738.

Theis, C. V., and Conover, C. 8., 1963, "Charts for
Determination of the Percentage of Pumped water Being
Diverted from a Stream or Drain", U.S. Geological Survey
Water Supply Paper 1545-C, p. C106-C109.

8.. and J- L- Wilson. 1980. WW

0
v-.!.< !

-_ .- ! i . ‘
Aggifgm31, Ralph M. Parsons Laboratory f

    

Townley, L.
:- I

             

or Water

Resources and ' Hydrodynamics, Report No. 252,
Massachusetts Institute of Technology , Cambridge ,
Massachusettes.

Trescott, P. C., G. F. Finder and S. P. Larson, 1976, "Finite
Difference Model for Aquifer Simulation in Two
Dimensions with Results of Numerical Experiments",
Techniques of Water-Resources Investigations of the USGS,
Chapter 1, Washington.

Walton, W. C., 1962, "Selected Analytical.Methods for Well and
Aquifer Evaluation", Illinois State Water Survey
Bulletin, 49.

Warrick, A. W., 1988, "Additional Solutions for Steady State
Evaporation from a Shallow Water Table", Soil Science,
Vol. 146, No. 2, p. 63-66.

William, B. A., J. B. Miller, and W. W. Wood, 1973,
"Availability of Water in.Kalamazoo County, Southwestern
Michigan", U.S. Geological Survey Water Supply Paper,
1972.

HICH

IGRN STATE UNIV. LIBRRRIES
I “WINW“NWI?"IIIWWII‘IWIWHM
1293008967352

 

3