- Email: [email protected]
Flow in a vertical porous column drained at its bottom at constant flux
Journal of Hydrology, 24 (1975) 135--153 © North-Holland Publishing Company, Amsterdam - - Printed in The Netherlands
FLOW IN A VERTICAL AT CONSTANT FLUX
POROUS
COLUMN DRAINED
AT ITS BOTTOM
U. KROSZYNSKI
Department of Civil Engineering, Technion, Israel Institute of Technology, Haifa (Israel) (Accepted for publication March 24, 1974)
ABSTRACT
Kroszynski, U., 1975. Flow in a vertical porous column drained at its b o t t o m at constant flux. J. Hydrol., 24: 135--153. A soil column, saturated at its lower portion and initially in condition of hydrostatical equilibrium is drained at a constant flux at its bottom. Due to the presence of the unsaturated zone above the water table, the motion of the saturation surface depends in a complex manner on time and on soil properties. After defining a characteristic length, the equation of unsaturated flow is solved numerically by a finite difference scheme for different values of drained flux, for two similar soils. Although moisture and head profiles are different for both soils, drawdown curves have been found to be in good agreement during a certain time interval which depends on the value of the drained flux. The same problem is solved analytically by a linearization procedure. A simple expression which predicts saturation surface drawdowns quite accurately for small values of the drained flux and for moderate time intervals, is derived. An effective porosity coefficient is defined and its variation with time is represented. The bearing of the results on flow in unconfined aquifers is discussed.
INTRODUCTION In unconfined hydrological systems operating under drainage conditions (e.g., a q u i f e r p u m p i n g , d r a i n a g e b y t i l e d r a i n s ) w a t e r is r e l e a s e d m a i n l y b e cause of the lowering of the water table and at a much lesser extent from elastic storage. I t is c u s t o m a r y in t h e t h e o r y o f g r o u n d w a t e r f l o w t o a s s u m e t h a t t h e v o l u m e o f w a t e r d r a i n e d ( p e r u n i t a r e a ) is e q u a l t o t h e d r o p o f t h e w a t e r t a b l e t i m e s t h e e f f e c t i v e p o r o s i t y c o e f f i c i e n t n, w h i c h is g e n e r a l l y t a k e n as c o n s t a n t . A s it h a s l o n g b e e n r e a l i z e d , d u e m a i n l y t o t h e e x i s t e n c e o f a n u n s a t u r a t e d z o n e a b o v e t h e w a t e r t a b l e , t h e r e is a " d e l a y e d y i e l d " a n d t h e c o n s t a n t d r a i n a b l e p o r o s i t y is p r e s u m a b l y a t t a i n e d o n l y a s y m p t o t i c a l l y , f o r l a r g e v a l u e s of time. For this reason, some inconsistencies arise when interpreting field m e a s u r e m e n t s o r p r e d i c t i n g f o r m a t i o n p a r a m e t e r s o n t h e b a s i s o f c o n s t a n t n. I n t h i s c o n t e x t , in t h e f r a m e w o r k o f a s t u d y o n t h e i n f l u e n c e o f t h e u n -
136 saturated zone u p o n the aquifer drainage process, we have considered as a first step the one dimensional vertical flow problem. The drainage of vertical porous columns has been already considered by different authors, especially for the case where drainage is caused by a sudden head d r o p at the b o t t o m of an initially saturated column (constant head t ype b o u n d a r y condition). Youngs (1960) c o m p u t e d the variation of the drained volume by assuming a single diameter capillary tube model of the porous medium, which is a constant porosity approach. Gardner (1962), Fujioka and Kitamura (1964) and Sw a r t zendr ube r (1969) have assumed that the interface between the saturated and unsaturated zones, or saturation surface (S.S.) drops instantaneously from its initial position at the column t op to its ultimate one near the column b o t t o m and solved the equation of unsaturated flow by linearizing it. Vachaud (1968) and Jackson and Whisler (1970) c o m p u t e d the variation of the drained volume by generalizing Youngs' model namely by considering a bundle of capillary tubes of different diameters. Numerical solutions of the drainage problem under constant head b o u n d a r y condition at the column b o t t o m were given by Day and Luthin (1956), Jensen and Hanks (1967), Watson (1967), Whisler and Watson (1968) and Hornberger and Remson (1970). Som ew ha t different b o u n d a r y conditions were also considered by Freeze (1969). All those works were based upon finite difference schemes. A numerical simulation of the physical process of drainage at constant head at the column b o t t o m was presented by Kastanek (1971). The present study differs f r om those m ent i oned before in two aspects. (1) We consider the case of drainage under constant flux at the column b o t t o m rather than constant head. Such a study is relevant to a few applications where the aquifer is drained at constant discharge as a b o u n d a r y condition. Moreover, in the case of a constant head drop, the drainage process depends on the initial length of the saturated zone while this length is immaterial in the case of constant flux. This way, general conclusions may be reached at the expense of fewer computations. (2) The analytical m e t h o d we suggest considers variable coefficients in the equation of unsaturated flow rather than constant ones, which may yield a more accurate solution. Childs and Poulovassilis (1962) solved a problem closer to this context. T h e y considered the case where bot h water table and moisture c o n t e n t profile translate at a constant dow nw a r d (upward) speed. Although exact, their solution might perhaps apply asymptotically after a large time interval. Even in this case t h e y clearly show that the moisture profile c a n n o t descend at constant speed in the zone of low moisture content. STATEMENT OF THE PROBLEM AND QUALITATIVE ANALYSIS OF THE DRAINAGE PROCESS We consider a vertical c ol um n of h o m o g e n e o u s and stable soil. The column is initially saturated with water up to a certain level z = 0 (Fig.l) with an un-
137
/~ q=O
kl
-,6--
~Z
\)\
t--ol
',\',t~o 0 \
"~o t~,e
S(t~4 -
e~s t=O t-O,
(a) (b) Fig.1. Drainage of a vertical c o l u m n : (a) the m o t i o n of the saturation surface, and (b) distribution of 0 and ~ initially and at time t (qualitative).
saturated portion above it in hydrostatical equilibrium. The initial moisture profile is assumed to correspond to the 0 (~) main drainage curve. Starting at t = 0 water is withdrawn from the column b o t t o m at a constant rate q0 and the flow is permanently downwards. There is no flux through the column top. The initial and boundary conditions are such that the soil undergoes drainage at every point and consequently the 0 ($) relationship remains univalued. By continuity, q0 will also be the flux at the saturation surface S.S. This fact enables us to deal with the unsaturated zone solely, provided a proper (constant flux) boundary condition is stated at the S.S., thus making the saturated length of the column immaterial. The column is, moreover, assumed to be sufficiently long so that in the time interval to be considered, its b o t t o m remains at saturation while its top is always at low moisture content. The aim of the present work is to determine the motion of the saturation surface S.S. (z = --s(t)), which is of paramount importance in hydrologic applications. Under the considered b o u n d a r y conditions this is almost identical to the motion of the phreatic surface. If the flow in the unsaturated zone is neglected, as usually done in existing approaches or alternatively, if it is assumed t h a t the moisture profile translates without change of shape after the S.S., the solution to the problem obviously becomes: SSAW(t) =
qo t / n
(1)
where n is the so called drainable porosity, n can therefore be defined from eq.1 as the ratio between the volume of water released from the soil and the volume of soil spanned by the S.S. (See also Appendix C.)
138
Actually, there is a delay in the release of water from the unsaturated zone and the S.S. moves downwards faster than SsnT(t ) predicted by eq.1. As a result, the effective porosity: t/s(t)
neff = qo
(2)
is smaller than the drainable porosity of eq.1. Moreover, neff depends in a complex manner on the unsaturated properties of the soil, the discharge q0 and time. Its computation is one of the aims of the present study. We define now the problem in precise mathematical terms. The equation of unsaturated flow is written in terms of the head $. Above the S.S. we have:
~ ( ~ ) Ot
~z
(~)
+~(~)~z
where a(~) = d0/d~ and ~($) = dK/d~. The boundary and initial conditions are: (i) Zero flux through the column top: .....
~z
1
(4)
(z=H,t>~O)
(ii) Continuity of head and flux at the S.S.: q o = K s -(~ z~ )+ 1
; ~=~e
(z=--s(t),t>O)
(5)
(iii) The initial hydrostatical equilibrium distribution: =--z
(H > z /> 0, t = 0)
+ ~c
(6)
Eq.3 with conditions 4, 5, 6 allow only in principle to determine ~ (z, t) within the unsaturated portion of the column and the S.S. location s(t). Two major difficulties make a closed analytical solution almost impossible: the nonlinearity of eq.3 and the fact that the lower boundary z --- -- s(t) is an u n k n o w n of the problem. On the other hand the solution within the saturated portion of the column is readily given by: =
--
--
x (z + s(t))
(z < - - s ( t ) )
(7)
It is worthwhile to clarify two points at this stage. First, when dealing with a saturation surface which is sufficiently deep below the top, it is mathematically convenient and physically sound to transfer the upper boundary from the column top to z = H ' (Fig.2) where H ' is selected such that the hydraulic conductivity is negligibly small for z/> H ' (we have selected H ' in the examples presented in the sequel such that K(O L )/Ks = 1 0 - 4 (see Fig.2)). It seems reasonable to assume that practically no flow will take
139
Zj
H.z~ COLUMNTOP
Fig.2. Location of the upper boundary and the characteristic length for the unsaturated zone.
place in th at region in the time interval of interest. Moreover, this transfer makes the initial moisture c o n t e n t at the upper b o u n d a r y 0 = 8 L (Fig.2) a well defined q u a n t i t y rather than the irreducible moisture c o n t e n t 0 r which is some° tim~s difficult to detect. In any case we can define a characteristic length of the problem L as the smaller of the t w o quantities H and H ' and take b o u n d a r y conditio.l 4 at z = L. Secondly, on purely dimensional grounds, we may write the solution of eq.3 for a given soil as: ~ -function L
z
q0
Kst
' L
1>~-~>-L L
'
(8)
reducing thus the n u m b e r of variables. Since for q0 = 0 there is no flow and t h er ef o r e $ = $0 = -- z + ~c and s ( t ) = 0, it is legitimate to assume that for small values o f e = q o / K s , eq.8 m ay be written as: ~/L = ~o/L
+ e$ I/L + e2~2/L
+ . . .
(9)
Similarly for $ = $c in eq.8: s(t)/L
= esl(t)/L
+ e2s2(t)/L
+...
(10)
where ~ o / L + e ~ I / L and e s l ( t ) / L represent the "first order solution", obtained by neglecting quadratic and higher order terms in eqs.9 and 10. We now proceed to solving the pr obl e m numerically for two soils and different values of the parameter e. A simple a p p r o x i m a t e analytical solution based on the a f o r e m e n t i o n e d points is derived subsequently. Soil 1 is a river sand whose r e t e n t i o n and hydraulic c onduc ti vi t y curves (definition curves) have been taken f r om the publication by Vachaud (1968). Soil 2 is a hypothetical sand with the same 0s, 0r, $c and K s of soil 1 but whose definition curves follow exact analytical expressions which serve as basis for the
140
analytical solution (see eqs.24, 25). Such a choice has a twofold purpose: to study the influence of the difference of shape of the definition curves between both soils other things being equal and to appraise the validity of the "first order solution" {derived analytically) by comparing it with the numerical solution for soil 2. For the sake of completeness the definition curves of both soils are presented in Fig.3. This study will be restricted to comparison between numerical and analytical solutions solely since no experimental work dealing with constant flux boundary conditions has been done so far. -~cmofwater
60
g ecru/rain)
I --
~s: 0.342
--"t -
81-: 0.075 Ks= 1.09cm/min U)c:-16.5crn.
/4
soil 2
/ /
1.0
/
I
a= oDgcrn-'
4O
08 \
!
,,
/
3O
0£
\\
12 04,
20
10
S
A
i I
02
0 0
10
20
30
eC%)
Fig.3. Retention and conductivity curves for soil 1 (after Vachaud, 1968) and soil 2.
NUMERICAL EXPERIMENTS The problem formulated in eqs.3, 4, 5, 6 has been solved numerically by a finite difference scheme for both soils. The scheme used is that of Hornberger and Remson (1970) adapted to our boundary conditions. The scheme proved to be accurate when dealing with the constant head boundary condition drainage problem. The flow domain -- s ( t ) <~ z <. L is first mapped by the transformation ~ = ( L - - z ) / ( L + s ( t ) ) onto the interval 0 ~< ~ ~< 1. The differential equation and
141 b o u n d a r y c o n d i t i o n s are r e w r i t t e n in t e r m s o f t h e n e w variables ~ a n d K s t / L and t h e algebraic s y s t e m resulting f r o m t h e finite d i f f e r e n c e a p p r o x i m a t i o n is solved b y t h e D o u g l a s - J o n e s ( 1 9 6 3 ) p r e d i c t o r c o r r e c t o r t e c h n i q u e . T h e locat i o n a n d v e l o c i t y o f t h e S.S. are d e t e r m i n e d iteratively, t h e i t e r a t i o n p r o c e d u r e i m p o s i n g t h e r e l a t i o n s h i p b e t w e e n t h e space a n d t i m e i n c r e m e n t s . T h e acc u r a c y o f t h e s o l u t i o n was c h e c k e d a m o n g o t h e r criteria b y t h e m a t e r i a l balance, i.e. b y c o m p a r i n g t h e d r a i n e d v o l u m e q0 A t w i t h t h e d i f f e r e n c e bet w e e n t h e m o i s t u r e profiles i n t e g r a t e d o v e r z at t and t + A t. T o ensure acc u r a c y it was f o u n d t h a t q u i t e a d e n s e space grid is n e e d e d . T h e r e a s o n f o r t h a t is t h e s t e e p v a r i a t i o n o f t h e $ (and 0) profiles n e a r t h e S.S. As t i m e elapses t h e profiles b e c o m e s t e e p e r in t h e n o r m a l i z e d variable ~, n e a r ~ = 1. On the o t h e r h a n d a sparser m e s h c o u l d be used in t h e u p p e r p o r t i o n o f the uns a t u r a t e d zone, w h e r e a " t a i l " o f a l m o s t c o n s t a n t l o w m o i s t u r e c o n t e n t is d e v e l o p e d . H o w e v e r , in o r d e r to a v o i d c o n s i d e r a b l e c o m p l i c a t i o n s , a fixed 60 n o d e m e s h was used f r o m t h e beginning. T h e t i m e interval was selected as 60 m i n , a n d t h e values o f e = q o / K s used were 0.05, 0.1, 0.2. A c c o r d i n g t o t h e c r i t e r i o n K ( O L ) / K s = 10 --4 (see p r e c e d ing section), L was t a k e n 25.6 cm. Curtailing t h e e x t e n t o f t h e u n s a t u r a t e d z o n e at z = L saved c o m p u t e r t i m e w i t h o u t a f f e c t i n g a c c u r a c y . A n u m e r i c a l e x p e r i m e n t w i t h L = 40 c m and a m e s h o f 80 n o d e s was p e r f o r m e d f o r soil 1 yielding t h e s a m e results as t h o s e s h o w n in full lines in Fig.4 w h i c h c o r r e s p o n d t o L = 25.6 c m a n d 60 nodes. In t h e case s h o w n in Fig.4 f o r soil 1 t h e I.B.M. 3 6 0 / 7 0 T e c h n i o n c o m p u t e r p e r f o r m e d t h e n e a r l y 1 4 0 0 t i m e steps n e e d e d in a b o u t 102 seconds, c o m p i l a t i o n included. Z(cm) I i~, ,, 'i',
'
,o r I
4::0.0 5
xL'"IL,'
IZ i
20
\ \x
\\\\ 10 ~
\
l
\\ \\\\\x \'\X\ ,x\\
t=o iI
\\
i 0'-
t=O t:~O mI'
t:3s.,.
_1o
[,
10
t=60~,~ 20 "8-(%)
(cm)
30
-40
\\
h" XNX} ~
""~,
0 ',
~-io \,t =E,O i \'~I
-30 -2~0 -16 5 ~J(cm orwater)
Fig.4. Moisture content and head distributions at different values of time for e = 0.05 (numerical solution).
142
For details on the numerical scheme the reader is referred to Hornberger and Remson (1970) and to Appendix A. The results, namely the $ (z, t) and 0 (z, t) profiles for both soils are pre* sented in Figs.4, 5 and 6. Although the trends are similar the profiles are quite different for both soils. However, it should be noted that, as time increases, the profiles have a tendency to translate with the S.S. except for the "tail" zone of low moistures. Moreover, the motion of the S.S. appears to be little influenced by the shape of the definition curves at least for the times represented.
20
-E=0.1
~
!20
o
t=
\ 10
.',XC °
20 @(%)
30
t = 29.~ -40
~ 30 -20 -16 5 $ (cm ofwater)
Moisture c o n t e n t and head distributions at different values o f time for • = 0 . 1 (numerical s o l u t i o n ) . Fig.5.
Z(cm) \ \
Z(cm)
20 ~
{=(]2 ~
,~
soil 1
20 \\
\\\ \~,~':
10
10 t=0
°
\
-10
--10
"X
~ ' ~ t=5 mitt
~ ' ~ ' 20 30 e (%)
t=11 rain. -40
",Ix:
\\ \ \~-10
-,.t0 -20 -16.5 ~ (¢m. of water )
Fig.6. Moisture c o n t e n t and head distributions at different values o f t i m e for • = 0.2 (numerical solution).
143 T h e m o s t i n t e r e s t i n g p a r a m e t e r , n a m e l y t h e S.S. d r a w d o w n s(t) is represented in t e r m s o f d i m e n s i o n l e s s variables in Fig.7. s(t)/L has b e e n divided b y e in o r d e r t o c h e c k t h e v a l i d i t y o f t h e first o r d e r a p p r o x i m a t i o n w h i c h implies t h a t sl (t)/L b e c o m e s i n d e p e n d e n t o f e (see e q . 1 0 ) . T h e d i m e n s i o n l e s s t i m e ~ = Kst/nL (n = 0 s - - 6r) selected as t h e i n d e p e n d e n t variable m a k e s possible t o c o m p a r e b o t h soils on a c o m m o n basis. M o r e o v e r , t h e c o n s t a n t p o r o s i t y solut i o n w h i c h is t h e usual o n e a d o p t e d in h y d r o l o g y yields t h e straight line SsAw(t)/cL = r (see eq.1).
J
:
A
0
10-J
'
2
.......
4
---s0il 2 . . . . (E q 29 or 31 )
~tur solut,on(Eq1/
6
8
%=Kst/nL
E-~:],2/ / ¢~[0,I / > ~ J ~ ~,/~"~ J " ~ ' " .~ ~ ~ ~.~.-r~
f-.......
0k,~rd/P'4"..... 0
2
4
6
8
T'=K s t / n L
B
Fig.7. The motion o£ the saturation surface z = -- s(t).
T h e closeness b e t w e e n t h e d r a w d o w n curves f o r e = 0.05 a n d e = 0.1 suggests t h a t t h e first o r d e r s o l u t i o n ( n a m e l y t h e first t e r m in eq.10) w o u l d p r o v i d e a fair a p p r o x i m a t i o n f o r soil 1, f o r e < 0.1 until r ~-- 5 while f o r soil 2 o n l y until r ~ 2. N o n - l i n e a r t e r m s b e c o m e r e l e v a n t o u t s i d e t h a t range, as s h o w n b y the d i f f e r e n c e w i t h t h e d r a w d o w n curve c o r r e s p o n d i n g t o e = 0.2 o v e r alm o s t t h e entire t i m e range. O n t h e o t h e r h a n d b o t h soils have f o r e = 0 . 0 5 v e r y similar d r a w d o w n curves o v e r 0 ~ r ~ 10 t h u s suggesting t h a t in t h a t case or f o r smaller values o f e, t h e s h a p e o f t h e d e f i n i t i o n curves d o e s n o t a f f e c t d r a w d o w n s . F u r t h e r i n s p e c t i o n o f Fig.7 s h o w s t h a t at t h e b e g i n n i n g o f t h e d r a i n a g e pro-
144
cess, say for 0 < r < 1, the S.S. descends m u ch faster, and at a variable speed, as co mp ar ed with the constant por os i t y drawdown. This is precisely the interval in which the unsaturated zone delays the release of water and as a result s ( t ) becomes larger than sSAT. Later, the speed of descent becomes nearly constant corresponding to the phase in which most of the unsaturated profile almost translates behind the S.S. These two regimes are illustrated even be t t er by Fig.8 where the variation of the effective porosity coefficient neff (eq.2) with time is represented. The fact that n eff does n o t reach the ultimate value n = 0 s - 0r during the time interval considered, is explained by the existence of a certain volume of water retained in the " t a i l " o f the moisture profiles. Finally, it is seen t hat interpretation of actual draw dow n curves in order to derive a constant hypot het i cal porosity n = qo t / s m ay lead to quite small values of n, especially when drawdowns are measured shortly after the beginning o f drainage. APPROXIMATE ANALYTICAL SOLUTION The numerical solution can be used in given particular cases and it seems that its extension t o two- and three-dimensional cases is difficult and expensive. F o r these reasons it was considered useful to derive an a p p r o x i m a t e analytical solution of the drainage problem. Such a solution is also motivated n : (126 7 c
A
0
2
[
-
[
:
4
6
~
~
.....
8
T=P,: b ' n ~
n=0.267'
(Eqs 2 a n d 2 9 o r 3 1 )
of B
Fig.& Variation of the effective porosity coefficient defined in eq.2.
145 by the fact that for e = q o / K s sufficiently small, as encountered in application, the "first order solution" in eq.10, represents quite accurately the numerical solution (Figs.7, 8). The proposed m e t h o d takes advantage of this property. We derive the linearized approximation by substituting expansions 9 and 10 in the exact equations of flow 3 to 6, expanding the a, K and ~ coefficients in Taylor series near C0 and neglect terms of order e2 or higher. Collecting terms independent of e we get the zero order approximation to C representing the initial state of rest. Co = - - z + Cc
(11)
Collecting now the terms of order e we obtain the following set for C l (z, t) and sl (t) (for details see Dagan and Kroszynski (1971)). ao (z) OC~ - a [K o (z)
at1 ]
(L >~z ~> 0)
(12)
- 0
(z = L)
(13)
- 1
(z = 0)
(14)
o
(t = o)
(151)
at
az
~zJ
aC~ az
at1 az ~
--
s~ ( t ) = - - ~ ~ ( O , t )
(161)
where in eq. 12 a0 (z) = a (C 0 ) and K0 (z) = K(C0 ) are known functions of z based on the initial state. In contrast with usual linearization procedures where one or both these coefficients are taken as constant, in the present case t h e y are variable. The m e t h o d used is known as m e t h o d of small perturbation expansions. It is widely used in different branches of physics and has been applied before to saturated groundwater flow problems (Dagan, 1964, 1966, 1968). Eqs.12 to 15, satisfied by C~, are much simpler than the exact original ones with two respects: (a) rather than the non-linear eq.3,we have now the linear eq.12; and (b) the boundary conditions (5) applying to the moving boundary z = - - s ( t ) a r e transferred to the fixed z = 0. Still c o m p u t a t i o n of the drawdown is possible by means of eq.16. Considering the function : 1--
g(z,t)=C~(z,t)+L
(17)
2
we obtain for Y by substitution in eqs.12--15 the following set: ao(Z) at - az
Ko (z)
+ 1 --
(L/> z >/0)
(18)
146
3Y 3z
-0
Y=-2
(z = 0; z = L)
(19)
(t = O)
(20)
1--
Separation of variables in the form: g =~
(21)
ui(z) G i ( t )
1
yields ui(z) as e igen-functions of the Sturm-Liouville equation: dz
Ko
dz
+ OiaoU i
0
(L ~> z ~> 0)
(22)
(z = 0;z
(23)
with b o u n d a r y conditions: dui dz
- 0
= L)
Eq. 22 can be solved approximatively to a desired accuracy by Galerkin's m e t h o d (this has been done in our previous work, Dagan and Kroszynski, 1971) However, a much simpler solution is obtained if we assume some analytical representation for the K(O) and 0 (~) relationships. We have eventually adopted the following: K _ [ 0 --Or ~3 |
Ks
\ 0s -- 0r
0 --0 r
- ea(•
(24)
! --
~c)
(25)
0s--0 r Eq.24 has been derived by Irmay (1954) and proved satisfactory for a large collection of soils. Eq.25 does fit measured retention curves quite well, at least in limited intervals of ~ variation, for a given value of the parameter a. Soil 2 definition curves were ehosen following those expressions and one can see in Fig.3 the agreement with soil 1. Substituting in eq.22 the following closed solutions are obtained after a few manipulations. The eigen-values: i 2 rr2 K s a °i = (0 s - 0L) e aL (e aL -- 1)
i = 0, 1, 2 , . . .
(26)
147 The eigen-functions: v0 = (0s-- 0L) - ' a
vi(z) =
eal~
- sin ilrv -- e az cos brv
0 s -- 0L)
(i = 1, 2, 3, ..)(27)
l~
where 0 L = O(z = L; t = 0) and v = (e az -- 1)/(e ai -- 1). The derivation of the Gi(t) functions is outlined in A p p e n d i x B. After s u m m a t i o n 21 is performed remembering from eq.25 that 0 s - 0 L = (0 s - Or)(eai -- 1)/e aL and the pertinent back substitutions are done, one finally obtains:
{ eaL (Z,t) = - z +@c- - eL Tea L_ I ~ 12 ~
i=l
~sininv \ i3~~
+
(eaL-- l ): I 2v 3 6aL
-
3V2--6v+2 --
3V:
+
eaL-1
+
e az cosinvl exp ( aLi2n2r I- I e e L - 1 i~ n~ ! - ( ~ ~ l)l) :] d
(0 <~z <.L)
(28)
and the drawdown:
{ s(t)=eL
e aL --
rea-~ff~l +-
eaL--1
3aL
[ 1--6
~ ~
i= 1
1 ~
( exp
ka_Li2n2r ]~ I (eaL -- 1)2]-J J
(29)
where:
r = Kst/nL Inspection of eqs.28 and 29 shows that, as time elapses, the influence of the time exponential terms decreases and @ and s become almost translatory solutions with linear dependence on time (which incidentally coincides with Childs and Poulovassilis' solution (1962}). In the case where a semi-infinite column is considered (aL+~}, we obtain from eqs.28 and 29 the simple expressions:
$(z,t) = - - z + ~ c - - e L
~ - ~ - ~al-,x/~ exp I t eaz-l~v' (r
s ( t ) =eL
{2
r+x/_~/r
}
(eaZ-1)2~4aLTJ + e ~ a L !)erfc (eaZ
1'1 (31)
Expressions 30 and 31 are independent of L though it appears because of the definition of r. By approximate curve fitting on semi-log paper we have found for soil 1:
148
a = 0.09 cm -1 . This value was intentionally selected from the beginning for soil 2 in order to apply the analytical solution only once for both soils. The corresponding definition curves based on eqs.24 and 25 are given (dashed lines) in Fig.3 (soil 2). With this value of a, O profiles have been computed and represented in Figs.9, 10 and 11 together with numerically computed profiles for soil 2. Although agreement is not satisfactory, the analytical solution (eq.28, eq.30) reproduces the general trends of the numerical one. The discrepancy is due mainly to the linearity of the solution below z = 0 (eqs.28 or 30 hold for
L>~z>~O). Z(cm), I !l/ Ili"~ 201-11\~,~
]
..... .......
'% ~
£=0.05 soil 2 numerK:al solulion analytical solu[ion Eqs 25and 2~5
\
-10 -
~ ' - t =60~o
10
20
30
34.2
0(%)
Volumetric moisture content profiles corresponding to numerical and analytical solutions, at different values of time for e = 0 . 0 5 . Fig.9.
Z(cm) 20
\ i~ i~ ii\~
~, =0.1 ..... soil 2numerical s o l u h o n . . . . . . . analytical solution Eqs. 2 5 a n d 2 8 ..... " " Eqs. 25and 30
~ ~,'~\~\
i 10
i 20
i 30
i 0(%)
Volumetric moisture content profiles corresponding to numerical and analytical solutions, at different values of time for e = 0 . 1 . Fig.10.
149
Z (:m) i
£=0.2
23~
, i
....
soil
--- •
ana!ybcal
2 numerical so:,uton so~Jhon
Eqs 2 5 a n d 2 8 Eqs 2 5 a P d 3 0
\ ~\\
s
i
\
-.
[ 10
t
-4 20
30
6(I;?'
Fig.11. Volumetric moisture content profiles corresponding to numerical and analytical solutions, at different values of time for e = 0.2. A mu ch b ette r agreement is obtained for the m o t i o n of the S.S. (Fig.7). Eqs.29 and 31 yield almost coincident values for r ~< 10. Here, eq.31 despite its simplicity proves to be very accurate in the time interval of interest, in which the unsaturated zone plays an i m p o r t a n t role. This result shows again t h a t the m o t i o n of the S.S. is n o t very sensitive to changes in L provided it is long enough, or to changes in a {this has been actually verified by selecting different values of a and comparing dr aw down curves). Since both, measured data and selection of L or a are subject to some u n c e r t a i n t y this p r o p e r t y of the solution is particularly useful. It should be kept in mind, however, t hat the behavior o f the first order solution fails to consider non-linear "tail effects" which b e c o m e apparent after some time in the presumably exact numerical solution. Th er ef or e, the suggested analytical expression 31 is restricted to small e and m o d e r a t e values of time. But, the major effects of the unsaturated zone are present precisely during this time interval. The analytical m e t h o d used here has been also applied previously (Dagan and Kroszynski, 1971, 1973) to the cases of vertical drainage by a sudden head drop at the c o l u m n b o t t o m and b y a head variation according to a prescribed f u n ctio n of time and the results have been also satisfactory. CONCLUSIONS The problem of drainage of a vertical porous column with constant flux at its b o t t o m has been solved b o t h numerically and in an a p p r o x i m a t e manner, analytically. In t he first stage of the drainage process, for say z ~< 1, the un° saturated zone has a delaying ef f ect on drainage. As a result, the water table drops faster than predicted by the usual constant porosity approach. After-
150
wards, drainage takes place approximately at constant effective porosity. This conclusion has immediate bearing on hydrological applications. It suggests, for instance, that in pumping tests in unconfined aquifers, the values of n determined by Theis' m e t h o d or any other m e t h o d underlain by the assumption of constant n is likely to yield values of n which change with time and distance from the well. Only for large times or distances the values are presumably close to the drainable porosity, but the interpretation of data becomes then difficult. Although unsatisfactory when predicting $ (z, t) profiles, the analytical method suggested provides a simple and reliable means of predicting the drawdown of the saturation surface (and therefore of the phreatic surface). We are presently engaged in the extension of our study to two- and threedimensional flows and hope to report results in the near future. ACKNOWLEDGEMENT
The present work is part of a Thesis to be submitted by U. Kroszynski as partial fulfilment of the requirements for a D.Sc. degree at the Technion, Israel Institute of Technology. The research has been partially sponsored by the U.S.D.A. through Grant No. FG-IS-287 under PL 480. APPENDIX A
Some details on the numerical scheme The boundary condition at the moving interface from eq.5 becomes:
- - (s(t),t) = e - - 1
(A-l)
az
instead of e q . l d in Hornberger and Remson (1970). This changes their eq.6d to:
Ou ~ ( 1 , ~ ) = R (1' - - e)
(h-2)
The iterative scheme in our case is given by the following expressions:
R*(m)
Kse
]+ 1
=-O,s
2
A~O,s
(
1
Uc-- Un-- l , j + l )
R}m++)l ~:e)A~--- ~(m-)Rj +1 J+ 1
+';)
(A-3)
(A4)
where R * = dR/dT other symbols being equal to those of Hornberger and Remson. For the first time step we use R~ m + 1) = 1 + fl~R,*(m) instead of (A-4). In any case an estimation of Un - 1,i + 1 is needed in (A-3) and it is taken as
151
U n - 1,1" The values obtained that way from (A-3) (A-4) are inserted into the predictor formula then in the corrector. With the corrected results of u one could apply the corrector formula again. This was actually done but no change in the results was observed. With the values ui, ] + 1 obtained from the corrector we actualize the values of 0 i 1 + 1, K~, J + 1 and from (A-3) (A-4) we get the actual values R j + 1, RT+ 1 using now Un - 1, j + 1 in (A-3). The procedure is then carried on forward in time with a different time step if desired. For the borders i = 0 and i = n a finite difference approximation to the boundary conditions is used leaving the predictor and corrector formulae as simple tridiagonal systems solved by Thomas' algorithm. Although the Predictor-corrector scheme is unconditionally stable with a truncation error of 0(A~ 2 + AT3/2) (Remson et al., 1971, p.89) for the iterative procedure to converge one needs Ar/A~ < i O's i R~/3. Although considered a multi (time) step method, in this case information is needed at the immediate previous time step solely. A Fortran IV version of the program is available upon request from the author. APPENDIX B
The derivation o f the Gi(t) function (eq. 21) The eigen-functions vi(z ) (eq.27) have the very important property of being orthonormal with respect to the weight function a0 (z), i.e.: L fo C~o(Z) vi(z) vy(z) dz = 6i1=
0 ~ (iCj) t (i = j ) 1
(B-l)
This fact enables us to express any function of z as a series of vi(z). In particular we can find the coefficients hi of the expansion:
ao(Z) dz
0(z) (1 --
=~o hj vj(z)
(B-2)
merely by multiplying both sides by a0 (z) x vi(z) and integrating over z:
hi:Lfo vi(Z) d z
Ko(z)
1--L
dz
(B-3)
Substitution of eq.21 into 18 and use of eqs.22 and B-2 yield: Oo
i= 0
oo,z, i,z,
+ oiGi(t) -- hi dt
)
- 0
(B-4)
This imposes for Gi(t) the following differential equation: dGi - - + oiGi(t) -- hi= 0 dt
(B-5)
152
T h e initial c o n d i t i o n t o (B-5) results f r o m e q . 2 0 b y use o f eq.21 and the orthon o r m a l i t y p r o p e r t y B-1 :
Gi(O) =-2 ~ ao(Z) 1 -o
(B-6)
vi(z) dz
T h e solution is:
Gi(t) =--hi + [ Gi(O)-- hi 1 e - a i t oi ~ii-
(i = 1, 2, 3 . . . )
(B-7)
while for i = 0 Go(t) = Go(0) + h o t
(B-S)
where ho = -- uoKs (see eq.B-3). APPENDIX C - NOTATION List of symbols a
Gi H K Ks L rl
r / e tf
q0
s(t) sl(t) SSAT
t V
Y(z, t) Z
c
0 0L 0r 0s
vi(z) ai T
~c
o(z) ~g(z,t)
empirical parameter (cm -1 ) time dependent coefficients height of the unsaturated portion of the column (cm) hydraulic conductivity (unsaturated) (cm/min) hydraulic conductivity at saturation (cm/min) characteristic length (cm) drainable porosity of porous medium effective porosity imposed constant flux drained from the column (cm/min) drawdown of the saturation surface (cm) linear approximation to s(t) (cm) drawdown for constant porosity conditions (cm) time (rain) non-dimensional variable involving space auxiliary function space variable (cm) small parameter (qJKs) volumetric moisture content volumetric moisture content at selected upper boundary residual volumetric moisture content volumetric moisture content at saturation eigenfunctions non-dimensional space variable eigenvalues (min 1) non-dimensional time head (negative) (cm of water) critical head (negative), air entry value head distribution for the initial state of rest first order term in the perturbation expansion of
Symbols used in Appendix A correspond to those used by Hornberger and Remson (1970).
153 REFERENCES Childs, E.C. and Poulovassilis, A., 1962. The moisture profile above a moving water table, J. Soil Sci., 13(2): 271--285. Dagan, G., 1964. Second order linearized theory of free surface flow in porous media. Houille Blanche, 8: 901--910. Dagan, G., 1966. The solution of the linearized equations of free surface flow in porous media. J. Mech., 5(2): 207--215. Dagan, G., 1968. A derivation of Dupuit theory of steady free-surface flow toward wells by matched asymptotic expansions. Water Resour. Res., 4(2): 403--412. Dagan, G. and Kroszynski, U., 1971. Drainage of a vertical column, USDA project No. A10-SWC-77. In: Development of Methods, Tools and Solutions for Unsaturated Flow. First Ann. Rept. (Part 2), Hyd. Lab. Technion, 71 pp. Dagan, G. and Kroszynski, U., 1973. Drainage of a vertical column. In: Physical Aspects of Soil Water and Salts in Ecosystems, Springer, Berlin, pp.17--28. Day, R.D. and Luthin, J.N., 1956. A numerical solution of the differential equation of flow for a vertical drainage problem. Soil Sci. Soc. Am., Proc., 20: 443--447. Douglas, J. and Jones, B.F., 1963. On predictor-corrector methods for non-linear parabolic differential equations, J. Soc. Ind. Appl. Math., 11: 195--204. Freeze, R.A., 1969. The mechanism of natural groundwater recharge and discharge 1 : Onedimensional, vertical, unsteady, unsaturated flow above a recharging or discharging ground-water flow system. Water Resour. Res., 5 : 153--171. Fujioka, Y. and Kitamura, T., 1964. Approximate solution of a vertical drainage problem. J. Geophys. Res., 69(24): 5249--5255. Gardner, W.R., 1962. Approximate solution of a nonsteady state drainage problem. Proc. Soil Sci. Soc. Am., 26: 129--132. Hornberger, G.M. and Remson, I., 1970. A moving boundary model of a one-dimensional saturated--unsaturated, transient porous flow system. Water Resour. Res., 6: 898--905. Irmay, S., 1954. On the hydraulic conductivity of unsaturated soils. Trans. Am. Geophys. Union, 35(3): 463--467. Jackson, R.D. and Whisler, F.D., 1970. Equations for approximating vertical nonsteadystate drainage of soil columns. Soil Sci. Soc. Am., Proc., 34: 715--718. Jensen, M.E. and Hanks, R.J., 1967. Nonsteady-state drainage from porous media. J. Irr. Drain. Div., ASCE, 93, No. IR3, pp.209--231. Kastanek, F., 1971. Numerical simulation technique for vertical drainage from a soil column. J. Hydrol., 14: 213--232. Philip, J.R., 1969. Theory of infiltration. In: Ven Te Chow (Editor), Advances in Hydroscience, Vol. 5. Academic Press, London, pp.215--296. Remson, I., Hornberger, G.M. and Molz, F.J., 1971. Numerical Methods in Subsurface Hydrology. Wiley-Interscience, New York, N.Y., p.89. Swartzendruber, D., 1969. The flow of water in unsaturated soils. In: R.J.M. de Wiest (Editor), Flow Through Porous Media. Academic Press, London, pp.215--292. Vachaud, G., 1968. Contribution ~ l'~tude des probl~mes d'~coulement en milieux poreux n o n satur~es. Th4se pour obtenir le grade de Docteur des Sciences physiques. Universit~ de Grenoble. Watson, K.K., 1967. Experimental and numerical study of column drainage. J. Hydrol. Div., ASCE, 93: 1--15. Whisler, F.D. and Watson, K.K., 1968. One-dimensional gravity drainage of uniform columns of porous materials. J. Hydrol., 6: 277--296. Youngs, E.G., 1960. The drainage of liquids from porous materials. J. Geophys. Res., 65(12): 4025--4030.
FLOW IN A VERTICAL AT CONSTANT FLUX
POROUS
COLUMN DRAINED
AT ITS BOTTOM
U. KROSZYNSKI
Department of Civil Engineering, Technion, Israel Institute of Technology, Haifa (Israel) (Accepted for publication March 24, 1974)
ABSTRACT
Kroszynski, U., 1975. Flow in a vertical porous column drained at its b o t t o m at constant flux. J. Hydrol., 24: 135--153. A soil column, saturated at its lower portion and initially in condition of hydrostatical equilibrium is drained at a constant flux at its bottom. Due to the presence of the unsaturated zone above the water table, the motion of the saturation surface depends in a complex manner on time and on soil properties. After defining a characteristic length, the equation of unsaturated flow is solved numerically by a finite difference scheme for different values of drained flux, for two similar soils. Although moisture and head profiles are different for both soils, drawdown curves have been found to be in good agreement during a certain time interval which depends on the value of the drained flux. The same problem is solved analytically by a linearization procedure. A simple expression which predicts saturation surface drawdowns quite accurately for small values of the drained flux and for moderate time intervals, is derived. An effective porosity coefficient is defined and its variation with time is represented. The bearing of the results on flow in unconfined aquifers is discussed.
INTRODUCTION In unconfined hydrological systems operating under drainage conditions (e.g., a q u i f e r p u m p i n g , d r a i n a g e b y t i l e d r a i n s ) w a t e r is r e l e a s e d m a i n l y b e cause of the lowering of the water table and at a much lesser extent from elastic storage. I t is c u s t o m a r y in t h e t h e o r y o f g r o u n d w a t e r f l o w t o a s s u m e t h a t t h e v o l u m e o f w a t e r d r a i n e d ( p e r u n i t a r e a ) is e q u a l t o t h e d r o p o f t h e w a t e r t a b l e t i m e s t h e e f f e c t i v e p o r o s i t y c o e f f i c i e n t n, w h i c h is g e n e r a l l y t a k e n as c o n s t a n t . A s it h a s l o n g b e e n r e a l i z e d , d u e m a i n l y t o t h e e x i s t e n c e o f a n u n s a t u r a t e d z o n e a b o v e t h e w a t e r t a b l e , t h e r e is a " d e l a y e d y i e l d " a n d t h e c o n s t a n t d r a i n a b l e p o r o s i t y is p r e s u m a b l y a t t a i n e d o n l y a s y m p t o t i c a l l y , f o r l a r g e v a l u e s of time. For this reason, some inconsistencies arise when interpreting field m e a s u r e m e n t s o r p r e d i c t i n g f o r m a t i o n p a r a m e t e r s o n t h e b a s i s o f c o n s t a n t n. I n t h i s c o n t e x t , in t h e f r a m e w o r k o f a s t u d y o n t h e i n f l u e n c e o f t h e u n -
136 saturated zone u p o n the aquifer drainage process, we have considered as a first step the one dimensional vertical flow problem. The drainage of vertical porous columns has been already considered by different authors, especially for the case where drainage is caused by a sudden head d r o p at the b o t t o m of an initially saturated column (constant head t ype b o u n d a r y condition). Youngs (1960) c o m p u t e d the variation of the drained volume by assuming a single diameter capillary tube model of the porous medium, which is a constant porosity approach. Gardner (1962), Fujioka and Kitamura (1964) and Sw a r t zendr ube r (1969) have assumed that the interface between the saturated and unsaturated zones, or saturation surface (S.S.) drops instantaneously from its initial position at the column t op to its ultimate one near the column b o t t o m and solved the equation of unsaturated flow by linearizing it. Vachaud (1968) and Jackson and Whisler (1970) c o m p u t e d the variation of the drained volume by generalizing Youngs' model namely by considering a bundle of capillary tubes of different diameters. Numerical solutions of the drainage problem under constant head b o u n d a r y condition at the column b o t t o m were given by Day and Luthin (1956), Jensen and Hanks (1967), Watson (1967), Whisler and Watson (1968) and Hornberger and Remson (1970). Som ew ha t different b o u n d a r y conditions were also considered by Freeze (1969). All those works were based upon finite difference schemes. A numerical simulation of the physical process of drainage at constant head at the column b o t t o m was presented by Kastanek (1971). The present study differs f r om those m ent i oned before in two aspects. (1) We consider the case of drainage under constant flux at the column b o t t o m rather than constant head. Such a study is relevant to a few applications where the aquifer is drained at constant discharge as a b o u n d a r y condition. Moreover, in the case of a constant head drop, the drainage process depends on the initial length of the saturated zone while this length is immaterial in the case of constant flux. This way, general conclusions may be reached at the expense of fewer computations. (2) The analytical m e t h o d we suggest considers variable coefficients in the equation of unsaturated flow rather than constant ones, which may yield a more accurate solution. Childs and Poulovassilis (1962) solved a problem closer to this context. T h e y considered the case where bot h water table and moisture c o n t e n t profile translate at a constant dow nw a r d (upward) speed. Although exact, their solution might perhaps apply asymptotically after a large time interval. Even in this case t h e y clearly show that the moisture profile c a n n o t descend at constant speed in the zone of low moisture content. STATEMENT OF THE PROBLEM AND QUALITATIVE ANALYSIS OF THE DRAINAGE PROCESS We consider a vertical c ol um n of h o m o g e n e o u s and stable soil. The column is initially saturated with water up to a certain level z = 0 (Fig.l) with an un-
137
/~ q=O
kl
-,6--
~Z
\)\
t--ol
',\',t~o 0 \
"~o t~,e
S(t~4 -
e~s t=O t-O,
(a) (b) Fig.1. Drainage of a vertical c o l u m n : (a) the m o t i o n of the saturation surface, and (b) distribution of 0 and ~ initially and at time t (qualitative).
saturated portion above it in hydrostatical equilibrium. The initial moisture profile is assumed to correspond to the 0 (~) main drainage curve. Starting at t = 0 water is withdrawn from the column b o t t o m at a constant rate q0 and the flow is permanently downwards. There is no flux through the column top. The initial and boundary conditions are such that the soil undergoes drainage at every point and consequently the 0 ($) relationship remains univalued. By continuity, q0 will also be the flux at the saturation surface S.S. This fact enables us to deal with the unsaturated zone solely, provided a proper (constant flux) boundary condition is stated at the S.S., thus making the saturated length of the column immaterial. The column is, moreover, assumed to be sufficiently long so that in the time interval to be considered, its b o t t o m remains at saturation while its top is always at low moisture content. The aim of the present work is to determine the motion of the saturation surface S.S. (z = --s(t)), which is of paramount importance in hydrologic applications. Under the considered b o u n d a r y conditions this is almost identical to the motion of the phreatic surface. If the flow in the unsaturated zone is neglected, as usually done in existing approaches or alternatively, if it is assumed t h a t the moisture profile translates without change of shape after the S.S., the solution to the problem obviously becomes: SSAW(t) =
qo t / n
(1)
where n is the so called drainable porosity, n can therefore be defined from eq.1 as the ratio between the volume of water released from the soil and the volume of soil spanned by the S.S. (See also Appendix C.)
138
Actually, there is a delay in the release of water from the unsaturated zone and the S.S. moves downwards faster than SsnT(t ) predicted by eq.1. As a result, the effective porosity: t/s(t)
neff = qo
(2)
is smaller than the drainable porosity of eq.1. Moreover, neff depends in a complex manner on the unsaturated properties of the soil, the discharge q0 and time. Its computation is one of the aims of the present study. We define now the problem in precise mathematical terms. The equation of unsaturated flow is written in terms of the head $. Above the S.S. we have:
~ ( ~ ) Ot
~z
(~)
+~(~)~z
where a(~) = d0/d~ and ~($) = dK/d~. The boundary and initial conditions are: (i) Zero flux through the column top: .....
~z
1
(4)
(z=H,t>~O)
(ii) Continuity of head and flux at the S.S.: q o = K s -(~ z~ )+ 1
; ~=~e
(z=--s(t),t>O)
(5)
(iii) The initial hydrostatical equilibrium distribution: =--z
(H > z /> 0, t = 0)
+ ~c
(6)
Eq.3 with conditions 4, 5, 6 allow only in principle to determine ~ (z, t) within the unsaturated portion of the column and the S.S. location s(t). Two major difficulties make a closed analytical solution almost impossible: the nonlinearity of eq.3 and the fact that the lower boundary z --- -- s(t) is an u n k n o w n of the problem. On the other hand the solution within the saturated portion of the column is readily given by: =
--
--
x (z + s(t))
(z < - - s ( t ) )
(7)
It is worthwhile to clarify two points at this stage. First, when dealing with a saturation surface which is sufficiently deep below the top, it is mathematically convenient and physically sound to transfer the upper boundary from the column top to z = H ' (Fig.2) where H ' is selected such that the hydraulic conductivity is negligibly small for z/> H ' (we have selected H ' in the examples presented in the sequel such that K(O L )/Ks = 1 0 - 4 (see Fig.2)). It seems reasonable to assume that practically no flow will take
139
Zj
H.z~ COLUMNTOP
Fig.2. Location of the upper boundary and the characteristic length for the unsaturated zone.
place in th at region in the time interval of interest. Moreover, this transfer makes the initial moisture c o n t e n t at the upper b o u n d a r y 0 = 8 L (Fig.2) a well defined q u a n t i t y rather than the irreducible moisture c o n t e n t 0 r which is some° tim~s difficult to detect. In any case we can define a characteristic length of the problem L as the smaller of the t w o quantities H and H ' and take b o u n d a r y conditio.l 4 at z = L. Secondly, on purely dimensional grounds, we may write the solution of eq.3 for a given soil as: ~ -function L
z
q0
Kst
' L
1>~-~>-L L
'
(8)
reducing thus the n u m b e r of variables. Since for q0 = 0 there is no flow and t h er ef o r e $ = $0 = -- z + ~c and s ( t ) = 0, it is legitimate to assume that for small values o f e = q o / K s , eq.8 m ay be written as: ~/L = ~o/L
+ e$ I/L + e2~2/L
+ . . .
(9)
Similarly for $ = $c in eq.8: s(t)/L
= esl(t)/L
+ e2s2(t)/L
+...
(10)
where ~ o / L + e ~ I / L and e s l ( t ) / L represent the "first order solution", obtained by neglecting quadratic and higher order terms in eqs.9 and 10. We now proceed to solving the pr obl e m numerically for two soils and different values of the parameter e. A simple a p p r o x i m a t e analytical solution based on the a f o r e m e n t i o n e d points is derived subsequently. Soil 1 is a river sand whose r e t e n t i o n and hydraulic c onduc ti vi t y curves (definition curves) have been taken f r om the publication by Vachaud (1968). Soil 2 is a hypothetical sand with the same 0s, 0r, $c and K s of soil 1 but whose definition curves follow exact analytical expressions which serve as basis for the
140
analytical solution (see eqs.24, 25). Such a choice has a twofold purpose: to study the influence of the difference of shape of the definition curves between both soils other things being equal and to appraise the validity of the "first order solution" {derived analytically) by comparing it with the numerical solution for soil 2. For the sake of completeness the definition curves of both soils are presented in Fig.3. This study will be restricted to comparison between numerical and analytical solutions solely since no experimental work dealing with constant flux boundary conditions has been done so far. -~cmofwater
60
g ecru/rain)
I --
~s: 0.342
--"t -
81-: 0.075 Ks= 1.09cm/min U)c:-16.5crn.
/4
soil 2
/ /
1.0
/
I
a= oDgcrn-'
4O
08 \
!
,,
/
3O
0£
\\
12 04,
20
10
S
A
i I
02
0 0
10
20
30
eC%)
Fig.3. Retention and conductivity curves for soil 1 (after Vachaud, 1968) and soil 2.
NUMERICAL EXPERIMENTS The problem formulated in eqs.3, 4, 5, 6 has been solved numerically by a finite difference scheme for both soils. The scheme used is that of Hornberger and Remson (1970) adapted to our boundary conditions. The scheme proved to be accurate when dealing with the constant head boundary condition drainage problem. The flow domain -- s ( t ) <~ z <. L is first mapped by the transformation ~ = ( L - - z ) / ( L + s ( t ) ) onto the interval 0 ~< ~ ~< 1. The differential equation and
141 b o u n d a r y c o n d i t i o n s are r e w r i t t e n in t e r m s o f t h e n e w variables ~ a n d K s t / L and t h e algebraic s y s t e m resulting f r o m t h e finite d i f f e r e n c e a p p r o x i m a t i o n is solved b y t h e D o u g l a s - J o n e s ( 1 9 6 3 ) p r e d i c t o r c o r r e c t o r t e c h n i q u e . T h e locat i o n a n d v e l o c i t y o f t h e S.S. are d e t e r m i n e d iteratively, t h e i t e r a t i o n p r o c e d u r e i m p o s i n g t h e r e l a t i o n s h i p b e t w e e n t h e space a n d t i m e i n c r e m e n t s . T h e acc u r a c y o f t h e s o l u t i o n was c h e c k e d a m o n g o t h e r criteria b y t h e m a t e r i a l balance, i.e. b y c o m p a r i n g t h e d r a i n e d v o l u m e q0 A t w i t h t h e d i f f e r e n c e bet w e e n t h e m o i s t u r e profiles i n t e g r a t e d o v e r z at t and t + A t. T o ensure acc u r a c y it was f o u n d t h a t q u i t e a d e n s e space grid is n e e d e d . T h e r e a s o n f o r t h a t is t h e s t e e p v a r i a t i o n o f t h e $ (and 0) profiles n e a r t h e S.S. As t i m e elapses t h e profiles b e c o m e s t e e p e r in t h e n o r m a l i z e d variable ~, n e a r ~ = 1. On the o t h e r h a n d a sparser m e s h c o u l d be used in t h e u p p e r p o r t i o n o f the uns a t u r a t e d zone, w h e r e a " t a i l " o f a l m o s t c o n s t a n t l o w m o i s t u r e c o n t e n t is d e v e l o p e d . H o w e v e r , in o r d e r to a v o i d c o n s i d e r a b l e c o m p l i c a t i o n s , a fixed 60 n o d e m e s h was used f r o m t h e beginning. T h e t i m e interval was selected as 60 m i n , a n d t h e values o f e = q o / K s used were 0.05, 0.1, 0.2. A c c o r d i n g t o t h e c r i t e r i o n K ( O L ) / K s = 10 --4 (see p r e c e d ing section), L was t a k e n 25.6 cm. Curtailing t h e e x t e n t o f t h e u n s a t u r a t e d z o n e at z = L saved c o m p u t e r t i m e w i t h o u t a f f e c t i n g a c c u r a c y . A n u m e r i c a l e x p e r i m e n t w i t h L = 40 c m and a m e s h o f 80 n o d e s was p e r f o r m e d f o r soil 1 yielding t h e s a m e results as t h o s e s h o w n in full lines in Fig.4 w h i c h c o r r e s p o n d t o L = 25.6 c m a n d 60 nodes. In t h e case s h o w n in Fig.4 f o r soil 1 t h e I.B.M. 3 6 0 / 7 0 T e c h n i o n c o m p u t e r p e r f o r m e d t h e n e a r l y 1 4 0 0 t i m e steps n e e d e d in a b o u t 102 seconds, c o m p i l a t i o n included. Z(cm) I i~, ,, 'i',
'
,o r I
4::0.0 5
xL'"IL,'
IZ i
20
\ \x
\\\\ 10 ~
\
l
\\ \\\\\x \'\X\ ,x\\
t=o iI
\\
i 0'-
t=O t:~O mI'
t:3s.,.
_1o
[,
10
t=60~,~ 20 "8-(%)
(cm)
30
-40
\\
h" XNX} ~
""~,
0 ',
~-io \,t =E,O i \'~I
-30 -2~0 -16 5 ~J(cm orwater)
Fig.4. Moisture content and head distributions at different values of time for e = 0.05 (numerical solution).
142
For details on the numerical scheme the reader is referred to Hornberger and Remson (1970) and to Appendix A. The results, namely the $ (z, t) and 0 (z, t) profiles for both soils are pre* sented in Figs.4, 5 and 6. Although the trends are similar the profiles are quite different for both soils. However, it should be noted that, as time increases, the profiles have a tendency to translate with the S.S. except for the "tail" zone of low moistures. Moreover, the motion of the S.S. appears to be little influenced by the shape of the definition curves at least for the times represented.
20
-E=0.1
~
!20
o
t=
\ 10
.',XC °
20 @(%)
30
t = 29.~ -40
~ 30 -20 -16 5 $ (cm ofwater)
Moisture c o n t e n t and head distributions at different values o f time for • = 0 . 1 (numerical s o l u t i o n ) . Fig.5.
Z(cm) \ \
Z(cm)
20 ~
{=(]2 ~
,~
soil 1
20 \\
\\\ \~,~':
10
10 t=0
°
\
-10
--10
"X
~ ' ~ t=5 mitt
~ ' ~ ' 20 30 e (%)
t=11 rain. -40
",Ix:
\\ \ \~-10
-,.t0 -20 -16.5 ~ (¢m. of water )
Fig.6. Moisture c o n t e n t and head distributions at different values o f t i m e for • = 0.2 (numerical solution).
143 T h e m o s t i n t e r e s t i n g p a r a m e t e r , n a m e l y t h e S.S. d r a w d o w n s(t) is represented in t e r m s o f d i m e n s i o n l e s s variables in Fig.7. s(t)/L has b e e n divided b y e in o r d e r t o c h e c k t h e v a l i d i t y o f t h e first o r d e r a p p r o x i m a t i o n w h i c h implies t h a t sl (t)/L b e c o m e s i n d e p e n d e n t o f e (see e q . 1 0 ) . T h e d i m e n s i o n l e s s t i m e ~ = Kst/nL (n = 0 s - - 6r) selected as t h e i n d e p e n d e n t variable m a k e s possible t o c o m p a r e b o t h soils on a c o m m o n basis. M o r e o v e r , t h e c o n s t a n t p o r o s i t y solut i o n w h i c h is t h e usual o n e a d o p t e d in h y d r o l o g y yields t h e straight line SsAw(t)/cL = r (see eq.1).
J
:
A
0
10-J
'
2
.......
4
---s0il 2 . . . . (E q 29 or 31 )
~tur solut,on(Eq1/
6
8
%=Kst/nL
E-~:],2/ / ¢~[0,I / > ~ J ~ ~,/~"~ J " ~ ' " .~ ~ ~ ~.~.-r~
f-.......
0k,~rd/P'4"..... 0
2
4
6
8
T'=K s t / n L
B
Fig.7. The motion o£ the saturation surface z = -- s(t).
T h e closeness b e t w e e n t h e d r a w d o w n curves f o r e = 0.05 a n d e = 0.1 suggests t h a t t h e first o r d e r s o l u t i o n ( n a m e l y t h e first t e r m in eq.10) w o u l d p r o v i d e a fair a p p r o x i m a t i o n f o r soil 1, f o r e < 0.1 until r ~-- 5 while f o r soil 2 o n l y until r ~ 2. N o n - l i n e a r t e r m s b e c o m e r e l e v a n t o u t s i d e t h a t range, as s h o w n b y the d i f f e r e n c e w i t h t h e d r a w d o w n curve c o r r e s p o n d i n g t o e = 0.2 o v e r alm o s t t h e entire t i m e range. O n t h e o t h e r h a n d b o t h soils have f o r e = 0 . 0 5 v e r y similar d r a w d o w n curves o v e r 0 ~ r ~ 10 t h u s suggesting t h a t in t h a t case or f o r smaller values o f e, t h e s h a p e o f t h e d e f i n i t i o n curves d o e s n o t a f f e c t d r a w d o w n s . F u r t h e r i n s p e c t i o n o f Fig.7 s h o w s t h a t at t h e b e g i n n i n g o f t h e d r a i n a g e pro-
144
cess, say for 0 < r < 1, the S.S. descends m u ch faster, and at a variable speed, as co mp ar ed with the constant por os i t y drawdown. This is precisely the interval in which the unsaturated zone delays the release of water and as a result s ( t ) becomes larger than sSAT. Later, the speed of descent becomes nearly constant corresponding to the phase in which most of the unsaturated profile almost translates behind the S.S. These two regimes are illustrated even be t t er by Fig.8 where the variation of the effective porosity coefficient neff (eq.2) with time is represented. The fact that n eff does n o t reach the ultimate value n = 0 s - 0r during the time interval considered, is explained by the existence of a certain volume of water retained in the " t a i l " o f the moisture profiles. Finally, it is seen t hat interpretation of actual draw dow n curves in order to derive a constant hypot het i cal porosity n = qo t / s m ay lead to quite small values of n, especially when drawdowns are measured shortly after the beginning o f drainage. APPROXIMATE ANALYTICAL SOLUTION The numerical solution can be used in given particular cases and it seems that its extension t o two- and three-dimensional cases is difficult and expensive. F o r these reasons it was considered useful to derive an a p p r o x i m a t e analytical solution of the drainage problem. Such a solution is also motivated n : (126 7 c
A
0
2
[
-
[
:
4
6
~
~
.....
8
T=P,: b ' n ~
n=0.267'
(Eqs 2 a n d 2 9 o r 3 1 )
of B
Fig.& Variation of the effective porosity coefficient defined in eq.2.
145 by the fact that for e = q o / K s sufficiently small, as encountered in application, the "first order solution" in eq.10, represents quite accurately the numerical solution (Figs.7, 8). The proposed m e t h o d takes advantage of this property. We derive the linearized approximation by substituting expansions 9 and 10 in the exact equations of flow 3 to 6, expanding the a, K and ~ coefficients in Taylor series near C0 and neglect terms of order e2 or higher. Collecting terms independent of e we get the zero order approximation to C representing the initial state of rest. Co = - - z + Cc
(11)
Collecting now the terms of order e we obtain the following set for C l (z, t) and sl (t) (for details see Dagan and Kroszynski (1971)). ao (z) OC~ - a [K o (z)
at1 ]
(L >~z ~> 0)
(12)
- 0
(z = L)
(13)
- 1
(z = 0)
(14)
o
(t = o)
(151)
at
az
~zJ
aC~ az
at1 az ~
--
s~ ( t ) = - - ~ ~ ( O , t )
(161)
where in eq. 12 a0 (z) = a (C 0 ) and K0 (z) = K(C0 ) are known functions of z based on the initial state. In contrast with usual linearization procedures where one or both these coefficients are taken as constant, in the present case t h e y are variable. The m e t h o d used is known as m e t h o d of small perturbation expansions. It is widely used in different branches of physics and has been applied before to saturated groundwater flow problems (Dagan, 1964, 1966, 1968). Eqs.12 to 15, satisfied by C~, are much simpler than the exact original ones with two respects: (a) rather than the non-linear eq.3,we have now the linear eq.12; and (b) the boundary conditions (5) applying to the moving boundary z = - - s ( t ) a r e transferred to the fixed z = 0. Still c o m p u t a t i o n of the drawdown is possible by means of eq.16. Considering the function : 1--
g(z,t)=C~(z,t)+L
(17)
2
we obtain for Y by substitution in eqs.12--15 the following set: ao(Z) at - az
Ko (z)
+ 1 --
(L/> z >/0)
(18)
146
3Y 3z
-0
Y=-2
(z = 0; z = L)
(19)
(t = O)
(20)
1--
Separation of variables in the form: g =~
(21)
ui(z) G i ( t )
1
yields ui(z) as e igen-functions of the Sturm-Liouville equation: dz
Ko
dz
+ OiaoU i
0
(L ~> z ~> 0)
(22)
(z = 0;z
(23)
with b o u n d a r y conditions: dui dz
- 0
= L)
Eq. 22 can be solved approximatively to a desired accuracy by Galerkin's m e t h o d (this has been done in our previous work, Dagan and Kroszynski, 1971) However, a much simpler solution is obtained if we assume some analytical representation for the K(O) and 0 (~) relationships. We have eventually adopted the following: K _ [ 0 --Or ~3 |
Ks
\ 0s -- 0r
0 --0 r
- ea(•
(24)
! --
~c)
(25)
0s--0 r Eq.24 has been derived by Irmay (1954) and proved satisfactory for a large collection of soils. Eq.25 does fit measured retention curves quite well, at least in limited intervals of ~ variation, for a given value of the parameter a. Soil 2 definition curves were ehosen following those expressions and one can see in Fig.3 the agreement with soil 1. Substituting in eq.22 the following closed solutions are obtained after a few manipulations. The eigen-values: i 2 rr2 K s a °i = (0 s - 0L) e aL (e aL -- 1)
i = 0, 1, 2 , . . .
(26)
147 The eigen-functions: v0 = (0s-- 0L) - ' a
vi(z) =
eal~
- sin ilrv -- e az cos brv
0 s -- 0L)
(i = 1, 2, 3, ..)(27)
l~
where 0 L = O(z = L; t = 0) and v = (e az -- 1)/(e ai -- 1). The derivation of the Gi(t) functions is outlined in A p p e n d i x B. After s u m m a t i o n 21 is performed remembering from eq.25 that 0 s - 0 L = (0 s - Or)(eai -- 1)/e aL and the pertinent back substitutions are done, one finally obtains:
{ eaL (Z,t) = - z +@c- - eL Tea L_ I ~ 12 ~
i=l
~sininv \ i3~~
+
(eaL-- l ): I 2v 3 6aL
-
3V2--6v+2 --
3V:
+
eaL-1
+
e az cosinvl exp ( aLi2n2r I- I e e L - 1 i~ n~ ! - ( ~ ~ l)l) :] d
(0 <~z <.L)
(28)
and the drawdown:
{ s(t)=eL
e aL --
rea-~ff~l +-
eaL--1
3aL
[ 1--6
~ ~
i= 1
1 ~
( exp
ka_Li2n2r ]~ I (eaL -- 1)2]-J J
(29)
where:
r = Kst/nL Inspection of eqs.28 and 29 shows that, as time elapses, the influence of the time exponential terms decreases and @ and s become almost translatory solutions with linear dependence on time (which incidentally coincides with Childs and Poulovassilis' solution (1962}). In the case where a semi-infinite column is considered (aL+~}, we obtain from eqs.28 and 29 the simple expressions:
$(z,t) = - - z + ~ c - - e L
~ - ~ - ~al-,x/~ exp I t eaz-l~v' (r
s ( t ) =eL
{2
r+x/_~/r
}
(eaZ-1)2~4aLTJ + e ~ a L !)erfc (eaZ
1'1 (31)
Expressions 30 and 31 are independent of L though it appears because of the definition of r. By approximate curve fitting on semi-log paper we have found for soil 1:
148
a = 0.09 cm -1 . This value was intentionally selected from the beginning for soil 2 in order to apply the analytical solution only once for both soils. The corresponding definition curves based on eqs.24 and 25 are given (dashed lines) in Fig.3 (soil 2). With this value of a, O profiles have been computed and represented in Figs.9, 10 and 11 together with numerically computed profiles for soil 2. Although agreement is not satisfactory, the analytical solution (eq.28, eq.30) reproduces the general trends of the numerical one. The discrepancy is due mainly to the linearity of the solution below z = 0 (eqs.28 or 30 hold for
L>~z>~O). Z(cm), I !l/ Ili"~ 201-11\~,~
]
..... .......
'% ~
£=0.05 soil 2 numerK:al solulion analytical solu[ion Eqs 25and 2~5
\
-10 -
~ ' - t =60~o
10
20
30
34.2
0(%)
Volumetric moisture content profiles corresponding to numerical and analytical solutions, at different values of time for e = 0 . 0 5 . Fig.9.
Z(cm) 20
\ i~ i~ ii\~
~, =0.1 ..... soil 2numerical s o l u h o n . . . . . . . analytical solution Eqs. 2 5 a n d 2 8 ..... " " Eqs. 25and 30
~ ~,'~\~\
i 10
i 20
i 30
i 0(%)
Volumetric moisture content profiles corresponding to numerical and analytical solutions, at different values of time for e = 0 . 1 . Fig.10.
149
Z (:m) i
£=0.2
23~
, i
....
soil
--- •
ana!ybcal
2 numerical so:,uton so~Jhon
Eqs 2 5 a n d 2 8 Eqs 2 5 a P d 3 0
\ ~\\
s
i
\
-.
[ 10
t
-4 20
30
6(I;?'
Fig.11. Volumetric moisture content profiles corresponding to numerical and analytical solutions, at different values of time for e = 0.2. A mu ch b ette r agreement is obtained for the m o t i o n of the S.S. (Fig.7). Eqs.29 and 31 yield almost coincident values for r ~< 10. Here, eq.31 despite its simplicity proves to be very accurate in the time interval of interest, in which the unsaturated zone plays an i m p o r t a n t role. This result shows again t h a t the m o t i o n of the S.S. is n o t very sensitive to changes in L provided it is long enough, or to changes in a {this has been actually verified by selecting different values of a and comparing dr aw down curves). Since both, measured data and selection of L or a are subject to some u n c e r t a i n t y this p r o p e r t y of the solution is particularly useful. It should be kept in mind, however, t hat the behavior o f the first order solution fails to consider non-linear "tail effects" which b e c o m e apparent after some time in the presumably exact numerical solution. Th er ef or e, the suggested analytical expression 31 is restricted to small e and m o d e r a t e values of time. But, the major effects of the unsaturated zone are present precisely during this time interval. The analytical m e t h o d used here has been also applied previously (Dagan and Kroszynski, 1971, 1973) to the cases of vertical drainage by a sudden head drop at the c o l u m n b o t t o m and b y a head variation according to a prescribed f u n ctio n of time and the results have been also satisfactory. CONCLUSIONS The problem of drainage of a vertical porous column with constant flux at its b o t t o m has been solved b o t h numerically and in an a p p r o x i m a t e manner, analytically. In t he first stage of the drainage process, for say z ~< 1, the un° saturated zone has a delaying ef f ect on drainage. As a result, the water table drops faster than predicted by the usual constant porosity approach. After-
150
wards, drainage takes place approximately at constant effective porosity. This conclusion has immediate bearing on hydrological applications. It suggests, for instance, that in pumping tests in unconfined aquifers, the values of n determined by Theis' m e t h o d or any other m e t h o d underlain by the assumption of constant n is likely to yield values of n which change with time and distance from the well. Only for large times or distances the values are presumably close to the drainable porosity, but the interpretation of data becomes then difficult. Although unsatisfactory when predicting $ (z, t) profiles, the analytical method suggested provides a simple and reliable means of predicting the drawdown of the saturation surface (and therefore of the phreatic surface). We are presently engaged in the extension of our study to two- and threedimensional flows and hope to report results in the near future. ACKNOWLEDGEMENT
The present work is part of a Thesis to be submitted by U. Kroszynski as partial fulfilment of the requirements for a D.Sc. degree at the Technion, Israel Institute of Technology. The research has been partially sponsored by the U.S.D.A. through Grant No. FG-IS-287 under PL 480. APPENDIX A
Some details on the numerical scheme The boundary condition at the moving interface from eq.5 becomes:
- - (s(t),t) = e - - 1
(A-l)
az
instead of e q . l d in Hornberger and Remson (1970). This changes their eq.6d to:
Ou ~ ( 1 , ~ ) = R (1' - - e)
(h-2)
The iterative scheme in our case is given by the following expressions:
R*(m)
Kse
]+ 1
=-O,s
2
A~O,s
(
1
Uc-- Un-- l , j + l )
R}m++)l ~:e)A~--- ~(m-)Rj +1 J+ 1
+';)
(A-3)
(A4)
where R * = dR/dT other symbols being equal to those of Hornberger and Remson. For the first time step we use R~ m + 1) = 1 + fl~R,*(m) instead of (A-4). In any case an estimation of Un - 1,i + 1 is needed in (A-3) and it is taken as
151
U n - 1,1" The values obtained that way from (A-3) (A-4) are inserted into the predictor formula then in the corrector. With the corrected results of u one could apply the corrector formula again. This was actually done but no change in the results was observed. With the values ui, ] + 1 obtained from the corrector we actualize the values of 0 i 1 + 1, K~, J + 1 and from (A-3) (A-4) we get the actual values R j + 1, RT+ 1 using now Un - 1, j + 1 in (A-3). The procedure is then carried on forward in time with a different time step if desired. For the borders i = 0 and i = n a finite difference approximation to the boundary conditions is used leaving the predictor and corrector formulae as simple tridiagonal systems solved by Thomas' algorithm. Although the Predictor-corrector scheme is unconditionally stable with a truncation error of 0(A~ 2 + AT3/2) (Remson et al., 1971, p.89) for the iterative procedure to converge one needs Ar/A~ < i O's i R~/3. Although considered a multi (time) step method, in this case information is needed at the immediate previous time step solely. A Fortran IV version of the program is available upon request from the author. APPENDIX B
The derivation o f the Gi(t) function (eq. 21) The eigen-functions vi(z ) (eq.27) have the very important property of being orthonormal with respect to the weight function a0 (z), i.e.: L fo C~o(Z) vi(z) vy(z) dz = 6i1=
0 ~ (iCj) t (i = j ) 1
(B-l)
This fact enables us to express any function of z as a series of vi(z). In particular we can find the coefficients hi of the expansion:
ao(Z) dz
0(z) (1 --
=~o hj vj(z)
(B-2)
merely by multiplying both sides by a0 (z) x vi(z) and integrating over z:
hi:Lfo vi(Z) d z
Ko(z)
1--L
dz
(B-3)
Substitution of eq.21 into 18 and use of eqs.22 and B-2 yield: Oo
i= 0
oo,z, i,z,
+ oiGi(t) -- hi dt
)
- 0
(B-4)
This imposes for Gi(t) the following differential equation: dGi - - + oiGi(t) -- hi= 0 dt
(B-5)
152
T h e initial c o n d i t i o n t o (B-5) results f r o m e q . 2 0 b y use o f eq.21 and the orthon o r m a l i t y p r o p e r t y B-1 :
Gi(O) =-2 ~ ao(Z) 1 -o
(B-6)
vi(z) dz
T h e solution is:
Gi(t) =--hi + [ Gi(O)-- hi 1 e - a i t oi ~ii-
(i = 1, 2, 3 . . . )
(B-7)
while for i = 0 Go(t) = Go(0) + h o t
(B-S)
where ho = -- uoKs (see eq.B-3). APPENDIX C - NOTATION List of symbols a
Gi H K Ks L rl
r / e tf
q0
s(t) sl(t) SSAT
t V
Y(z, t) Z
c
0 0L 0r 0s
vi(z) ai T
~c
o(z) ~g(z,t)
empirical parameter (cm -1 ) time dependent coefficients height of the unsaturated portion of the column (cm) hydraulic conductivity (unsaturated) (cm/min) hydraulic conductivity at saturation (cm/min) characteristic length (cm) drainable porosity of porous medium effective porosity imposed constant flux drained from the column (cm/min) drawdown of the saturation surface (cm) linear approximation to s(t) (cm) drawdown for constant porosity conditions (cm) time (rain) non-dimensional variable involving space auxiliary function space variable (cm) small parameter (qJKs) volumetric moisture content volumetric moisture content at selected upper boundary residual volumetric moisture content volumetric moisture content at saturation eigenfunctions non-dimensional space variable eigenvalues (min 1) non-dimensional time head (negative) (cm of water) critical head (negative), air entry value head distribution for the initial state of rest first order term in the perturbation expansion of
Symbols used in Appendix A correspond to those used by Hornberger and Remson (1970).
153 REFERENCES Childs, E.C. and Poulovassilis, A., 1962. The moisture profile above a moving water table, J. Soil Sci., 13(2): 271--285. Dagan, G., 1964. Second order linearized theory of free surface flow in porous media. Houille Blanche, 8: 901--910. Dagan, G., 1966. The solution of the linearized equations of free surface flow in porous media. J. Mech., 5(2): 207--215. Dagan, G., 1968. A derivation of Dupuit theory of steady free-surface flow toward wells by matched asymptotic expansions. Water Resour. Res., 4(2): 403--412. Dagan, G. and Kroszynski, U., 1971. Drainage of a vertical column, USDA project No. A10-SWC-77. In: Development of Methods, Tools and Solutions for Unsaturated Flow. First Ann. Rept. (Part 2), Hyd. Lab. Technion, 71 pp. Dagan, G. and Kroszynski, U., 1973. Drainage of a vertical column. In: Physical Aspects of Soil Water and Salts in Ecosystems, Springer, Berlin, pp.17--28. Day, R.D. and Luthin, J.N., 1956. A numerical solution of the differential equation of flow for a vertical drainage problem. Soil Sci. Soc. Am., Proc., 20: 443--447. Douglas, J. and Jones, B.F., 1963. On predictor-corrector methods for non-linear parabolic differential equations, J. Soc. Ind. Appl. Math., 11: 195--204. Freeze, R.A., 1969. The mechanism of natural groundwater recharge and discharge 1 : Onedimensional, vertical, unsteady, unsaturated flow above a recharging or discharging ground-water flow system. Water Resour. Res., 5 : 153--171. Fujioka, Y. and Kitamura, T., 1964. Approximate solution of a vertical drainage problem. J. Geophys. Res., 69(24): 5249--5255. Gardner, W.R., 1962. Approximate solution of a nonsteady state drainage problem. Proc. Soil Sci. Soc. Am., 26: 129--132. Hornberger, G.M. and Remson, I., 1970. A moving boundary model of a one-dimensional saturated--unsaturated, transient porous flow system. Water Resour. Res., 6: 898--905. Irmay, S., 1954. On the hydraulic conductivity of unsaturated soils. Trans. Am. Geophys. Union, 35(3): 463--467. Jackson, R.D. and Whisler, F.D., 1970. Equations for approximating vertical nonsteadystate drainage of soil columns. Soil Sci. Soc. Am., Proc., 34: 715--718. Jensen, M.E. and Hanks, R.J., 1967. Nonsteady-state drainage from porous media. J. Irr. Drain. Div., ASCE, 93, No. IR3, pp.209--231. Kastanek, F., 1971. Numerical simulation technique for vertical drainage from a soil column. J. Hydrol., 14: 213--232. Philip, J.R., 1969. Theory of infiltration. In: Ven Te Chow (Editor), Advances in Hydroscience, Vol. 5. Academic Press, London, pp.215--296. Remson, I., Hornberger, G.M. and Molz, F.J., 1971. Numerical Methods in Subsurface Hydrology. Wiley-Interscience, New York, N.Y., p.89. Swartzendruber, D., 1969. The flow of water in unsaturated soils. In: R.J.M. de Wiest (Editor), Flow Through Porous Media. Academic Press, London, pp.215--292. Vachaud, G., 1968. Contribution ~ l'~tude des probl~mes d'~coulement en milieux poreux n o n satur~es. Th4se pour obtenir le grade de Docteur des Sciences physiques. Universit~ de Grenoble. Watson, K.K., 1967. Experimental and numerical study of column drainage. J. Hydrol. Div., ASCE, 93: 1--15. Whisler, F.D. and Watson, K.K., 1968. One-dimensional gravity drainage of uniform columns of porous materials. J. Hydrol., 6: 277--296. Youngs, E.G., 1960. The drainage of liquids from porous materials. J. Geophys. Res., 65(12): 4025--4030.