- Email: [email protected]
A comparison of modeling approaches for steady-state unconfined flow
Journal
ELSEVIER
Journal of Hydrology 181 (1996) 189-209
A comparison of modeling approaches for steady-state unconfined flow T.P. Clementa, William R. Wiseb.*, Fred J. MolzC, Menghong Wend aBattelle, Pacific Northwest Laboratories, Battelle Boulevard, P.O. Box 999, Richland, WA 99352, USA bDepartment of Environmental Engineering Sciences, University of Florida, A.P. Black Hall, Gainesville, FL 32611-6450, USA ‘Environmental Systems Engineering Department, Clemson University, Rich Environmental Research Laboratories, Clemson Research Park, Clemson. SC 29634-0919, USA ‘Department of Civil Engineering, Auburn University, 238 Harbert Engineering Center, Auburn, AL 36849-5337, USA Received 14 October 1994; revision accepted 20 August 1995
Abstract The Dupuit-Forchheimer, the fully saturated flow, and the variably saturated flow models, are compared for problems involving steady-state, unconfined flow through porous media. The variably saturated flow model is the most comprehensive of the three and requires more parameters. The performances of the three models are compared for different soil properties, problem dimensions, and flow geometries. There are certain types of problems where the simpler models may yield satisfactory results. For soils with large pores and/or broad poresize-density functions, the variably saturated flow model solutions to steady-state problems approach those of the fully saturated flow model, owing to the manner in which the soil-water retention curve and relative permeability function, respectively, affect the variably saturated flow model solutions. For problems of significant size, the fully saturated flow model may be sufficient, as the effects of the vadose zone are relatively diminished. For problems with radial symmetry (e.g. steady flow to a well), the fully saturated flow model performs well because the variably saturated flow model is relatively insensitive to the parameters describing the soil properties, as the amount of vadose zone flow, compared with the total flow, is relatively insignificant in such problems.
Corresponding author.
??
0022-1694/96/!§15.000 1996 - Elsevier Science B.V. All rights reserved SSDZ 0022-l 694(95)02904-4
190
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
1. Introduction In the literature, three distinct models (listed in increasing level of sophistication and accuracy), the Dupuit-Forchheimer, the fully saturated flow, and the variably saturated flow models, are used to solve unconfined flow problems. Each model is based upon a unique set of assumptions and involves model-specific formulation of boundary conditions. In this paper, results from these three models are compared to shed insight into how two-dimensional flow near the exit boundary and flow through the vadose zone are treated by these models. The bases for comparison of these models are the predicted phreatic surface locations, seepage face heights, and system discharges. This analysis is limited to steady-state, two-dimensional problems. Some of the features presented herein have been addressed by Todsen (1973), Bear (1979), and Shamsai and Narasimhan (1991). However, the existing body of literature contains gaps, and in some cases misunderstandings, in the descriptions of how the three models perform for the types of problems posed within this work. At issue are what modeled behaviors are attributable to accounting for the seepage face versus those attributable to accounting for the vadose zone. This study is intended to help fill the gaps and develop a more comprehensive understanding of how these models perform relative to one another.
2. Mathematical models of steady-state unconfined flow This section briefly reviews the various mathematical formulations used to describe steady-state flow through unconfined aquifers. Only the governing equations are presented. The applicable boundary conditions for each model-problem combination discussed may be found in the cited literature. Of central importance to the discussion are the contrasting sets of assumptions made within the three models. These assumptions are numbered as follows: (1) the Dupuit assumptions and (2) no flow through the vadose zone. It should be noted that Assumption 1 may be invoked only after Assumption 2 has been asserted. The Dupuit-Forchheimer model is predicated upon both assumptions, the fully saturated flow model is based only upon the second, and both are relaxed for the variably saturated flow model. 2.1. The Dupuit-Forchheimer
model
Often, steady-state, two-dimensional hydraulics of unconfined flow is analyzed by invoking the Dupuit assumptions: equipotential surfaces are vertical and flow is essentially horizontal (Bear, 1979, pp. 74-78). Strack (1984) fine-tuned this description to suggest that vertical flow is allowed, but that resistance to such flow is ignored. The Dupuit assumptions preclude both seepage faces and flow through the vadose zone. The performance of the Dupuit-Forchheimer equations in modeling unconfined-flow problems has been investigated by Muskat (1937), Babbitt and Caldwell (1948), Boulton (1951), and Hantush (1964). Their research indicates that the Dupuit-Forchheimer discharge formula estimates the flow within l-2% of
T.P. Clement et al. 1 Journal of Hydrology 181 (1996) 189-209
191
experimental values. By contrast, the performance of the Dupuit-Forchheimer formula for predicting the location of phreatic surfaces has been historically found to be less satisfactory (Wyckoff et al., 1932; Muskat, 1937; Bear, 1979, p. 80). Analytic solutions to the unconfined-flow equations can be developed using the Dupuit assumptions. Using these assumptions, a two-dimensional flow problem in a homogeneous domain is approximated by a one-dimensional equation, which can be solved analytically subject to the appropriate boundary conditions. For rectangularflow systems, the well-known equations for the discharge and phreatic-surface are (Bear, 1979, pp. 78-80) Q’
=
K
s
(hi- hi) 2L
and h(x)
=
1 h;
_
(hi--hi)x] l’*
L
(2)
A
respectively, where Q’ [L* T-‘1 is the discharge per unit width of the aquifer, KS [L T-‘1 is the saturated hydraulic conductivity, ho [L] is the hydraulic head or the height of the phreatic surface at the position x = 0, and hL [L] is the head at x = L. Similarly, for radial-flow systems, the equations are (Bear, 1979, pp. 308-311)
and h(r) =
{
(h’,+ (h’R- hi)
[1~~;1,1}
I’*
(4)
where Q [L3 T-l] is the total discharge through the system, h, [L] is the elevation of the water table at the radius of the well r = rW, and hR [L] is the height of the water table at the radius of influence, r = rR [L]. The appropriate boundary conditions for both the rectangular and radial problems are well defined and illustrated in the cited reference. 2.2. Two-dimensional,
fully saturatedflow
model
In the fully saturated flow modeling approach, Laplace’s equation is solved for a homogeneous, isotropic, saturated flow domain neglecting the influence of the vadose zone. That is, the first assumption (the Dupuit assumptions) is relaxed, whereas the second is maintained. Laplace’s equation for a homogeneous rectangular-flow domain is written as
where x [L] and z [L] are the horizontal and vertical coordinates, respectively, and h
192
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
[L] is the hydraulic head. Although Laplace’s equation is linear, it is difficult to solve for phreatic-surface problems because of the nonlinear nature of phreatic-surface ‘boundaries’; that is, their positions are unknown before solution of Eq. (5). In this study, a numerical code based on the boundary-integral-equation method of Liggett and Liu (1983) is developed, verified against the analytic solutions of PolubarinovaKochina (1962), and used to develop the fully saturated flow estimates of phreatic surfaces, seepage faces, and discharges. The boundary-integral-equation method is also used to solve radial-flow problems, which are modeled using the threedimensional version of the Laplace equation:
#h $h #h (6) =+2+g=o 8Y The appropriate boundary conditions for both rectangular and radial problems have been presented by Liggett and Liu (1983, p. 66 and pp. 99-100, respectively). For each symmetry, these conditions include a seepage face. The phreatic surface is treated as a no-flow boundary; that is, no water flows through the vadose zone (which is outside of the domain of the fully saturated flow model). 2.3. Two-dimensional, variably saturatedflow model The transport of water within a two-dimensional, variably saturated domain can be written as (Richards, 1931) (7) where Q [L] is the pressure head, K(8) [L T-l] is the hydraulic conductivity function, 0 is the moisture content, and x [L] and z [L] are the Cartesian coordinates in the horizontal and vertical directions, respectively. Eq. (7) describes the variably saturated flow model of two-dimensional, steady-state, transport of water in a homogeneous, isotropic, unconfined soil matrix. For radial-flow systems, the steady-state form of the variably saturated flow equation is written in radial coordinates as
(8) where r [L] and z [L] are the radial and vertical coordinates, respectively. Owing to their nonlinear natures, Eqs. (7) and (8) cannot, in general, be solved analytically. The variably saturated lIow model is solved using a two-dimensional, finitedifference, variably saturated flow algorithm (Clement et al., 1994). For all of the simulations reported in this paper, the domains are discretized using grids of 100 nodes x 100 nodes. The boundary conditions required for the variably saturated formulation of both rectangular and radial problems have been addressed by Clement et al. (1994). Like those for the fully saturated flow model discussed above, these include a seepage face which is found iteratively (Clement et al., 1994). However, the domain for the variably
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
193
saturated flow model includes the vadose zone. As will be seen below, the phreatic surface is not a boundary to flow (not a streamline as suggested by Bear (1979, p. 76)), but merely the locus of all points within the porous medium where the water pressure is atmospheric. In fact, flow is observed to cross the phreatic surface in both directions (at different regions along the surface) for the variably saturated solutions to example problems presented later in this work. In this study, Van Genuchten’s (1980) closed-form equation for the soil-water retention curve and Mualem’s (1976) unsaturated hydraulic conductivity function are used to describe the soil properties for the variably saturated flow model. The relationship between water content and pressure head is given by (Van Genuchten, 1980)
Q= [l +(d,*,J”
(9
where CY~, n,, and m, = 1 - (1 /n,) are the Van Genuchten parameters, whose values depend upon the soil properties, and 8 is the effective saturation given by the relation
e- e,
0 =es- 8,
(10)
where 0, is the saturated water content and 0, is the residual water content of the soil. Of course, for nonnegative Q’,0 = 1, that is, the medium is saturated. The parameter (Y, [L-*1 is a measure of the first moment of the pore-size-density function (as (Y, increases, so does the first moment) and n, is an inverse measure of the second moment of the pore-size-density function (as n, increases, the pore-size-density function becomes narrower) (Wise, 1991; Wise et al., 1994). Based on Mualem’s (1976) model, the relation between water content and hydraulic conductivity is given by (Van Genuchten, 1980) K(Q) = KS{1 - [l - Qi’Q]mu}*Q”*
(11)
where KS [LT-l]is the saturated hydraulic conductivity. It should be noted that the relative permeability, K(Q)/K,, does not depend upon the value of (Y,,.The values for K(8)in Eqs. (7) and (8) are computed using Eqs. (9)-(1 l), simultaneously. The nature of the sensitivities of the soil-water retention curve and the relative permeability function to the parameters a, and n, have been presented by Van Genuchten and Nielsen (1985). Although this study considers the Van Genuchten (1980) and Mualem (1976) formulations for the water retention and transmission characteristics, respectively, results similar to those that follow would hold for other constitutive relations (e.g. the Brooks and Corey (1966) and Burdine (1953) models) where the underlying physical description (pore size information) is embedded in the pertinent parameters.
3. Steady-state unconfined flow in rectangular-flow systems A series of numerical experiments are performed to analyze the details of simulated
194
T.P. Clement et al. / Journal of Hydrology I81 (1996) 189-209
unconfined flow in various rectangular domains. Steady-state flow through a square embankment (see Fig. 1) is examined first. The 10 m x 10 m square embankment has no-flow boundary conditions on the base and at the top. The water levels at the left and right boundaries are 10 m and 2 m, respectively. The soil has a saturated hydraulic conductivity of 1.Om day-’ . As this is a steady-state flow problem, the head field is not a function of the saturated hydraulic conductivity, although the discharge is. The values for the Van Genuchten soil parameters are (Y, = 0.64 m-l and n, = 4.65, corresponding to a uniform, fine sand. This steady-state problem is solved using the Dupuit-Forchheimer, the fully saturated flow, and the variably saturated flow models. Of course, the saturated permeability, KS, is the only soil property required by the first two models. The phreatic surfaces computed using these three models are shown in Fig. 1, which illustrates that the height of the seepage face computed by the fully saturated flow model is smaller than the height predicted by the variably saturated flow model. As the fully saturated flow model neglects flow in the vadose zone, it underestimates the height of the seepage face and the discharge. The net discharge per unit width computed by the fully saturated flow model is 4.80 m2 day-‘, which is exactly equal to the value calculated by the Dupuit-Forchheimer model (Eq. (1)). (Bear (1988,‘~. 367) noted that this agreement is mandated.) The discharge computed by the more comprehensive, variably saturated flow model is 5.25 m2 day-‘. The additional flow computed by the variably saturated flow model is attributed to water transport through the vadose zone. As observed by Todsen (1973) towards the upgradient boundary, the fully saturated flow model overestimates the location of the phreatic surface (which is a
‘8. t $
*
r
-
-
-
Variably-Saturated Model I
0 0
Model
-Fully-Saturated
2
.
.
*
I 4
*
.
.
1
*
6
*
*
1. 8
”
10
X Coordinate (In) Fig. 1. Comparison of Dupuit-Forchheimer model, fully saturated model, and variably saturated model estimates of the phreatic surface and seepage face for the 10 m x 10 m square embankment problem (KS = 1 mday-‘, (Y, = 0.64 m-‘, and n, = 4.65).
195
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
no-flow boundary for this formulation). However, in the variably saturated solution, there is a substantial contribution to flow occurring above the phreatic surface, which can be seen from the velocity vectors in the vadose zone shown in Fig. 2. This figure indicates fairly horizontal flow near the upgradient boundary, which crosses the phreatic surface, staying predominantly in the capillary fringe. In the vicinity of the exit boundary, vadose-zone water re-enters the phreatic zone. In this region, the twodimensional flow converges towards the external open-water boundary. Fig. 2 shows that the inflow horizontal velocities are relatively uniform with respect to vertical position and that the outflow velocities vary considerably and approach their maximum value near the elevation of the downstream open-water level. Below this open-water level, outflow velocities decrease rapidly and approach a constant value. Above the open-water level, i.e. through the seepage face, horizontal flow components decrease rapidly with elevation and become zero at the top of the seepage face, where vertical flow components are dominant. Peterson (1957) has experimentally observed this type of outflow velocity distribution in a radial-flow system toward a gravity well. The fully saturated flow model also predicts similar types of inflow and outflow velocity distributions, given its less accurate predictions for both the seepage-face height and the discharge. As the Dupuit-Forchheimer model neglects resistance to vertical flow, the position of the phreatic surface and the base-pressure-head distribution are exactly the same for this model. By contrast, for the two-dimensional, variably saturated flow model they are different. From Fig. 1, it can be observed that the Dupuit-Forchheimer model gives a poor estimate for the location of the phreatic surface for this problem. Clement et al. (1994) illustrated that the base-pressure-head distributions predicted by
0
2
4
6
8
10
X Coordinate (m) Fig.2. Flowvectors Fig. 1.
corresponding
to the variably saturated model solution to the problem
illustratedin
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
196
both models are reasonably close. Consequently, the Dupuit-Forchheimer model gives a much better approximation for base-pressure-head distribution than it does for the location of the phreatic surface. Therefore, the Dupuit-Forchheimer model probably yields satisfactory estimates of the hydraulics just above the lower (impermeable) boundary of this unconfined system. 3.1. Rectangular flow:
influence of problem scale
The sensitivity of seepage face height to the scale of the problem was observed by Cooley (1983). He solved a 10 m x 10 m problem and a 100 m x 100 m problem using a two-dimensional, finite-element, variably saturated flow model, and observed that the seepage face did not scale to the problem dimensions. To shed further light onto this observation, several square-domain problems are solved using both the variably saturated and the fully saturated flow models. Simulation results are tabulated in Table 1. The seepage-face heights and discharges computed by the fully saturated flow model scale exactly with the problem dimensions, probably because the fully saturated flow equation, Eq. (5), is linear inside the flow domain (even though the bounding phreatic surface must be found iteratively). In the variably saturated flow model, the flow is governed by a highly nonlinear partial differential equation (Eq. (7)), and hence the problem does not scale simply with the problem dimensions. Table 1 documents that the effect of vadose-zone flow has more influence on smaller-scale problems. This trend is present because, in comparison with the total water transmitted through a system, more water flows through the vadose zone (where the nonlinearity arises) in a smaller-scale problem than in a larger-scale problem. Hence, the disparities between the fully saturated flow model and variably saturated flow model results are more pronounced in smaller-scale problems. This effect should be of particular interest to those trying to model laboratory-scale experiments. 3.2. RectangularJow:
influence of aspect ratio
In this section, the sensitivity of the seepage face determined
by the variably
Table 1 Discharge and seepage face sensitivity to problem scale for the various models Problem dimensions (m x m)
IX1
10 x 25 x 50 x 100 x
10 25 50 100
Tailwater level (m)
DupuitForchheimer model discharge (m’ day-‘)
Fully saturated model discharge (m* day-‘)
Variably saturated model discharge (m* day-‘)
Vadosezone flow contribution (% of flow)
Fully saturated model seepage face (m)
Variably saturated model seepage (m)
0.2 2.0 5.0 10.0 20.0
0.48 4.80 12.00 24.00 48.00
0.48 4.80 12.00 24.00 48.00
0.55 5.25 12.56 24.56 48.53
12.7 8.6 4.5 2.3 1.1
0.40 4.0 10.0 20.0 40.0
0.52 4.8 11.0 21.0 40.0
197
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
2
0
4
6
8
10
X Coordinate (m) Fig. 3. Variation of the phreatic surface predicted by the variably saturated model owing to changing downstream open-water levels for the square embankment problem (KJ = 1 mday-‘, Q, = 0.64 m-‘, and n, = 4.65).
t
t
I
0
.
.
..I
.
.
.
1
.
1
.
.
.
I
4
3
X &ordinate
.
.
.
.
.
5
(m)
Fig. 4. Variation of the phreatic surface predicted by the variably saturated model owing to changing downstream open-water levels for the ‘tall’ 5 m x 25 m, rectangular embankment problem (KS = 1 mday-‘, a, = 0.64 m-l, and n, = 4.65).
T.P. Clementet al. / Journalof Hydrology 181 (19%) 189-209
198
saturated flow model to the aspect ratio of the flow domain is examined. Three rectangular-flow regions with dimensions of 10 m x 10 m, 5 m x 25 m, and 25 m x 5 m are considered. The length-to-height ratio of the problems are 1,0.2, and 5, respectively. The soil properties used are the same as for the first example. For all three cases, four variably saturated flow simulations are performed, using tail-water levels of 20%, 40%, 60%, and 80% of the total height of the flow domains. Results of these simulations are shown in Figs. 3-5. For the square problem (Fig. 3) with an aspect ratio of unity, the absolute length of the seepage face decreases significantly with an increase in the tail-water level. The decrease in the seepage-face height appears to be due to a decrease in the discharge through the system at lower hydraulic gradients coupled with the concurrent increase in cross-sectional area for flow at higher tail-water levels. For the tall problem with an aspect ratio of 0.2, Fig. 4 shows that the position of the phreatic surface is relatively insensitive to tail-water levels. For the long problem, with an aspect ratio of five, no seepage faces are discernible for the various downgradient conditions (see Fig. 5) indicating that the effects of seepage faces are diminished for long, thin problems. Fig. 5 also demonstrates how the Dupuit-Forchheimer model gives fairly reasonable estimates for the positions of phreatic surfaces when seepage faces are negligible. 3.3. RectangularJlow: sensitivity to soil properties Steady-state water tables and seepage faces are insensitive to changes in the hydraulic conductivity. The resulting discharges are directly proportional to the
. Variably-Saturated Model: solid line
: Dupuit-Forchheimer Model: dashedline I,, ., I * * . . I . . . . I. 0
5
10
15
20
. . . 25
X Coordinate (m) Fig. 5. Variations of the phreatic surfaces predicted by the variably saturated model and the DupuitForchheimer model owing to changing downstream open-water levels for the ‘elongated’, 25 m x 5 m rectangular embankment problem (K, = 1 mday-‘, (Y,,= 0.64 m-‘, and n, = 4.65).
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
199
saturated hydraulic conductivity. These relations hold for all three models and are determinate from perusal of Eqs. (l)-(8), and (11). The hydraulic properties of an unsaturated soil are described by soil-water retention and relative-permeability functions. In this study, as described earlier, the Van Genuchten (1980) and Mualem (1976) models are used to describe the retention and transmission of water through porous media, respectively. Fig. 6 shows the sensitivity of the variably saturated flow model to variations in the parameter o,. In this figure, the value of (Y,is varied from 0.2 to 4.0 m-l, holding the parameter n, constant at a value of 2.5 and KS constant at 1 m day-‘. The Van Genuchten parameter a, describes a measure of the mean pore size of a porous medium. A simple sensitivity analysis for the effect of varying a, at a constant value of n, reveals that as the value of CY,increases, the pore sizes increase; bigger pores relate to a smaller capillary fringe (owing to less capillary action), and a relatively less extensive vadose zone available for the transmission of water. Fig. 6 shows that for (a relatively high value of) (Y, = 4.0 m-l, the position of the phreatic surface predicted by the variably saturated flow model almost relaxes to the position predicted by the fully saturated flow model. This result is expected, because as o, increases, the size of the capillary fringe above the phreatic surface decreases. For large values of ov, the height of the capillary fringe is negligible and the water content above the phreatic surface diminishes. Under these conditions, almost all of water movement (and retention) occurs below the phreatic surface. This situation is similar to that in the fully saturated flow model, which can be conceptually visualized as the variably saturated flow model with the following assumption: the hydraulic conductivity and/or water content of the porous medium is negligible above the
4
_
-
. ---
-
-
Alpha = 0.2
Alpha=OJ
Fully-Satumted Model
0 0
8
2
10
X C400rdinate6(m)
Fig. 6. Sensitivity analysis of the variably saturated flow model to the Van Genuchten (K, = 1 mday-’ and n,, = 2.5) for the square embankment problem.
parameter
cu,
200
T.P. Clement et al. / Journal OjHydrology 181 (1996) 189-209
Table 2 Variably saturated discharge sensitivity to a, (10 m x 10 m domain; n, = 2.5; Dupui-Forchheimer chargea is 4.80 m2 day-‘)
(m-9
Variably saturated model discharge (m2 day-‘)
0.2 0.5 1.0 4.0
5.38 5.20 5.04 4.85
cu,
dis-
a Exactly equal to the fully saturated model discharge.
phreatic surface. Table 2 gives the values of total discharge transmitted through the system for different values of a, (n, = 2.5). As a, increases, the discharge estimates reduce and approach the value predicted by the fully saturated flow model (and the Dupuit-Forchheimer model). Previous investigations by Muskat (1937), Boulton (1951) and Hantush (1964) suggest that the deviation from the Dupuit discharge is within l-2%. Recently, Shamsai and Narasimhan (1991) concluded that the deviation may range up to lo- 12%. The numerical experiments performed in this study clearly indicate that the deviation from Dupuit-Forchheimer flow is a function of soil properties and problem dimensions. For the 10 m x 10 m rectangular-flow domain studied in this section, the deviations are as high as 12% for o, = 0.2 m-l and as low as 1% for q, = 4.0 m-t
.
. . . II = 2.0
w
-- --- n = 1.25
- -
Fully-SammtedModel I
0 0
2
*
.
.
I
6
.
.
.
1..
6
.
I. a
.
. 10
X Coordinate (m) Fig. 7. Sensitivity analysis of the variably saturated flow model to the Van Genuchten parameter n, (K, = 1 mday-’ and a, = 0.5 m-‘) for the square embankment problem.
T.P. Clement et al. 1 Journal of Hydrology 181 (1996) 189-209
201
(n, = 2.5). Consequently, it is impossible and may be misleading to attempt to define a typical range of this disparity in discharge estimates. Shamsai and Narasimhan (1991) also suggested that the Dupuit-Forchheimer model will tend to underestimate the discharge because the flux through the seepage face is ignored. However, Bear (1988, p. 367) clearly demonstrated that the fully saturated flow model, which includes the seepage face, must exactly agree with the Dupuit-Forchheimer model, as far as discharge is concerned. This study shows that the underestimation of discharge by the Dupuit-Forchheimer model is instead a result of ignoring flow through the vadose zone. Accounting for the seepage face alone (by only relaxing Assumption l), as in the fully saturated flow model, does not improve the discharge estimate. Fig. 7 shows the sensitivity of the variably saturated flow model to variations in the soil parameter n,. The value of n, is varied from 1.25 to 7.0, holding constant the values of a, = 0.5 m-l and K, = 1 mday-i. Table 3 gives the values of total discharge transmitted through the system for different values of n,. As n, decreases, the discharge estimates reduce and approach the value predicted by the fully saturated flow model (and the Dupuit-Forchheimer model). The parameter n, isan inverse measure of the breadth of the pore-size-density function; as n, decreases the width of the pore-size-density function increases and vice versa. Mualem’s (1976) relative permeability model, Eq. (1 l), is a function of n,. When the value of n, is reduced, the relative abundance of smaller pores increases in a porous medium. The larger pores, which have less resistance to fluid transmission, tend to be drained first, constraining water to flow through the remaining (smaller) wetted pores. Hence, as the value of n, decreases, the relative permeability of the unsaturated zone decreases more rapidly with decreasing water content. As illustrated in Fig. 7 for the (relatively low) value of n, = 1.25, the variably saturated flow model predictions are almost identical to the fully saturated flow model predictions. It should be noted that for small values of n, (poorly sorted media), the low permeability of the unsaturated zone forces the variably saturated flow model to behave like a fully saturated flow model. By contrast, for high values of a, the low water content in the vadose zone forces the variably saturated flow model predictions to approach those of the fully saturated flow model. Although these differences are
Table 3 Variably saturated discharge sensitivity to n, (10 m x 10 m domain; ay, = 0.5 m-‘; Dupui-Forchheimer dischargea is 4.80 mz day-‘)
(m-l)
Variably saturated model discharge (m* day-‘)
7.0 4.0 2.0 1.25
5.36 5.29 5.12 4.91
n”
a Exactly equal to the fully saturated model discharge.
202
T.P. Clement et al. / Journal of Hydrology I81 (1996) 189-209
subtle, they are important in understanding the manner in which the variably saturated flow model responds to changes in the parameters relating to soil properties. An approximate upper bound for the effects of soil properties upon the discharge estimate, given a selected problem geometry, may be determined by choosing a relatively small value for Q, and a relatively large value for n,. For the 10 m x 10 m geometry, using a, = 0.2 m-l and n, = 7.0 results in a discharge of 5.51 m2 day-‘, which varies by almost 15% from the Dupuit-Forchheimer estimate (4.80 m2 day-‘). The interested reader is directed to the Appendix to explore some subtle modeling details involved in performing the types of sensitivity analyses addressed above.
4. Steady-state unconfined flow in radial-flow systems The dimensions of radial-flow problems considered in this study are similar to the rectangular-flow problems discussed above. As a first example, radial flow toward a pumping well in a finite, unconfined, cylindrical domain is analyzed. The radius of the flow domain (radius of influence) is 10 m, the height of the domain is 10 m and the radius of the well is 0.1 m. The saturated hydraulic conductivity of the soil is set to 1.0 m day-‘, and the values of the Van Genuchten parameters are o, = 1.O m-’ and n, = 2.5. The steady-state, radial-flow problem is solved using all three models. Fig. 8 shows the location of phreatic surfaces predicted by the different models when the pumped water level is 2 m. Reasonable agreement between the phreatic surfaces computed by
-
- -
2
4
variably-samrated Fhreadc swface Fully-saturated Model -Rxchhehner Model y-Saturated Base Ressure
6
8
10
Radial Coordinate (m) Fig. 8. Comparison of Dupuit-Forchheimer model, fully saturated model, and variably saturated model estimates of the phreatic surface and seepage face for a 10m x 10 m radial-flow problem (K, = 1 mday-‘, (Y, = 1.0 m-l, and n, = 2.5).
203
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
the fully saturated flow model and the variably saturated flow model suggests that the effect of transport through the vadose zone is relatively less important in radial than in rectangular problems. This conclusion is also verified by comparing the discharge estimates computed by the models. Using Eq. (3), the Dupuit-Forchheimer discharge (the same as the fully saturated flow discharge) for this problem is computed as Q = 65.5 m3day-‘, whereas the discharge computed using the variably saturated flow model is Q = 66.9 m3day-‘. The Dupuit-Forchheimer discharge formula for radial flow underestimates the discharge by only 2.1%. (Contrast this value with the 4.8% deviation for the similarly sized rectangular-flow problem with the same soil properties.) Owing to this small deviation, the effect of flow through the vadose zone does not have much influence on the extent of the seepage face. However, when compared with a similarly sized planar-flow problem (see Fig. l), the absolute length of the seepage face is more pronounced in the radial case (as observed by Shamsai and Narasimhan (1991)). This disparity is due to the convergent nature of radial flow, which requires more seepage area to accommodate the induced flow. Fig. 8 also shows that, similarly to its application to rectangular-flow problems, the DupuitForchheimer model approximates the base-pressure-head distribution reasonably well for this radial-flow problem. In the Dupuit-Forchheimer model for radial problems, the volume available for flow to a well is always restricted to the saturated zone below the well-water level. In reality, owing to the existence of the seepage face, a larger volume is available for flow. Hence, the Dupuit-Forchheimer model should not be used to estimate the velocity distribution near a well. The one-dimensional velocities computed by the Dupuit-Forchheimer model may have large deviations from the true two-
2
4
6
8
10
Radial Coordinate (m) Fig. 9. Variation of the phreatic surface predicted by the variably saturated model owing to changing wellwater levels for the 10 m x 10 m radial-flow problem (K, = 1 mday-‘, CY,= 1.0 m-l, and n, = 2.5).
204
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
dimensional velocity distribution. Moreover, as seen from Fig. 8, the cone of depression predicted by the variably saturated flow model is much smaller than that predicted by the Dupuit-Forchheimer model. This difference can play a significant role when designing remediation systems for contaminated aquifers using pump-and-treat systems. Seepage faces and the associated rises in phreatic surfaces are generally advantageous to pump-and-treat schemes because they expose larger aquifer volumes to flow, facilitating the removal of contaminants during cleanup efforts. At higher pumping rates, the water level inside a well further decreases, but owing to the seepage face, the volume exposed to the flow remains relatively unchanged over a significant variation of drawdowns (see Fig. 9). This fact is particularly important as remediation wells are often installed in the most contaminated regions at a site. Dupuit-Forchheimer-based phreatic surfaces should not be used when designing remediation schemes, as they are inherently inaccurate, and may force excessively low flow rates to be used for reasons of contacting contaminants which, in actuality, will be within the treatment zone at significantly higher flow rates. 4.1. Radialflow: influence of aspect ratio In this section, the sensitivity of phreatic surfaces predicted by the variably saturated flow model to variations in the aspect ratio of radial-flow problems is examined. Two flow domains with dimensions of 10 m x 10 m and 50 m x 5 m are considered. The length-to-height ratios of these problems are unity and ten, respectively. The soil properties used are the same as those used in the previous example. Similar to the rectangular analysis, variably saturated flow simulations are repeated four times using tail-water levels of 20%, 40%, 60%, and 80% of the
Variably&turated Model: solid line
0
10
40 R8~UaYYCoordinZ.e
50
(m)
Fig. 10. Variation of the phreatic surface predicted by the variably saturated model owing to changing wellwater levels for the 50 m x 5 m radial-flow problem (KS = 1 mday-‘, CZ,= 1.0 m-l, and n, = 2.5).
‘T.P. Clement et al. 1 Journal of Hydrology 181 (19%)
189-209
20s
total height of the flow domain. Results of these simulations are summarized in Figs. 9 and 10. For the problem with an aspect ratio of unity, changes in the well-water level have a minor influence on the phreatic surface. By comparing Figs. 3 and 9, changes in the height of the downstream boundary are shown to have a more significant influence on planar problems than on radial problems. This fact suggests that seepage faces are more dominant and persistent in radial-flow problems (as observed by Shamsai and Narasimhan (1991)). Fig. 10 illustrates both the Dupuit-Forchheimer model and the variably saturated flow model predictions for the elongated problem with an aspect ratio of ten, for the indicated drawdowns. The influence of the seepage face is still significant and, hence, the Dupuit-Forchheimer model is a poor approximation, especially near the well. Away from the tail-water boundary in long, field-scale problems, the Dupuit-Forchheimer model overestimates the height of the phreatic surface in both radial- and rectangular-flow domains (see Figs. 5 and 10). The Dupuit-Forchheimer approach is somewhat more successful at describing elongated rectangular-flow problems compared with elongated radial-flow problems. 4.2. Radial flow: sensitivity to soil parameters Clement et al. (1994) showed that the phreatic surfaces in radial problems are relatively insensitive to these soil properties (but the same trends toward the fully saturated flow model exist as CX,, increases and n, decreases). This relative insensitivity is one reason why, in spite of using different types of soil characteristics, several workers, including Shamsai and Narasimhan (1991), Cooley (1983) and Taylor and Luthin (1969), were able to model Hall’s (1955) experimental (radial flow toward a gravity well in a sand-filled, wedge-shaped box) phreatic-surface data reasonably well, despite the fact that Hall did not report any capillary properties of his soil. This apparent insensitivity also helps to elucidate why the fully saturated flow model approximates the position of the phreatic surface in a radial-flow problem better than in planar-flow problems. 4.3. Radialjlow: sensitivity to well radius
Clement et al. (1994) also discussed the sensitivities of the seepage-face and phreatic-surface positions to changes in the well radius. As the well radius increases, the height of the seepage face decreases. Even though the absolute changes are not pronounced, Clement et al. (1994) observed that variably saturated, radial-flow problems seem to be relatively more sensitive to the well radius than to the soil properties.
5. Conclusions The numerical experiments indicate that the major contribution to a seepage face arises primarily from the two-dimensional nature of the flow near the exit boundary. Ignoring the transport through the vadose zone and applying a two-dimensional,
206
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
fully saturated flow model generally results in an underestimation of the seepage face and the net discharge through the system. Hence, modeling unconfined flow using a fully saturated flow model may not be appropriate, especially in smaller-scale problems (e.g. simulating a laboratory-scale apparatus with rectangular symmetry). The height of the seepage face and the vadose-zone flow contribution predicted by the variably saturated flow model are functions of the parameters describing the soil properties, the problem dimensions, and the symmetry of the flow domain. Sensitivity analyses suggest that solutions to rectangular-flow problems are more sensitive to soil parameters than are those for radial-flow problems. The natures of these sensitivities are that the fully saturated flow modeling approach becomes more valid as the porous medium of interest becomes either coarser or more poorly sorted, owing to the effects of the retention and transmission characteristics of the vadose zone, respectively. Simulations show that in limiting conditions the location of the phreatic surface and the flow predicted by the variably saturated flow model relax to those predicted by the fully saturated flow model. Even though the absolute changes are not pronounced, solutions to the variably saturated, radial-flow model are generally more sensitive to changes in the well radius than to variations in the parameters describing the soil properties. Simulation results show that the effect of seepage faces may be ignored in some (notably elongated) rectangular field-scale problems. For similarly sized problems, deviations from the Dupuit-Forchheimer discharge are more significant in rectangular problems than in radial problems. However, owing to the convergent nature of radial flow, seepage faces are more pronounced in radial problems. Thus, much caution should be used when applying a Dupuit-Forchheimer analysis to evaluate or design remediation systems because flow velocities into a well will be overestimated (if the water level in the well is used to describe the vertical extent of the flow at the well), as will the depth of the cone of depression. Seepage faces in radial-flow systems are generally advantageous to remediation efforts because they expose larger aquifer volumes to flow, thereby facilitating the removal of contaminants.
Acknowledgments
This work was supported, in part, by the US Environmental Protection Agency (Contract CR-8 187 17-O1-O) 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.
Appendix: a brief caveat
Previous findings by the first three authors (Wise et al., 1994) indicate that, to represent physical reality,, it is more appropriate to perform sensitivity analyses for
207
.T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
Table A 1 Variably saturated discharge sensitivity to a, (10 m x 10 m domain; n, = 2.5)
0.2 0.5 1.0 4.0
Saturated hydraulic conductivity (m day-‘)
Variably saturated model discharge (m’ day-‘)
Dupuit-Forchheimer dischargeb (m’ day-‘)
0.011 0.070 0.278 4.451
0.059 0.364 1.401 21.587
0.053 0.336 1.334 21.365
a Determined by Wise et al. (1994, Eqs. (16) and (17)), assuming the difference between the saturated and residual water contents is 0.287. b Exactly equal to the fully saturated model discharge.
unconfined flow problems using two degrees of freedom (varying (Y, and n, and inferring K,) rather than three degrees of freedom (simultaneously varying Q,, nu, and K,), as done above. The authors have used three degrees of freedom (even though one, K,, was not exercised) purposely because this paper is aimed at highlighting differences in modeling approaches; this objective is achieved more easily with the course followed above. Moreover, as the head distributions (and consequently, the phreatic surfaces and seepage faces) for steady-state unconfined flow are independent of the saturated hydraulic conductivity, only the discharge estimates are affected. Use of this approach does not affect the comparisons between the three modeling approaches discussed in this paper, as the parameterized saturated hydraulic conductivities (Wise et al., 1994) would be used for all three models in a more exhaustive comparison. It does, however, affect the absolute discharge estimates of these models, as the discharge is directly proportional to the saturated hydraulic conductivity for all of the models. Were the relation between the saturated hydraulic conductivity, KS, and the Van Genuchten (1980) parameters, (Y,and nv, presented by Wise et al. (1994, their Eqs. (16) and (17)) used, the relation between o,, and discharge (predicted by all models) identified above would typically be reversed (a higher o, would result in a higher discharge); the relation between n, and discharge (predicted Table A2 Variably saturated discharge sensitivity to n, (10 m x 10 m domain; o, = 0.5) a,
Saturated hydraulic conductivitya (m day-‘)
Variably saturated model discharge (m* day-‘)
Dupuit-Forchheimer dischargeb (m’ day-‘)
7.0 4.0 2.0 1.25
2.555 0.360 0.032 0.006
13.695 1.904 0.164 0.029
12.264 1.728 0.154 0.029
a Determined by Wise et al. (1994, Eqs. (16) and (17)), assuming the difference between the saturated and residual water contents is 0.287. b Exactly equal to the fully saturated model discharge.
208
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
by all models) identified above would be exaggerated (a higher n, would result in an even higher discharge). These effects are illustrated in Tables Al and A2, which present the same sensitivity analyses of Tables 2 and 3, respectively, given the empirical coupling between KS, and the Van Genuchten (1980) parameters, q, and n,, presented by Wise et al. (1994). The difference between the saturated and residual water contents is set to 0.287, as this is the value which brings the first sample problem of this paper (KS = 1.0 mday-‘, a, = 0.64 m-l, and n, = 4.65) into concordance with the desired relation. Although the discharges presented in Tables Al and A2 are different from those presented in Tables 2 and 3, respectively, the ratios of variably saturated flow model discharge estimates to Dupuit-Forchheimer (and fully saturated flow model) discharge estimates are exactly preserved.
References Babbitt, H.E. and Caldwell, D.H., 1948. The free surface around and interface between gravity wells. Ill. Univ. Eng. Exp. Stn Bull. 374. Bear, J., 1979. Hydraulics of Groundwater. McGraw-Hill, New York. Bear, J., 1988. Dynamics of Fluids in Porous Media. Dover, New York, 764 pp. Boulton, N.S., 1951. The flow pattern near a gravity well in a uniform water bearing medium. J. Inst. Civ. Eng., 36(10): 543-550. Brooks, R.H. and Corey, A.T., 1966. Properties of porous media affecting fluid flow. J. Brig. Drainage Div. ASCE, 92(IR2): 61-88. Burdine, N.T., 1953. Relative permeability calculations from pore sire distribution data. Pet. Trans. AIME, 198(T.P. 3519): 71-78. Clement, T.P., Wise, W.R. and Molz, F.J., 1994. A physically based, two-dimensional, finite-difference algorithm for modeling variably saturated flow. J. Hydrol., 161: 71-90. Cooley, R.L., 1983. Some new procedures for numerical solution of variably saturated flow problems. Water Resour. Res., 19: 1271-1285. Hall, H.P., 1955. An investigation of steady state flow toward a gravity well. Houille Blanche, 10; 8-35. Hantush, M.S., 1964. Hydraulics of wells. In: Advances in Hydroscience, Vol. 1. Academic Press, San Diego, CA, pp. 281-432. Liggett, J.A. and Liu, P.L-F., 1983. The Boundary Integral Equation Method for Porous Media Flow. George Allen and Unwin, London. Mualem, Y., 1976. A new model for predicting hydraulic conductivity of unsaturated porous media. Water Resour. Res., 12: 513-522. Muskat, M., 1937. The Flow of Homogeneous Fluids through Porous Media. McGraw-Hill, New York. Peterson, D.F., 1957. Hydraulics of wells. Trans. ASCE, 122: 502-517. Polubarinova-Kochina, P. Ya., 1962:Theory of Groundwater Movement (translated from Russian by R.J.M. de Wiest). Princeton University Press, Princeton, NJ. Richards, L.A., 1931. Capillary conduction of liquids through porous mediums. Physics, 1: 318-333. Shamsai, A. and Narasimhan, T.N., 1991. A numerical investigation of free surface-seepage face relationship under steady state flow conditions. Water Resour. Res., 27(3): 409-421. Strack, O.D.L., 1984. Three-dimensional streamlines in Dupuit-Forchheimer models. Water Resour. Res., 20(7): 8 12-822. Taylor, G.S. and Luthin, J.N., 1969. Computer methods for transient analysis of water table aquifers. Water Resour. Res., 5(l): 144152. Todsen, M., 1973. Numerical studies of two-dimensional saturated/unsaturated drainage models. J. Hydrol., 20: 31 l-326. Van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Sot. Am. J., 44: 892-898.
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
209
Van Genuchten, M.Th. and Nielsen, D.R., 1985. On describing and predicting the hydraulic properties of unsaturated soils. Ann. Geophys., 3(5): 615-628. 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, v. 28, no. 5, pp. 115-711. Ground Water, 29(2): 212-213. Wise, W.R., Clement, T.P. and Mob, F.J., 1994. Variably saturated modeling of transient drainage: sensitivity to soil properties. J. Hydrol., 161: 91-108. Wyckoff, R.D., Botset, H.G. and Muskat, M., 1932. Flow of liquids through porous media under the action of gravity. Physics, 2: 90-l 13.
ELSEVIER
Journal of Hydrology 181 (1996) 189-209
A comparison of modeling approaches for steady-state unconfined flow T.P. Clementa, William R. Wiseb.*, Fred J. MolzC, Menghong Wend aBattelle, Pacific Northwest Laboratories, Battelle Boulevard, P.O. Box 999, Richland, WA 99352, USA bDepartment of Environmental Engineering Sciences, University of Florida, A.P. Black Hall, Gainesville, FL 32611-6450, USA ‘Environmental Systems Engineering Department, Clemson University, Rich Environmental Research Laboratories, Clemson Research Park, Clemson. SC 29634-0919, USA ‘Department of Civil Engineering, Auburn University, 238 Harbert Engineering Center, Auburn, AL 36849-5337, USA Received 14 October 1994; revision accepted 20 August 1995
Abstract The Dupuit-Forchheimer, the fully saturated flow, and the variably saturated flow models, are compared for problems involving steady-state, unconfined flow through porous media. The variably saturated flow model is the most comprehensive of the three and requires more parameters. The performances of the three models are compared for different soil properties, problem dimensions, and flow geometries. There are certain types of problems where the simpler models may yield satisfactory results. For soils with large pores and/or broad poresize-density functions, the variably saturated flow model solutions to steady-state problems approach those of the fully saturated flow model, owing to the manner in which the soil-water retention curve and relative permeability function, respectively, affect the variably saturated flow model solutions. For problems of significant size, the fully saturated flow model may be sufficient, as the effects of the vadose zone are relatively diminished. For problems with radial symmetry (e.g. steady flow to a well), the fully saturated flow model performs well because the variably saturated flow model is relatively insensitive to the parameters describing the soil properties, as the amount of vadose zone flow, compared with the total flow, is relatively insignificant in such problems.
Corresponding author.
??
0022-1694/96/!§15.000 1996 - Elsevier Science B.V. All rights reserved SSDZ 0022-l 694(95)02904-4
190
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
1. Introduction In the literature, three distinct models (listed in increasing level of sophistication and accuracy), the Dupuit-Forchheimer, the fully saturated flow, and the variably saturated flow models, are used to solve unconfined flow problems. Each model is based upon a unique set of assumptions and involves model-specific formulation of boundary conditions. In this paper, results from these three models are compared to shed insight into how two-dimensional flow near the exit boundary and flow through the vadose zone are treated by these models. The bases for comparison of these models are the predicted phreatic surface locations, seepage face heights, and system discharges. This analysis is limited to steady-state, two-dimensional problems. Some of the features presented herein have been addressed by Todsen (1973), Bear (1979), and Shamsai and Narasimhan (1991). However, the existing body of literature contains gaps, and in some cases misunderstandings, in the descriptions of how the three models perform for the types of problems posed within this work. At issue are what modeled behaviors are attributable to accounting for the seepage face versus those attributable to accounting for the vadose zone. This study is intended to help fill the gaps and develop a more comprehensive understanding of how these models perform relative to one another.
2. Mathematical models of steady-state unconfined flow This section briefly reviews the various mathematical formulations used to describe steady-state flow through unconfined aquifers. Only the governing equations are presented. The applicable boundary conditions for each model-problem combination discussed may be found in the cited literature. Of central importance to the discussion are the contrasting sets of assumptions made within the three models. These assumptions are numbered as follows: (1) the Dupuit assumptions and (2) no flow through the vadose zone. It should be noted that Assumption 1 may be invoked only after Assumption 2 has been asserted. The Dupuit-Forchheimer model is predicated upon both assumptions, the fully saturated flow model is based only upon the second, and both are relaxed for the variably saturated flow model. 2.1. The Dupuit-Forchheimer
model
Often, steady-state, two-dimensional hydraulics of unconfined flow is analyzed by invoking the Dupuit assumptions: equipotential surfaces are vertical and flow is essentially horizontal (Bear, 1979, pp. 74-78). Strack (1984) fine-tuned this description to suggest that vertical flow is allowed, but that resistance to such flow is ignored. The Dupuit assumptions preclude both seepage faces and flow through the vadose zone. The performance of the Dupuit-Forchheimer equations in modeling unconfined-flow problems has been investigated by Muskat (1937), Babbitt and Caldwell (1948), Boulton (1951), and Hantush (1964). Their research indicates that the Dupuit-Forchheimer discharge formula estimates the flow within l-2% of
T.P. Clement et al. 1 Journal of Hydrology 181 (1996) 189-209
191
experimental values. By contrast, the performance of the Dupuit-Forchheimer formula for predicting the location of phreatic surfaces has been historically found to be less satisfactory (Wyckoff et al., 1932; Muskat, 1937; Bear, 1979, p. 80). Analytic solutions to the unconfined-flow equations can be developed using the Dupuit assumptions. Using these assumptions, a two-dimensional flow problem in a homogeneous domain is approximated by a one-dimensional equation, which can be solved analytically subject to the appropriate boundary conditions. For rectangularflow systems, the well-known equations for the discharge and phreatic-surface are (Bear, 1979, pp. 78-80) Q’
=
K
s
(hi- hi) 2L
and h(x)
=
1 h;
_
(hi--hi)x] l’*
L
(2)
A
respectively, where Q’ [L* T-‘1 is the discharge per unit width of the aquifer, KS [L T-‘1 is the saturated hydraulic conductivity, ho [L] is the hydraulic head or the height of the phreatic surface at the position x = 0, and hL [L] is the head at x = L. Similarly, for radial-flow systems, the equations are (Bear, 1979, pp. 308-311)
and h(r) =
{
(h’,+ (h’R- hi)
[1~~;1,1}
I’*
(4)
where Q [L3 T-l] is the total discharge through the system, h, [L] is the elevation of the water table at the radius of the well r = rW, and hR [L] is the height of the water table at the radius of influence, r = rR [L]. The appropriate boundary conditions for both the rectangular and radial problems are well defined and illustrated in the cited reference. 2.2. Two-dimensional,
fully saturatedflow
model
In the fully saturated flow modeling approach, Laplace’s equation is solved for a homogeneous, isotropic, saturated flow domain neglecting the influence of the vadose zone. That is, the first assumption (the Dupuit assumptions) is relaxed, whereas the second is maintained. Laplace’s equation for a homogeneous rectangular-flow domain is written as
where x [L] and z [L] are the horizontal and vertical coordinates, respectively, and h
192
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
[L] is the hydraulic head. Although Laplace’s equation is linear, it is difficult to solve for phreatic-surface problems because of the nonlinear nature of phreatic-surface ‘boundaries’; that is, their positions are unknown before solution of Eq. (5). In this study, a numerical code based on the boundary-integral-equation method of Liggett and Liu (1983) is developed, verified against the analytic solutions of PolubarinovaKochina (1962), and used to develop the fully saturated flow estimates of phreatic surfaces, seepage faces, and discharges. The boundary-integral-equation method is also used to solve radial-flow problems, which are modeled using the threedimensional version of the Laplace equation:
#h $h #h (6) =+2+g=o 8Y The appropriate boundary conditions for both rectangular and radial problems have been presented by Liggett and Liu (1983, p. 66 and pp. 99-100, respectively). For each symmetry, these conditions include a seepage face. The phreatic surface is treated as a no-flow boundary; that is, no water flows through the vadose zone (which is outside of the domain of the fully saturated flow model). 2.3. Two-dimensional, variably saturatedflow model The transport of water within a two-dimensional, variably saturated domain can be written as (Richards, 1931) (7) where Q [L] is the pressure head, K(8) [L T-l] is the hydraulic conductivity function, 0 is the moisture content, and x [L] and z [L] are the Cartesian coordinates in the horizontal and vertical directions, respectively. Eq. (7) describes the variably saturated flow model of two-dimensional, steady-state, transport of water in a homogeneous, isotropic, unconfined soil matrix. For radial-flow systems, the steady-state form of the variably saturated flow equation is written in radial coordinates as
(8) where r [L] and z [L] are the radial and vertical coordinates, respectively. Owing to their nonlinear natures, Eqs. (7) and (8) cannot, in general, be solved analytically. The variably saturated lIow model is solved using a two-dimensional, finitedifference, variably saturated flow algorithm (Clement et al., 1994). For all of the simulations reported in this paper, the domains are discretized using grids of 100 nodes x 100 nodes. The boundary conditions required for the variably saturated formulation of both rectangular and radial problems have been addressed by Clement et al. (1994). Like those for the fully saturated flow model discussed above, these include a seepage face which is found iteratively (Clement et al., 1994). However, the domain for the variably
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
193
saturated flow model includes the vadose zone. As will be seen below, the phreatic surface is not a boundary to flow (not a streamline as suggested by Bear (1979, p. 76)), but merely the locus of all points within the porous medium where the water pressure is atmospheric. In fact, flow is observed to cross the phreatic surface in both directions (at different regions along the surface) for the variably saturated solutions to example problems presented later in this work. In this study, Van Genuchten’s (1980) closed-form equation for the soil-water retention curve and Mualem’s (1976) unsaturated hydraulic conductivity function are used to describe the soil properties for the variably saturated flow model. The relationship between water content and pressure head is given by (Van Genuchten, 1980)
Q= [l +(d,*,J”
(9
where CY~, n,, and m, = 1 - (1 /n,) are the Van Genuchten parameters, whose values depend upon the soil properties, and 8 is the effective saturation given by the relation
e- e,
0 =es- 8,
(10)
where 0, is the saturated water content and 0, is the residual water content of the soil. Of course, for nonnegative Q’,0 = 1, that is, the medium is saturated. The parameter (Y, [L-*1 is a measure of the first moment of the pore-size-density function (as (Y, increases, so does the first moment) and n, is an inverse measure of the second moment of the pore-size-density function (as n, increases, the pore-size-density function becomes narrower) (Wise, 1991; Wise et al., 1994). Based on Mualem’s (1976) model, the relation between water content and hydraulic conductivity is given by (Van Genuchten, 1980) K(Q) = KS{1 - [l - Qi’Q]mu}*Q”*
(11)
where KS [LT-l]is the saturated hydraulic conductivity. It should be noted that the relative permeability, K(Q)/K,, does not depend upon the value of (Y,,.The values for K(8)in Eqs. (7) and (8) are computed using Eqs. (9)-(1 l), simultaneously. The nature of the sensitivities of the soil-water retention curve and the relative permeability function to the parameters a, and n, have been presented by Van Genuchten and Nielsen (1985). Although this study considers the Van Genuchten (1980) and Mualem (1976) formulations for the water retention and transmission characteristics, respectively, results similar to those that follow would hold for other constitutive relations (e.g. the Brooks and Corey (1966) and Burdine (1953) models) where the underlying physical description (pore size information) is embedded in the pertinent parameters.
3. Steady-state unconfined flow in rectangular-flow systems A series of numerical experiments are performed to analyze the details of simulated
194
T.P. Clement et al. / Journal of Hydrology I81 (1996) 189-209
unconfined flow in various rectangular domains. Steady-state flow through a square embankment (see Fig. 1) is examined first. The 10 m x 10 m square embankment has no-flow boundary conditions on the base and at the top. The water levels at the left and right boundaries are 10 m and 2 m, respectively. The soil has a saturated hydraulic conductivity of 1.Om day-’ . As this is a steady-state flow problem, the head field is not a function of the saturated hydraulic conductivity, although the discharge is. The values for the Van Genuchten soil parameters are (Y, = 0.64 m-l and n, = 4.65, corresponding to a uniform, fine sand. This steady-state problem is solved using the Dupuit-Forchheimer, the fully saturated flow, and the variably saturated flow models. Of course, the saturated permeability, KS, is the only soil property required by the first two models. The phreatic surfaces computed using these three models are shown in Fig. 1, which illustrates that the height of the seepage face computed by the fully saturated flow model is smaller than the height predicted by the variably saturated flow model. As the fully saturated flow model neglects flow in the vadose zone, it underestimates the height of the seepage face and the discharge. The net discharge per unit width computed by the fully saturated flow model is 4.80 m2 day-‘, which is exactly equal to the value calculated by the Dupuit-Forchheimer model (Eq. (1)). (Bear (1988,‘~. 367) noted that this agreement is mandated.) The discharge computed by the more comprehensive, variably saturated flow model is 5.25 m2 day-‘. The additional flow computed by the variably saturated flow model is attributed to water transport through the vadose zone. As observed by Todsen (1973) towards the upgradient boundary, the fully saturated flow model overestimates the location of the phreatic surface (which is a
‘8. t $
*
r
-
-
-
Variably-Saturated Model I
0 0
Model
-Fully-Saturated
2
.
.
*
I 4
*
.
.
1
*
6
*
*
1. 8
”
10
X Coordinate (In) Fig. 1. Comparison of Dupuit-Forchheimer model, fully saturated model, and variably saturated model estimates of the phreatic surface and seepage face for the 10 m x 10 m square embankment problem (KS = 1 mday-‘, (Y, = 0.64 m-‘, and n, = 4.65).
195
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
no-flow boundary for this formulation). However, in the variably saturated solution, there is a substantial contribution to flow occurring above the phreatic surface, which can be seen from the velocity vectors in the vadose zone shown in Fig. 2. This figure indicates fairly horizontal flow near the upgradient boundary, which crosses the phreatic surface, staying predominantly in the capillary fringe. In the vicinity of the exit boundary, vadose-zone water re-enters the phreatic zone. In this region, the twodimensional flow converges towards the external open-water boundary. Fig. 2 shows that the inflow horizontal velocities are relatively uniform with respect to vertical position and that the outflow velocities vary considerably and approach their maximum value near the elevation of the downstream open-water level. Below this open-water level, outflow velocities decrease rapidly and approach a constant value. Above the open-water level, i.e. through the seepage face, horizontal flow components decrease rapidly with elevation and become zero at the top of the seepage face, where vertical flow components are dominant. Peterson (1957) has experimentally observed this type of outflow velocity distribution in a radial-flow system toward a gravity well. The fully saturated flow model also predicts similar types of inflow and outflow velocity distributions, given its less accurate predictions for both the seepage-face height and the discharge. As the Dupuit-Forchheimer model neglects resistance to vertical flow, the position of the phreatic surface and the base-pressure-head distribution are exactly the same for this model. By contrast, for the two-dimensional, variably saturated flow model they are different. From Fig. 1, it can be observed that the Dupuit-Forchheimer model gives a poor estimate for the location of the phreatic surface for this problem. Clement et al. (1994) illustrated that the base-pressure-head distributions predicted by
0
2
4
6
8
10
X Coordinate (m) Fig.2. Flowvectors Fig. 1.
corresponding
to the variably saturated model solution to the problem
illustratedin
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
196
both models are reasonably close. Consequently, the Dupuit-Forchheimer model gives a much better approximation for base-pressure-head distribution than it does for the location of the phreatic surface. Therefore, the Dupuit-Forchheimer model probably yields satisfactory estimates of the hydraulics just above the lower (impermeable) boundary of this unconfined system. 3.1. Rectangular flow:
influence of problem scale
The sensitivity of seepage face height to the scale of the problem was observed by Cooley (1983). He solved a 10 m x 10 m problem and a 100 m x 100 m problem using a two-dimensional, finite-element, variably saturated flow model, and observed that the seepage face did not scale to the problem dimensions. To shed further light onto this observation, several square-domain problems are solved using both the variably saturated and the fully saturated flow models. Simulation results are tabulated in Table 1. The seepage-face heights and discharges computed by the fully saturated flow model scale exactly with the problem dimensions, probably because the fully saturated flow equation, Eq. (5), is linear inside the flow domain (even though the bounding phreatic surface must be found iteratively). In the variably saturated flow model, the flow is governed by a highly nonlinear partial differential equation (Eq. (7)), and hence the problem does not scale simply with the problem dimensions. Table 1 documents that the effect of vadose-zone flow has more influence on smaller-scale problems. This trend is present because, in comparison with the total water transmitted through a system, more water flows through the vadose zone (where the nonlinearity arises) in a smaller-scale problem than in a larger-scale problem. Hence, the disparities between the fully saturated flow model and variably saturated flow model results are more pronounced in smaller-scale problems. This effect should be of particular interest to those trying to model laboratory-scale experiments. 3.2. RectangularJow:
influence of aspect ratio
In this section, the sensitivity of the seepage face determined
by the variably
Table 1 Discharge and seepage face sensitivity to problem scale for the various models Problem dimensions (m x m)
IX1
10 x 25 x 50 x 100 x
10 25 50 100
Tailwater level (m)
DupuitForchheimer model discharge (m’ day-‘)
Fully saturated model discharge (m* day-‘)
Variably saturated model discharge (m* day-‘)
Vadosezone flow contribution (% of flow)
Fully saturated model seepage face (m)
Variably saturated model seepage (m)
0.2 2.0 5.0 10.0 20.0
0.48 4.80 12.00 24.00 48.00
0.48 4.80 12.00 24.00 48.00
0.55 5.25 12.56 24.56 48.53
12.7 8.6 4.5 2.3 1.1
0.40 4.0 10.0 20.0 40.0
0.52 4.8 11.0 21.0 40.0
197
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
2
0
4
6
8
10
X Coordinate (m) Fig. 3. Variation of the phreatic surface predicted by the variably saturated model owing to changing downstream open-water levels for the square embankment problem (KJ = 1 mday-‘, Q, = 0.64 m-‘, and n, = 4.65).
t
t
I
0
.
.
..I
.
.
.
1
.
1
.
.
.
I
4
3
X &ordinate
.
.
.
.
.
5
(m)
Fig. 4. Variation of the phreatic surface predicted by the variably saturated model owing to changing downstream open-water levels for the ‘tall’ 5 m x 25 m, rectangular embankment problem (KS = 1 mday-‘, a, = 0.64 m-l, and n, = 4.65).
T.P. Clementet al. / Journalof Hydrology 181 (19%) 189-209
198
saturated flow model to the aspect ratio of the flow domain is examined. Three rectangular-flow regions with dimensions of 10 m x 10 m, 5 m x 25 m, and 25 m x 5 m are considered. The length-to-height ratio of the problems are 1,0.2, and 5, respectively. The soil properties used are the same as for the first example. For all three cases, four variably saturated flow simulations are performed, using tail-water levels of 20%, 40%, 60%, and 80% of the total height of the flow domains. Results of these simulations are shown in Figs. 3-5. For the square problem (Fig. 3) with an aspect ratio of unity, the absolute length of the seepage face decreases significantly with an increase in the tail-water level. The decrease in the seepage-face height appears to be due to a decrease in the discharge through the system at lower hydraulic gradients coupled with the concurrent increase in cross-sectional area for flow at higher tail-water levels. For the tall problem with an aspect ratio of 0.2, Fig. 4 shows that the position of the phreatic surface is relatively insensitive to tail-water levels. For the long problem, with an aspect ratio of five, no seepage faces are discernible for the various downgradient conditions (see Fig. 5) indicating that the effects of seepage faces are diminished for long, thin problems. Fig. 5 also demonstrates how the Dupuit-Forchheimer model gives fairly reasonable estimates for the positions of phreatic surfaces when seepage faces are negligible. 3.3. RectangularJlow: sensitivity to soil properties Steady-state water tables and seepage faces are insensitive to changes in the hydraulic conductivity. The resulting discharges are directly proportional to the
. Variably-Saturated Model: solid line
: Dupuit-Forchheimer Model: dashedline I,, ., I * * . . I . . . . I. 0
5
10
15
20
. . . 25
X Coordinate (m) Fig. 5. Variations of the phreatic surfaces predicted by the variably saturated model and the DupuitForchheimer model owing to changing downstream open-water levels for the ‘elongated’, 25 m x 5 m rectangular embankment problem (K, = 1 mday-‘, (Y,,= 0.64 m-‘, and n, = 4.65).
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
199
saturated hydraulic conductivity. These relations hold for all three models and are determinate from perusal of Eqs. (l)-(8), and (11). The hydraulic properties of an unsaturated soil are described by soil-water retention and relative-permeability functions. In this study, as described earlier, the Van Genuchten (1980) and Mualem (1976) models are used to describe the retention and transmission of water through porous media, respectively. Fig. 6 shows the sensitivity of the variably saturated flow model to variations in the parameter o,. In this figure, the value of (Y,is varied from 0.2 to 4.0 m-l, holding the parameter n, constant at a value of 2.5 and KS constant at 1 m day-‘. The Van Genuchten parameter a, describes a measure of the mean pore size of a porous medium. A simple sensitivity analysis for the effect of varying a, at a constant value of n, reveals that as the value of CY,increases, the pore sizes increase; bigger pores relate to a smaller capillary fringe (owing to less capillary action), and a relatively less extensive vadose zone available for the transmission of water. Fig. 6 shows that for (a relatively high value of) (Y, = 4.0 m-l, the position of the phreatic surface predicted by the variably saturated flow model almost relaxes to the position predicted by the fully saturated flow model. This result is expected, because as o, increases, the size of the capillary fringe above the phreatic surface decreases. For large values of ov, the height of the capillary fringe is negligible and the water content above the phreatic surface diminishes. Under these conditions, almost all of water movement (and retention) occurs below the phreatic surface. This situation is similar to that in the fully saturated flow model, which can be conceptually visualized as the variably saturated flow model with the following assumption: the hydraulic conductivity and/or water content of the porous medium is negligible above the
4
_
-
. ---
-
-
Alpha = 0.2
Alpha=OJ
Fully-Satumted Model
0 0
8
2
10
X C400rdinate6(m)
Fig. 6. Sensitivity analysis of the variably saturated flow model to the Van Genuchten (K, = 1 mday-’ and n,, = 2.5) for the square embankment problem.
parameter
cu,
200
T.P. Clement et al. / Journal OjHydrology 181 (1996) 189-209
Table 2 Variably saturated discharge sensitivity to a, (10 m x 10 m domain; n, = 2.5; Dupui-Forchheimer chargea is 4.80 m2 day-‘)
(m-9
Variably saturated model discharge (m2 day-‘)
0.2 0.5 1.0 4.0
5.38 5.20 5.04 4.85
cu,
dis-
a Exactly equal to the fully saturated model discharge.
phreatic surface. Table 2 gives the values of total discharge transmitted through the system for different values of a, (n, = 2.5). As a, increases, the discharge estimates reduce and approach the value predicted by the fully saturated flow model (and the Dupuit-Forchheimer model). Previous investigations by Muskat (1937), Boulton (1951) and Hantush (1964) suggest that the deviation from the Dupuit discharge is within l-2%. Recently, Shamsai and Narasimhan (1991) concluded that the deviation may range up to lo- 12%. The numerical experiments performed in this study clearly indicate that the deviation from Dupuit-Forchheimer flow is a function of soil properties and problem dimensions. For the 10 m x 10 m rectangular-flow domain studied in this section, the deviations are as high as 12% for o, = 0.2 m-l and as low as 1% for q, = 4.0 m-t
.
. . . II = 2.0
w
-- --- n = 1.25
- -
Fully-SammtedModel I
0 0
2
*
.
.
I
6
.
.
.
1..
6
.
I. a
.
. 10
X Coordinate (m) Fig. 7. Sensitivity analysis of the variably saturated flow model to the Van Genuchten parameter n, (K, = 1 mday-’ and a, = 0.5 m-‘) for the square embankment problem.
T.P. Clement et al. 1 Journal of Hydrology 181 (1996) 189-209
201
(n, = 2.5). Consequently, it is impossible and may be misleading to attempt to define a typical range of this disparity in discharge estimates. Shamsai and Narasimhan (1991) also suggested that the Dupuit-Forchheimer model will tend to underestimate the discharge because the flux through the seepage face is ignored. However, Bear (1988, p. 367) clearly demonstrated that the fully saturated flow model, which includes the seepage face, must exactly agree with the Dupuit-Forchheimer model, as far as discharge is concerned. This study shows that the underestimation of discharge by the Dupuit-Forchheimer model is instead a result of ignoring flow through the vadose zone. Accounting for the seepage face alone (by only relaxing Assumption l), as in the fully saturated flow model, does not improve the discharge estimate. Fig. 7 shows the sensitivity of the variably saturated flow model to variations in the soil parameter n,. The value of n, is varied from 1.25 to 7.0, holding constant the values of a, = 0.5 m-l and K, = 1 mday-i. Table 3 gives the values of total discharge transmitted through the system for different values of n,. As n, decreases, the discharge estimates reduce and approach the value predicted by the fully saturated flow model (and the Dupuit-Forchheimer model). The parameter n, isan inverse measure of the breadth of the pore-size-density function; as n, decreases the width of the pore-size-density function increases and vice versa. Mualem’s (1976) relative permeability model, Eq. (1 l), is a function of n,. When the value of n, is reduced, the relative abundance of smaller pores increases in a porous medium. The larger pores, which have less resistance to fluid transmission, tend to be drained first, constraining water to flow through the remaining (smaller) wetted pores. Hence, as the value of n, decreases, the relative permeability of the unsaturated zone decreases more rapidly with decreasing water content. As illustrated in Fig. 7 for the (relatively low) value of n, = 1.25, the variably saturated flow model predictions are almost identical to the fully saturated flow model predictions. It should be noted that for small values of n, (poorly sorted media), the low permeability of the unsaturated zone forces the variably saturated flow model to behave like a fully saturated flow model. By contrast, for high values of a, the low water content in the vadose zone forces the variably saturated flow model predictions to approach those of the fully saturated flow model. Although these differences are
Table 3 Variably saturated discharge sensitivity to n, (10 m x 10 m domain; ay, = 0.5 m-‘; Dupui-Forchheimer dischargea is 4.80 mz day-‘)
(m-l)
Variably saturated model discharge (m* day-‘)
7.0 4.0 2.0 1.25
5.36 5.29 5.12 4.91
n”
a Exactly equal to the fully saturated model discharge.
202
T.P. Clement et al. / Journal of Hydrology I81 (1996) 189-209
subtle, they are important in understanding the manner in which the variably saturated flow model responds to changes in the parameters relating to soil properties. An approximate upper bound for the effects of soil properties upon the discharge estimate, given a selected problem geometry, may be determined by choosing a relatively small value for Q, and a relatively large value for n,. For the 10 m x 10 m geometry, using a, = 0.2 m-l and n, = 7.0 results in a discharge of 5.51 m2 day-‘, which varies by almost 15% from the Dupuit-Forchheimer estimate (4.80 m2 day-‘). The interested reader is directed to the Appendix to explore some subtle modeling details involved in performing the types of sensitivity analyses addressed above.
4. Steady-state unconfined flow in radial-flow systems The dimensions of radial-flow problems considered in this study are similar to the rectangular-flow problems discussed above. As a first example, radial flow toward a pumping well in a finite, unconfined, cylindrical domain is analyzed. The radius of the flow domain (radius of influence) is 10 m, the height of the domain is 10 m and the radius of the well is 0.1 m. The saturated hydraulic conductivity of the soil is set to 1.0 m day-‘, and the values of the Van Genuchten parameters are o, = 1.O m-’ and n, = 2.5. The steady-state, radial-flow problem is solved using all three models. Fig. 8 shows the location of phreatic surfaces predicted by the different models when the pumped water level is 2 m. Reasonable agreement between the phreatic surfaces computed by
-
- -
2
4
variably-samrated Fhreadc swface Fully-saturated Model -Rxchhehner Model y-Saturated Base Ressure
6
8
10
Radial Coordinate (m) Fig. 8. Comparison of Dupuit-Forchheimer model, fully saturated model, and variably saturated model estimates of the phreatic surface and seepage face for a 10m x 10 m radial-flow problem (K, = 1 mday-‘, (Y, = 1.0 m-l, and n, = 2.5).
203
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
the fully saturated flow model and the variably saturated flow model suggests that the effect of transport through the vadose zone is relatively less important in radial than in rectangular problems. This conclusion is also verified by comparing the discharge estimates computed by the models. Using Eq. (3), the Dupuit-Forchheimer discharge (the same as the fully saturated flow discharge) for this problem is computed as Q = 65.5 m3day-‘, whereas the discharge computed using the variably saturated flow model is Q = 66.9 m3day-‘. The Dupuit-Forchheimer discharge formula for radial flow underestimates the discharge by only 2.1%. (Contrast this value with the 4.8% deviation for the similarly sized rectangular-flow problem with the same soil properties.) Owing to this small deviation, the effect of flow through the vadose zone does not have much influence on the extent of the seepage face. However, when compared with a similarly sized planar-flow problem (see Fig. l), the absolute length of the seepage face is more pronounced in the radial case (as observed by Shamsai and Narasimhan (1991)). This disparity is due to the convergent nature of radial flow, which requires more seepage area to accommodate the induced flow. Fig. 8 also shows that, similarly to its application to rectangular-flow problems, the DupuitForchheimer model approximates the base-pressure-head distribution reasonably well for this radial-flow problem. In the Dupuit-Forchheimer model for radial problems, the volume available for flow to a well is always restricted to the saturated zone below the well-water level. In reality, owing to the existence of the seepage face, a larger volume is available for flow. Hence, the Dupuit-Forchheimer model should not be used to estimate the velocity distribution near a well. The one-dimensional velocities computed by the Dupuit-Forchheimer model may have large deviations from the true two-
2
4
6
8
10
Radial Coordinate (m) Fig. 9. Variation of the phreatic surface predicted by the variably saturated model owing to changing wellwater levels for the 10 m x 10 m radial-flow problem (K, = 1 mday-‘, CY,= 1.0 m-l, and n, = 2.5).
204
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
dimensional velocity distribution. Moreover, as seen from Fig. 8, the cone of depression predicted by the variably saturated flow model is much smaller than that predicted by the Dupuit-Forchheimer model. This difference can play a significant role when designing remediation systems for contaminated aquifers using pump-and-treat systems. Seepage faces and the associated rises in phreatic surfaces are generally advantageous to pump-and-treat schemes because they expose larger aquifer volumes to flow, facilitating the removal of contaminants during cleanup efforts. At higher pumping rates, the water level inside a well further decreases, but owing to the seepage face, the volume exposed to the flow remains relatively unchanged over a significant variation of drawdowns (see Fig. 9). This fact is particularly important as remediation wells are often installed in the most contaminated regions at a site. Dupuit-Forchheimer-based phreatic surfaces should not be used when designing remediation schemes, as they are inherently inaccurate, and may force excessively low flow rates to be used for reasons of contacting contaminants which, in actuality, will be within the treatment zone at significantly higher flow rates. 4.1. Radialflow: influence of aspect ratio In this section, the sensitivity of phreatic surfaces predicted by the variably saturated flow model to variations in the aspect ratio of radial-flow problems is examined. Two flow domains with dimensions of 10 m x 10 m and 50 m x 5 m are considered. The length-to-height ratios of these problems are unity and ten, respectively. The soil properties used are the same as those used in the previous example. Similar to the rectangular analysis, variably saturated flow simulations are repeated four times using tail-water levels of 20%, 40%, 60%, and 80% of the
Variably&turated Model: solid line
0
10
40 R8~UaYYCoordinZ.e
50
(m)
Fig. 10. Variation of the phreatic surface predicted by the variably saturated model owing to changing wellwater levels for the 50 m x 5 m radial-flow problem (KS = 1 mday-‘, CZ,= 1.0 m-l, and n, = 2.5).
‘T.P. Clement et al. 1 Journal of Hydrology 181 (19%)
189-209
20s
total height of the flow domain. Results of these simulations are summarized in Figs. 9 and 10. For the problem with an aspect ratio of unity, changes in the well-water level have a minor influence on the phreatic surface. By comparing Figs. 3 and 9, changes in the height of the downstream boundary are shown to have a more significant influence on planar problems than on radial problems. This fact suggests that seepage faces are more dominant and persistent in radial-flow problems (as observed by Shamsai and Narasimhan (1991)). Fig. 10 illustrates both the Dupuit-Forchheimer model and the variably saturated flow model predictions for the elongated problem with an aspect ratio of ten, for the indicated drawdowns. The influence of the seepage face is still significant and, hence, the Dupuit-Forchheimer model is a poor approximation, especially near the well. Away from the tail-water boundary in long, field-scale problems, the Dupuit-Forchheimer model overestimates the height of the phreatic surface in both radial- and rectangular-flow domains (see Figs. 5 and 10). The Dupuit-Forchheimer approach is somewhat more successful at describing elongated rectangular-flow problems compared with elongated radial-flow problems. 4.2. Radial flow: sensitivity to soil parameters Clement et al. (1994) showed that the phreatic surfaces in radial problems are relatively insensitive to these soil properties (but the same trends toward the fully saturated flow model exist as CX,, increases and n, decreases). This relative insensitivity is one reason why, in spite of using different types of soil characteristics, several workers, including Shamsai and Narasimhan (1991), Cooley (1983) and Taylor and Luthin (1969), were able to model Hall’s (1955) experimental (radial flow toward a gravity well in a sand-filled, wedge-shaped box) phreatic-surface data reasonably well, despite the fact that Hall did not report any capillary properties of his soil. This apparent insensitivity also helps to elucidate why the fully saturated flow model approximates the position of the phreatic surface in a radial-flow problem better than in planar-flow problems. 4.3. Radialjlow: sensitivity to well radius
Clement et al. (1994) also discussed the sensitivities of the seepage-face and phreatic-surface positions to changes in the well radius. As the well radius increases, the height of the seepage face decreases. Even though the absolute changes are not pronounced, Clement et al. (1994) observed that variably saturated, radial-flow problems seem to be relatively more sensitive to the well radius than to the soil properties.
5. Conclusions The numerical experiments indicate that the major contribution to a seepage face arises primarily from the two-dimensional nature of the flow near the exit boundary. Ignoring the transport through the vadose zone and applying a two-dimensional,
206
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
fully saturated flow model generally results in an underestimation of the seepage face and the net discharge through the system. Hence, modeling unconfined flow using a fully saturated flow model may not be appropriate, especially in smaller-scale problems (e.g. simulating a laboratory-scale apparatus with rectangular symmetry). The height of the seepage face and the vadose-zone flow contribution predicted by the variably saturated flow model are functions of the parameters describing the soil properties, the problem dimensions, and the symmetry of the flow domain. Sensitivity analyses suggest that solutions to rectangular-flow problems are more sensitive to soil parameters than are those for radial-flow problems. The natures of these sensitivities are that the fully saturated flow modeling approach becomes more valid as the porous medium of interest becomes either coarser or more poorly sorted, owing to the effects of the retention and transmission characteristics of the vadose zone, respectively. Simulations show that in limiting conditions the location of the phreatic surface and the flow predicted by the variably saturated flow model relax to those predicted by the fully saturated flow model. Even though the absolute changes are not pronounced, solutions to the variably saturated, radial-flow model are generally more sensitive to changes in the well radius than to variations in the parameters describing the soil properties. Simulation results show that the effect of seepage faces may be ignored in some (notably elongated) rectangular field-scale problems. For similarly sized problems, deviations from the Dupuit-Forchheimer discharge are more significant in rectangular problems than in radial problems. However, owing to the convergent nature of radial flow, seepage faces are more pronounced in radial problems. Thus, much caution should be used when applying a Dupuit-Forchheimer analysis to evaluate or design remediation systems because flow velocities into a well will be overestimated (if the water level in the well is used to describe the vertical extent of the flow at the well), as will the depth of the cone of depression. Seepage faces in radial-flow systems are generally advantageous to remediation efforts because they expose larger aquifer volumes to flow, thereby facilitating the removal of contaminants.
Acknowledgments
This work was supported, in part, by the US Environmental Protection Agency (Contract CR-8 187 17-O1-O) 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.
Appendix: a brief caveat
Previous findings by the first three authors (Wise et al., 1994) indicate that, to represent physical reality,, it is more appropriate to perform sensitivity analyses for
207
.T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
Table A 1 Variably saturated discharge sensitivity to a, (10 m x 10 m domain; n, = 2.5)
0.2 0.5 1.0 4.0
Saturated hydraulic conductivity (m day-‘)
Variably saturated model discharge (m’ day-‘)
Dupuit-Forchheimer dischargeb (m’ day-‘)
0.011 0.070 0.278 4.451
0.059 0.364 1.401 21.587
0.053 0.336 1.334 21.365
a Determined by Wise et al. (1994, Eqs. (16) and (17)), assuming the difference between the saturated and residual water contents is 0.287. b Exactly equal to the fully saturated model discharge.
unconfined flow problems using two degrees of freedom (varying (Y, and n, and inferring K,) rather than three degrees of freedom (simultaneously varying Q,, nu, and K,), as done above. The authors have used three degrees of freedom (even though one, K,, was not exercised) purposely because this paper is aimed at highlighting differences in modeling approaches; this objective is achieved more easily with the course followed above. Moreover, as the head distributions (and consequently, the phreatic surfaces and seepage faces) for steady-state unconfined flow are independent of the saturated hydraulic conductivity, only the discharge estimates are affected. Use of this approach does not affect the comparisons between the three modeling approaches discussed in this paper, as the parameterized saturated hydraulic conductivities (Wise et al., 1994) would be used for all three models in a more exhaustive comparison. It does, however, affect the absolute discharge estimates of these models, as the discharge is directly proportional to the saturated hydraulic conductivity for all of the models. Were the relation between the saturated hydraulic conductivity, KS, and the Van Genuchten (1980) parameters, (Y,and nv, presented by Wise et al. (1994, their Eqs. (16) and (17)) used, the relation between o,, and discharge (predicted by all models) identified above would typically be reversed (a higher o, would result in a higher discharge); the relation between n, and discharge (predicted Table A2 Variably saturated discharge sensitivity to n, (10 m x 10 m domain; o, = 0.5) a,
Saturated hydraulic conductivitya (m day-‘)
Variably saturated model discharge (m* day-‘)
Dupuit-Forchheimer dischargeb (m’ day-‘)
7.0 4.0 2.0 1.25
2.555 0.360 0.032 0.006
13.695 1.904 0.164 0.029
12.264 1.728 0.154 0.029
a Determined by Wise et al. (1994, Eqs. (16) and (17)), assuming the difference between the saturated and residual water contents is 0.287. b Exactly equal to the fully saturated model discharge.
208
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
by all models) identified above would be exaggerated (a higher n, would result in an even higher discharge). These effects are illustrated in Tables Al and A2, which present the same sensitivity analyses of Tables 2 and 3, respectively, given the empirical coupling between KS, and the Van Genuchten (1980) parameters, q, and n,, presented by Wise et al. (1994). The difference between the saturated and residual water contents is set to 0.287, as this is the value which brings the first sample problem of this paper (KS = 1.0 mday-‘, a, = 0.64 m-l, and n, = 4.65) into concordance with the desired relation. Although the discharges presented in Tables Al and A2 are different from those presented in Tables 2 and 3, respectively, the ratios of variably saturated flow model discharge estimates to Dupuit-Forchheimer (and fully saturated flow model) discharge estimates are exactly preserved.
References Babbitt, H.E. and Caldwell, D.H., 1948. The free surface around and interface between gravity wells. Ill. Univ. Eng. Exp. Stn Bull. 374. Bear, J., 1979. Hydraulics of Groundwater. McGraw-Hill, New York. Bear, J., 1988. Dynamics of Fluids in Porous Media. Dover, New York, 764 pp. Boulton, N.S., 1951. The flow pattern near a gravity well in a uniform water bearing medium. J. Inst. Civ. Eng., 36(10): 543-550. Brooks, R.H. and Corey, A.T., 1966. Properties of porous media affecting fluid flow. J. Brig. Drainage Div. ASCE, 92(IR2): 61-88. Burdine, N.T., 1953. Relative permeability calculations from pore sire distribution data. Pet. Trans. AIME, 198(T.P. 3519): 71-78. Clement, T.P., Wise, W.R. and Molz, F.J., 1994. A physically based, two-dimensional, finite-difference algorithm for modeling variably saturated flow. J. Hydrol., 161: 71-90. Cooley, R.L., 1983. Some new procedures for numerical solution of variably saturated flow problems. Water Resour. Res., 19: 1271-1285. Hall, H.P., 1955. An investigation of steady state flow toward a gravity well. Houille Blanche, 10; 8-35. Hantush, M.S., 1964. Hydraulics of wells. In: Advances in Hydroscience, Vol. 1. Academic Press, San Diego, CA, pp. 281-432. Liggett, J.A. and Liu, P.L-F., 1983. The Boundary Integral Equation Method for Porous Media Flow. George Allen and Unwin, London. Mualem, Y., 1976. A new model for predicting hydraulic conductivity of unsaturated porous media. Water Resour. Res., 12: 513-522. Muskat, M., 1937. The Flow of Homogeneous Fluids through Porous Media. McGraw-Hill, New York. Peterson, D.F., 1957. Hydraulics of wells. Trans. ASCE, 122: 502-517. Polubarinova-Kochina, P. Ya., 1962:Theory of Groundwater Movement (translated from Russian by R.J.M. de Wiest). Princeton University Press, Princeton, NJ. Richards, L.A., 1931. Capillary conduction of liquids through porous mediums. Physics, 1: 318-333. Shamsai, A. and Narasimhan, T.N., 1991. A numerical investigation of free surface-seepage face relationship under steady state flow conditions. Water Resour. Res., 27(3): 409-421. Strack, O.D.L., 1984. Three-dimensional streamlines in Dupuit-Forchheimer models. Water Resour. Res., 20(7): 8 12-822. Taylor, G.S. and Luthin, J.N., 1969. Computer methods for transient analysis of water table aquifers. Water Resour. Res., 5(l): 144152. Todsen, M., 1973. Numerical studies of two-dimensional saturated/unsaturated drainage models. J. Hydrol., 20: 31 l-326. Van Genuchten, M.Th., 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Sot. Am. J., 44: 892-898.
T.P. Clement et al. / Journal of Hydrology 181 (1996) 189-209
209
Van Genuchten, M.Th. and Nielsen, D.R., 1985. On describing and predicting the hydraulic properties of unsaturated soils. Ann. Geophys., 3(5): 615-628. 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, v. 28, no. 5, pp. 115-711. Ground Water, 29(2): 212-213. Wise, W.R., Clement, T.P. and Mob, F.J., 1994. Variably saturated modeling of transient drainage: sensitivity to soil properties. J. Hydrol., 161: 91-108. Wyckoff, R.D., Botset, H.G. and Muskat, M., 1932. Flow of liquids through porous media under the action of gravity. Physics, 2: 90-l 13.