


In the finite difference modeling of unsaturated infiltration by Richards' equation, standard interblock conductivity means like the arithmetic, geometric and harmonic have been shown to produce flows different from true Darcian by up to orders of magnitude (Warrick, 1991; Baker, 1995a, 1998). This may have significant impact on models used for contaminant transport of nuclear and other hazardous waste. Previously, it has been impractical to calculate Darcian means in a running model. This paper introduces the analytic solution for Darcian means for the case of an exponential relative conductivity relation, and one possible approximation for general use, a piecewise BrooksCorey Darcian (PBCD) mean. It compares the PBCD mean to the arithmetic mean in a modern massconservative numerical model (Celia, et al., 1990) for vertical space steps from 0.8 to 10^{5} cm, using soil parameters from Haverkamp, et al., (1977). The results demonstrate that even in a model with mass balance good to one part in 10^{8}, nonDarcian flow errors are alive and well, costing up to 50% of the total mass infiltrated with the arithmetic mean and large model space steps. The new PBCD mean reduces that figure to 0.5% for this example soil.
Although Warrick (1991) demonstrated that standard interblock conductivity means, such as the arithmetic, geometric and upstream, produce nonDarcian flow, sometimes to orders of magnitude in error, many still use standard means in models of vertical unsaturated flow. The reason is as simple as standard means themselves  they do not take much time or power to calculate on a computer.
Baker (1993, 1995b) demonstrated that combining methods to conserve mass and reduce time step discretization errors can produce increases in computational accuracy and efficiency in an unsaturated flow model amounting to orders of magnitude. The stability of such models can depend on the type of interblock conductivity mean chosen. Baker (1995a, 1998) graphically demonstrated, with a new approach to Warrick’s concept of the Darcian mean, that Darcian means even in simple, homogeneous soils possess a much greater complexity than previously recognized. They vary radically with porous medium properties and model vertical space step size, while standard means are valid only in special cases.
Baker (1998) further developed the concept of Darcian means, discovering an analytic solution for the general exponential unsaturated conductivity relation [1], and suggesting a piecewise general approximation using the BrooksCorey conductivity relation. Until otherwise demonstrated, it currently seems that an analytic solution exists only for the exponential conductivity relation, making reasonable approximations imperative if Darcian means are to be used in running models. This paper will briefly review Darcian means and illustrate such an approximation, demonstrating the relative errors of this approximation and the arithmetic mean.
[1] , where
kr = unsaturated relative conductivity (dimensionless), l = fitting parameter (1/m), y = matric suction (m), and y_{d} = displacement pressure (m).
Darcian means come from the modeling assumption that Darcy’s Law for vertical unsaturated flow [2] holds constant between grid cells (Figure 1) in space [3] during one time step. It is usually applied as the difference equation [4], where one chooses a value for the vertical interblock conductivity mean, Kv. In a fully implicit model, Kv is held constant at the level of the next discrete time, usually improved to assure mass balance by iteration.
[2] , where rq = relative flow, Q/Ks, (dimensionless), Q = flow (m/s), Ks = saturated conductivity (m/s), H = total head (m), and x = vertical displacement (m).
[3]
[4] , where Kv = the Darcian relative conductivity mean for vertical flow (dimensionless).


