draft of

Some Analytical and Approximate Darcian Means

Abstract

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 Brooks-Corey Darcian (PBCD) mean. It compares the PBCD mean to the arithmetic mean in a modern mass-conservative numerical model (Celia, et al., 1990) for vertical space steps from 0.8 to 10⁵ 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⁸, non-Darcian 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.

Introduction

Although Warrick (1991) demonstrated that standard interblock conductivity means, such as the arithmetic, geometric and upstream, produce non-Darcian 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 Brooks-Corey 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

Brief Review of Darcian Means

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).

Figure 1: Two adjacent grid cells in a finite difference model of infiltration

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 right-hand-side 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 non-physical. 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⁵ or 10⁶, we have Kv from [4].

The Darcian Mean of the Exponential Form

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.

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 Brooks-Corey and van-Genuchten-style forms.

A Piecewise Brooks-Corey Darcian (PBCD) Mean Approximation

It currently appears that all of the more complicated relations of kr(y), such as the Brooks-Corey [15] or van-Genuchten 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 Brooks-Corey 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 Brooks-Corey 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]

A Test Case - Haverkamp's Soil

Consider now the infiltration problem in [19], using soil parameters and relations from Haverkamp, et al., (1977). The conductivity relation can be called a van-Genuchten-style form, which differs from the Brooks-Corey form in having a rounded knee where the diagonal and horizontal asymptotes cross in a log-log 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.

First we need to know the true Darcian means for this soil, derived from the solutions to [7]. Figures 2a-e 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³-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², 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 wet-over-dry on the left, and dry-over-wet 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.

Figure 2a: Kv for Dx = 0.2 cm

Figure 2b: Kv for Dx = 2 cm

Figure 2c: Kv for Dx = 20 cm

Figure 2d: Kv for Dx = 200 cm

Figure 2e: Kv for Dx = 2000 cm

Figure 3: Ka = (kr1 + kr2)/2

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 wet-over-dry and vastly underestimate flow for dry-over-wet. 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 log-space is log(kr1/kr2), keeping wet-over-dry and dry-over-wet 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 4b-d 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 wet-over-dry conditions. In Fig. 4d, the underestimation for wet-over-dry and gross overestimation for dry-over-wet conditions is clearly visible. This kind of plot is useful in selecting candidates for approximate Darcian means.

Figure 4a: Log(Ka/Kv) vs Log(kr1/kr2) for Dx = 0.2 cm	Figure 4b: Log(Ka/Kv) vs Log(kr1/kr2) for Dx = 20 cm
Figure 4c: Log(Ka/Kv) vs Log(kr1/kr2) for Dx = 200 cm	Figure 4d: Log(Ka/Kv) vs Log(kr1/kr2) for Dx = 2000 cm

If we apply the suggested PBCD approximation to this soil, we get log ratios of Kpbcd/Kv as plotted in Figures 5a-d. 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 2a-e, they would show the surface plot bending down near saturation, instead of up. This is the source of points in Fig. 5a-d 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.

Figure 5a: Log(Kpbcd/Kv) vs Log(kr1/kr2) for Dx = 0.2 cm	Figure 5b: Log(Kpbcd/Kv) vs Log(kr1/kr2) for Dx = 20 cm
Figure 5c: Log(Kpbcd/Kv) vs Log(kr1/kr2) for Dx = 200 cm	Figure 5d: Log(Kpbcd/Kv) vs Log(kr1/kr2) for Dx = 2000 cm

Results

Let us first solve [19] numerically with the h-based and the mass-conservative 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⁸. 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⁵ 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 mass-conservative 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 mass-conservative method falls of until it merges with the head-based 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 6: Ratio of total infiltrated mass, for a given vertical space step size, Dx, to the estimated mass infiltrated for dx = 0. Head-base line passes through circled points. Mass-conservative, theta-based line starts closer to ideal for small Dx.

Figure 7a shows the same curve for the mass-conservative method with the arithmetic mean extended up to Dx = 10⁵ 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 non-physical oscillations that rise above the infiltration water content, increasing with Dx, suggesting the reason why some runs failed.

Figure 7a: Relative mass infiltrated versus vertical space step, Dx (cm)	Figure 7b: Water content profiles for near the wetting front for Dx = 0.8 to 1500 cm

Figures 8a-b show similar plots for the same model, with mass balance again good to one part in 10⁸, using the piecewise Brooks-Corey 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.

Figure 8a: Relative mass infiltrated versus vertical space step, Dx (cm)	Figure 8b: Water content profiles for near the wetting front for Dx = 0.8 to 1500 cm

There is not room to show it here, but the "upstream" mean for infiltration, replacing Kv with kr2, produces plots similar to Figures 8a-b, 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 2a-e show that it is not appropriate for dry-over-wet conditions in soil. However, since it is a much faster to calculate, it is a viable candidate.

Conclusions

These results demonstrate that even in a modern, mass-conservative infiltration model with mass balance good to one part in 10⁸, a standard arithmetic interblock conductivity mean can produce both non-physical 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 mass-conservative 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.

Acknowledgments

This work was partly supported by a grant from the US Department of Energy, DE-FG02-97ER82329. 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 pre-grant moral support and employment given by Dr. H. Don Scott of the Agronomy Department of the University of Arkansas, Fayetteville.

References