- Email: [email protected]
Variably saturated modeling of transient drainage: sensitivity to soil properties
Journal of
Hydrology ELSEVIER
Journal of Hydrology 161 (1994) 91-108
[1]
Variably saturated modeling of transient drainage: sensitivity to soil properties W i l l i a m R. W i s e * , T . P . C l e m e n t l, F r e d J. M o l z Department of Civil Engineering, 238 Harbert Engineering Center, Auburn University, Auburn, AL 36849-5337, USA
Received 13 September 1993; revision accepted 16 March 1994
Abstract
Numerical analysis of a transient, two-dimensional (rectangular symmetry), unconfineddrainage problem, using a variably saturated flow model, shows that the location of the phreatic surface and the height of the seepage face are functions of the capillary forces exerted in the vadose zone. A sensitivity analysis investigating the effects of independent variations in the saturated hydraulic conductivity and the two Van Genuchten parameters is compared to one performed where the saturated hydraulic conductivity is functionally related to the Van Genuchten parameters. The relationship between these descriptive soil parameters, which may be used to describe the pore-size-density function of a porous medium, and the saturated hydraulic conductivity is surmised based upon limited empirical evidence. Expressing the permeability as a function of these soil properties reverses the sensitivities of the variably saturated model to the Van Genuchten parameters. From the preliminary work performed herein, consideration of the number of degrees-of-freedom used when performing a sensitivity analysis is shown to warrant great care, if the predicted sensitivities are to have meaningful physical interpretations.
I. Introduction
Numerical models are frequently employed to predict the transient behavior o f subsurface water. These models are developed based u p o n mathematical descriptions o f the physical processes considered. Consequently, the quality o f any particular numerical model directly depends u p o n the accuracy and completeness o f the physical understanding upon which the model is built. This linkage is of great significance when modeling * Corresponding author. l Present address: Battelle, Pacific Northwest Laboratories, Battelle Boulevard, P.O. Box 999, Richland, WA 99352, USA. 0022-1694/94/$07.00 © 1994 - Elsevier Science B.V. All rights reserved SSDI 0022-1694(94)02509-A
92
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
transient unconfined flow through the subsurface, where the mathematical descriptions are highly nonlinear. This communication is aimed at shedding insight into how a sensitivity analysis should be properly conducted for variably saturated models in order to maintain an alliance with the underlying physics of flow. While many studies of steady-state seepage faces have been performed (see Clement et al. (1994) for an extended discusssion), fewer transient problems have been investigated. Most of the transient unconfined drainage problems are restricted to onedimensional analyses (e.g. Young and Aggelides, 1976; Sisson et al., 1980; Watson and Awadalla, 1985). A few studies have investigated two-dimensional drainage processes using variably saturated flow models (e.g. Rubin, 1968; Verma and Brutsaert, 1970; Todsen, 1973; Vauclin et al., 1979; Ahmad et al., 1993). These studies reveal the importance of considering transient flow processes in the unsaturated zone and demonstrate the failure of fully saturated flow models (such as that posed by Liggett (1977), for example) in simulating these types of infiltration curedrainage problems. Common failings of fully saturated flow models for describing unconfined drainage problems include the observation that the fully saturated flow models predict slower responses in water tables than actually occur. Vauclin et al. (1975) present an experimental investigation of transient seepage-face and phreatic-surface positions in a two-dimensional (rectangular) flow problem. A 6 m × 2 m saturated porous medium was allowed to drain after a sudden drop in the external water table. This sudden drop problem was investigated by measuring the water content and pressure head distribution during the drainage phase. The variably saturated algorithm used in this work is verified successfully by simulating Vauclin et al.'s (1975) experimental results (Clement et al., 1994). Within this paper, certain terms are used as precisely as possible. The terms phreatic surface and water table are used interchangeably to describe the surface connecting the locus of points at atmospheric pressure, excluding those along the seepage face. The term vadose zone is used to describe every part of the soil system residing above the phreatic surface. Thus the capillary fringe, which is saturated, is included in the vadose zone. The term unsaturated zone refers to that part of the vadose zone residing above the capillary fringe. As the name implies, the saturation within the unsaturated zone is less than unity; within the vadose zone, by contrast, it is less than or equal to unity. The term permeability is used to refer to the ability of a porous medium to transmit a fluid which saturates the medium. By contrast, the term hydraulic conductivity is used to describe the transmission of water through the medium at a given water content under a unit hydraulic gradient. Precise use of these terms is helpful below, when unraveling some of the more subtle effects involved with modeling variably saturated flow. The objectives of the present work are to perform and contrast sensitivity analyses of the variably saturated model to the soil properties (the two Van Genuchten parameters and the saturated hydraulic conductivity) using two methods, one employing three degrees-of-freedom, the other, two. In order to accomplish the second part of this objective, a new empirical way to represent the saturated hydraulic conductivity as a function of the Van Genuchten parameters is introduced. To set the tone for the modeling exercises reported, the problem at hand is developed below.
93
W.R. Wise et al. / Journal o f Hydrology 161 (1994) 91 108 V B
C !
p(x, t)
V h2
Z
±
A
L Fig. 1. Schematic of the flow domain used in all simulations presented in this paper (L = 10m, hI = 10m and h 2 = 3 m for all simulations).
2. Problem description A schematic of the problem used to pursue the stated objective is presented in Fig. 1. The 10 m x 10 m square e m b a n k m e n t has no-flow boundaries on the base and on the top of the flow domain. The porous medium is assumed to be homogeneous and isotropic. Values of the soil parameters used in the test problem are: a saturated hydraulic conductivity of Ks = 5.1 m day -1, a porosity of 77 = 0.46, a residual water content of Or = 0.01, and a specific storage of Ss = 5 x 10 -5 m - I. In the base case, the values of the Van Genuchten soil parameters used are: C~v = 2.0 m -1 and nv = 2.8, which represent a poorly sorted, fine sand. The water table is initially at the top of the block and the flow system is at static equilibrium. The water table on the left boundary is maintained at a constant level of z = 10 m, while that on the right face of the block is lowered instantaneously (at t = 0) from z = 10 m to z = 3 m and is maintained there, subsequently. Transient locations of the phreatic surface and the seepage face are determined numerically. Due to the instantaneous drop in the external water level on the right boundary, a rapid drop in pressure occurs along the right edge of the block, resulting in large pressure gradients. Pressure changes inside the block outpace the drainage of the pores. Consequently, due to the vertical extent of the problem at hand, specific storage effects are important, especially at early times. (Unfortunately, no experimental data is found in the literature for analyzing the specific storage effects.)
3. Modeling two-dimensional transient flow with a seepage-face boundary The governing equation for water flow in a two-dimensional (rectangular
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
94
symmetry), variably saturated porous medium is written (Clement et al., 1994):
00~
--[ O0
Ss- --g? + ot
I
0 K(o) O~;
ox
0 KO)~
+--
Oz
(1)
where Ss is the specific storage of the medium (L-l), 0 is the water content, ~/is the porosity, K(O) is the hydraulic conductivity (LT-1), t is time (T), • = h - z, where is the pressure head (L), h is the hydraulic head (L), and x and z are the horizontal and vertical coordinates respectively (L). The soil properties 0(q), and K(O) are herein described using Van Genuchten's (1980) soil-water-retention function and Mualem's (1976) variably saturated hydraulic conductivity function, respectively. The relationship between water content and pressure head (under tension) is given by (Van Genuchten, 1980): 1
mv
where c~v, nv and mv = 1 - 1/nv are the Van Genuchten parameters, whose values depend upon the soil properties. The parameter c~v is, in part, a measure of the first moment of the pore-size-density function (L -1) (as ctv increases, so does the first moment); nv is an inverse measure of the second moment of the pore size-density function (as nv increases, the pore-size-density function becomes narrower) (Wise, 1991). O is the effective saturation given by the relation: 0 - 0r O- - 0s - Or
(3)
where 0s is the saturated water content (often approximated by the porosity 7/) and Or is the residual water content of the soil. Of course, for non-negative ffJ, 6) = 1, that is, the medium is saturated. Based on Mualem's (1976) model, the relation between water content and hydraulic conductivity is given by (Van Genuchten, 1980): K(O) = Ks{1 - [1 - O(1/mv)]mv}2(~ 1/2
(4)
where K s is the saturated hydraulic conductivity (LT l). Note that the relative hydraulic conductivity, K(O)/Ks, does not depend upon the value of c~v. The initial and boundary conditions for the variably saturated, transient drainage problem generically shown in Fig. 1 are as follows: h =-h i
Vx, Vz, t = 0 (initial condition; hydrostatic equilibrium),
(5)
h=hl
along AB, (constant-head boundary),
(6)
along BC, (no-flow boundary),
(7)
along CD, (no-flow boundary),
(8)
along DE, (seepage-face boundary),
(9)
Oh Oz Oh --~0 Ox
--~0
h=z
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
h = h2
Oh
0~ = 0
95
along EF, (constant-head boundary), and
(10)
along AF, (no-flow boundary).
(11)
Numerical solution to the variably saturated problem discussed earlier is accomplished by using the numerical algorithm presented by Clement et al. (1994). The nodal spacing in the x and z directions used is equal, Ax = Az = 0.2 m. A variable time step, At, ranging from 0.0001 day to 0.1 day, is used in all simulations. The details of the finite-difference solution techniques, implementation of the boundary conditions, and verification of the algorithm are discussed by Clement et al. (1994).
3.1. Results of the variably saturated model for the sample problem Shown in Fig. 2 are temporal profiles of the phreatic surface predicted by the variably saturated model for the base case. As seen from the figure, the decline of the phreatic surface is quite rapid during earlier times and slows at later times. Fig. 3 shows the distribution of the drainage rates due to the various processes in the early stages of the simulation. This figure illustrates that the early drainage (0 < t < 0.00012 day) is dominated by the specific storage effects which accommodate the swift pressure variations within the porous medium and the resulting rapid decline in the phreatic surface. However, as time progresses, the specific storage effects decrease rapidly and the pore drainage begins to have a more pronounced influence on the flow. The specific storage effects are negligible after about 0.0004 day (35 s). Because the specific storage effects are so short lived, there is no compelling reason to
•
"::
:22.2 .........
- ..................
o ¢I
...... t = O.OOff2d ..... t=O.O2d ----t=O.14d - - - t =0.60d SteadyStar, -
-
i
i
i
I
2
,
,
,
I
4
,
,
,
I
L
6
,
,
I
8
,
,
,
I0
X Coordinate (m) Fig. 2. T e m p o r a l p h r e a t i c - s u r f a c e profiles predicted by the v a r i a b l y s a t u r a t e d m o d e l (av = 2 . 0 m - I , nv = 2.8 a n d Ks = 5.1 m d a y - 1 ) .
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
96
250
"~
150
loo
50
o 0.0000
0.0001
0.0002
0.0003
0.0004
Elapsed Time (d) Fig. 3. F l o w histories d u e to the v a r i o u s c o n t r i b u t o r y p r o c e s s e s - - e a r l y t i m e (av = 2.0 m -1 , n v = 2.8 a n d K s = 5.1 m d a y - l ) .
explore the sensitivity of the model with respect to the specific storage, Ss. Fig. 4 illustrates the inflow, outflow and pore-drainage flow components at later times.
3.2. Observations of the drainage behavior The soil moisture distribution is always continuous, and, depending on the soil 50
,
,
,
,
<
40
"~
30
'-
20
Outflow Inflow
10
o
Pore Drainage i
0 0.0
i
i
I
0.2
i
i
i
I
0.4
l
i
i
I
i
0.6
i
l
[
0.8
,
l
l
1.0
Elapsed Time (d) Fig. 4. F l o w h i s t o r i e s d u e to the v a r i o u s c o n t r i b u t o r y p r o c e s s e s - - l a t e r t i m e (av = 2 . 0 m - l , nv = 2.8 a n d Ks = 5.1 m d a y - l ) .
W.R. Wise et al. / Journal o[Hydrology 161 (1994) 91-108
PhreaticSurface/4 " ~
97
100 %
¢
0
2
4
6
8
10
X Coordinate (m) Fig. 5. C o n t o u r d i a g r a m o f the m o i s t u r e c o n t e n t p r e d i c t e d by the v a r i a b l y s a t u r a t e d m o d e l at t = 0.14 d a y (C~v = 2 . 0 m - l , nv = 2.8 a n d Ks = 5.1 m d a y - l ) .
properties, the vadose zone above the phreatic surface can store a significant amount of water. Fig. 5 is a contour diagram of the effective saturation at time t = 0.14 day, predicted by the variably saturated model. The location of the phreatic surface predicted by the model at this time is also shown in the figure. This figure clearly illustrates that the porous medium retains significant moisture well above the phreatic surface. Similar observations can also be made from Fig. 6, which shows superposed 10
.... -
.
,._.~ -.
- ..............
- -..--....-:.~..
.......
-
E e~
.2 >
...... ..... ~-----0
0.0
t=0.14 d t = 0.27 d t = 0.40 d t = 0.60 d -t=l.0d
t = 8.0d Steady State I
I
I
i
0.2
0.4
0.6
0.8
Effective Saturation
(o)
i
i
1.0
at x = 8 m
Fig. 6. T e m p o r a l effective s a t u r a t i o n profiles w i t h d e p t h p r e d i c t e d by the v a r i a b l y s a t u r a t e d m o d e l at x = 8 m (c~v = 2 . 0 m - t , nv = 2.8 a n d K s = 5 . 1 t o d a y 1).
98
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
temporal snapshots of the vertical soil-moisture profile at the distance x = 8 m. It can be seen from Fig. 6 that at earlier times there is a substantial amount of water in the vadose zone. As time progresses, the soil drains and slowly relaxes toward the steadystate profile, which corresponds to the soil-water retention curve for the medium.
4. Capillary properties, pore sizes and permeability The manners in which the variably saturated model is influenced by its fundamental model parameters are addressed later in this communication, with particular attention given to the actual number of degrees-of-freedom involved in an informed sensitivity analysis. In order to lay the foundation for such a discussion, some basic principles relating capillarity, pore structure and permeability are addressed in this section. The capillary properties of a porous medium arise from the structure of the drainable porosity, the pore space between Or and 0s. In the present study, the capillary properties, and hence, descriptions of the pore structure, are parameterized by Van Genuchten's (1980) relation, and specifically by c~v and nv. At this point it is instructive to illustrate the sensitivites of the (drainable) pore-size-density function to variations in these parameters. Noting that the effective saturation, O, is the nondimensionalized cumulative measure of (drainable) pores filled at any given moisture content 0 (see Eq. (3)) and that smaller pores are filled first, the relation between O and the filled pore sizes can be expressed as: (12)
6) = F ( r )
where r is the largest pore radius wetted (L) at the effective saturation, 6), and F is the cumulative (drainable) pore-size-distribution function. Consequently, the (drainable) pore-size-density function, f (r), is directly obtained by differentiating equation (12): f(r) -
dF(r) dr
dO dr
(13)
The pressure head corresponds to a pore radius through the well known relation: -
2or cos/3 pgr
(14)
where a is the interfacial tension between air and water phases ( M T - 2 ) (here taken as 0.0725Nm-1), /3 is the contact angle formed in the air-water-solid system (here taken as zero), p is the density of water (ML -3) (here taken as 1000 kgm-3), and g is the gravitational constant (LT -2) (9.81ms-Z). Eqs. (2), (13), and (14) may be combined to obtain an analytical expression for the (drainable) pore-size-density function:
(15)
99
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
1.0
I
'
' ''|'1'I
'
I ' '''"I
'
~ (alpha= 0.5, n = 4.0) = ~" t - - - (alpha= 5.0, n=4.0) " ~" 0.8 1 ,~: ...... (alpha= 0.5
~,~
.~.Z~'t~ ~ll ~" 0.6
~,~ ~
0.4
' ' "]I'"I
'
I ' '''"I
] I /
I
.,.."":'
r" ..." ,,,..,ii, ..
0.0 -''°''" 104
'
10.~
10-6
10.5
10.4
10.3
Pore Radius, r (m)
Fig. 7. Normalized (drainable) pore-size density functions for three media: (c~v,nv)=(0.5m i, 4.0), (5.0m 1, 4.0) and (0.5m-1, 2.5). (Normalized by dividing by the maximum value of f(r) for the (0.5m ~, 4.0) soil.) Fig. 7 illustrates the (drainable) pore-size-density functions for soils parameterized by the following three (av, nv) pairs: (0.5m -l, 4.0), (5.0m -1, 4.0), and (0.5m -l , 2.5). These functions have all been non-dimensionalized by dividingf(r) for each medium by the peak value o f f ( r ) for the (0.5 m -1, 4.0) soil, in order to facilitate presentation in a unit-scaled graph. Note that the peak for the (5.0 m- l, 4.0) soil is 0.1 and is moved up along the abscissa exactly one order of magnitude from that of the (0.5 m -l, 4.0) soil. In log (r) space, these two density functions have the same width; naturally, therefore, the one on the right, (5.0m -l, 4.0), has a broader distribution when plotted in linear space. Consequently, the effects of increasing ~v by one order of magnitude are: a shift up of one order of magnitude in absolute pore size and a corresponding one order-of-magnitude decrease in the peak height of the (drainable) pore-size-density function to preserve the property of unit total probability. Comparing the (drainable) pore-size-density functions of the soils parameterized by (0.5 m 1,4.0) and (0.5 m -i, 2.5), one sees that the lower value ofnv corresponds to the broader density function with a somewhat lower, left-shifted peak. One remarkable disparity between these media is the relative abundance of smaller pores which accompanies the smaller value of nv. These smaller pores can serve as chokes, or throttles, within a porous medium, and thus restrict the flow of fluid through it. After all, it is the interconnectivity of larger pores which facilitates permeability for a given medium; the smallest pore through which water must flow when traversing a porous medium limits the permeability (Wise, 1992). Given the simplistic understanding of pore-size-density functions developed above, it is desirable to begin to extend this understanding to the realms of permeability. Most pore-structure based permeability models relate the contribution of a pore to the permeability of the medium of interest by a factor of the square of the pore radius
W.R. Wise et al. / Journal o f Hydrology 161 (1994) 91 108
100
10
'''
......
- -
=~:~.
I
•
'
'
'
'''"1
'
'
'''"I
,
,
i,
.... It
(alpha = 0.5, n = 4.0)
- - - (alpha = 5.0, n = 4.0) . . . . . . (alpha = 0.5, n = 2.5)
z,do_.
I l
~
6
~
I f
4
I
tI tL
2
fl
e,I '~"
A._I ,
0 lO-S
,
, , , , , d
lO"~
J
t
i
.....
I
IIl
I
10 .6
10.5
10.4
10 .3
Pore Radius, r (m) Fig. 8. N o r m a l i z e d potential permeability m e a s u r e s for three media: ( a v , n v ) = ( 0 . 5 m -1, 4.0), ( 5 . 0 m -1, 4.0), and ( 0 . 5 m -1, 2.5). (Normalized by dividing by the m a x i m u m vaule o f f ( r ) r 2 for the ( 0 . 5 m I, 4.0) soil.)
based upon a conjectured Poiseuille-type flow through the pore (see Dullien (1979)). Consequently, a measure of the potential contribution to permeability for a (drainable) pore-size-density function, f(r), is defined as f(r)r 2. Fig. 8 illustrates the permeability measures so defined for the same three soils whose (drainable) poresize-density functions are illustrated in Fig. 7. Again, these three functions are nondimensionalized by dividing by the peak value off(r)r 2 for the (0.5 m -1, 4.0) soil. Although the (drainable) pore-size-density function, f(r), for the larger value of av is one-tenth that for the smaller (same nv), the permeability measure, f(r)r 2, for the former medium is ten times that for the latter, owing to the r 2 factor and the one order-of-magnitude shift to the right for the larger value of av. Note that these curves do not represent actual 'permeability-density functions', but rather theoretical potential pore contributions to the permeability structure of these media. In fact, for the (0.5m -1, 2.5) soil, the smaller pores (sub-micron level) have basically no relative ability to transmit water and, hence, their potential contributions to the permeability, here measured byf(r)r 2, are negligible in Fig. 8. However, their potential for choking flow must be kept in mind, so neither Fig. 7 nor Fig. 8 are, in and of themselves, adequate to develop a comprehensive understanding of the permeability structure for these media. Wise (1992) identifies a methodology for identifying the smallest pore through which water must flow when traversing a given medium, the extreme path value radius, repv. Essentially, repv serves as the lower limit of the (drainable) pore-size-density function through which water is transmitted. Identification of the repv for the hypothetical media discussed in this work is not considered herein. However, the reader is reminded that there is probably at least one more parameter, in addition to av and nv, which influences the permeability for a given
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
101
porous medium. Detailed discussion of this third parameter is beyond the scope of this work. In the next section, the sensitivities of the variably saturated model to Ks, av and n~ are explored, first by varying these parameters independently (three degrees-of-freedom) and second using the above discussion to tie the capillary and permeability characteristics together, resulting in a system with two degrees-of-freedom. 5. Sensitivity analysis - - the traditional approach
In this section, the sensitivities of the location of the phreatic surface predicted by the variably saturated model to the variations in the saturated hydraulic conductivity and the Van Genuchten parameters are analyzed. The parameters Ks, ctv and nv are varied independently; there are three degrees-of-freedom within this sensitivity analysis.
5.1. Sensitivity to the saturated hydraulic conductivity, Ks Fig. 9 shows the locations of the transient phreatic surface at time t = 0.02 day, for the values of saturated hydraulic conductivity Ks = 0.051, 0.51, 5.1 and 51 mday 1, respectively. The values of c~v and nv are held constant at 1.0 m -1 and 2.8, respectively. The figure shows that as the value of the saturated hydraulic conductivity is increased, the soil tends to drain at a faster rate, and hence, the location of the phreatic surface is lowered more rapidly. This result is quite intuitive.
5.2. Sensitivity to the Van Genuchten parameter, %, Fig. 10 shows the sensitivity of the phreatic surface to variations in the Van 10
....... ==::~::~21122"----72~i~
....... i
e~
...... K =0.051 m/d . . . . . K =0.51 m/d - - - - K =5.1 m/d - - K = 5 1 . ngd 0
,
0
2
,
I
4
,
L
,
1
L
6
,
.
I
8
,
,
,
10
X Coordinate (m) Fig. 9. Sensitivity o f t h e v a r i a b l y s a t u r a t e d m o d e l s o l u t i o n (t = 0 . 0 2 d a y ) to v a r i a t i o n s in the s a t u r a t e d h y d r a u l i c c o n d u c t i v i t y , K s (c~v = 1 . 0 m -1 a n d nv = 2.8).
102
W.R. Wise et al. / Journal o f Hydrology 161 (1994) 91-108
°~
O~
...... alpha = 0,5 (l/m) . . . . . alpha = 1.0 (l/m) - - - - alpha = 2.0 (l/m) --
0
alpha = 4.0 (l/m)
2
4
6
8
10
X Coordinate (m) F i g . 10. S e n s i t i v i t y o f t h e v a r i a b l y s a t u r a t e d m o d e l s o l u t i o n G e n u c h t e n p a r a m e t e r a v (nv = 2.8 a n d Ks = 5.1 m d a y - l ) .
(t = 0 . 0 2 d a y )
to v a r i a t i o n s
in t h e v a n
Genuchten parameter av. The figure shows the location of the transient phreatic surface at time t = 0.02day, for the value of C~v = 0.5, 1.0, 2.0 and 4.0m -1 at constant values of K s = 5 . 1 m d a y -1 and nv = 2.8. As illustrated in Fig. 7, the Van Genuchten parameter c~v is a measure of absolute pore size. A porous medium with small pores has a small value of c~v and a large capillary fringe height. Soils with small values of av tend to retain considerable amounts of moisture within the vadose zone due to capillary forces. As seen in Fig. 10, as av decreases, the location of the phreatic surface is lowered more rapidly, since at small values of av more water is retained under tension, allowing the movement of the surface of atmospheric pressure well below the region of significant moisture. Verma and Brutsaert (1970) also observed that the simulated water table falls more rapidly (for a similar type of drainage problem) when examining the sensitivity to media composed of smaller pores. 5.3. Sensitivity to the Van Genuchten parameter, n~, Fig. 11 shows the sensitivity of the phreatic surface position to variations in the soil parameter nv at time t = 0.02day. The value of nv is varied from 1.4 to 4.5, holding constant the values of K s = 5.1 m d a y I and c~v = 1.0m --1. As illustrated in Fig. 7, the Van Genuchten soil parameter nv is an inverse measure of the breadth of the pore-size-density function; as nv decreases, the width of the pore-size-density function increases and vice versa. When the value of nv is reduced, the relative abundance of smaller pores (in comparison to the mean pore size) increases in a porous medium. These smaller pores are difficult to drain, due to their larger viscous effects and their higher affinity for water. Consequently, in Fig. 11, for smaller values of nv, the phreatic surface is lower, since more water is retained in the vadose zone.
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
103
2 ...... n=l.4 ..... n=l,8 ----n=3.0 --n=4.5 i
0
i
,
I
2
,
,
I
4
L
,
,
I
,
6
,
i
I
8
i
i
10
X Coordinate (m) Fig. 11. Sensitivity o f the v a r i a b l y s a t u r a t e d m o d e l solution (t = 0.02 day) to v a r i a t i o n s in the Van Genu c h t e n p a r a m e t e r nv (av = 1.0m -I a n d K s = 5.1 m d a y - ] ) .
Moreover, Mualem's (1976) relative hydraulic conductivity, subject to Van Genuchten's (1980) parameterization, Eq. (4), is a function of nv. The larger pores, which have less resistance to fluid transmission, tend to be drained first, constraining water to flow through the remaining (smaller) pores. Hence, as the value of nv decreases, the relative hydraulic conductivity of the unsaturated zone decreases more rapidly with decreasing water content, which decreases the net drainage rate in the unsaturated zone. Todsen (1973) and Verma and Brutsaert (1970) observed this trend that the simulated water table falls more rapidly when the pore-size-density function is broad, rather than narrow. Note that the above sensitivity analysis is based upon the assumption that the saturated hydraulic conductivity, Ks, is independent of the pore sizes parameterized by av and nv. The sensitivities observed by Todsen (1973) and Verma and Brutsaert (1970) were also the result of invoking this assumption. As mentioned above, common sense dictates that the hydraulic conductivity of a porous medium should be a direct function of the pore-size-density function, which, of course, also governs the capillaric properties of that medium. Below, an attempt is made to perform a more sophisticated sensitivity analysis, treating the hydraulic conductivity as a function of the Van Genuchten parameters, C~vand nv. The exact nature of this function is, however, speculative; the relation is based on very limited empirical evidence.
6. Sensitivity analysis - - a more realistic approach
Previous studies have indicated that the saturated hydraulic conductivity of a porous medium is a function of the pore structure of the medium (see Dullien (1979) for an extended discussion on pore structure models of permeability). Mishra
104
W.R. Wise et al./JournalofHydrology 161 (1994) 91-108 102 [] oL
m
10 l desired fit
E
10°
L.._a .~_
10-1
"o
n
e~ 10-z [] a
Mishra & Parker's (1990) Model Revised Model I
10"s 0
2
I
I
4
6
8
10
12
V a n G e n u c h t e n P a r a m e t e r nv Fig. 12. Comparison of Mishra and Parker's (1990) model, eq. (16), and its revised form, eq. (17), for predicting the saturated hydraulic conductivity of the media presented by Van Genuchten (1980). Note that these models are non-dimensionalized against the measured saturated hydraulic conductivities for these media.
and Parker (1990) proposed a model, based upon Mualem's (1976) expression for the relative hydraulic conductivity, which relates the saturated hydraulic conductivity to av, given by: KMp = 9.3(0s
- 0r)2"50~ 2
(16)
where KMp is the estimated value of saturated hydraulic conductivity, in metres per day, (0s - Or) is the drainable porosity, c~v is measured per metre and the factor of 9.3 lumps various information concerning the (assumed) geometry of porous media and the physical properties of water. This model was used to describe the saturated conductivity of 48 soils over a limited range of nv values, varying from 1.11 to 2.28 (Mishra and Parker, 1990). Wise (1991) suggested that this model should be refined to include nv, because he noted that there seemed to be a systematic scattering of Mishra and Parker's (1990) data, with respect to this parameter. This scattering effect is even more pronounced when examining data from Van Genuchten (1980), where nv varies from 1.17 to 10.4• In this study, Mishra and Parker's (1990) model is modified to account for the effects of nv, with the implied impact on the pore-size-density function• This refinement is empirical and is developed by fitting the saturated hydraulic conductivity to the Ks, c~v and nv data reported by Van Genuchten (1980). The modified model is written as: Ks ~ ~ 1 vlxMpnv3.5
(17)
where the symbol ~ reminds the reader of the empirical nature of the relation. Fig. 12 compares the performance of both these models (Eqs. (16) and (17)) in predicting the saturated hydraulic conductivity data presented by Van Genuchten (1980). The figure
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108 '
*
I
'
'
'
I
'
'
'
I
*
,
,
I
,
'
'
I
'
I
,
'
105
'
¢I
4
...... alpha: 0.5(l/m) ..... alpha= 1.0(l/m) ---- alpha: 2.0(I/m) -,
alpha=4.0(I/m) ,
,
I
0
,
.
,
2
[
,
4 6 X Coordinate (m)
8
10
Fig. 13. Sensitivity of the variably saturated model solution (t = 0.02 day) to variations in the Van Genuchten parameter c~v using the revised model, Eq. (17), for predicting the saturated hydraulic conductivity (nv = 2.8). s h o w s t h a t t h e r e v i s e d e m p i r i c a l e q u a t i o n i m p r o v e s t h e fit s i g n i f i c a n t l y w h e n c o m p a r e d to M i s h r a a n d P a r k e r ' s m o d e l f o r t h o s e d a t a . I n a g r e e m e n t w i t h t h e a b o v e d i s c u s s i o n r e g a r d i n g Figs. 7 a n d 8, E q . (17) r e l a t e s t h e s a t u r a t e d h y d r a u l i c c o n d u c t i v i t y to b o t h c~v a n d n v. T o t h e a u t h o r s ' k n o w l e d g e , this is t h e first t i m e t h a t this issue has b e e n explicitly a d d r e s s e d w h e n p e r f o r m i n g a sensitivity analysis. T h i s d i s c u s s i o n is a n initial effort to 10
,
,
[
'
,
,
[
'
,
,
I
'
i
i
[
'
,
i
'.,:::-:.,
...... n=l.4 ..... n=l.8 ----n:3.0 --n=4.5 L
0
I
'
L
2
4
6
8
10
X Coordinate (m) Fig. 14. Sensitivity of the variably saturated model solution (t = 0,02day) to variations in the Van Genuchten parameter nv using the revised model, Eq, (17), for predicting the saturated hydraulic conductivity (c~v = 1.0m-l).
106
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
come up with an approximate empirical formulation relating the properties of capillarity and permeability for porous media. More experimental data are needed in this area to verify and develop more robust models. In fact, the above mentioned extreme path value radius, repv, or a similarly defined parameter, should be included in a comprehensive model of saturated permeability. Such level of detail, however, is forgone in the present analysis for reasons of tractability. The sensitivity analysis now discussed is directed toward the influences ofa~ and nv (both with their intrinsic impacts upon Ks). Each component of the sensitivity analysis is carried out by evaluating the saturated hydraulic conductivity through Eq. (17) for different values of av and nv, thus restricting the analysis to two degrees-of-freedom. Figs. 13 and 14 show the locations of the phreatic surface at time t = 0.02 day, for different values of the soil parameters av and nv, respectively. The values of the saturated hydraulic conductivity computed from Eq. (17) for the values Ofav = 0.5, 1.0, 2.0, 4.0m -1 are Ks = 0.32, 1.3, 5.1 and 20.3mday -l, respectively (nv = 2.8). Similarly for the values ofnv = 1.4, 1.8, 3.0, 4.5 the computed values of hydraulic conductivity are, Ks = 0.11, 0.27, 1.6 and 6 . 7 m d a y -1, respectively (av = 1.0m-l). In both cases, as the values of the soil parameters av and nv increase, the saturated hydraulic conductivity increases. In the case of av, the larger values of av correspond to larger, more permeable pores. For nv, this effect is caused by the uniformity of the pore-size-density function achieved at larger values of nv; this uniformity eliminates smaller, less permeable pores from the pore structure of media so parameterized. These increased saturated hydraulic conductivity values lead to faster drainage rates and hence more rapid falls in the phreatic surface. The effects of these changes in the saturated hydraulic conductivity values seem to dominate the capillaric effects discussed in the previous section (contrast Figs. 13 and 14 with Figs. 10 and 11, respectively). Consequently, it appears that K s, av and nv should not be treated as three independent variables when performing a sensitivity analysis for a variably saturated flow model. Their interrelation significantly affects the predicted behavior and precludes the use of these three degrees-of-freedom in sensitivity investigations. This coupling of capillary and permeability effects may also make it less appropriate to average soil parameters within identified soil classes, as the parameters are not independently distributed (as assumed during the traditionally used averaging schemes), but rather, are jointly distributed as a result of the fact that pore-size-density functions govern both capillary and permeability phenomena.
7. Conclusions
It is important to understand how variations in the soil parameters of variably saturated flow models affect simulations. Heuristic arguments support a claim that the capillaric properties of a porous medium should be closely related to the permeability of that medium. A sensitivity analysis based upon this principle reveals that simulation results are more sensitive to the permeability than to the capillarity of porous media. In fact, when the capillary properties are tied to the saturated hydraulic conductivity through a preliminary model, the sensitivities to the two capillary
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
107
parameters are, in general, reversed from that obtained with three degrees-offreedom, due to the manner in which variations in these parameters affect the permeability. The relationship between the descriptive soil parameters, which indirectly partially parameterize the pore structure, and the saturated hydraulic conductivity is still an open area in soil physics; more experimental investigations are needed to develop robust models. However, from the preliminary work performed herein, consideration of the number of degrees-of-freedom used when performing a sensitivity analysis is shown to warrant great care, if the predicted sensitivities are to have meaningful physical interpretations. Hopefully, the results of this research will provide a stepping stone to a more realistic understanding of the transient behavior involved with variably saturated flow problems.
Acknowledgment This work was supported, in part, by the US Environmental Protection Agency (contrast CR-818717-01-0) through the Environmental Research Laboratory in Athens, Georgia. However, it has not been through the agency's required peer and administrative review. Therefore, it does not necessarily reflect the views of the agency, and no official endorsement should be inferred.
References Ahmad, S., Kashyap, D. and Mathur, B.S., 1993. Mathematical modeling of saturated-unsaturated flow to drains. J. Irrig. Drain. Eng., 119: 18-33. Clement, T.P., Wise, W.R. and Molz, F.J., 1994. A physically based, two-dimensional finite-difference algorithm for variably saturated flow. J. Hydrol., 161 :, 71-90. Clement, T.P., Wise, W.R., Molz, F.J. and Wen, M., 1994. A comprehensive study of the modeling of steady-state seepage faces. Water Resour. Res., in press. Dullien, F.A.L., 1979. Porous Media Fluid Transport and Pore Structure. Academic Press, San Diego, 396 pp. Liggett, J.A., 1977. Location of free surface in porous media. J. Hydraul. Div., 103(HY4): 353-365. Mishra, S. and Parker, J.C., 1990. On the relation between saturated conductivity and capillary retention characteristics. Ground Water, 28: 775-777. M ualem, Y., 1976. A new model for predicting hydraulic conductivity of unsaturated porous media. Water Resour. Res., 12: 513-522. Rubin, J., 1968. Theoretical analysis of two-dimensional, transient flow of water in unsaturated and partly saturated soils. Soil Sci. Soc. Am. Proc,, 32: 607-615. Sisson, J.B., Ferguson, A.H. and Van Genuchten, M.Th., 1980. Simple method for predicting drainage from field plots. Soil Sci. Soc. Am. J., 44:1147 1152. Todsen, M., 1973. Numerical studies of two-dimensional saturated/unsaturated drainage models. J. Hydrol., 20:311-326. Van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J., 44: 892-898. Vauclin, M., Vachaud, G. and Khanji, J., 1975. Two dimensional numerical analysis of transient water transfer in saturated-unsaturated soils. In: G.C. Vansteenkiste (Editor), Modeling and Simulation of Water Resources Systems, North-Holland, Amsterdam, pp. 299-323.
108
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
Vauclin, M., Khanji, D. and Vachaud, G., 1979. Experimental and numerical study of a transient, twodimensional unsaturated-saturated water table recharge problem. Water Resour. Res., 15:1089 1101. Verma, R.D. and Brutsaert, W., 1970. Unconfined aquifer seepage by capillary flow theory. J. Hydraul. Div., 96(HY 6): 1331-1344. Watson, K.K. and Awadalla, S.A., 1985. Comparative study of Green Ampt analysis for a falling water table in a homogeneous sand profile. Water Resour. Res., 21:1157-1164. Wise, W.R., 1991. Discussion of 'On the relation between saturated conductivity and capillary retention characteristics,' by S. Mishra and J.C. Parker, September-October 1990 issue, vol. 28, no. 5, pp. 775777. Ground Water, 29: 272-273. Wise, W.R., 1992. A new insight on pore structure and permeability. Water Resour. Res., 28: 189-198. Young, E.G. and Aggelides, S., 1976. Drainage to a water table analysed by the Green-Ampt approach. J. Hydrol., 31:67 79.
Hydrology ELSEVIER
Journal of Hydrology 161 (1994) 91-108
[1]
Variably saturated modeling of transient drainage: sensitivity to soil properties W i l l i a m R. W i s e * , T . P . C l e m e n t l, F r e d J. M o l z Department of Civil Engineering, 238 Harbert Engineering Center, Auburn University, Auburn, AL 36849-5337, USA
Received 13 September 1993; revision accepted 16 March 1994
Abstract
Numerical analysis of a transient, two-dimensional (rectangular symmetry), unconfineddrainage problem, using a variably saturated flow model, shows that the location of the phreatic surface and the height of the seepage face are functions of the capillary forces exerted in the vadose zone. A sensitivity analysis investigating the effects of independent variations in the saturated hydraulic conductivity and the two Van Genuchten parameters is compared to one performed where the saturated hydraulic conductivity is functionally related to the Van Genuchten parameters. The relationship between these descriptive soil parameters, which may be used to describe the pore-size-density function of a porous medium, and the saturated hydraulic conductivity is surmised based upon limited empirical evidence. Expressing the permeability as a function of these soil properties reverses the sensitivities of the variably saturated model to the Van Genuchten parameters. From the preliminary work performed herein, consideration of the number of degrees-of-freedom used when performing a sensitivity analysis is shown to warrant great care, if the predicted sensitivities are to have meaningful physical interpretations.
I. Introduction
Numerical models are frequently employed to predict the transient behavior o f subsurface water. These models are developed based u p o n mathematical descriptions o f the physical processes considered. Consequently, the quality o f any particular numerical model directly depends u p o n the accuracy and completeness o f the physical understanding upon which the model is built. This linkage is of great significance when modeling * Corresponding author. l Present address: Battelle, Pacific Northwest Laboratories, Battelle Boulevard, P.O. Box 999, Richland, WA 99352, USA. 0022-1694/94/$07.00 © 1994 - Elsevier Science B.V. All rights reserved SSDI 0022-1694(94)02509-A
92
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
transient unconfined flow through the subsurface, where the mathematical descriptions are highly nonlinear. This communication is aimed at shedding insight into how a sensitivity analysis should be properly conducted for variably saturated models in order to maintain an alliance with the underlying physics of flow. While many studies of steady-state seepage faces have been performed (see Clement et al. (1994) for an extended discusssion), fewer transient problems have been investigated. Most of the transient unconfined drainage problems are restricted to onedimensional analyses (e.g. Young and Aggelides, 1976; Sisson et al., 1980; Watson and Awadalla, 1985). A few studies have investigated two-dimensional drainage processes using variably saturated flow models (e.g. Rubin, 1968; Verma and Brutsaert, 1970; Todsen, 1973; Vauclin et al., 1979; Ahmad et al., 1993). These studies reveal the importance of considering transient flow processes in the unsaturated zone and demonstrate the failure of fully saturated flow models (such as that posed by Liggett (1977), for example) in simulating these types of infiltration curedrainage problems. Common failings of fully saturated flow models for describing unconfined drainage problems include the observation that the fully saturated flow models predict slower responses in water tables than actually occur. Vauclin et al. (1975) present an experimental investigation of transient seepage-face and phreatic-surface positions in a two-dimensional (rectangular) flow problem. A 6 m × 2 m saturated porous medium was allowed to drain after a sudden drop in the external water table. This sudden drop problem was investigated by measuring the water content and pressure head distribution during the drainage phase. The variably saturated algorithm used in this work is verified successfully by simulating Vauclin et al.'s (1975) experimental results (Clement et al., 1994). Within this paper, certain terms are used as precisely as possible. The terms phreatic surface and water table are used interchangeably to describe the surface connecting the locus of points at atmospheric pressure, excluding those along the seepage face. The term vadose zone is used to describe every part of the soil system residing above the phreatic surface. Thus the capillary fringe, which is saturated, is included in the vadose zone. The term unsaturated zone refers to that part of the vadose zone residing above the capillary fringe. As the name implies, the saturation within the unsaturated zone is less than unity; within the vadose zone, by contrast, it is less than or equal to unity. The term permeability is used to refer to the ability of a porous medium to transmit a fluid which saturates the medium. By contrast, the term hydraulic conductivity is used to describe the transmission of water through the medium at a given water content under a unit hydraulic gradient. Precise use of these terms is helpful below, when unraveling some of the more subtle effects involved with modeling variably saturated flow. The objectives of the present work are to perform and contrast sensitivity analyses of the variably saturated model to the soil properties (the two Van Genuchten parameters and the saturated hydraulic conductivity) using two methods, one employing three degrees-of-freedom, the other, two. In order to accomplish the second part of this objective, a new empirical way to represent the saturated hydraulic conductivity as a function of the Van Genuchten parameters is introduced. To set the tone for the modeling exercises reported, the problem at hand is developed below.
93
W.R. Wise et al. / Journal o f Hydrology 161 (1994) 91 108 V B
C !
p(x, t)
V h2
Z
±
A
L Fig. 1. Schematic of the flow domain used in all simulations presented in this paper (L = 10m, hI = 10m and h 2 = 3 m for all simulations).
2. Problem description A schematic of the problem used to pursue the stated objective is presented in Fig. 1. The 10 m x 10 m square e m b a n k m e n t has no-flow boundaries on the base and on the top of the flow domain. The porous medium is assumed to be homogeneous and isotropic. Values of the soil parameters used in the test problem are: a saturated hydraulic conductivity of Ks = 5.1 m day -1, a porosity of 77 = 0.46, a residual water content of Or = 0.01, and a specific storage of Ss = 5 x 10 -5 m - I. In the base case, the values of the Van Genuchten soil parameters used are: C~v = 2.0 m -1 and nv = 2.8, which represent a poorly sorted, fine sand. The water table is initially at the top of the block and the flow system is at static equilibrium. The water table on the left boundary is maintained at a constant level of z = 10 m, while that on the right face of the block is lowered instantaneously (at t = 0) from z = 10 m to z = 3 m and is maintained there, subsequently. Transient locations of the phreatic surface and the seepage face are determined numerically. Due to the instantaneous drop in the external water level on the right boundary, a rapid drop in pressure occurs along the right edge of the block, resulting in large pressure gradients. Pressure changes inside the block outpace the drainage of the pores. Consequently, due to the vertical extent of the problem at hand, specific storage effects are important, especially at early times. (Unfortunately, no experimental data is found in the literature for analyzing the specific storage effects.)
3. Modeling two-dimensional transient flow with a seepage-face boundary The governing equation for water flow in a two-dimensional (rectangular
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
94
symmetry), variably saturated porous medium is written (Clement et al., 1994):
00~
--[ O0
Ss- --g? + ot
I
0 K(o) O~;
ox
0 KO)~
+--
Oz
(1)
where Ss is the specific storage of the medium (L-l), 0 is the water content, ~/is the porosity, K(O) is the hydraulic conductivity (LT-1), t is time (T), • = h - z, where is the pressure head (L), h is the hydraulic head (L), and x and z are the horizontal and vertical coordinates respectively (L). The soil properties 0(q), and K(O) are herein described using Van Genuchten's (1980) soil-water-retention function and Mualem's (1976) variably saturated hydraulic conductivity function, respectively. The relationship between water content and pressure head (under tension) is given by (Van Genuchten, 1980): 1
mv
where c~v, nv and mv = 1 - 1/nv are the Van Genuchten parameters, whose values depend upon the soil properties. The parameter c~v is, in part, a measure of the first moment of the pore-size-density function (L -1) (as ctv increases, so does the first moment); nv is an inverse measure of the second moment of the pore size-density function (as nv increases, the pore-size-density function becomes narrower) (Wise, 1991). O is the effective saturation given by the relation: 0 - 0r O- - 0s - Or
(3)
where 0s is the saturated water content (often approximated by the porosity 7/) and Or is the residual water content of the soil. Of course, for non-negative ffJ, 6) = 1, that is, the medium is saturated. Based on Mualem's (1976) model, the relation between water content and hydraulic conductivity is given by (Van Genuchten, 1980): K(O) = Ks{1 - [1 - O(1/mv)]mv}2(~ 1/2
(4)
where K s is the saturated hydraulic conductivity (LT l). Note that the relative hydraulic conductivity, K(O)/Ks, does not depend upon the value of c~v. The initial and boundary conditions for the variably saturated, transient drainage problem generically shown in Fig. 1 are as follows: h =-h i
Vx, Vz, t = 0 (initial condition; hydrostatic equilibrium),
(5)
h=hl
along AB, (constant-head boundary),
(6)
along BC, (no-flow boundary),
(7)
along CD, (no-flow boundary),
(8)
along DE, (seepage-face boundary),
(9)
Oh Oz Oh --~0 Ox
--~0
h=z
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
h = h2
Oh
0~ = 0
95
along EF, (constant-head boundary), and
(10)
along AF, (no-flow boundary).
(11)
Numerical solution to the variably saturated problem discussed earlier is accomplished by using the numerical algorithm presented by Clement et al. (1994). The nodal spacing in the x and z directions used is equal, Ax = Az = 0.2 m. A variable time step, At, ranging from 0.0001 day to 0.1 day, is used in all simulations. The details of the finite-difference solution techniques, implementation of the boundary conditions, and verification of the algorithm are discussed by Clement et al. (1994).
3.1. Results of the variably saturated model for the sample problem Shown in Fig. 2 are temporal profiles of the phreatic surface predicted by the variably saturated model for the base case. As seen from the figure, the decline of the phreatic surface is quite rapid during earlier times and slows at later times. Fig. 3 shows the distribution of the drainage rates due to the various processes in the early stages of the simulation. This figure illustrates that the early drainage (0 < t < 0.00012 day) is dominated by the specific storage effects which accommodate the swift pressure variations within the porous medium and the resulting rapid decline in the phreatic surface. However, as time progresses, the specific storage effects decrease rapidly and the pore drainage begins to have a more pronounced influence on the flow. The specific storage effects are negligible after about 0.0004 day (35 s). Because the specific storage effects are so short lived, there is no compelling reason to
•
"::
:22.2 .........
- ..................
o ¢I
...... t = O.OOff2d ..... t=O.O2d ----t=O.14d - - - t =0.60d SteadyStar, -
-
i
i
i
I
2
,
,
,
I
4
,
,
,
I
L
6
,
,
I
8
,
,
,
I0
X Coordinate (m) Fig. 2. T e m p o r a l p h r e a t i c - s u r f a c e profiles predicted by the v a r i a b l y s a t u r a t e d m o d e l (av = 2 . 0 m - I , nv = 2.8 a n d Ks = 5.1 m d a y - 1 ) .
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
96
250
"~
150
loo
50
o 0.0000
0.0001
0.0002
0.0003
0.0004
Elapsed Time (d) Fig. 3. F l o w histories d u e to the v a r i o u s c o n t r i b u t o r y p r o c e s s e s - - e a r l y t i m e (av = 2.0 m -1 , n v = 2.8 a n d K s = 5.1 m d a y - l ) .
explore the sensitivity of the model with respect to the specific storage, Ss. Fig. 4 illustrates the inflow, outflow and pore-drainage flow components at later times.
3.2. Observations of the drainage behavior The soil moisture distribution is always continuous, and, depending on the soil 50
,
,
,
,
<
40
"~
30
'-
20
Outflow Inflow
10
o
Pore Drainage i
0 0.0
i
i
I
0.2
i
i
i
I
0.4
l
i
i
I
i
0.6
i
l
[
0.8
,
l
l
1.0
Elapsed Time (d) Fig. 4. F l o w h i s t o r i e s d u e to the v a r i o u s c o n t r i b u t o r y p r o c e s s e s - - l a t e r t i m e (av = 2 . 0 m - l , nv = 2.8 a n d Ks = 5.1 m d a y - l ) .
W.R. Wise et al. / Journal o[Hydrology 161 (1994) 91-108
PhreaticSurface/4 " ~
97
100 %
¢
0
2
4
6
8
10
X Coordinate (m) Fig. 5. C o n t o u r d i a g r a m o f the m o i s t u r e c o n t e n t p r e d i c t e d by the v a r i a b l y s a t u r a t e d m o d e l at t = 0.14 d a y (C~v = 2 . 0 m - l , nv = 2.8 a n d Ks = 5.1 m d a y - l ) .
properties, the vadose zone above the phreatic surface can store a significant amount of water. Fig. 5 is a contour diagram of the effective saturation at time t = 0.14 day, predicted by the variably saturated model. The location of the phreatic surface predicted by the model at this time is also shown in the figure. This figure clearly illustrates that the porous medium retains significant moisture well above the phreatic surface. Similar observations can also be made from Fig. 6, which shows superposed 10
.... -
.
,._.~ -.
- ..............
- -..--....-:.~..
.......
-
E e~
.2 >
...... ..... ~-----0
0.0
t=0.14 d t = 0.27 d t = 0.40 d t = 0.60 d -t=l.0d
t = 8.0d Steady State I
I
I
i
0.2
0.4
0.6
0.8
Effective Saturation
(o)
i
i
1.0
at x = 8 m
Fig. 6. T e m p o r a l effective s a t u r a t i o n profiles w i t h d e p t h p r e d i c t e d by the v a r i a b l y s a t u r a t e d m o d e l at x = 8 m (c~v = 2 . 0 m - t , nv = 2.8 a n d K s = 5 . 1 t o d a y 1).
98
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
temporal snapshots of the vertical soil-moisture profile at the distance x = 8 m. It can be seen from Fig. 6 that at earlier times there is a substantial amount of water in the vadose zone. As time progresses, the soil drains and slowly relaxes toward the steadystate profile, which corresponds to the soil-water retention curve for the medium.
4. Capillary properties, pore sizes and permeability The manners in which the variably saturated model is influenced by its fundamental model parameters are addressed later in this communication, with particular attention given to the actual number of degrees-of-freedom involved in an informed sensitivity analysis. In order to lay the foundation for such a discussion, some basic principles relating capillarity, pore structure and permeability are addressed in this section. The capillary properties of a porous medium arise from the structure of the drainable porosity, the pore space between Or and 0s. In the present study, the capillary properties, and hence, descriptions of the pore structure, are parameterized by Van Genuchten's (1980) relation, and specifically by c~v and nv. At this point it is instructive to illustrate the sensitivites of the (drainable) pore-size-density function to variations in these parameters. Noting that the effective saturation, O, is the nondimensionalized cumulative measure of (drainable) pores filled at any given moisture content 0 (see Eq. (3)) and that smaller pores are filled first, the relation between O and the filled pore sizes can be expressed as: (12)
6) = F ( r )
where r is the largest pore radius wetted (L) at the effective saturation, 6), and F is the cumulative (drainable) pore-size-distribution function. Consequently, the (drainable) pore-size-density function, f (r), is directly obtained by differentiating equation (12): f(r) -
dF(r) dr
dO dr
(13)
The pressure head corresponds to a pore radius through the well known relation: -
2or cos/3 pgr
(14)
where a is the interfacial tension between air and water phases ( M T - 2 ) (here taken as 0.0725Nm-1), /3 is the contact angle formed in the air-water-solid system (here taken as zero), p is the density of water (ML -3) (here taken as 1000 kgm-3), and g is the gravitational constant (LT -2) (9.81ms-Z). Eqs. (2), (13), and (14) may be combined to obtain an analytical expression for the (drainable) pore-size-density function:
(15)
99
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
1.0
I
'
' ''|'1'I
'
I ' '''"I
'
~ (alpha= 0.5, n = 4.0) = ~" t - - - (alpha= 5.0, n=4.0) " ~" 0.8 1 ,~: ...... (alpha= 0.5
~,~
.~.Z~'t~ ~ll ~" 0.6
~,~ ~
0.4
' ' "]I'"I
'
I ' '''"I
] I /
I
.,.."":'
r" ..." ,,,..,ii, ..
0.0 -''°''" 104
'
10.~
10-6
10.5
10.4
10.3
Pore Radius, r (m)
Fig. 7. Normalized (drainable) pore-size density functions for three media: (c~v,nv)=(0.5m i, 4.0), (5.0m 1, 4.0) and (0.5m-1, 2.5). (Normalized by dividing by the maximum value of f(r) for the (0.5m ~, 4.0) soil.) Fig. 7 illustrates the (drainable) pore-size-density functions for soils parameterized by the following three (av, nv) pairs: (0.5m -l, 4.0), (5.0m -1, 4.0), and (0.5m -l , 2.5). These functions have all been non-dimensionalized by dividingf(r) for each medium by the peak value o f f ( r ) for the (0.5 m -1, 4.0) soil, in order to facilitate presentation in a unit-scaled graph. Note that the peak for the (5.0 m- l, 4.0) soil is 0.1 and is moved up along the abscissa exactly one order of magnitude from that of the (0.5 m -l, 4.0) soil. In log (r) space, these two density functions have the same width; naturally, therefore, the one on the right, (5.0m -l, 4.0), has a broader distribution when plotted in linear space. Consequently, the effects of increasing ~v by one order of magnitude are: a shift up of one order of magnitude in absolute pore size and a corresponding one order-of-magnitude decrease in the peak height of the (drainable) pore-size-density function to preserve the property of unit total probability. Comparing the (drainable) pore-size-density functions of the soils parameterized by (0.5 m 1,4.0) and (0.5 m -i, 2.5), one sees that the lower value ofnv corresponds to the broader density function with a somewhat lower, left-shifted peak. One remarkable disparity between these media is the relative abundance of smaller pores which accompanies the smaller value of nv. These smaller pores can serve as chokes, or throttles, within a porous medium, and thus restrict the flow of fluid through it. After all, it is the interconnectivity of larger pores which facilitates permeability for a given medium; the smallest pore through which water must flow when traversing a porous medium limits the permeability (Wise, 1992). Given the simplistic understanding of pore-size-density functions developed above, it is desirable to begin to extend this understanding to the realms of permeability. Most pore-structure based permeability models relate the contribution of a pore to the permeability of the medium of interest by a factor of the square of the pore radius
W.R. Wise et al. / Journal o f Hydrology 161 (1994) 91 108
100
10
'''
......
- -
=~:~.
I
•
'
'
'
'''"1
'
'
'''"I
,
,
i,
.... It
(alpha = 0.5, n = 4.0)
- - - (alpha = 5.0, n = 4.0) . . . . . . (alpha = 0.5, n = 2.5)
z,do_.
I l
~
6
~
I f
4
I
tI tL
2
fl
e,I '~"
A._I ,
0 lO-S
,
, , , , , d
lO"~
J
t
i
.....
I
IIl
I
10 .6
10.5
10.4
10 .3
Pore Radius, r (m) Fig. 8. N o r m a l i z e d potential permeability m e a s u r e s for three media: ( a v , n v ) = ( 0 . 5 m -1, 4.0), ( 5 . 0 m -1, 4.0), and ( 0 . 5 m -1, 2.5). (Normalized by dividing by the m a x i m u m vaule o f f ( r ) r 2 for the ( 0 . 5 m I, 4.0) soil.)
based upon a conjectured Poiseuille-type flow through the pore (see Dullien (1979)). Consequently, a measure of the potential contribution to permeability for a (drainable) pore-size-density function, f(r), is defined as f(r)r 2. Fig. 8 illustrates the permeability measures so defined for the same three soils whose (drainable) poresize-density functions are illustrated in Fig. 7. Again, these three functions are nondimensionalized by dividing by the peak value off(r)r 2 for the (0.5 m -1, 4.0) soil. Although the (drainable) pore-size-density function, f(r), for the larger value of av is one-tenth that for the smaller (same nv), the permeability measure, f(r)r 2, for the former medium is ten times that for the latter, owing to the r 2 factor and the one order-of-magnitude shift to the right for the larger value of av. Note that these curves do not represent actual 'permeability-density functions', but rather theoretical potential pore contributions to the permeability structure of these media. In fact, for the (0.5m -1, 2.5) soil, the smaller pores (sub-micron level) have basically no relative ability to transmit water and, hence, their potential contributions to the permeability, here measured byf(r)r 2, are negligible in Fig. 8. However, their potential for choking flow must be kept in mind, so neither Fig. 7 nor Fig. 8 are, in and of themselves, adequate to develop a comprehensive understanding of the permeability structure for these media. Wise (1992) identifies a methodology for identifying the smallest pore through which water must flow when traversing a given medium, the extreme path value radius, repv. Essentially, repv serves as the lower limit of the (drainable) pore-size-density function through which water is transmitted. Identification of the repv for the hypothetical media discussed in this work is not considered herein. However, the reader is reminded that there is probably at least one more parameter, in addition to av and nv, which influences the permeability for a given
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
101
porous medium. Detailed discussion of this third parameter is beyond the scope of this work. In the next section, the sensitivities of the variably saturated model to Ks, av and n~ are explored, first by varying these parameters independently (three degrees-of-freedom) and second using the above discussion to tie the capillary and permeability characteristics together, resulting in a system with two degrees-of-freedom. 5. Sensitivity analysis - - the traditional approach
In this section, the sensitivities of the location of the phreatic surface predicted by the variably saturated model to the variations in the saturated hydraulic conductivity and the Van Genuchten parameters are analyzed. The parameters Ks, ctv and nv are varied independently; there are three degrees-of-freedom within this sensitivity analysis.
5.1. Sensitivity to the saturated hydraulic conductivity, Ks Fig. 9 shows the locations of the transient phreatic surface at time t = 0.02 day, for the values of saturated hydraulic conductivity Ks = 0.051, 0.51, 5.1 and 51 mday 1, respectively. The values of c~v and nv are held constant at 1.0 m -1 and 2.8, respectively. The figure shows that as the value of the saturated hydraulic conductivity is increased, the soil tends to drain at a faster rate, and hence, the location of the phreatic surface is lowered more rapidly. This result is quite intuitive.
5.2. Sensitivity to the Van Genuchten parameter, %, Fig. 10 shows the sensitivity of the phreatic surface to variations in the Van 10
....... ==::~::~21122"----72~i~
....... i
e~
...... K =0.051 m/d . . . . . K =0.51 m/d - - - - K =5.1 m/d - - K = 5 1 . ngd 0
,
0
2
,
I
4
,
L
,
1
L
6
,
.
I
8
,
,
,
10
X Coordinate (m) Fig. 9. Sensitivity o f t h e v a r i a b l y s a t u r a t e d m o d e l s o l u t i o n (t = 0 . 0 2 d a y ) to v a r i a t i o n s in the s a t u r a t e d h y d r a u l i c c o n d u c t i v i t y , K s (c~v = 1 . 0 m -1 a n d nv = 2.8).
102
W.R. Wise et al. / Journal o f Hydrology 161 (1994) 91-108
°~
O~
...... alpha = 0,5 (l/m) . . . . . alpha = 1.0 (l/m) - - - - alpha = 2.0 (l/m) --
0
alpha = 4.0 (l/m)
2
4
6
8
10
X Coordinate (m) F i g . 10. S e n s i t i v i t y o f t h e v a r i a b l y s a t u r a t e d m o d e l s o l u t i o n G e n u c h t e n p a r a m e t e r a v (nv = 2.8 a n d Ks = 5.1 m d a y - l ) .
(t = 0 . 0 2 d a y )
to v a r i a t i o n s
in t h e v a n
Genuchten parameter av. The figure shows the location of the transient phreatic surface at time t = 0.02day, for the value of C~v = 0.5, 1.0, 2.0 and 4.0m -1 at constant values of K s = 5 . 1 m d a y -1 and nv = 2.8. As illustrated in Fig. 7, the Van Genuchten parameter c~v is a measure of absolute pore size. A porous medium with small pores has a small value of c~v and a large capillary fringe height. Soils with small values of av tend to retain considerable amounts of moisture within the vadose zone due to capillary forces. As seen in Fig. 10, as av decreases, the location of the phreatic surface is lowered more rapidly, since at small values of av more water is retained under tension, allowing the movement of the surface of atmospheric pressure well below the region of significant moisture. Verma and Brutsaert (1970) also observed that the simulated water table falls more rapidly (for a similar type of drainage problem) when examining the sensitivity to media composed of smaller pores. 5.3. Sensitivity to the Van Genuchten parameter, n~, Fig. 11 shows the sensitivity of the phreatic surface position to variations in the soil parameter nv at time t = 0.02day. The value of nv is varied from 1.4 to 4.5, holding constant the values of K s = 5.1 m d a y I and c~v = 1.0m --1. As illustrated in Fig. 7, the Van Genuchten soil parameter nv is an inverse measure of the breadth of the pore-size-density function; as nv decreases, the width of the pore-size-density function increases and vice versa. When the value of nv is reduced, the relative abundance of smaller pores (in comparison to the mean pore size) increases in a porous medium. These smaller pores are difficult to drain, due to their larger viscous effects and their higher affinity for water. Consequently, in Fig. 11, for smaller values of nv, the phreatic surface is lower, since more water is retained in the vadose zone.
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
103
2 ...... n=l.4 ..... n=l,8 ----n=3.0 --n=4.5 i
0
i
,
I
2
,
,
I
4
L
,
,
I
,
6
,
i
I
8
i
i
10
X Coordinate (m) Fig. 11. Sensitivity o f the v a r i a b l y s a t u r a t e d m o d e l solution (t = 0.02 day) to v a r i a t i o n s in the Van Genu c h t e n p a r a m e t e r nv (av = 1.0m -I a n d K s = 5.1 m d a y - ] ) .
Moreover, Mualem's (1976) relative hydraulic conductivity, subject to Van Genuchten's (1980) parameterization, Eq. (4), is a function of nv. The larger pores, which have less resistance to fluid transmission, tend to be drained first, constraining water to flow through the remaining (smaller) pores. Hence, as the value of nv decreases, the relative hydraulic conductivity of the unsaturated zone decreases more rapidly with decreasing water content, which decreases the net drainage rate in the unsaturated zone. Todsen (1973) and Verma and Brutsaert (1970) observed this trend that the simulated water table falls more rapidly when the pore-size-density function is broad, rather than narrow. Note that the above sensitivity analysis is based upon the assumption that the saturated hydraulic conductivity, Ks, is independent of the pore sizes parameterized by av and nv. The sensitivities observed by Todsen (1973) and Verma and Brutsaert (1970) were also the result of invoking this assumption. As mentioned above, common sense dictates that the hydraulic conductivity of a porous medium should be a direct function of the pore-size-density function, which, of course, also governs the capillaric properties of that medium. Below, an attempt is made to perform a more sophisticated sensitivity analysis, treating the hydraulic conductivity as a function of the Van Genuchten parameters, C~vand nv. The exact nature of this function is, however, speculative; the relation is based on very limited empirical evidence.
6. Sensitivity analysis - - a more realistic approach
Previous studies have indicated that the saturated hydraulic conductivity of a porous medium is a function of the pore structure of the medium (see Dullien (1979) for an extended discussion on pore structure models of permeability). Mishra
104
W.R. Wise et al./JournalofHydrology 161 (1994) 91-108 102 [] oL
m
10 l desired fit
E
10°
L.._a .~_
10-1
"o
n
e~ 10-z [] a
Mishra & Parker's (1990) Model Revised Model I
10"s 0
2
I
I
4
6
8
10
12
V a n G e n u c h t e n P a r a m e t e r nv Fig. 12. Comparison of Mishra and Parker's (1990) model, eq. (16), and its revised form, eq. (17), for predicting the saturated hydraulic conductivity of the media presented by Van Genuchten (1980). Note that these models are non-dimensionalized against the measured saturated hydraulic conductivities for these media.
and Parker (1990) proposed a model, based upon Mualem's (1976) expression for the relative hydraulic conductivity, which relates the saturated hydraulic conductivity to av, given by: KMp = 9.3(0s
- 0r)2"50~ 2
(16)
where KMp is the estimated value of saturated hydraulic conductivity, in metres per day, (0s - Or) is the drainable porosity, c~v is measured per metre and the factor of 9.3 lumps various information concerning the (assumed) geometry of porous media and the physical properties of water. This model was used to describe the saturated conductivity of 48 soils over a limited range of nv values, varying from 1.11 to 2.28 (Mishra and Parker, 1990). Wise (1991) suggested that this model should be refined to include nv, because he noted that there seemed to be a systematic scattering of Mishra and Parker's (1990) data, with respect to this parameter. This scattering effect is even more pronounced when examining data from Van Genuchten (1980), where nv varies from 1.17 to 10.4• In this study, Mishra and Parker's (1990) model is modified to account for the effects of nv, with the implied impact on the pore-size-density function• This refinement is empirical and is developed by fitting the saturated hydraulic conductivity to the Ks, c~v and nv data reported by Van Genuchten (1980). The modified model is written as: Ks ~ ~ 1 vlxMpnv3.5
(17)
where the symbol ~ reminds the reader of the empirical nature of the relation. Fig. 12 compares the performance of both these models (Eqs. (16) and (17)) in predicting the saturated hydraulic conductivity data presented by Van Genuchten (1980). The figure
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108 '
*
I
'
'
'
I
'
'
'
I
*
,
,
I
,
'
'
I
'
I
,
'
105
'
¢I
4
...... alpha: 0.5(l/m) ..... alpha= 1.0(l/m) ---- alpha: 2.0(I/m) -,
alpha=4.0(I/m) ,
,
I
0
,
.
,
2
[
,
4 6 X Coordinate (m)
8
10
Fig. 13. Sensitivity of the variably saturated model solution (t = 0.02 day) to variations in the Van Genuchten parameter c~v using the revised model, Eq. (17), for predicting the saturated hydraulic conductivity (nv = 2.8). s h o w s t h a t t h e r e v i s e d e m p i r i c a l e q u a t i o n i m p r o v e s t h e fit s i g n i f i c a n t l y w h e n c o m p a r e d to M i s h r a a n d P a r k e r ' s m o d e l f o r t h o s e d a t a . I n a g r e e m e n t w i t h t h e a b o v e d i s c u s s i o n r e g a r d i n g Figs. 7 a n d 8, E q . (17) r e l a t e s t h e s a t u r a t e d h y d r a u l i c c o n d u c t i v i t y to b o t h c~v a n d n v. T o t h e a u t h o r s ' k n o w l e d g e , this is t h e first t i m e t h a t this issue has b e e n explicitly a d d r e s s e d w h e n p e r f o r m i n g a sensitivity analysis. T h i s d i s c u s s i o n is a n initial effort to 10
,
,
[
'
,
,
[
'
,
,
I
'
i
i
[
'
,
i
'.,:::-:.,
...... n=l.4 ..... n=l.8 ----n:3.0 --n=4.5 L
0
I
'
L
2
4
6
8
10
X Coordinate (m) Fig. 14. Sensitivity of the variably saturated model solution (t = 0,02day) to variations in the Van Genuchten parameter nv using the revised model, Eq, (17), for predicting the saturated hydraulic conductivity (c~v = 1.0m-l).
106
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
come up with an approximate empirical formulation relating the properties of capillarity and permeability for porous media. More experimental data are needed in this area to verify and develop more robust models. In fact, the above mentioned extreme path value radius, repv, or a similarly defined parameter, should be included in a comprehensive model of saturated permeability. Such level of detail, however, is forgone in the present analysis for reasons of tractability. The sensitivity analysis now discussed is directed toward the influences ofa~ and nv (both with their intrinsic impacts upon Ks). Each component of the sensitivity analysis is carried out by evaluating the saturated hydraulic conductivity through Eq. (17) for different values of av and nv, thus restricting the analysis to two degrees-of-freedom. Figs. 13 and 14 show the locations of the phreatic surface at time t = 0.02 day, for different values of the soil parameters av and nv, respectively. The values of the saturated hydraulic conductivity computed from Eq. (17) for the values Ofav = 0.5, 1.0, 2.0, 4.0m -1 are Ks = 0.32, 1.3, 5.1 and 20.3mday -l, respectively (nv = 2.8). Similarly for the values ofnv = 1.4, 1.8, 3.0, 4.5 the computed values of hydraulic conductivity are, Ks = 0.11, 0.27, 1.6 and 6 . 7 m d a y -1, respectively (av = 1.0m-l). In both cases, as the values of the soil parameters av and nv increase, the saturated hydraulic conductivity increases. In the case of av, the larger values of av correspond to larger, more permeable pores. For nv, this effect is caused by the uniformity of the pore-size-density function achieved at larger values of nv; this uniformity eliminates smaller, less permeable pores from the pore structure of media so parameterized. These increased saturated hydraulic conductivity values lead to faster drainage rates and hence more rapid falls in the phreatic surface. The effects of these changes in the saturated hydraulic conductivity values seem to dominate the capillaric effects discussed in the previous section (contrast Figs. 13 and 14 with Figs. 10 and 11, respectively). Consequently, it appears that K s, av and nv should not be treated as three independent variables when performing a sensitivity analysis for a variably saturated flow model. Their interrelation significantly affects the predicted behavior and precludes the use of these three degrees-of-freedom in sensitivity investigations. This coupling of capillary and permeability effects may also make it less appropriate to average soil parameters within identified soil classes, as the parameters are not independently distributed (as assumed during the traditionally used averaging schemes), but rather, are jointly distributed as a result of the fact that pore-size-density functions govern both capillary and permeability phenomena.
7. Conclusions
It is important to understand how variations in the soil parameters of variably saturated flow models affect simulations. Heuristic arguments support a claim that the capillaric properties of a porous medium should be closely related to the permeability of that medium. A sensitivity analysis based upon this principle reveals that simulation results are more sensitive to the permeability than to the capillarity of porous media. In fact, when the capillary properties are tied to the saturated hydraulic conductivity through a preliminary model, the sensitivities to the two capillary
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
107
parameters are, in general, reversed from that obtained with three degrees-offreedom, due to the manner in which variations in these parameters affect the permeability. The relationship between the descriptive soil parameters, which indirectly partially parameterize the pore structure, and the saturated hydraulic conductivity is still an open area in soil physics; more experimental investigations are needed to develop robust models. However, from the preliminary work performed herein, consideration of the number of degrees-of-freedom used when performing a sensitivity analysis is shown to warrant great care, if the predicted sensitivities are to have meaningful physical interpretations. Hopefully, the results of this research will provide a stepping stone to a more realistic understanding of the transient behavior involved with variably saturated flow problems.
Acknowledgment This work was supported, in part, by the US Environmental Protection Agency (contrast CR-818717-01-0) through the Environmental Research Laboratory in Athens, Georgia. However, it has not been through the agency's required peer and administrative review. Therefore, it does not necessarily reflect the views of the agency, and no official endorsement should be inferred.
References Ahmad, S., Kashyap, D. and Mathur, B.S., 1993. Mathematical modeling of saturated-unsaturated flow to drains. J. Irrig. Drain. Eng., 119: 18-33. Clement, T.P., Wise, W.R. and Molz, F.J., 1994. A physically based, two-dimensional finite-difference algorithm for variably saturated flow. J. Hydrol., 161 :, 71-90. Clement, T.P., Wise, W.R., Molz, F.J. and Wen, M., 1994. A comprehensive study of the modeling of steady-state seepage faces. Water Resour. Res., in press. Dullien, F.A.L., 1979. Porous Media Fluid Transport and Pore Structure. Academic Press, San Diego, 396 pp. Liggett, J.A., 1977. Location of free surface in porous media. J. Hydraul. Div., 103(HY4): 353-365. Mishra, S. and Parker, J.C., 1990. On the relation between saturated conductivity and capillary retention characteristics. Ground Water, 28: 775-777. M ualem, Y., 1976. A new model for predicting hydraulic conductivity of unsaturated porous media. Water Resour. Res., 12: 513-522. Rubin, J., 1968. Theoretical analysis of two-dimensional, transient flow of water in unsaturated and partly saturated soils. Soil Sci. Soc. Am. Proc,, 32: 607-615. Sisson, J.B., Ferguson, A.H. and Van Genuchten, M.Th., 1980. Simple method for predicting drainage from field plots. Soil Sci. Soc. Am. J., 44:1147 1152. Todsen, M., 1973. Numerical studies of two-dimensional saturated/unsaturated drainage models. J. Hydrol., 20:311-326. Van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J., 44: 892-898. Vauclin, M., Vachaud, G. and Khanji, J., 1975. Two dimensional numerical analysis of transient water transfer in saturated-unsaturated soils. In: G.C. Vansteenkiste (Editor), Modeling and Simulation of Water Resources Systems, North-Holland, Amsterdam, pp. 299-323.
108
W.R. Wise et al. / Journal of Hydrology 161 (1994) 91-108
Vauclin, M., Khanji, D. and Vachaud, G., 1979. Experimental and numerical study of a transient, twodimensional unsaturated-saturated water table recharge problem. Water Resour. Res., 15:1089 1101. Verma, R.D. and Brutsaert, W., 1970. Unconfined aquifer seepage by capillary flow theory. J. Hydraul. Div., 96(HY 6): 1331-1344. Watson, K.K. and Awadalla, S.A., 1985. Comparative study of Green Ampt analysis for a falling water table in a homogeneous sand profile. Water Resour. Res., 21:1157-1164. Wise, W.R., 1991. Discussion of 'On the relation between saturated conductivity and capillary retention characteristics,' by S. Mishra and J.C. Parker, September-October 1990 issue, vol. 28, no. 5, pp. 775777. Ground Water, 29: 272-273. Wise, W.R., 1992. A new insight on pore structure and permeability. Water Resour. Res., 28: 189-198. Young, E.G. and Aggelides, S., 1976. Drainage to a water table analysed by the Green-Ampt approach. J. Hydrol., 31:67 79.