If one integrates [2] between x1 and x2 (Warrick, 1991; Baker, 1995a, 1998) and solves for Kv in [4], one gets [5]. Kx is the mean Darcian conductivity for gravity or advective flow and Kh is the mean Darcian conductivity for horizontal and diffusive flow. When finite differences are used to solve Richards’ equation for unsaturated flow [6], some may express the righthandside using separately estimated Kx and Kh components. In many cases, Kh is analytically integrable for the chosen form of kr(Y). Note that [5] shows Kv ® Kh, in the limit as Dx ® 0, and that Kv ® Kx as Dx ® ¥ . But since Kx itself is changing from a prototypical form to kr2 as 0 ® Dx ® ¥ , the second property is much less useful.
But there is considerable difficulty in knowing the distribution of kr(x) between grid centers. It is nearly impossible (Baker, 1998) to balance any approximation of Kx with Kh in [5] to produce a mean that is in general Darcian. The errors can reach orders of magnitude near hydrostatic conditions, Dx d Dy, and even produce Kv < 0, which is very nonphysical. It is better to use the total Darcian mean Kv in [4], or an approximation.
[5] , where and
[6]
We solve numerically for the true Darcian mean, Kv, by substituting [2] into [3] to get [7]. We solve [7] by iteration until [3] is satisfied between each pair of fine grid cells to an arbitrarily small amount, expressed as the extreme spread in rq between grid cells, divided by the mean rq. Since this can involve about 1000 or more fine grid cells for each value of Kv between the grid centers in Figure 1, it is not practical to use in a running model. When the numerical solution converges to a value of rq to one part in, say, 10^{5} or 10^{6}, we have Kv from [4].
[7] ,
where kr(x1) = kr1, kr(x2) = kr2, Y (x1) = Y(kr1), and Y (x1) = Y(kr2), given the relationship kr(Y ).
However, the exponential form of kr (Y ), [1], has an analytic solution [8] to [7], from which approximations for other forms can be developed. The reader can verify this solution by direct substitution of [1] and [8] into [7]. (The lead author has found the symbolic computation package, Maple V, to be very useful.) One can use equation [5] to develop Kh [9], Kx [10] and Kv [11] for this form. Note that as Dx or u d 0, Kx d (kr1+kr2)/2 and Kv d Kh. When kr1 = kr2, Kx = Kv = Kh = kr2, demonstrating that the means properly pass through the kr1=kr2 diagonal in (kr1,kr2)space. As Dx or u d +¥ , Kv d Kx d kr2.
[8] ,
where a = kr1, , u = l× Dx , and Dx = x2x1
[9]
[10]
[11]
Many investigators construct an estimate of Kv from a weighted sum of kr1 and kr2. But we see from this and previous work that Kv really transforms from Kh to kr2, so if a weighted sum is used, it should really be of the form [12]. The analytic solution to the weight, wv, is quite complicated and will not be presented here. But the value of wv on the kr1=kr2 diagonal, wvo [13], is much simpler, and more useful for approximations to other forms of kr(y). Because of the exponential nature of [1], the scaling parameter, u, can be defined as [14]. Baker (1998) found that this equation tends to hold with [13] on the kr1=kr2 diagonal for other forms of kr(y), such as the BrooksCorey and vanGenuchtenstyle forms.
[12]
[13]
[14]
It currently appears that all of the more complicated relations of kr(y), such as the BrooksCorey [15] or vanGenuchten style, have no explicit analytic solutions for Kv. Therefore, an approximation is necessary. It is tempting to substitute an approximation of Kh into [11] and create a version of Kv(kr1,kr2,Kh). Unfortunately, it suffers from the same discontinuity about hydrostatic conditions as [5] in some cases. It is better to put any approximation or calculation of Kh into [12], substituting wvo [13] for wv. Since wvo is bounded in the range of [0,1], this at least guarantees that Kv will be bounded by the range [Kh,kr2].
[15] ,
where yd = "displacement pressure" (cm), and h is the "pore size distribution index".
Suppose that we have a arbitrary kr(y) with at least one curve in the (logy, log(kr)) plane. Consider two points, (y1,kr1) and (y2,kr2) on that relation from Figure 1. If they are close enough, we might approximate that part of the curve with a straight line, fitting that line to [15] to get h and yd [16]. It happens that for a BrooksCorey conductivity relation Kh has an analytic expression [17]. Let us use [14] and [15] to make the approximation for the scaling parameter, u [18]. If we substitute wvo [13] for wv and the fitted BrooksCorey Kh in [12], we have an approximation for Kv. Call it Kpbcd. Note that while we are working with relative conductivities here, one can use the full conductivity in most of these equations and still maintain proper dimensions, but [18] still requires relative conductivities.
[16]
[17]
[18]
Consider now the infiltration problem in [19], using soil parameters and relations from Haverkamp, et al., (1977). The conductivity relation can be called a vanGenuchtenstyle form, which differs from the BrooksCorey form in having a rounded knee where the diagonal and horizontal asymptotes cross in a loglog plot of kr and y (y = h). The parameters a and A indicate "displacement pressures" for each relation of 36.936 cm and 19.081cm, respectively, indicating a moderately to coarsely textured soil. If [13] and [14] hold for this soil on the kr1=kr2 diagonal, then they may be the basis for an approximation to Darcian means in a running model. The initial condition of h = 929.8 cm corresponds to a value of kr(y) = 10^{8} cm/s.
[19]
First we need to know the true Darcian means for this soil, derived from the solutions to [7]. Figures 2ae show the solutions to [7] for Kv using the kr(y) relation in [19] for a range of kr1 and kr1 from 1(10^{8}) to 1 in log^{3}space, (log(kr1), log(kr2), log(Kv)), on a 37 by 37 grid. The bottom right axis is log(kr1), marked in ratio steps of 1:10^{2}, matching the steps on the bottom left, log(kr2), and vertical, log(Kv), axes. The kr1=kr2 diagonal is the line from the bottom center to the top center. Conditions are wetoverdry on the left, and dryoverwet on the right. Recall that kr1 and kr2 are the grid point conductivities in Figure 1, and Kv is the true Darcian mean for those conductivities at vertical spacing, Dx. For comparison, Figure 3 is the arithmetic mean, Ka, on the same scale.










We can see from the plots that Kh (Figure 2a, approximately) is nothing like Ka. Near saturation, the rounded knee of kr(y) produces a sharp rise in the Kv » Kh surface, which is still far below the surface for Ka. So for small space steps, Ka might overestimate the flow if used for Kv in [4]. For large space steps (Figure 2e), Kv has almost completely transformed to kr2, where Ka would underestimate flow for wetoverdry and vastly underestimate flow for dryoverwet. But for small Dx (Dx < yd), the adjacent grid point conductivities are likely to be close together, and thus near the kr1=kr2 diagonal, where all valid interblock conductivity means approach each other, ameliorating the overestimation of flow. For large Dx (Dx > yd), the adjacent grid point conductivities are more likely to be farther apart, especially for sharp fronts. So for infiltration, as in [19], we should expect using Ka for Kv to produce an underestimation of flow by up to a factor of two.
While Figures 2 & 3 are good for an intuitive understanding of the process, they do not show precisely the magnitudes of error over these large dynamic ranges. It works better to look at the log ratio of the approximation to Kv, for example log(Ka/Kv). If one looks directly up the kr1=kr2 diagonal in Figures 2 & 3, from bottom to top center, the normal coordinate in logspace is log(kr1/kr2), keeping wetoverdry and dryoverwet on the left and right, respectively.
If the values of Ka in Figure 3 are divided by the values of Kv in Figure 2a at each point, and plotted separately against log(kr1/kr2), we have Figure 4a, for Dx = 0.2 cm. Figures 4bd similarly correspond to Dx = 20, 200 and 2000 cm, respectively. The center (0,0) point is the end view of the kr1=kr2 diagonal. Here we can see in Fig. 4a that Ka overestimates Kv by an order of magnitude at log(kr1/kr2) = ± 4. The pair of slightly lower lines of circles come from the grid points with kr1 or kr2 = 1, at saturation. We can see in Fig. 4c that Ka and Kv intersect in wetoverdry conditions. In Fig. 4d, the underestimation for wetoverdry and gross overestimation for dryoverwet conditions is clearly visible. This kind of plot is useful in selecting candidates for approximate Darcian means.








If we apply the suggested PBCD approximation to this soil, we get log ratios of Kpbcd/Kv as plotted in Figures 5ad. It turns out that this approximation does not do well where (y1,kr1) and (y2,kr2) are widely separated across rounded knee of kr(y) in [19], where one is near saturation and the other is dry. If there were room here for plots of Kpbcd(kr1,kr2) as in Figures 2ae, they would show the surface plot bending down near saturation, instead of up. This is the source of points in Fig. 5ad that loop below the main body. Fortunately, we have chosen the inflow saturation in [19] below the knee, which avoids that difficulty and still makes necessary points.
Notice that the ratio of Kpbcd/Kv stays closer to the ideal value of 1 farther from the kr1=kr2 diagonal (the (0,0) center point in the plots) than does Ka/Kv. This should mean less flow error in an infiltration model. In particular, for large vertical space steps (Figures 4d and 5d), notice how the Ka ratio drops more sharply through the kr1=kr2 diagonal while the Kpbcd ratio is much flatter.








Let us first solve [19] numerically with the hbased and the massconservative modified Picard methods detailed in Celia, et al., (1990), to establish the difference between those two methods, both using the arithmetic mean, Ka, for a variety of vertical model step sizes, Dx. This example differs from the first in Celia, et al., with a much drier initial condition, and larger space steps and infiltration times. Also mass conservation in the modified Picard method is maintained here to one part in 10^{8}. Here, we chose the length of the problem, xl (cm) to be as large as necessary to contain the modeled wetting front at the end time, tl = 10^{5} seconds.
Figure 6 shows the total infiltrated mass, for Dx = 0.8 cm to 2000 cm, divided by the estimated total infiltrated mass for Dx = 0. The estimate is obtained by polynomial extrapolation from the values for the massconservative method. We can see that mass conservation reduces the error in total mass infiltrated from about 2.5% to an abitrarily small amount from Dx = 0 to about 20 cm. Above that, the amount for the massconservative method falls of until it merges with the headbased method at 150 cm. They both lose, or misplace, equally large amounts of mass at Dx = 1500 cm, nearly 50%. Note that these losses occur for Dx significantly larger than the yd parameter in the conductivity relation.


Figure 7a shows the same curve for the massconservative method with the arithmetic mean extended up to Dx = 10^{5} cm. Between Dx = 150 and 1500 cm, the method failed to converge for the values tried. It shows that the infiltration does indeed lose 50% of the expected mass as Dx gets large. Figure 7b shows the resulting final moisture profiles at t = tl for x = 1000 to 3000 cm. The profile for Dx = 1500 cm shows as a single line segment, because the spreadsheet package used will not plot to points outside the local plot area. For several values of Dx below 200 cm, the plots show nonphysical oscillations that rise above the infiltration water content, increasing with Dx, suggesting the reason why some runs failed.




Figures 8ab show similar plots for the same model, with mass balance again good to one part in 10^{8}, using the piecewise BrooksCorey Darcian mean approximation, Kpbcd. Note that the model ran for all the values of Dx tried, and that the oscillations have disappeared. Where using the arithmetic mean lost 50% of the total mass infiltrated, the PBCD mean reduces that error to less than 0.5%, an improvement of two orders of magnitude.




There is not room to show it here, but the "upstream" mean for infiltration, replacing Kv with kr2, produces plots similar to Figures 8ab, with about twice the error. This might explain the popularity of such means in petroleum geophysics (Oldenberg and Pruess, 1993), but it should be used with caution because Figures 2ae show that it is not appropriate for dryoverwet conditions in soil. However, since it is a much faster to calculate, it is a viable candidate.
These results demonstrate that even in a modern, massconservative infiltration model with mass balance good to one part in 10^{8}, a standard arithmetic interblock conductivity mean can produce both nonphysical oscillations and a loss or misplacement of total mass infiltrated of up to 50% with large vertical model space steps. For one example of infiltration into a moderate to coarse soil, an imperfect approximation to true Darcian means demonstrates significant reduction in oscillatory behavior and an improvement of up to two orders of magnitude in the loss of mass infiltrated. This demonstrates that the common perception that modern massconservative models are "perfect" (Celia, et al., 1990), and thus no longer subject to improvement, may be somewhat overstated. This may have implications for fracture flow models at such places as the Yucca Mountain Project, where the model space step sizes can be orders of magnitude greater then the displacement pressure head length. It at least suggests that the measurement of unsaturated flow parameters may be more critical in such projects than previously thought.
This work was partly supported by a grant from the US Department of Energy, DEFG0297ER82329. Although none of the contents of this paper represent any official position or endorsement of the Department of Energy, the lead author was very grateful to receive support to do this work. The lead author also wishes to express his profound appreciation for the pregrant moral support and employment given by Dr. H. Don Scott of the Agronomy Department of the University of Arkansas, Fayetteville.
Baker, D.L., 1993, A secondorder, diagonallyimplicit, RungeKutta timestepping method, GROUND WATER, NovDec 1993, 31(6):890895.
Baker, D.L., 1995a, Darcian weighted interblock conductivity means for vertical unsaturated flow, GROUND WATER, MayJun 1995, 33(3):385390.
Baker, D.L., 1995b, Applying higherorder DIRK time steps to the "modified Picard" method, GROUND WATER, MarApr 1995, 33(2):259263.
Baker, D.L., 1998, Developing Darcian means in application to Topopah Spring welded volcanic tuff, Report DOE/ER/823292 to the U.S. Department of Energy, June 1998, 118p.
Celia, M.A., E.T. Bouloutas and R.L. Zarba, 1990, A general massconservative numerical solution for the unsaturated flow equation, Water Resour. Res., 26(7):14831496.
Haverkamp, R., M. Vauclin, J. Touma, P. Wierenga and G. Vachaud, Comparison of numerical simulation models for onedimensional infiltration, Soil Sci. Soc. Am. J., 41:285294.
Oldenberg, C.M. and K. Pruess, 1993, On numerical modeling of capillary barriers, Water Resour. Res., 29(4):10451056.
Warrick, A.W., 1991, Numerical approximation of Darcian flow through unsaturated soil, Water Resour. Res., 27(6):12151222.